In drizopoulos/JMbayes2: Extended Joint Models for Longitudinal and Time-to-Event Data

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library("JMbayes2")

Competing Risks

Prepare data

The first step to fit a joint model for competing events in JMbayes2 is to prepare the data for the event process. If there are $K$ competing events, then each subject needs to have $K$ rows, one for each possible cause. The observed event time $T_i$ of each subject is repeated $K$ times, and there are two indicator variables, namely one identifying the cause, and one indicating whether the corresponding event type is the one that occurred. Standard survival datasets that included a single row per patient, can be easily transformed to the competing risks long format using function crisk_setup(). This function accepts as main arguments the survival data in the standard format that has a single row per patient, the name of the status variable, and the level in this status variable that corresponds to censoring. We illustrate the use of this function in the PBC data, in which we treat as competing risks transplantation and death:

pbc2.id[pbc2.id$id %in% c(1, 2, 5), c("id", "years", "status")]

pbc2.idCR <- crisk_setup(pbc2.id, statusVar = "status", censLevel = "alive", 
                         nameStrata = "CR")

pbc2.idCR[pbc2.idCR$id %in% c(1, 2, 5), 
          c("id", "years", "status", "status2", "CR")]

Note that each patient is now represented by two rows (we have two possible causes of discontinuation from the study, death and transplantation), the event time variable years is identical in both rows of each patient, variable CR denotes the cause for the specific line of the long dataset, and variable status2 equals 1 if the corresponding event occurred.

Fit models

For the event process, we specify cause-specific relative risks models. Using dataset pbc2.idCR, we fit the corresponding cause-specific Cox regressions by including the interaction terms of age and treatment with variable CR, which is treated as a stratification variable using the strata() function:

CoxFit_CR <- coxph(Surv(years, status2) ~ (age + drug) * strata(CR),
                     data = pbc2.idCR)

For the longitudinal process, we include two longitudinal outcomes, namely serum bilirubin and the prothrombin time. For the former we use quadratic orthogonal polynomials in the fixed- and random-effects parts, and for the latter linear evolutions:

fm1 <- lme(log(serBilir) ~ poly(year, 2) * drug, data = pbc2, 
           random = ~ poly(year, 2) | id)
fm2 <- lme(prothrombin ~ year * drug, data = pbc2, random = ~ year | id)

To specify that each longitudinal outcome has a separate association coefficient per competing risk, we define the corresponding functional forms:

CR_forms <- list(
    "log(serBilir)" = ~ value(log(serBilir)):CR,
    "prothrombin" = ~ value(prothrombin):CR
)

Finally, the competing risks joint model is fitted with the following call to jm() (due to the complexity of the model, we have increased the number of MCMC iterations and the burn-in period per chain):

jFit_CR <- jm(CoxFit_CR, list(fm1, fm2), time_var = "year", 
              functional_forms = CR_forms, 
              n_iter = 25000L, n_burnin = 5000L, n_thin = 5L)

summary(jFit_CR)

Dynamic predictions

Based on the fitted competing risks joint model, we will illustrate how (dynamic) predictions for the cause-specific cumulative risk probabilities can be calculated. As an example, we will show these calculations for Patient 81 from the PBC dataset. First, we extract the data of this subject.

ND_long <- pbc2[pbc2$id == 81, ]
ND_event <- pbc2.idCR[pbc2.idCR$id == 81, ]
ND_event$status2 <- 0
ND <- list(newdataL = ND_long, newdataE = ND_event)

The first line extracts the longitudinal measurements, and the second line extracts the event times per cause (i.e., death and transplantation). This particular patient died at 6.95 years, but to make the calculation of cause-specific cumulative risk more relevant, we presume that she did not have the event, and we set the event status variable status2 to zero. The last line combines the two datasets in a list. Note: this last step is a prerequisite from the predict() method for competing risks joint model. That is, the datasets provided in the arguments newdata and newdata2 need to be named lists with two components. The first component needs to be named newdataL and contain the dataset with the longitudinal measurements, and the second component needs to be named newdataE and contain the dataset with the event information.

The predictions are calculated using the predict() method. The first call to this function calculates the prediction for the longitudinal outcomes at the times provided in the times argument, and the second call calculates the cause-specific cumulative risk probabilities. By setting the argument return_newdata to TRUE in both calls, we can use the corresponding plot() method to depict the predictions:

predLong <- predict(jFit_CR, newdata = ND, return_newdata = TRUE,
                    times = seq(6.5, 15, length = 25))

predEvent <- predict(jFit_CR, newdata = ND, return_newdata = TRUE,
                     process = "event")

plot(predLong, predEvent, outcomes = 1:2, ylim_long_outcome_range = FALSE,
     col_line_event = c("#03BF3D", "#FF0000"), 
     fill_CI_event = c("#03BF3D4D", "#FF00004D"), pos_ylab_long = c(1.5, 11.5))
legend(x = 8.1, y = 0.45, legend = levels(pbc2.idCR$CR), 
       lty = 1, lwd = 2, col = c("#03BF3D", "#FF0000"), bty = "n", cex = 0.8)

drizopoulos/JMbayes2 documentation built on April 20, 2024, 8:52 a.m.