Is there a way to obtain event probability from a coxme model?

user_11102244 0
Viewed 601

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.

lsmeans survival r
1 Answers

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.

user_127287480
Related