We’re rewarding the question askers & reputations are being recalculated! Read more.
3 of 5 added a graph showing the outcome and some tests for over-dispersion

Poisson regression in a survival setting on a simulated set

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_factor_1 <- list("frequency"=.1, 
                      "event" = base_risk$event*.5, 
                      "cmprsk" = base_risk$cmprsk*2)
risk_factor_2 <- list("frequency"=.05, 
                      "event" = base_risk$event*1, 
                      "cmprsk" = base_risk$cmprsk*1)
risk_factor_3 <- 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=3, nrow=n)
time_2_event <- rep(0, n)
time_2_comprsk <- rep(0, n)
for(i in 1:n){
    x[i, 1] <- rbinom(1, 1, risk_factor_1$frequency)[1]
    x[i, 2] <- rbinom(1, 1, risk_factor_2$frequency)[1]
    x[i, 3] <- rbinom(1, 1, risk_factor_3$frequency)[1]

    # Add risk factors
    event_risk <- base_risk$event + 
            x[i, 1]*risk_factor_1$event +
            x[i, 2]*risk_factor_2$event +
            x[i, 3]*risk_factor_3$event

    # Add risk factors
    comp_risk <- base_risk$cmprsk + 
            x[i, 1]*risk_factor_1$cmprsk +
            x[i, 2]*risk_factor_2$cmprsk +
            x[i, 3]*risk_factor_3$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]
}

colnames(x) <- c("RF 1", "RF 2", "RF 3")

# 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:

A forestplot comparing the different methods - poisson: 1.152  1.509  0.794, crr: 1.151 1.524 0.812, cox ph: 1.897 1.931 0.798

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?