对于我的解释变量的给定值和给定时间步,我想从 coxme 模型中获得对事件概率的预测。我的数据结构如下:
# 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)
对于每个处理:起源组合,我想在时间 = 3 处获得预测的生存可能性。如果我有一个 coxph 模型(即没有随机效应),这可以通过生存包中的 survfit 轻松完成:
# 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
但是,survfit
不适用于coxme object
. 我知道 apredict.coxme
中存在coxme package
但是当我尝试使用它时,我总是得到错误:“找不到函数”predict.coxme”。我使用的是 2.2-10 版本的 coxme 包,所以 predict.coxme应包含功能(请参阅https://cran.r-project.org/web/packages/coxme/news.html)。
我已经看到coxme
包支持对象emmeans
,lsmeans
但我不确定这些包是否可以用于预测生存。在此先感谢您的帮助。