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?