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?
 A: 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)

## Follow example for cch
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.
