knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library("JMbayes2")
It is often the case that a subject may transition between multiple states and we are interested to assess the association of longitudinal marker(s) with each of these transitions. In this vignette we will illustrate how to achieve this using JMbayes2.
We will consider a simple case with one longitudinal outcome and a three-state (illness-death) model, but this application can be extended for the cases of multiple longitudinal markers and more than three states.
First we will simulate data from a joint model with a single linear mixed effects model and a multi-state process with three possible states. The multi-state process can be visualized as:
library("ggplot2") d <- data.frame("xmin" = c(0, 45, 22.5), "xmax" = c(15, 60, 37.5), "ymin" = c(50, 50, 0), "ymax" = c(60, 60, 15)) dline <- data.frame("x1" = c(15), "y1" = c(55), "xend" = c(45), "yend" = c(55)) dcurve <- data.frame("x1" = c(7.5, 52.5), "y1" = c(50, 50), "xend" = c(22.4, 37.4), "yend" = c(7.5, 7.5), "curvature" = c(1, -1)) ggplot() + geom_rect(data = d, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax), fill = "#ffffff", color = 'black', size = 1) + geom_text(aes(x = 7.5, y = 55, label = "Healthy"), size = 6.5) + geom_text(aes(x = 52.5, y = 55, label = "Illness"), size = 6.5) + geom_text(aes(x = 30, y = 7, label = "Death"), size = 6.5) + geom_text(aes(x = 30, y = 52.5, label = 'h[1][2](t)'), size = 4, parse = TRUE) + geom_text(aes(x = 12.5, y = 30, label = "h[1][3](t)"), size = 4, parse = TRUE) + geom_text(aes(x = 47.5, y = 30, label = "h[2][3](t)"), size = 4, parse = TRUE) + geom_segment(data = dline, aes(x = x1, y = y1, xend = xend, yend = yend), size = 1, arrow = arrow(length = unit(0.02, "npc"))) + geom_curve(data = dcurve[1, ], aes(x = x1, y = y1, xend = xend, yend = yend), size = 1, curvature = 0.3, arrow = arrow(length = unit(0.02, "npc"))) + geom_curve(data = dcurve[2, ], aes(x = x1, y = y1, xend = xend, yend = yend), size = 1, curvature = -0.3, arrow = arrow(length = unit(0.02, "npc"))) + ylim(0, 60) + xlim(0, 60) + theme_void()
where all subjects start from state "Healthy" and then can transition to either state "Illness" and then state "Death" or directly to state "Death". In this case, states "Healthy" and "Illness" are transient states as the subject, when occupying these states, can still transition to other states whereas "Death" is an absorbing state as when a subject reaches this state then no further transitions can occur. This means that three transitions are possible: $1 \rightarrow 2$, $1 \rightarrow 3$ and $2 \rightarrow 3$ with transition intensities $h_{12}\left(t\right)$, $h_{13}\left(t\right)$ and $h_{23}\left(t\right)$ respectively.
For our example the default functional form is assumed, i.e., that the linear predictor $\eta(t)$ of the mixed model is associated with the each transition intensity at time $t$. The following piece of code simulates the data:
set.seed(1234) # number of subjects N <- 500 # number of measurements per subject n <- 20 # vector of ids id <- rep(1:N, each = n) # minimum and maximum follow-up times min.t <- 0.01 max.t <- 16 # sample time-points time <- replicate(N, c(0, sort(runif(n - 1, min = min.t, max = max.t))), simplify = FALSE) time <- do.call(c, time) # sample continuous covariate values Xcov.s <- rnorm(N, mean = 4.763, sd = 2.8) # wide version Xcov <- rep(Xcov.s, each = n) # long # initiate data frame to store results DF <- data.frame("id" = id, "time" = time, "X" = Xcov) # design matrices for fixed and random effects X <- model.matrix(~ 1 + time + X, data = DF) Z <- model.matrix(~ 1 + time, data = DF) D11 <- 1.0 # variance of random intercepts D22 <- 0.5 # variance of random slopes # we simulate random effects b <- cbind(rnorm(N, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22))) # fixed effects coefficients true.betas <- c(-0.482, 0.243, 1.52) # linear predictor eta.y <- as.vector(X %*% true.betas + rowSums(Z * b[id, ])) # residual standard error sigma.e <- 1.242 # sample longitudinal outcome values DF$y <- rnorm(N * n, eta.y, sigma.e) # values for the association parameter per transition alpha <- c("alpha.01" = 0.8, "alpha.02" = 0.55, "alpha.12" = 1.25) # shape of Weibull for each transition phi <- c("phi.01" = 12.325, "phi.02" = 8.216, "phi.12" = 3.243) # regression coefficients for transition intensities gammas <- c("(Intercept)1" = -22.25, "X" = 1.2, "(Intercept)2" = -18.25, "X" = 1.2, "(Intercept)3" = -19.25, "X" = 1.2) # design matrix transition intensities W <- cbind("(Intercept)1"= rep(1, N), Xcov[seq(1, by = n, N*n)], "(Intercept)2"= rep(1, N), Xcov[seq(1, by = n, N*n)], "(Intercept)3"= rep(1, N), Xcov[seq(1, by = n, N*n)]) ## linear predictor for transition: 0 -> 1 eta.t1 <- as.vector(W[, c(1, 2), drop = FALSE] %*% gammas[1:2]) ## linear predictor for transition: 0 -> 2 eta.t2 <- as.vector(W[, c(3, 4), drop = FALSE] %*% gammas[3:4]) ## linear predictor for transition: 1 -> 2 eta.t3 <- as.vector(W[, c(5, 6), drop = FALSE] %*% gammas[5:6]) # we simulate event times using inverse transform sampling invS01 <- function(t, u, i) { h <- function(s) { XX <- cbind(1, s, Xcov[i]) ZZ <- cbind(1, s) f1 <- as.vector(XX %*% true.betas + rowSums(ZZ * b[rep(i, nrow(ZZ)), ])) exp(log(phi["phi.01"]) + (phi["phi.01"] - 1)*log(s) + eta.t1[i] + f1*alpha["alpha.01"]) } integrate(h, lower = 0, upper = t, subdivisions = 10000L)$value + log(u) } invS02 <- function(t, u, i) { h <- function(s) { XX <- cbind(1, s, Xcov[i]) ZZ <- cbind(1, s) f1 <- as.vector(XX %*% true.betas + rowSums(ZZ * b[rep(i, nrow(ZZ)), ])) exp(log(phi["phi.02"]) + (phi["phi.02"] - 1)*log(s) + eta.t2[i] + f1*alpha["alpha.02"]) } integrate(h, lower = 0, upper = t, subdivisions = 10000)$value + log(u) } invS12 <- function (t, u, i) { h <- function (s) { XX <- cbind(1, s, Xcov[i]) ZZ <- cbind(1, s) f1 <- as.vector(XX %*% true.betas + rowSums(ZZ * b[rep(i, nrow(ZZ)), ])) exp(log(phi["phi.12"]) + (phi["phi.12"] - 1) * log(s) + eta.t3[i] + f1 * alpha["alpha.12"]) } integrate(h, lower = 0, upper = t, subdivisions = 10000)$value + log(u) } # Probability for each transition u01 <- runif(N, 0, 1) u02 <- runif(N, 0, 1) u12 <- runif(N, 0, 1) # initiate vectors to save true event times trueT01 <- numeric(N) trueT02 <- numeric(N) trueT12 <- numeric(N) # sample censoring times mean.Cens <- 9 Ctimes <- runif(N, 0, 2 * mean.Cens) # simulate time-to-event data for (i in 1:N) { Root01 <- NULL Root02 <- NULL Root12 <- NULL Up <- 50 tries <- 5 # Transition 0->1 Up <- 200 Root01 <- try(uniroot(invS01, interval = c(1e-05, Up), u = u01[i], i = i)$root, TRUE) trueT01[i] <- if (!inherits(Root01, "try-error")) Root01 else 500 # Transition 0->2 Up <- 200 Root02 <- try(uniroot(invS02, interval = c(1e-05, Up), u = u02[i], i = i)$root, TRUE) trueT02[i] <- if (!inherits(Root02, "try-error")) Root02 else 500 # Transition 1->2 if(as.numeric(trueT01[i]) < as.numeric(trueT02[i]) && as.numeric(trueT01[i]) < Ctimes[i]) { Up <- Up + 200 Root12 <- try(uniroot(invS12, interval = c(as.numeric(trueT01[i]), Up), u = u12[i], i = i)$root, TRUE) } else {Root12 <- 500} trueT12[i] <- if (!inherits(Root12, "try-error")) Root12 else 500 } # initiate multi-state dataset in wide format data_mstate <- data.frame('id' = 1:N, 'trueT01' = trueT01, 'trueT02' = trueT02, 'trueT12' = trueT12, 'Ctimes' = Ctimes, 'X' = Xcov.s) # split by id data_mstate_split.by.id <- split(data_mstate, data_mstate$id) # function to pass to lapply to prepare multi-state data per id ms_arrange <- function (x) { if (x$Ctimes < min(x$trueT01, x$trueT02)) { x_new <- data.frame('id' = rep(x$id, 2), 'from_state' = rep(1, 2), 'to_state' = 2:3, 'transition' = 1:2, 'Tstart' = rep(0, 2), 'Tstop' = x$Ctimes, 'status' = rep(0, 2), 'X' = x$X) } else { if (x$trueT02 < x$trueT01) { x_new <- data.frame('id' = rep(x$id, 2), 'from_state' = rep(1, 2), 'to_state' = 2:3, 'transition' = 1:2, 'Tstart' = rep(0, 2), 'Tstop' = x$trueT02, 'status' = c(0, 1), 'X' = x$X) } else { if (x$Ctimes < x$trueT12) { x_new <- data.frame('id' = rep(x$id, 3), 'from_state' = c(1, 1, 2), 'to_state' = c(2:3, 3), 'transition' = 1:3, 'Tstart' = c(rep(0, 2), x$trueT01), 'Tstop' = c(rep(x$trueT01, 2), x$Ctimes), 'status' = c(1, 0, 0), 'X' = x$X) } else { x_new <- data.frame('id' = rep(x$id, 3), 'from_state' = c(1, 1, 2), 'to_state' = c(2:3, 3), 'transition' = 1:3, 'Tstart' = c(rep(0, 2), x$trueT01), 'Tstop' = c(rep(x$trueT01, 2), x$trueT12), 'status' = c(1, 0, 1), 'X' = x$X) } } } } data_mstate_split.by.id <- lapply(data_mstate_split.by.id, ms_arrange) data_mstate <- do.call(rbind, data_mstate_split.by.id) data_mstate$transition <- factor(data_mstate$transition) Tstop <- tapply(data_mstate$Tstop, data_mstate$id, max) Tstop <- Tstop[id] DF <- DF[DF$time <= Tstop, ]
The data for the multi-state process need to be in the appropriate long format:
head(data_mstate, n = 5L)
for example subject 1 experienced the following transition: $1 \rightarrow 2$ and therefore is represented in 3 rows, one for each transition, because all of these transitions were plausible. On the other hand subject 2 is only represented by two rows, only for transitions $1 \rightarrow 2$ and $1 \rightarrow 3$ since these are the only transitions that are possible from state 1. That is, since subject 2 never actually transitioned to state 2, transition $2 \rightarrow 3$ was never possible and therefore no row for this transition is in the dataset. It is also important to note that the time in the dataset follows the counting process formulation with intervals specified by Tstart
and Tstop
and that there is a variable (in this case transition
) which indicates to which transition the row corresponds to.
As soon data in the appropriate format are available, fitting the model is very straightforward. First we fit a linear mixed model using the lme()
function from package nlme:
mixedmodel <- lme(y ~ time + X, random = ~ time | id, data = DF)
then we fit a multi-state model using function coxph()
from package survival making sure we use the counting process specification and that we add strata(transition)
to stratify by the transition indicator variable in the dataset. Furthermore we add an interaction between covariate X
and each transition to allow the effect of this covariate to vary across transitions.
msmodel <- coxph(Surv(Tstart, Tstop, status) ~ X * strata(transition), data = data_mstate)
Finally, to fit the joint model we simply run:
jm_ms_model <- jm(msmodel, mixedmodel, time_var = "time", functional_forms = ~ value(y):transition, n_iter = 10000L) summary(jm_ms_model)
which differs from a default call to jm()
by the addition of the functional_forms
argument by which we specify that we want an "interaction" between the value of the marker and each transition which translates into a separate association parameter for the longitudinal marker and each transition.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.