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
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.