Reputation: 21
I fitted a mixed effects cox model using the package coxme. The model includes a two-level variable as fixed effect and a random intercept.
I want to predict survival probability and 95% confidence interval over time for each level of the categorical variable based on the fixed effects only (similar to re.form=NA in predict.merMod).
Constructing predictions is straightforward when there is no random effect using predict.coxph in the survival package, but predict.coxme does not provide standard deviations and new datasets are not yet supported.
How can I compute predicted survivor curve and 95%CI based on the coxme model (with a new dataset and fixed effects only)?
Sample code using lung dataset in the survival package:
library(survival)
library(coxme)
# format data
df <- na.omit(lung[, c("inst", "time", "status","sex")])
df$sex <- factor(df$sex)
# fit model
cme <- coxme(Surv(time, status) ~ sex + (1|inst), data = df)
summary(cme)
# build dummy dataset
pred_data = data.frame(sex=unique(lung$sex))
# get predicted values and se (does not work)
pred <- predict(cme, newdata, type= "lp", se.fit =T)
Related questions and answers
The code here seems close to what is needed, but I was not able to adapt it to get predicted values and 95%CI for a new dataset and using only fixed effects: R coxme: How to get study-specific treatment effect and 95% confidence interval from a mixed effect model?
This question is also similar: Plotting adjusted survival curve for a mixed-effects cox regression and/or time interaction? but there are convergence issues when I fit my original data set with a frailty coxph model and I would like to know if there is a way to get predicted survival probability and confidence bands for a coxme model.
Upvotes: 2
Views: 347
Reputation: 72828
You could try ehahelper::predict_coxme
.
> # remotes::install_github('junkka/ehahelper') ## install package
> ehahelper::predict_coxme(cme, pred_data, se.fit=TRUE)
$fit
1 2
-0.3191101 -0.8478545
$se.fit
[1] 0.1010445 0.2684685
Upvotes: 1