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](http://i.imgur.com/kEF4a.png)

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?