3

I would like to obtain prediction of an event probability from a coxme model, for a given values of my explanatory variables, and at a given timestep. My data are structured like this: 

# Generate data
set.seed(123)
mydata <- data.frame(Site = as.factor(sample(c("SiteA", "SiteB", "SiteC"), 1000, replace = TRUE)), 
                     Block = as.factor(sample(c("A", "B", "C", "D", "E", "F"), 1000, replace = TRUE)),
                     Treatment = as.factor(sample(c("Treat.A", "Treat.B"), 1000, replace = TRUE)),
                     Origin = as.factor(sample(c("Origin.A", "Origin.B"), 1000, replace = TRUE)),
                     Time = sample(seq(3), 1000, replace = TRUE), 
                     Surv = sample(c(0, 1), 1000, replace = TRUE)) # Alive is 0, death is 1

# Coxme model
mymodel <- coxme(Surv(Time , Surv) ~ Treatment*Origin + 
                   (1|Site/Block), data = mydata)

I want to obtain the predicted survival likelihood at Time = 3, for each Treatment:Origin combination. If I had a coxph model (i.e. without random effects), this could easily be done with survfit from the survival package: 

# use expand.grid to get a table with all possible combinations of Site and Treatment
newdata.surv <- with(mydata, expand.grid(Site = unique(Origin), Treatment = unique(Treatment)))

# run survfit to predict the new values
fitted <- survival::survfit(mymodel, newdata = newdata.surv)

# extract the fitted values for the time slice of interest: 3
newdata.surv$fit <- fitted$surv[3,]
newdata.surv$lower <- fitted$lower[3,] # Lower confidence interval
newdata.surv$upper <- fitted$upper[3,] # Upper confidence interval

However, survfit doesn't work with coxme object. I am aware that a predict.coxme exists in the coxme package but when I try to use it, I always get the error: "could not find function "predict.coxme". I am using the 2.2-10 version of the coxme package so the predict.coxme function should be included (see https://cran.r-project.org/web/packages/coxme/news.html).

I have seen that coxme objects are supported by the packages emmeans and lsmeans but I am not sure that these packages could be use for predicting survival. Thanks in advance for your help.

Julien Barrère
  • 51
  • 4

1 Answers1

3

The coxme.predict function is not exported per se, but you can invoke it using predict(mymodel) which will then invoke this method (or you could directly call coxme:::predict.coxme(mymodel) (with 3 colons)). See ?coxme:::predict.coxme for a brief description. It appears it currently does not support a newdata argument, so I am not sure how useful it is for your use case.

user12728748
  • 8,106
  • 2
  • 9
  • 14
  • Indeed, I was only using two colons that's why I coundn't make the function work (```coxme::predict.coxme```) so now it works. However, the only type of prediction that the function function return are either the linear predictor, either the absolute risk, and none of these variables can be used to compute survival probability (as far as I know). Any other ideas ? – Julien Barrère Mar 01 '20 at 15:56
  • See https://stats.stackexchange.com/questions/45538/how-to-generate-predicted-survivor-curves-from-frailty-models-using-r-coxph for a discussion and potential workaround. It seems nothing has changed since then. – user12728748 Mar 01 '20 at 22:05