1
$\begingroup$

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?

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

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – Nixuz
    Mar 2, 2014 at 1:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.