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

Question

I've read that the best way to report the effect of spline terms on an outcome is to plot the predicted probabilities. I've been able to do this using the frailty coxph model in the survival package, by creating a new dataset with a dummy individual where the only differing variable is the spline variable, and then plotting the predictions for that dataset.

The code I've written looks like this:

### 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)

However, in reality, my model has two levels of random effects, which is only possible to model using the coxme package. Although there is a predict.coxme function, it doesn't work on a new dataset:

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")

This results in the following error: Error in predict.coxme(fit2, newdata, se.fit = TRUE, type = "risk") : newdata argument not yet supported

The package author has stated that he would like to implement this but it is unlikely to occur any time soon. So my previous method won't work, and I can't seem to find any packages that will support this.

I'm left wondering if there's a way of visualising the spline effects from a coxme model. Does anyone have any ideas or recommendations?

If not then... well, adding the levels does not significantly change the coefficients of the equivalent coxph model, so I'm tempted to plot the visualisations using predictions from the coxph model and report the more accurate coefficients in an accompanying table. Does that sound acceptable or am I heading into dangerous waters?

Just viewing this it appears you gave an example of a successful effort but you did not provide an example of a situation where you had difficulty. Am I interpreting this correctly? I ask because I had an idea of something to try but then when I "reached" for the MCVE, it appeared to have not been constructed. — IRTFM, Aug 16 '21 at 20:36
@IRTFM apologies, it's my first stackoverflow post so I might have failed to pose the question correctly. I'm not 100% certain what you mean by reaching for an MCVE, but I've updated my post assuming you mean I hadn't provided the code I would like to work and the accompanying error...? — Ee Cheng Ooi, Aug 17 '21 at 23:58
I am not familiar with this kind of model, but there must be a way to "manually" extract the parameters from the fitted model object so that you can calculate it yourself. You might have to read the source for `predict.coxme` to see how it works internally, and write the `newdata` functionality yourself. — shadowtalker, Aug 17 '21 at 23:59
@IRTFM I believe the MCVE is to replace the `fit <-` line with the `fit2 <-` line, and to replace `predict(fit, ...)` with `predict(fit2, ...)` acordingly. — shadowtalker, Aug 18 '21 at 00:00
Based on [the vignette](https://cran.r-project.org/web/packages/coxme/vignettes/laplace.pdf), it looks like you need to extract variables that correspond to: the baseline hazard function `lambda_0`, the fixed effect coefficients `beta`, and the variance of the random effect coefficients `theta`. If you can't figure it out from the source code, you might want to ask a separate question like "how do I make predictions from a fitted Multilevel Cox model?" on [Cross Validated](https://stats.stackexchange.com/). — shadowtalker, Aug 18 '21 at 00:12
@shadowtalker Oh, wow! Thanks very much. I think that's given me a few ideas... if I work it out, I'll post the code here. — Ee Cheng Ooi, Aug 18 '21 at 00:39

Ee Cheng Ooi · Accepted Answer · 2021-09-01T05:32:47.903

For anyone who has come to this, the ehahelper package provides a predict_coxme function (as well as some other useful tools for working with coxme objects). It's definitely worth reading the documentation though, to see how the developer deals with the random effects and calculates the denominator for the hazard ratio.

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)

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

1 Answers1