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

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.

• Cox proportional hazard model is commonly used for medical studies (e.g. spc.univ-lyon1.fr/lecture-critique/TP/travail_sur_article/…) . The plot generally mentions the number of cases and events at different time points of follow-up. The number of cases reduces as time of follow-up increases.
– rnso
May 13 '15 at 1:09

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.

• Just to get it clear. You are saying that if the times of censoring for different classes are not all the same then we should not assume that the censoring mechanism is independent of event time? May 13 '15 at 11:54
• I've edited the answer, hope the explanation is more clear now. I'm saying that if we included class as a covariate we can still assume independent right-censoring - explanation is found above.
– swmo
May 13 '15 at 12:17
• Thanks @swmo. I've included a picture of a model I've run with R to clear my question. Coefficients are negative and decrease because of shorter and shorter follow-up times! May 13 '15 at 15:41
• I've added some simulations and a guess at what's happening.
– swmo
May 13 '15 at 16:48
• Thanks again, and thank you too @Joshua. There is indeed a cohort effect. The exam got tougher as the years went by, so the negative coefficients are plausible, I was just scared that the censoring mechanism could in some way inflate them. May 13 '15 at 17:14

In survival analysis, the "time of study" isn't very relevant. What is relevant is how long the participant was in the study. As @rnso alludes and provides a link for, in medical studies, it is extremely uncommon for all the participants to start at the same time. Enrollment is on a rolling basis, meaning patients are enrolled for a certain period of time, say 6 months. Then they are followed for some period of time, say 2 years. That means that patients enrolled on Day 1 can be in the study for up to 2.5 years while patients enrolled on Day 180 can be in the study for up to 2 years. If anyone is still in the study at the end of the time, they are censored.

Translating to your situation, what is input into the survival analysis is the time the student entered the study until either 1) they pass the exam [event happened where $time = ExamDate - EntranceTime$], 2) they drop out [censored at $time = ExitTime - EntranceTime$], or 3) 2013 happens and everyone still in the study leaves it [censored at $time = 2013 - EntranceTime$]. So your students from 2009 will give you the most information since they can be on the study the longest, but the 2013 students can still be useful. In survival analysis, those 2013 students will either complete the exam before a year is up or they are censored after 1 year, since the study ends.

To answer your question directly, the time of study being different for each student is factored into the time you calculate for each student as they are on the study. The 2013 students only provide you a year's worth of information before they are censored, which is accounted for by none of those students having a time value longer than 1 year.