Survfit equivalent for coxme in R?

Question

I am struggling with the coxme package in R. I would like to use a function like survfit() - the way it would ordinarily be used for a coxph() model - to plot adjusted survival curves and find the median survival at different parameter values.

If I fit the model using coxph without random effects I can do the following:

library(KMsurv)
data(burn)

my.surv <- with(burn, Surv(T1, D1))

cox_nr = coxph(my.surv ~ Z1 , data = burn)

survfit(cox_nr, newdata = data.frame(Z1 =1))

This provides survival estimates. But if I fit the same model with coxme:

library(coxme)
cox_r = coxme(my.surv ~ Z1 + (1|Z11), data = burn)

survfit(cox_r, newdata = data.frame(Z1 = 1))

Error in UseMethod("survfit", formula) : no applicable method for 'survfit' applied to an object of class "coxme"

So survfit.coxme doesn't seem to exist and from reading the coxme package documentation, I don't see an equivalent. Is there something fundamentally wrong about what I am attempting to do? If not, how can I get these estimates?

Answer 1

I think the reason why there's no survfit methods for coxme is because of the frailty model. The log-rank or wilcoxon tests rely on the one-to-one correspondence between failure/censoring outcomes to individuals in the risk sets. This allows you to consistently estimate their survival curves using the non-parametric kaplan meier curves, which are monotonic and non-increasing always. That's not the case if an individual can have more than one outcome, which is what the coxme(frailty) is handling. In the case of herpes outbreaks, as an example, if individuals can have more than one outbreak, or if you can have any number of outbreaks in a cluster, then you can't estimate the survival curve with a KMcurve and you can't perform the log-rank test.

However, the inference on the Cox model using the summary command is asymptotically equivalent to the log-rank test for basic univariate linear Cox models. You can argue that taking a summary of the frailty model serves as a stratified equivalent test handling multiple endpoints and that the p-value represents a scientifically interesting component. For a graphical way of depicting failures within clusters, consider using cumulative incidence curves instead.