# Predicting new data using a Case-cohort Proportional Hazards Model

I am fitting a Case-Cohort Proportional Harzards (CCH) model using the Survival (version 2.37-4) package in R 2.15.3. Normally with a Cox Proportional Harzards (coxPH) model, I can use the survfit function to predict the probability of a new patient being survival free at a specific time point; however, there is no survfit function for the CCH model in the survival package. How do I generate survival predictions on new data using a CCH model?

• Did you choose a particular method in the cch function or did you leave as the default Prentice? Feb 26, 2014 at 9:39
• I've left it as default. Feb 26, 2014 at 21:59
• I am interested and will investigate. No promises that I will be of any help though! Feb 26, 2014 at 22:02
• Can you post the cch command you are using? Feb 27, 2014 at 14:03

OK, here is an attempt. According to Self and Prentice (1988) there is a natural baseline cumulative hazard function:

$$\tilde{\Lambda}(t)=\tilde{n}n^{-1}\int_0^t{\frac{d\bar{N}(u)}{\sum_{i\in\tilde{C}}Y_i(u)r\{\beta^\prime Z_i(u)\}}}$$

Where

• $\tilde{n}$: size of the subcohort
• $n$: size of total cohort
• $\tilde{C}$ is the subcohort
• $\bar{N}(u)$: counting process for number of failures up to and including time $u$
• $Y_i(u)$: indicator for whether subject $i$ is at risk at time $u$
• $r$: fixed function
• $\beta$: coefficients
• $Z_i(u)$: covariates for subject $i$ at time $u$

This can be simplified if the covariates are constant for each subject, i.e., $Z_i(u)=Z_i$ since then $r\{\beta^\prime Z_i(u)\}=v_i$. I assume that $r\{\cdot\}$ is the exponential function.

To simplify we then index all the failure times $t_j$, define $\bar{N}_j$ as the number of failures up to and including time $t_j$ and let $\bar{N}_0=0$. We let $Y_{ij}$ denote whether subject $i$ is at risk up to time $t_j$ such that $Y_{ij}=1$ indicates that subject $i$ could be observed to fail at time $t_j$.

We now obtain the simplified formula:

$$\tilde{\Lambda}(t) = \tilde{n}n^{-1}\sum_{\{j:\,t_j\le t\}}{\frac{\bar{N}_j-\bar{N}_{j-1}}{\sum_{i\in\tilde{C}}{Y_{ij}v_i}}}$$

I am not a regular R user but I have attempted to give an example, please excuse the no-doubt risible code (note I edited to fix $\tilde{C}$ in the calculation of $Y_{ij}$):

library(survival)

subcoh <- nwtco$in.subcohort selccoh <- with(nwtco, rel==1|subcoh==1) ccoh.data <- nwtco[selccoh,] ccoh.data$subcohort <- subcoh[selccoh]
ccoh.data$histol <- factor(ccoh.data$histol,labels=c("FH","UH"))
ccoh.data$stage <- factor(ccoh.data$stage,labels=c("I","II","III","IV"))
ccoh.data$age <- ccoh.data$age/12 # Age in years
fit.ccP <- cch(Surv(edrel, rel) ~ stage + histol + age, data =ccoh.data,
subcoh = ~subcohort, id=~seqno, cohort.size=4028)

## Extract the times at which any in whole cohort relapsed
t_rel <- c(0,unique(sort(ccoh.data$edrel[ccoh.data$rel==1])))

## Count the number relapsed up to each time
N <- numeric(length=length(t_rel))
for (j in 1:length(t_rel)) {
N[j] <- sum(ccoh.data$rel * (ccoh.data$edrel <= t_rel[j]))
}

## Compute individual risk factors v_i
v <- numeric(length=dim(ccoh.data)[1])
for (i in 1:dim(ccoh.data)[1]) {
coef <- fit.ccP$coef['XstageII'] * (ccoh.data[i,'stage'] == 'II') + fit.ccP$coef['XstageIII'] * (ccoh.data[i,'stage'] == 'III')
+ fit.ccP$coef['XstageIV'] * (ccoh.data[i,'stage'] == 'IV') + fit.ccP$coef['XhistolUH'] * (ccoh.data[i,'histol'] == 'UH')
+ fit.ccP$coef['Xage'] * ccoh.data[i,'age'] v[i] <- exp(coef) } ## Compute observation matrix Y_ij Y <- array(0, c(dim(ccoh.data)[1], length(t_rel))) for (i in 1:dim(ccoh.data)[1]) { for (j in 1:length(t_rel)) { Y[i,j] <- ccoh.data$subcohort[i] & (t_rel[j] <= ccoh.data$edrel[i]) # EDITED FROM ORIGINAL POST } } ## Compute baseline hazard baseline_cumhaz <- numeric(length=length(t_rel)) n_cohort <- 4028 n_subcohort <- sum(ccoh.data$subcohort)
n_ratio <- n_subcohort / n_cohort
baseline_cumhaz[1] <- 0
for (j in 2:length(t_rel)) {
baseline_cumhaz[j] <- baseline_cumhaz[j-1] + n_ratio * (N[j] - N[j-1]) / sum(Y[,j]*v)
}

## Compute baseline survival
survival <- numeric(length=length(t_rel))
prev_S <- 1
prev_H <- 0
for (i in 1:length(t_rel)) {
dH <- baseline_cumhaz[i] - prev_H
survival[i] <- prev_S * (1 - dH)
prev_H <- baseline_cumhaz[i]
prev_S <- survival[i]
}

## Compare survival across stages assuming FH histology and age 40
patient.stageI.v <- exp(fit.ccP$coef['Xage']*40.0) patient.stageII.v <- exp(fit.ccP$coef['Xage']*40.0 + fit.ccP$coef['XstageII']) patient.stageIII.v <- exp(fit.ccP$coef['Xage']*40.0 + fit.ccP$coef['XstageIII']) patient.stageIV.v <- exp(fit.ccP$coef['Xage']*40.0 + fit.ccP\$coef['XstageIV'])
patient.stageI.survival <- numeric(length=length(t_rel))
patient.stageII.survival <- numeric(length=length(t_rel))
patient.stageIII.survival <- numeric(length=length(t_rel))
patient.stageIV.survival <- numeric(length=length(t_rel))
patient.stageI.survival[1] <- 1.0
patient.stageII.survival[1] <- 1.0
patient.stageIII.survival[1] <- 1.0
patient.stageIV.survival[1] <- 1.0
for (i in 2:length(t_rel)) {
dH <- baseline_cumhaz[i] - baseline_cumhaz[i-1]
patient.stageI.survival[i] <- patient.stageI.survival[i-1] * (1 - patient.stageI.v * dH)
patient.stageII.survival[i] <- patient.stageII.survival[i-1] * (1 - patient.stageII.v * dH)
patient.stageIII.survival[i] <- patient.stageIII.survival[i-1] * (1 - patient.stageIII.v * dH)
patient.stageIV.survival[i] <- patient.stageIV.survival[i-1] * (1 - patient.stageIV.v * dH)
}

plot(t_rel, patient.stageI.survival, type='s', ylim=c(0, 1))
lines(t_rel, patient.stageII.survival, type='s', lty=2)
lines(t_rel, patient.stageIII.survival, type='s', lty=3)
lines(t_rel, patient.stageIV.survival, type='s', lty=4)


The same formula is used in Prentice (1986) so there shouldn't be any changes to make on that basis. Also be aware the code has not undergone rigorous checking, it would be worth at least running it on simulated data as a sanity check.

• Thank you for the answer, I'm going to need a little time to go through it and try it. Once I do, if everything checks out I'll mark it as the correct answer. Mar 2, 2014 at 1:24