Plotting adjusted survival curve for a mixed-effects cox regression and/or time interaction?

Question

I have created a mixed effects cox regression using coxme.

Q1) I would like to plot the coefficients of the fixed effects in an adjusted survival curve. However, it seems this functionality in packages like survminer is only possible for coxph objects without a frailty term.

Is there anyway that this could be calculated and plotted manually in R instead?

I have the below model for demonstration purposes:

library(survival)
library(coxme)
kidney <-  data.frame(kidney)

coxme <- coxme(Surv(time, status) ~ age + sex + (1|disease), 
                  data = kidney)

Q2) Additionally, is there anyway that this could be plotted for a time interaction (tt() in coxph)?

For instance with the model:

library(survival)
kidney <-  data.frame(kidney)

coxph.me <- coxph(Surv(time, status) ~ age + sex + tt(sex) + frailty(disease), 
                  data = kidney,
                  tt=function(x,t,...) x*t)

Thanks in advance

Answer 1

The most popular way to obtain adjusted survival curves from a model is g-computation, which requires a function to make predictions for survival probabilities at t given a vector of covariates and t. Unfortunately, the coxme function does not naturally support such predictions, e.g. there is no predict.coxme function included in the package.

You can, however, use the adjustedCurves package in conjunction with the riskRegression package to get what you want if you switch from coxme to a standard coxph model with a frailty() term. Below I give an example on how this could be done. It is a little hacky because of a bug in riskRegression but it works just fine.

library(riskRegression)
library(adjustedCurves)
library(survival)

set.seed(42)

# simulate some example data
sim_dat <- sim_confounded_surv(n=500, max_t=1.2)
sim_dat$group <- as.factor(sim_dat$group)
sim_dat$cluster <- sample(1:10, size=nrow(sim_dat), replace=TRUE)

# outcome model
# NOTE: your frailty term needs to be named "cluster" due to some
#       sub-optimal coding in the riskRegression package
cox_mod <- coxph(Surv(time, event) ~ x1 + x2 + x4 + x5 + group + 
                   frailty(cluster),
                 data=sim_dat, x=TRUE)

predict_fun <- function(...) {
  1 - predictRisk(...)
}

# using direct adjustment 
adjsurv <- adjustedsurv(data=sim_dat,
                        variable="group",
                        ev_time="time",
                        event="event",
                        method="direct",
                        outcome_model=cox_mod,
                        predict_fun=predict_fun)

plot(adjsurv)

If you need confidence intervals as well, you can obtain those using bootstrapping by changing the code a little bit (note that this may take a while, especially if you have lots of data).

# using direct adjustment with bootstrap confidence intervals
adjsurv <- adjustedsurv(data=sim_dat,
                        variable="group",
                        ev_time="time",
                        event="event",
                        method="direct",
                        bootstrap=TRUE,
                        n_boot=500,
                        outcome_model=cox_mod,
                        predict_fun=predict_fun)

plot(adjsurv, conf_int=TRUE, use_boot=TRUE)

None of this works with a tt() term in the model formula though.