1 of 5

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!