Solved – Different “end of study” times for different cohorts – Cox PH model in survival analysis

cox-modelstratificationsurvival

I have a dataset with 4 cohorts of about the same size (~700 people each). I'm trying to apply a Cox PH model using the time needed to pass a very difficult exam as my "time" variable. The cohorts differ because they are different classes (class of 2009, class of 2010, 2011 and 2012). These are clearly also the time they enter in the study.

All times are censored after 2013. Is there any way of accounting for the fact that the "time of study" is different for each cohort? I was thinking of stratifying for cohorts, but the coefficients would clearly be negative and would decrease faster and faster. This however would be due not to a real decrease in the hazard, but to the fact that the latter cohorts are monitored for less time.

To clear up my last point, I will include some code I've just run in R.

![R code]http://i58.tinypic.com/25txbfa.jpg

cox2 <- coxph(ml ~ as.factor(Cohort), data = data)
summary(cox2)
# Call:
# coxph(formula = m1 ~ as.factor(Cohort), data = data)
#   n= 1865, number of events= 621
#    (63237 observations deleted due to missingness)
# 
#                           coef exp(coef) se(coef)     z Pr(>lzl)
# as.factor(Cohort)2010 -0.23715  0.78887  0.09601 -2.470  0.01351 *
# as.factor(Cohort)2011 -0.35811  0.69899  0.12345 -2.901  0.00372 **
# as.factor(Cohort)2012 -0.76813  0.46388  0.17111 -4.489 7.15e-06 ***
# ---
# Signif. codes: 0 ‘***' 0.001 ‘**' 0.01 ‘*' 0.05 ‘.' 0.1 ‘ ' 1
# 
#                         exp(coef) exp(-coef) lower .95 upper .95
# as.factor(Cohort)2010      0.7889      1.268    0.6536    0.9522
# as.factor(Cohort)2011      0.6990      1.431    0.5488    0.8903
# as.factor(Cohort)2012      0.4639      2.156    0.3317    0.6487
# 
# Concordance= 0.568 (se = 0.012 )
# Rsquare= 0.013 (max possible= 0.99 )
# Likelihood ratio test= 24.5 on 3 df, p=1.968e-05
# Wald test            = 23.31 on 3 df, p=3.482e-05
# Score (logrank) test = 23.85 on 3 df, p=2.683e-05

As I expected, the coefficients are all negative and fastly decreasing. I'm sure that this is due not to a real decrease in hazard, but to the fact that follow-up times for each cohort is decreasing as the class covariate increases.

Best Answer

Basically, we're assuming a Cox model with four groups. I'm not sure what you mean by stratifying for class as this wouldn't help you compare the four classes.

The question is whether or not we can use a standard Cox model when we're applying different censoring mechanisms to different groups of observations. And actually, we can.

An assumption of the model is independent right-censoring. What is actually meant by this, is independence of the censoring times and the event times conditional on the covariates of the model. In your example, we have to determine if it's reasonable to assume independent right-censoring. When we condition on the covariates (class being one of them), the censoring times are constants, thus independent of the event times. Therefore, the censoring times and event times are actually independent conditional on the covariates. When you know the value of the class covariate, the dependence between censoring time and event time vanishes, so to speak. So in your example, the assumption of independent right-censoring is still reasonable.

This is of course only true when the class covariate that influences the censoring is actually included in the model. Imagine that you didn't include class and furthermore assume that class and event time are not independent. Then (heuristically) information on the censoring time would give you information on the class, which in turn would give you information on the event time. Thus, event times and censoring times are not independent any longer and an assumption of the model just broke down.

A reference on this would be PK Andersen et al., Statistical Models based on Counting Processes, 1997, pp. 139-146.

A small simulation experiment might shed some light on the situation. The following R-code generates data from four groups, all with the same hazard.

data <- data.frame(grp = c(rep(1, 1000), rep(2, 1000), rep(3, 1000), rep(4, 1000)))

data$x <- rweibull(n = 4000, shape = 1, scale = 1)


data$event <- data$x < (data$grp == 1)*1 + (data$grp == 2)*.7 + 
  (data$grp == 3)*.5 + (data$grp == 4)*.3
data$x <- pmin(data$x, (data$grp == 1)*1 + (data$grp == 2)*.7 + 
                 (data$grp == 3)*.5 + (data$grp == 4)*.3)


coxph(Surv(time = x, event = event) ~ as.factor(grp), data = data)

If you run the code, you'll see that the estimates you get are good, we get coefficients close to one.

Without the data it's difficult to say why you're seeing this artefact of decreasing coefficients for increasing class. Off the top of my head:

What you're seeing might be due to a kind of selection bias. For each year some of the 700 people will not ever pass the exam (I'm guessing). If we imagine that a student fails the exam in a early year, he or she might be removed from the cohort all together. For instance, because the student quits the programme and therefore is no longer in the process of passing the exam. This would clearly constitute a selection bias as students that are passing the exam slowly (or never passing it) are removed from the early cohorts whereas in the later cohorts they have less time to quit the programme. That would give the estimates you report, even if there's no real differences between classes. See the following simulation experiment. It's identical to the above, but I'm just removing slow students, and more so in early years (this corresponds to drop-outs).

data <- data.frame(grp = c(rep(1, 1000), rep(2, 1000), rep(3, 1000), rep(4, 1000)))

data$x <- rweibull(n = 4000, shape = 1, scale = 1)

data <- subset(data, x < (data$grp == 1)*1 + (data$grp == 2)*1.3 + 
                 (data$grp == 3)*1.5 + (data$grp == 4)*1.7)

data$event <- data$x < (data$grp == 1)*1 + (data$grp == 2)*.7 + 
  (data$grp == 3)*.5 + (data$grp == 4)*.3
data$x <- pmin(data$x, (data$grp == 1)*1 + (data$grp == 2)*.7 + 
                 (data$grp == 3)*.5 + (data$grp == 4)*.3)


coxph(Surv(time = x, event = event) ~ as.factor(grp), data = data)

It goes without saying that this is just a guess of what could be happening with your data. Assuming that the negative coefficients are an artefact of some sort, it's not caused by the censoring mechanism, though.

Related Question