# Retrieving the variance partition coefficient of the frailty() term (mixed effects) for cox regression in R

I am new to running survival analyses. I am attempting to construct a cox regression using the coxph() function. I want to include a frailty() term for mixed effects, but I also need to include a time interaction (tt()) for one of my categorical predictors, as this variable breaks the proportional hazards assumption. Therefore, I have opted to use the coxph() function from the survival package instead of coxme, as it appears coxme does not currently have the functionality to include a time interaction. The frailty term is to simulate a random intercept equivalent for in cox regressions, for instance (1|disease).

Below is an example model, for the sake of demonstration:

library(survival)
kidney <-  data.frame(kidney)

coxph.me <- coxph(Surv(time, status) ~ age + sex + tt(sex) + frailty(disease),
data = kidney,
tt=function(x,t,...) x*t)


I would appreciate any answers for the two questions below:

1. Now that I have included the frailty term, I want to be able to extract the variance partition coefficients or intra class correlations of the random effects, so that I can determine how much variance the random effect explains in the outcome.

2. I am aware that in cox regressions there is no such thing as an intercept and that the notation (1|disease) or frailty(disease) reflects the baseline hazard function. Would I be able to produce caterpillar plot of the random baseline hazard function/ hazard ratios as would be possible with the random intercepts in lme4? Or am I conceptualising this incorrectly?

See T.A. Balan and H. Putter, A Tutorial on frailty models, Statistical Methods in Medical Research 29: 3424–3454 (2020) for an overview of different types of heterogeneity in survival models, including those due to frailty.

First, be sure that the apparent proportional hazards (PH) violation isn't due to an omitted covariate. In particular, Section 2.1.2 of the above reference shows how omitting one outcome-associated predictor from a Cox model can lead to an apparent PH violation by an included predictor even though PH would be met for both in a model that included both. The Cox regression coefficient estimate for the included predictor then is a type of time-averaged value that is biased in magnitude toward 0. Newcomers to survival analysis are often surprised that, unlike with omitted-variable bias in linear regression, this happens even if the omitted covariate is uncorrelated with the included covariate.

Second, if the PH violation is relatively small you might choose to ignore it, accept the coefficient estimate as a type of time-averaged value, and show a plot of scaled Schoenfeld residuals over time to document the magnitude of the apparent violation.

Third, make sure that you really want a conditional survival model with a random-effect/frailty term. A marginal model using a cluster() term can account for simple intra-group correlations. That provides a robust coefficient-covariance matrix estimate similar to sandwich estimates in linear regression. That robust estimate then allows for inference on the (time-averaged) coefficient estimates even under PH violation. See this page. That doesn't, of course, directly model the variance among groups.

Fourth, if you do want a random-effect/frailty model, ask yourself whether things like intra-class correlation are really what you care about. Survival models typically have a very large inherent variance in survival times. My sense is that while readers might well be interested in how much including a frailty term improves a Cox model, they might not care so much about the fraction of variation that's explained by it.

Fifth, if you do want an estimate of variation explained by the frailty term, recall that a Cox model doesn't directly estimate event times. It only models the rank-ordering of event times. Thus calculations for things like "variation explained" typically start with a monotone transformation of event times $$g(T)$$ that effectively converts the semi-parametric Cox model to a Weibull or exponential form. G. Honerkamp-Smitha and R. Xu, Three measures of explained variation for correlated survival data under the proportional hazards mixed-effects model,* Statist. Med. 35: 4153–4165 (2016), write this in a form that incorporates both fixed effects $$Z$$ and random effects $$W$$ and corresponding coefficient vectors $$\beta$$ and $$b$$:

$$g(T_{ij})= -(\beta' Z_{ij} + b'W_{ij}) + \epsilon_{ij}$$

for the $$j^{th}$$observation in cluster $$i$$. In that form, the $$\epsilon_{ij}$$ have a minimum extreme value distribution with variance $$\pi^2/6$$ . With only a frailty term for the baseline hazard, you then could simply apply a standard intraclass correlation coefficient formula with $$\pi^2/6$$ for the "unobserved noise."

Sixth, there are arguments that the above formulation for variation explained isn't the best approach for Cox models with censored event times. J.T. Kent and J. O'Quigley, in Measures of Dependence for Censored Data, Biometrika 75: 525-34 (1988), start with that $$\pi^2/6$$ variance but abandon it for a "stronger notion of dependence... defined using the concept of information gain." P. Royston, in Explained Variation for Survival Models, Stata Journal 6:83-96 (2006), discusses some more recent literature and presents an alternative.

Finally, a few technical points. Yes, you can extract the per-group frailty values from a coxph model. They are in the frail slot of the model, for example coxph.me\$frail in your example. Also, if coxme() can work with counting-process data (I'm not sure that it can) you might be able to use it instead of coxph(). To model a time-varying coefficient, the tt() function internally calculates a time-varying covariate value, based on the supplied function, for each individual at risk at each event time in the data set. You could generate a corresponding extended data set yourself, in the Surv(startTime, stopTime, event) counting-process format. See Section 5 of the time-dependence vignette. "The disadvantage of this is the very large data set, of course, but it is no larger than the internal one created by the tt call. An advantage is that further processing of the model is available, such as residuals, the anova function, or survival curves."

*The author manuscript is freely available from PubMed Central, but it includes some errors that were corrected in the published, copy-edited version to which I link. In particular, Equation 5 should be the formula for $$\Omega^2$$, not for $$\Omega$$ as in the author manuscript.

• Thanks EdM for such a comprehensive answer. This was a huge help! Commented Dec 13, 2023 at 15:03