title: Is there a way to obtain event probability from a coxme model?
body: <p>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: </p>
  
  <pre><code># Generate data
  set.seed(123)
  mydata &lt;- 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 &lt;- coxme(Surv(Time , Surv) ~ Treatment*Origin + 
                     (1|Site/Block), data = mydata)
  </code></pre>
  
  <p>I want to obtain the predicted survival likelihood at Time = 3, for each Treatment:Origin combination. If I had a coxph model (<em>i.e.</em> without random effects), this could easily be done with survfit from the survival package: </p>
  
  <pre><code># use expand.grid to get a table with all possible combinations of Site and Treatment
  newdata.surv &lt;- with(mydata, expand.grid(Site = unique(Origin), Treatment = unique(Treatment)))
  
  # run survfit to predict the new values
  fitted &lt;- survival::survfit(mymodel, newdata = newdata.surv)
  
  # extract the fitted values for the time slice of interest: 3
  newdata.surv$fit &lt;- fitted$surv[3,]
  newdata.surv$lower &lt;- fitted$lower[3,] # Lower confidence interval
  newdata.surv$upper &lt;- fitted$upper[3,] # Upper confidence interval
  </code></pre>
  
  <p>However, <code>survfit</code> doesn't work with <code>coxme object</code>. I am aware that a <code>predict.coxme</code> exists in the <code>coxme package</code> 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 <a href="https://cran.r-project.org/web/packages/coxme/news.html" rel="nofollow noreferrer">https://cran.r-project.org/web/packages/coxme/news.html</a>). </p>
  
  <p>I have seen that <code>coxme</code> objects are supported by the packages <code>emmeans</code> and <code>lsmeans</code> but I am not sure that these packages could be use for predicting survival. Thanks in advance for your help. </p>