如何使用 R 中的 coxme 模型从样条项中获得预测？

Question

我读到报告样条项对结果影响的最佳方法是绘制预测概率图。我已经能够使用 survival 包中的 frailty 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

软件包作者已声明他想实现此功能，但不太可能很快实现。所以我以前的方法行不通，而且我似乎找不到任何支持这个的包。

我想知道是否有一种方法可以可视化 coxme 模型的样条效应。有没有人有什么想法或建议？

如果不是那么......好吧，添加水平不会显着改变等效 coxph 模型的系数，所以我很想使用 coxph 模型的预测来绘制可视化效果并报告更准确的系数伴随 table。这听起来是否可以接受table 还是我正在进入危险水域？

Answer 1

对于任何接触过这个的人来说，ehahelper 包提供了一个 predict_coxme 函数（以及其他一些用于处理 coxme 对象的有用工具）。不过，这绝对值得 reading the documentation，看看开发人员如何处理随机效应并计算风险比的分母。

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)

如何使用 R 中的 coxme 模型从样条项中获得预测？

How can I get predictions from spline terms with a coxme model in R?

r

spline

predict

survival