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")]
#>   id     years       status
#> 1  1  1.095170         dead
#> 2  2 14.152338        alive
#> 5  5  4.120578 transplanted

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")]
#>     id     years       status status2           CR
#> 1    1  1.095170         dead       1         dead
#> 1.1  1  1.095170         dead       0 transplanted
#> 2    2 14.152338        alive       0         dead
#> 2.1  2 14.152338        alive       0 transplanted
#> 5    5  4.120578 transplanted       0         dead
#> 5.1  5  4.120578 transplanted       1 transplanted

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)
#> 
#> Call:
#> jm(Surv_object = CoxFit_CR, Mixed_objects = list(fm1, fm2), time_var = "year", 
#>     functional_forms = CR_forms, n_iter = 25000L, n_burnin = 5000L, 
#>     n_thin = 5L)
#> 
#> Data Descriptives:
#> Number of Groups: 312        Number of events: 169 (27.1%)
#> Number of Observations:
#>   log(serBilir): 1945
#>   prothrombin: 1945
#> 
#>                  DIC     WAIC      LPML
#> marginal    10638.85 10843.25 -5736.199
#> conditional 15754.08 15437.73 -8218.645
#> 
#> Random-effects covariance matrix:
#>                                              
#>        StdDev    Corr                        
#> (Intr) 1.3331  (Intr)  p(,2)1  p(,2)2  (Intr)
#> p(,2)1 22.5073 0.6871                        
#> p(,2)2 12.4989 -0.2906 -0.1586               
#> (Intr) 0.7867  0.6315  0.4254  -0.3413       
#> year   0.3284  0.4328  0.3264  -0.0771 0.0307
#> 
#> Survival Outcome:
#>                                        Mean  StDev    2.5%   97.5%      P
#> age                                 -0.0480 0.0215 -0.0909 -0.0077 0.0177
#> drugD-penicil                       -0.3071 0.3958 -1.0971  0.4627 0.4378
#> age:strata(CR)dead                   0.1066 0.0205  0.0674  0.1455 0.0000
#> drugD-penicil:strata(CR)dead         0.2979 0.4289 -0.5265  1.1570 0.4948
#> value(log(serBilir)):CRtransplanted  1.0586 0.2124  0.6694  1.5045 0.0000
#> value(log(serBilir)):CRdead          1.4241 0.1096  1.2156  1.6460 0.0000
#> value(prothrombin):CRtransplanted    0.0792 0.1158 -0.1625  0.2915 0.4618
#> value(prothrombin):CRdead            0.1339 0.0477  0.0311  0.2199 0.0162
#>                                       Rhat
#> age                                 1.0876
#> drugD-penicil                       1.0319
#> age:strata(CR)dead                  1.1269
#> drugD-penicil:strata(CR)dead        1.0333
#> value(log(serBilir)):CRtransplanted 1.0353
#> value(log(serBilir)):CRdead         1.0023
#> value(prothrombin):CRtransplanted   1.3023
#> value(prothrombin):CRdead           1.0112
#> 
#> Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
#>                   Mean  StDev    2.5%   97.5%      P   Rhat
#> (Intercept)     1.1541 0.0978  0.9626  1.3484 0.0000 0.9999
#> poly(year, 2)1 26.0863 1.9731 22.2498 30.0239 0.0000 1.0014
#> poly(year, 2)2  0.8121 1.3953 -1.9108  3.5203 0.5580 1.0000
#> drugD-penicil  -0.1552 0.1346 -0.4175  0.1103 0.2462 1.0005
#> p(,2)1         -1.9547 2.2371 -6.2728  2.4660 0.3787 1.0004
#> p(,2)2         -0.9302 1.7152 -4.2955  2.4247 0.5847 1.0011
#> sigma           0.3022 0.0061  0.2905  0.3144 0.0000 1.0008
#> 
#> Longitudinal Outcome: prothrombin (family = gaussian, link = identity)
#>                       Mean  StDev    2.5%   97.5%      P   Rhat
#> (Intercept)        10.6228 0.0803 10.4671 10.7797 0.0000 1.0018
#> year                0.2904 0.0382  0.2157  0.3663 0.0000 1.0030
#> drugD-penicil      -0.0880 0.1121 -0.3080  0.1261 0.4367 1.0001
#> year:drugD-penicil -0.0192 0.0503 -0.1172  0.0803 0.7062 1.0004
#> sigma               1.0547 0.0204  1.0153  1.0958 0.0000 1.0011
#> 
#> MCMC summary:
#> chains: 3 
#> iterations per chain: 25000 
#> burn-in per chain: 5000 
#> thinning: 5 
#> time: 6.8 min

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)