How to generate predicted survivor curves from frailty models (using R coxph)?

Question

I want to compute predicted survivor function for a Cox proportional hazards model with frailty terms [using survival package]. It appears that when frailty terms are in the model, the predicted survivor function cannot be computed.

## Example 
require(survival)
data(rats)

## Create fake weight
set.seed(90989)
rats$weight<-runif(nrow(rats),0.2,0.9)

## Cox model with gamma frailty on litter
fit <- coxph(Surv(time, status) ~ rx+weight+frailty(litter,dist="gamma"),
data = rats) 

## Compute survival curve from the cox model for rx=0 and weight=0.5 kg
plot(survfit(fit, newdata=data.frame(rx=0,weight=0.5)),xlab = "time",
ylab="Survival") 

## Running this line, I get following error message:
Error in survfit.coxph(fit, newdata = data.frame(rx = 0, weight = 0.5)) : 
Newdata cannot be used when a model has sparse frailty terms

I tried using both sparse and non-sparse computation methods by using sparse=TRUE, Sparse =FALSE, sparse =0, sparse=5 options. However, none worked.

How do I compute predicted survivor curve based on my frailty model?

Answer 1

The problem here is the same as would be obtained trying to predict outcomes from a linear mixed effects model. Since the survival curve is non-collapsible, each litter in your example has a litter-specific survival curve according to the model you fit. A frailty as you may know is the same as a random intercept indicating common levels of confounding and prognostic variables endemic to each litter, presumably vis-à-vis genetic traits. Therefore the linear predictor for the hazard ratio is a mix of the observed fixed effects and random litter effects. Unlike mixed models, the Cox model fits the frailty term with penalized regression, the fitted object is of class coxph-penal and there is no method for survreg.coxph-penal, so the attempts to create the linear predictor fail. There are a couple workarounds.

1. Just fit the marginal model with centered covariates.

2. Center the covariates, fit 1, then fit the random effects model using coxme and extract the random effects, add them to the linear predictor with an offset to calculate the stratum specific survival curve for each litter.

3. Perform 2 and marginalize them by averaging all survival curves together, a separate approach to fitting the marginal model.

4. Use fixed effects or strata in a marginal Cox model to predict different survival curves for each litter.