The Stack Overflow podcast is back! Listen to an interview with our new CEO.

 5 Added a better graph and updated the data simulation script to a more flexible version edited Sep 4 '11 at 17:18 Max Gordon 3,49644 gold badges2626 silver badges4444 bronze badges library("cmprsk") # The time for the study accrual_time <- 10 followup_time <- 1 base_risk <- list("event" = .015, "cmprsk" = .1) risk_factor_1risk_factors <- list(list("frequency"=.1, "event" = base_risk$$event*.5, "cmprsk" = base_risk$$$$event*.5, "cmprsk" = base_risk$$cmprsk*2), risk_factor_2 <- list("frequency"=.05, "event" = base_risk$$event*1, "cmprsk" = base_risk$$$$event*1, "cmprsk" = base_risk$$cmprsk*1), risk_factor_3 <- list("frequency"=.05, "event" = base_risk$$event*-.5, "cmprsk" = 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=3ncol=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){ x[i,# 1]Set base risk event_risk <- rbinombase_risk$$event comp_risk <- base_risk$$cmprsk for(1,j 1,in risk_factor_1$$frequency) x[i, 2] <- rbinom(1, 1, risk_factor_2$$frequency1:length(risk_factors)){ x[i, 3]j] <- rbinom(1, 1, risk_factor_3$frequencyrisk_factors[[j]]$frequency) # Add risk factors # event_riskIf <-there base_risk$$event + x[i, 1]*risk_factor_1$$eventis + a risk factor defined if (x[i, 2]*risk_factor_2$$event + x[i, 3]*risk_factor_3$$event j] > 0){ # Add risk factors comp_risk event_risk <- base_risk$$cmprsk + x[i, 1]*risk_factor_1$$cmprskevent_risk + x[i, 2]*risk_factor_2$$cmprsk + x[i, 3]*risk_factor_3$$ 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) time_2_comprsk[i] <- rexp(1, rate=comp_risk) } colnamescn <- c(x) for(i in 1:length(risk_factors)){ ev_rsk <- crisk_factors[[i]]$$event/base_risk$$event+1 cmp_rsk <- risk_factors[[i]]$$cmprsk/base_risk$$cmprsk+1 name <- paste("RF"Risk 1"factor no: ", "RFi, 2""\n * ev=", "RFev_rsk, 3"" 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)  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) x[i, 2] <- rbinom(1, 1, risk_factor_2$$frequency) x[i, 3] <- rbinom(1, 1, risk_factor_3$frequency) # 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) time_2_comprsk[i] <- rexp(1, rate=comp_risk) } 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) 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) # 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) time_2_comprsk[i] <- rexp(1, rate=comp_risk) } 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) 4 added quasipoisson edited Sep 4 '11 at 11:06 Max Gordon 3,49644 gold badges2626 silver badges4444 bronze badges I've been suggested to do a poissonPoisson 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: Is the glm() code correct or should I somehow transform my data? Does the poissonPoisson output make any sense and how should if so interpret it? What are the benefits/pitfalls in using poissonPoisson regression for survival data?  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 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: 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? 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? 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: 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? 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 3 added a graph showing the outcome and some tests for over-dispersion edited Sep 4 '11 at 10:54 Max Gordon 3,49644 gold badges2626 silver badges4444 bronze badges 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:  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: 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: 2 added a graph showing the outcome and some tests for over-dispersion edited Sep 4 '11 at 10:41 Max Gordon 3,49644 gold badges2626 silver badges4444 bronze badges 1 asked Sep 3 '11 at 17:49 Max Gordon 3,49644 gold badges2626 silver badges4444 bronze badges