我如何从样条项与一个comme模型在R中得到预测?

Question

我读到过，报告样条项对结果的影响的最好方法是绘制预测概率。我已经能够使用survival包中的脆弱性coxph模型做到这一点，通过创建一个具有虚拟个体的新数据集，其中唯一不同的变量是样条变量，然后绘制该数据集的预测。

我写的代码是这样的:

### Example
library(tidyverse); library(survival)
# Create dataset
dat <- tibble(time = sample(seq(21), 1000, replace = TRUE),
status = sample(c(0,1), 1000, replace = TRUE), 
var1 = sample(seq(100), 1000, replace = TRUE), 
var2 = sample(seq(100), 1000, replace = TRUE), 
group_id = seq(1, 1000, by = 1))
# Run model
fit <- coxph(Surv(time, status) ~ var1 + ns(var2, df = 3) + frailty(group_id), data = dat)
# Create dummy dataset
spline_values <- tibble(var1 = seq(1, 1000, by = 1))
newdata <- dat %>% 
head(n = 1) %>% 
select(time, status, var1, var2, group_id) %>%
slice(rep(1:n(), each = nrow(spline_values))
newdata <- newdata %>% mutate(var1 = spline_values$var1)
pred <- predict(fit, newdata = newdata, se.fit = TRUE, type = "risk")
newdata <- newdata %>% add_column(rr = pred$fit, se.fit = pred$se.fit)
newdata <- newdata %>% mutate(
low = exp(log(rr) - 1.96 * se.fit), 
high = exp(log(rr) + 1.96 * se.fit))
ggplot(newdata, aes(var1)) +
geom_line(aes(y = rr)) +
geom_ribbon(aes(ymin = low, ymax = high), alpha = 0.1, linetype = 0)

然而，在现实中，我的模型有两个层次的随机效应，这只能使用coxme包来建模。虽然有一个predict.coxme函数，它不工作在一个新的数据集:

library(coxme)
# Run coxme model
fit2 <- coxme(Surv(time, status) ~ var1 + ns(var2, df = 3) + (1 | group_id), data = dat)
pred <- predict(fit2, newdata = newdata, se.fit = TRUE, type = "risk")

这会导致以下错误:Error in predict.coxme(fit2, newdata, se.fit = TRUE, type = "risk") : newdata argument not yet supported

包的作者已经声明他想要实现这个，但是它不太可能很快发生。所以我之前的方法不起作用，而且我似乎找不到任何支持这个的包。

我想知道是否有一种方法可以从comme模型中可视化样条效果。有人有什么想法或建议吗?

如果不是，那么…嗯，添加级别不会显著改变等效coxph模型的系数，所以我想使用coxph模型的预测来绘制可视化图，并在随附的表中报告更准确的系数。这听起来可以接受吗，还是我进入了危险的水域?

Answer 1

对于了解这一点的人来说，ehahelper包提供了一个predict_comme函数(以及其他一些用于处理comme对象的有用工具)。当然，还是值得一读文档，看看开发者是如何处理随机效应和计算风险比的分母的。

devtools::install_github('junkka/ehahelper')
library(ehahelper); library(coxme); library(tidyverse)
# Create dataset
dat <- tibble(time = sample(seq(21), 1000, replace = TRUE),
status = sample(c(0,1), 1000, replace = TRUE), 
var1 = sample(seq(100), 1000, replace = TRUE), 
var2 = sample(seq(100), 1000, replace = TRUE), 
group_id = seq(1, 1000, by = 1))
# Run model
fit2 <- coxme(Surv(time, status) ~ var1 + ns(var2, df = 3) + (1 | group_id), data = dat)
# Create dummy dataset
spline_values <- tibble(var1 = seq(1, 1000, by = 1))
newdata <- dat %>% 
head(n = 1) %>%
select(time, status, var1, var2, group_id) %>% 
slice(rep(1:n(), each = nrow(spline_values)))

newdata <- newdata %>% mutate(var1 = spline_values$var1)
# Ehahelper predict_coxme has difficulty with splines in a new dataset. 
# To get around this, append the new dataset to the old one, apply the predict_coxme function, then subset the predicted values of interest
dat <- bind_rows(dat, newdata)
pred <- ehahelper::predict_coxme(fit2, newdata = dat, se.fit = TRUE, type = "risk", strata_ref = FALSE)

dat <- dat %>% add_column(rr = pred$fit[,1], se.fit = pred$se.fit[,1])
# Subset values of interest
dat <- dat %>% tail(nrow(newdata))

dat <- dat %>% mutate(
low = exp(log(rr) - 1.96 * se.fit), 
high = exp(log(rr) + 1.96 * se.fit))

ggplot(dat, aes(var1)) +
geom_line(aes(y = rr)) +
geom_ribbon(aes(ymin = low, ymax = high), alpha = 0.1, linetype = 0)