Poisson regression in a survival setting on a simulated set

Question

I have a large dataset with patients and I'm studying a rare outcome (~ 2%) and death is a competing risk (mean age ~69 years). I've used the R "cmprsk" package for my statistics and it seems that competing risks and the Cox regression are performing similarly although the competing risk analysis is more conservative giving hazard ratios closer to 1.

I've been suggested to do a Poisson regression on the data but the results don't make any sense and I would be really grateful to get some input on the benefits of doing this kind of analysis on survival data. I've created this simulation for creating a dataset with similar risk factors:

library("cmprsk")
# The time for the study
accrual_time <- 10
followup_time <- 1

base_risk <- list("event" = .015, "cmprsk" = .1)

risk_factors <- list(list("frequency"=.1, 
                "event" = base_risk$event*.5, 
                "cmprsk" = base_risk$cmprsk*2),
        list("frequency"=.05, 
                "event" = base_risk$event*1, 
                "cmprsk" = base_risk$cmprsk*1),
        list("frequency"=.05, 
                "event" = base_risk$event*-.5, 
                "cmprsk" = base_risk$cmprsk*0))

# Number of subjects
n <- 5000

# Create base time, sequential inclusion
time_in_study <- rep(c(1:n)/n*accrual_time + followup_time, 1)

set.seed(100)

# Create empty sets
x <- matrix(0, ncol=length(risk_factors), nrow=n)
time_2_event <- rep(0, n)
time_2_comprsk <- rep(0, n)

# Create each studied observation and outcome
for(i in 1:n){
    # Set base risk
    event_risk <- base_risk$event 
    comp_risk <- base_risk$cmprsk

    for(j in 1:length(risk_factors)){
        x[i, j] <- rbinom(1, 1, risk_factors[[j]]$frequency)[1]

        # If there is a risk factor defined
        if (x[i, j] > 0){
            event_risk <- event_risk +
                    risk_factors[[j]]$event
            comp_risk <- comp_risk + 
                    risk_factors[[j]]$cmprsk
        }
    }

    # Time 2 event/risk is 1/rate meaning that higher number -> shorter time
    time_2_event[i] <- rexp(1, rate=event_risk)[1]
    time_2_comprsk[i] <- rexp(1, rate=comp_risk)[1]
}

cn <- c()
for(i in 1:length(risk_factors)){
    ev_rsk <- risk_factors[[i]]$event/base_risk$event+1
    cmp_rsk <- risk_factors[[i]]$cmprsk/base_risk$cmprsk+1
    name <- paste("Risk factor no: ", i, "\n * ev=", ev_rsk, " cr=", cmp_rsk, " *", sep="")
    cn <- c(cn, name)
}
colnames(x) <- cn

# Select the event that happens first: study ends, evenent occurs, a competing event occurs
time <- apply(cbind(time_in_study, time_2_event, time_2_comprsk), 1, min)

# Outcome identifiers
event <- (time_2_event == time) + 0
comprsk <- (time_2_comprsk == time) + 0
cens <- event+2*(event==0 & comprsk==1)

out.cox_ev <- coxph(Surv(time, event)~x)
summary(out.cox_ev)

out.crr_ev <- crr(time, cens, x, failcode=1)
summary(out.crr_ev)

out.cox_cmprsk <- coxph(Surv(time, comprsk)~x)
summary(out.cox_cmprsk)

out.crr_cmprsk <- crr(time, cens, x, failcode=2)
summary(out.crr_cmprsk)

The output makes sense but when I do a:

out.glm_pr <- glm(event ~ x, family="poisson")
summary(out.glm_pr)

It gives estimates of:

• RF 1 ~ .14
• RF 2 ~ .41
• RF 3 ~ -.23

My questions:

• Is the glm() code correct or should I somehow transform my data?
• Does the Poisson output make any sense and how should if so interpret it?
• What are the benefits/pitfalls in using Poisson regression for survival data?

Thanks!

---

UPDATE

After adding exp(out.glm_pr$coefficients) the results are almost identical to the competing risk regression, here's a forest plot that compares the three:

The x-axis is perhaps not entirely valid (should be "incident rate ratios" for the Poisson regression) but why are the outcomes for CRR & poisson almost identical?

As for testing over-dispersion I've found these two methods:

> library(qcc)
> qcc.overdispersion.test(event)

Overdispersion test Obs.Var/Theor.Var Statistic p-value
       poisson data         0.9391878      4695 0.99902
> 
> library(pscl)
> out.glm_nb <- glm.nb(event ~ x)
Warning messages:
1: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace >  :
  iteration limit reached
2: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace >  :
  iteration limit reached
> odTest(out.glm_nb)
Likelihood ratio test of H0: Poisson, as restricted NB model:
n.b., the distribution of the test-statistic under H0 is non-standard
e.g., see help(odTest) for details/references

Critical value of test statistic at the alpha= 0.05 level: 2.7055 
Chi-Square Test Statistic =  -0.0139 p-value = 0.5 

I conclude that there isn't any evidence of over-dispersion or are there other methods better suited for testing over-dispersion in this kind of survival data?

The quasipoisson analysis gives similar values:

> out.glm_quasi_pr <- glm(event ~ x, family=quasipoisson(link="log"))
> round(exp(out.glm_quasi_pr$coefficients), 3)
(Intercept)       xRF 1       xRF 2       xRF 3 
      0.059       1.152       1.509       0.794 

Answer 1

Answering two of your three questions, because I'm not comfortable enough in R to diagnose coding errors - though your code looks right to me.

For your estimates, as you didn't specify a link function, R is using the default log link function. In order for your estimates to make sense, you need to exp(estimate) - this will give you 1.15, 1.51 and 2.94 respectively. These should be close to the HRs coming off your Cox model. These numbers are "incident rate ratios" - the name is pretty suggestive of what they are. They can be interpreted very similarly to hazard ratios, and indeed should equal hazard ratios under certain assumption.

As for the benefits (and drawbacks) of Poisson regression. Poisson survival analysis is a fully parametric, maximum likelihood method of estimating differences in the survival between groups, which has some nice properties for some uses. It estimates the baseline hazard, which if you intend to use the baseline hazard in further analysis (I often do), is something a Cox model expressly does not estimate. The incidence density (# of cases / time at risk) is also vastly more intuitive than the hazard.

Now the drawbacks, of which there are many. The Poisson model is vulnerable to overdispersion - I'd rerun your model using quasipoisson or a negative binomial model to check and see if your results are sensitive to overdispersion. More important, IMO, is the assumptions the poisson model makes about the underlying survival function. The Cox Proportional Hazards model, as the name suggests, assumes the hazard function is proportional between the two groups over time. The Poisson model assumes not only are the hazards proportional, but constant. This is often a pretty major assumption, and should be checked. Adding a term or several in the model for time can help relax this assumption somewhat.

As for the similarity of your results: I'm not sure what kind of model your competing risk package is assuming, but I'm guessing it's estimating a parametric model of the survival function - the Cox model, which is "semi-parametric", doesn't estimate this baseline hazard. If the estimated parametric survival function is close to a exponential model (the distribution of the survival function assumed with a Poisson model), you may get very similar results if your data isn't sensitive to the assumption of no competing risks. Your results don't seem terribly vulnerable to this assumption.

What you do seem to be sensitive to is the violation of a constant hazard assumption, hence the dramatic difference in your estimates between Poisson and Cox models. If the hazard function was perfectly constant, the two models should actually produce the same estimate. I would try one of two things:

1. Add a time term or several to the Poisson model. Something like time and time*time, and see if your Poisson estimate moves closer to the Cox result.

2. You should be able to visualize the hazard function. The most common way to do this is to plot the log(-log(survival function)) versus the log of survival time for each of your variables, stratified by the groups. For the Cox model to be valid, they should be parallel. For the Poisson model to be valid, they should be straight and parallel. The survival functions in R must be able to do this, though I don't know how.

I'd find it very odd to see this in a simulated data set introduced essentially by accident, but I'd try looking at it anyway.