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