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
source | link
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)[1]
    x[i, 2] <- rbinom(1, 1, risk_factor_2$frequency1:length(risk_factors)[1]){
        x[i, 3]j] <- rbinom(1, 1, risk_factor_3$frequencyrisk_factors[[j]]$frequency)[1]

    # 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)[1]
    time_2_comprsk[i] <- rexp(1, rate=comp_risk)[1]
}

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)

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.798A 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

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)

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

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)

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

4 added quasipoisson
source | link

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?

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.798A 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

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?

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

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?

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

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
source | link

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: CRR, Cox PH & Poisson regression

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

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: CRR, Cox PH & Poisson regression

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

2 added a graph showing the outcome and some tests for over-dispersion
source | link
1
source | link