### Examples Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/coda_methods.html Examples of using coda methods for joint models. ```R # linear mixed model fits fit_lme1 <- lme(log(serBilir) ~ year:sex + age, random = ~ year | id, data = pbc2) fit_lme2 <- lme(prothrombin ~ sex, random = ~ year | id, data = pbc2) # cox model fit fit_cox <- coxph(Surv(years, status2) ~ age, data = pbc2.id) # joint model fit fit_jm <- jm(fit_cox, list(fit_lme1, fit_lme2), time_var = "year", n_chains = 1L) # trace plot for the fixed effects in the linear mixed submodels traceplot(fit_jm, parm = "betas") # density plot for the fixed effects in the linear mixed submodels densplot(fit_jm, parm = "betas") # cumulative quantile plot for the fixed effects in the linear mixed submodels cumuplot(fit_jm, parm = "betas") # trace plot for the fixed effects in the linear mixed submodels ggtraceplot(fit_jm, parm = "betas") ggtraceplot(fit_jm, parm = "betas", grid = TRUE) ggtraceplot(fit_jm, parm = "betas", custom_theme = c('1' = 'black')) # trace plot for the fixed effects in the linear mixed submodels ggdensityplot(fit_jm, parm = "betas") ggdensityplot(fit_jm, parm = "betas", grid = TRUE) ggdensityplot(fit_jm, parm = "betas", custom_theme = c('1' = 'black')) ``` -------------------------------- ### Example Usage Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/cr_setup.md An example demonstrating how to use the `crisk_setup` function with sample data and displaying the first few rows of the transformed data. ```R head(crisk_setup(pbc2.id, "status", "alive")) ``` -------------------------------- ### Example Survival Data Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/articles/Multi_State_Processes.md Displays the first 5 rows of the survival data, showing subject IDs, transitions, start and stop times, status, and a covariate X. ```r head(df_surv, n = 5L) #> id transition Tstart Tstop status X #> 1 1 1 0.000000 1.014016 1 0 #> 2 1 2 0.000000 1.014016 0 0 #> 3 1 3 1.014016 1.971082 1 0 #> 4 2 1 0.000000 7.231688 0 1 #> 5 2 2 0.000000 7.231688 0 1 ``` -------------------------------- ### Examples Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/coda_methods.md Example of fitting a joint model and generating trace, density, and cumulative quantile plots for the fixed effects. ```r # \donttest{ # linear mixed model fits fit_lme1 <- lme(log(serBilir) ~ year:sex + age, random = ~ year | id, data = pbc2) fit_lme2 <- lme(prothrombin ~ sex, random = ~ year | id, data = pbc2) # cox model fit fit_cox <- coxph(Surv(years, status2) ~ age, data = pbc2.id) # joint model fit fit_jm <- jm(fit_cox, list(fit_lme1, fit_lme2), time_var = "year", n_chains = 1L) # trace plot for the fixed effects in the linear mixed submodels traceplot(fit_jm, parm = "betas") # density plot for the fixed effects in the linear mixed submodels densplot(fit_jm, parm = "betas") # cumulative quantile plot for the fixed effects in the linear mixed submodels cumuplot(fit_jm, parm = "betas") ``` -------------------------------- ### Variogram Example Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/variogram.md An example demonstrating the usage of the variogram function with sample data. ```r ind <- aids$patient == 2 yy <- aids$CD4[ind] tt <- aids$obstime[ind] ids <- aids$patient[ind] variogram(yy, tt, ids) #> $svar #> time_lag diffs2 #> [1,] 6 1.61906969 #> [2,] 12 1.51724651 #> [3,] 18 0.87722340 #> [4,] 6 6.27097906 #> [5,] 12 4.87980798 #> [6,] 6 0.08712153 #> #> $sigma2 #> [1] NA #> #> attr(,"class") #> [1] "vrgm" ``` -------------------------------- ### Example Usage Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/accuracy.md This snippet shows an example of how to use the functionality, indicated by \ donttest{ which suggests it's a testable example. ```r # \donttest{ ``` -------------------------------- ### Examples Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/ppcheck.html Example of setting up data for Cox model and posterior predictive checks. ```R # Cox model for the composite event death or transplantation pbc2.id$status2 <- as.numeric(pbc2.id$status != 'alive') pbc2$status2 <- as.numeric(pbc2$status != 'alive') ``` -------------------------------- ### Examples Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/sliced_model_generics.html Example demonstrating the use of slice-aware model fitting generics with `slicer`. ```R if (FALSE) { # \dontrun{ slc <- slicer(n_slices = 2, id_var = "id", data_long = pbc2, data_surv = pbc2.id) n_cores <- max(parallel::detectCores() - 1L, 1L) lme_fit <- lme(fixed = log(serBilir) ~ year * sex, data = slc$long, random = ~ year | id, cores = n_cores) cox_fit <- coxph(formula = Surv(years, status2) ~ sex, data = slc$surv, cores = n_cores) mxm_fit <- mixed_model(fixed = ascites ~ year + sex, data = slc$long, random = ~ year | id, family = binomial(), cores = n_cores) } ``` -------------------------------- ### Example Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/cr_setup.html This example shows how to use the `crisk_setup` function with the `pbc2.id` dataset, specifying 'status' as the status variable and 'alive' as the censoring level. The output displays the first few rows of the transformed data frame, which is now in a long format suitable for competing risks analysis. ```r head(crisk_setup(pbc2.id, "status", "alive")) #> id years status drug age sex year ascites hepatomegaly #> 1 1 1.095170 dead D-penicil 58.76684 female 0 Yes Yes #> 1.1 1 1.095170 dead D-penicil 58.76684 female 0 Yes Yes #> 2 2 14.152338 alive D-penicil 56.44782 female 0 No Yes #> 2.1 2 14.152338 alive D-penicil 56.44782 female 0 No Yes #> 3 3 2.770781 dead D-penicil 70.07447 male 0 No No #> 3.1 3 2.770781 dead D-penicil 70.07447 male 0 No No #> spiders edema serBilir serChol albumin alkaline SGOT #> 1 Yes edema despite diuretics 14.5 261 2.60 1718 138.0 #> 1.1 Yes edema despite diuretics 14.5 261 2.60 1718 138.0 #> 2 Yes No edema 1.1 302 4.14 7395 113.5 #> 2.1 Yes No edema 1.1 302 4.14 7395 113.5 #> 3 No edema no diuretics 1.4 176 3.48 516 96.1 #> 3.1 No edema no diuretics 1.4 176 3.48 516 96.1 #> platelets prothrombin histologic status2 strata #> 1 190 12.2 4 1 dead #> 1.1 190 12.2 4 0 transplanted #> 2 221 10.6 3 0 dead #> 2.1 221 10.6 3 0 transplanted #> 3 151 12.0 4 1 dead #> 3.1 151 12.0 4 0 transplanted ``` -------------------------------- ### Examples Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/sliced_model_generics.md Example of using the slice-aware model fitting generics with sliced data. ```r if (FALSE) { # \dontrun{ slc <- slicer(n_slices = 2, id_var = "id", data_long = pbc2, data_surv = pbc2.id) n_cores <- max(parallel::detectCores() - 1L, 1L) lme_fit <- lme(fixed = log(serBilir) ~ year * sex, data = slc$long, random = ~ year | id, cores = n_cores) cox_fit <- coxph(formula = Surv(years, status2) ~ sex, data = slc$surv, cores = n_cores) mxm_fit <- mixed_model(fixed = ascites ~ year + sex, data = slc$long, random = ~ year | id, family = binomial(), cores = n_cores) } # } ``` -------------------------------- ### Example Usage of jmbayes2 Functions Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/methods.html Demonstrates fitting linear mixed models, a Cox model, and a joint model, followed by examples of extracting coefficients, fixed effects, random effects, terms, model frames, design matrices, family objects, and comparing joint models. ```R # linear mixed model fits fit_lme1 <- lme(log(serBilir) ~ year:sex + age, random = ~ year | id, data = pbc2) fit_lme2 <- lme(prothrombin ~ sex, random = ~ year | id, data = pbc2) # cox model fit fit_cox <- coxph(Surv(years, status2) ~ age, data = pbc2.id) # joint model fit fit_jm <- jm(fit_cox, list(fit_lme1, fit_lme2), time_var = "year", n_chains = 1L, n_iter = 11000L, n_burnin = 1000L) # coef(): fixed effects for the event process coef(fit_jm) #> $gammas #> Mean #> 0.06070451 #> #> $association #> value(log(serBilir)) value(prothrombin) #> 1.3604613 0.1063424 #> # fixef(): fixed effects for the first linear mixed submodel fixef(fit_jm, outcome = 1) #> (Intercept) age year:sexmale year:sexfemale #> 0.66594783 -0.00445036 0.23584373 0.16253524 # ranef(): random effects from all linear mixed submodels head(ranef(fit_jm)) #> [,1] [,2] [,3] [,4] #> [1,] 2.23319313 0.20276852 1.05194350 0.126098411 #> [2,] -0.36290322 0.00382308 -0.03130737 0.089480915 #> [3,] -0.19671871 0.06974686 0.51040201 0.210441867 #> [4,] 0.03344513 0.10239115 0.85942217 0.574548016 #> [5,] 0.32900876 0.22520488 -0.04018424 0.442201332 #> [6,] -0.62705384 -0.16161656 -0.11058573 -0.003544499 # terms(): random effects terms for the first linear mixed submodel terms(fit_jm, process = "longitudinal", type = "random")[[1]] #> ~year #> attr(,"variables") #> list(year) #> attr(,"factors") #> year #> year 1 #> attr(,"term.labels") #> [1] "year" #> attr(,"order") #> [1] 1 #> attr(,"intercept") #> [1] 1 #> attr(,"response") #> [1] 0 #> attr(," .Environment") #> #> attr(,"predvars") #> list(year) #> attr(,"dataClasses") #> year #> "numeric" # mode.frame(): model frame for the fixed effects in the second # linear mixed submodel head(model.frame(fit_jm, process = "longitudinal", type = "fixed")[[2]]) #> prothrombin sex #> 1 12.2 female #> 2 11.2 female #> 3 10.6 female #> 4 11.0 female #> 5 11.6 female #> 6 10.6 female # model.matrix(): fixed effects design matrix for the first linear # mixed submodel head(model.matrix(fit_jm)[[1]]) #> (Intercept) age year:sexmale year:sexfemale #> 1 1 58.76684 0 0.0000000 #> 2 1 58.76684 0 0.5256817 #> 3 1 56.44782 0 0.0000000 #> 4 1 56.44782 0 0.4983025 #> 5 1 56.44782 0 0.9993429 #> 6 1 56.44782 0 2.1027270 # family(): family objects from both linear mixed submodels family(fit_jm) #> [[1]] #> #> #> Family: gaussian #> Link function: identity #> #> #> #> [[2]] #> #> #> Family: gaussian #> Link function: identity #> #> #> # compare_jm(): compare two fitted joint models fit_lme1b <- lme(log(serBilir) ~ 1, random = ~ year | id, data = pbc2) fit_jm2 <- jm(fit_cox, list(fit_lme1b, fit_lme2), time_var = "year", n_chains = 1L, n_iter = 11000L, n_burnin = 1000L) compare_jm(fit_jm, fit_jm2) #> #> #> DIC WAIC LPML #> fit_jm2 10512.30 10540.46 -5268.228 #> fit_jm 10665.69 11146.24 -6097.345 #> #> The criteria are calculated based on the marginal log-likelihood. ``` -------------------------------- ### Examples Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/ppcheck.md Examples of using the ppcheck function. ```R # \donttest{ # Cox model for the composite event death or transplantation pbc2.id$status2 <- as.numeric(pbc2.id$status != 'alive') pbc2$status2 <- as.numeric(pbc2$status != 'alive') CoxFit <- coxph(Surv(years, status2) ~ sex, data = pbc2.id) # a linear mixed model for log serum bilirubin fm1 <- lme(log(serBilir) ~ ns(year, 3) * sex, data = pbc2, random = list(id = pdDiag(~ ns(year, 3)))) # the joint model jointFit <- jm(CoxFit, fm1, time_var = "year", save_random_effects = TRUE) ppcheck(jointFit) FF <- function (t, betas, bi, data) { sex <- as.numeric(data$sex == "female") NS <- ns(t, k = c(0.9911, 3.9863), B = c(0, 14.10579)) X <- cbind(1, NS, sex, NS * sex) Z <- cbind(1, NS) eta <- c(X %*% betas[[1]]) + rowSums(Z * bi) cbind(eta) } ppcheck(jointFit, process = "event", Fforms_fun = FF) # } ``` -------------------------------- ### EPCE Weights Calculation Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/articles/Super_Learning.md Example of calculating EPCE weights using the tvEPCE function. ```r EPCE_weights <- tvEPCE(Models_folds, newdata = CVdats$testing, Tstart = tstr, Thoriz = thor) EPCE_weights #> #> Cross-Validated Expected Predictive Cross-Entropy using the Library of Joint Models 'Models_folds' #> #> Super Learning Estimated EPCE: 0.3109 #> In the time interval: [6, 8) #> For the 166 subjects at risk at time 6 #> Number of subjects with an event in [6, 8): 18 #> Number of subjects with a censored time in [6, 8): 44 #> #> EPCE per model: 0.3568 0.3599 0.3639 0.3589 0.4157 #> Weights per model: 0.0013 0.5387 0.4599 0 1e-04 #> Number of folds: 5 ``` ``` -------------------------------- ### Examples Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/variogram.html An example demonstrating how to use the variogram function with sample data. ```R ind <- aids$patient == 2 yy <- aids$CD4[ind] tt <- aids$obstime[ind] ids <- aids$patient[ind] variogram(yy, tt, ids) #> $svar #> time_lag diffs2 #> [1,] 6 1.61906969 #> [2,] 12 1.51724651 #> [3,] 18 0.87722340 #> [4,] 6 6.27097906 #> [5,] 12 4.87980798 #> [6,] 6 0.08712153 #> #> $sigma2 #> [1] NA #> #> attr(,"class") #> [1] "vrgm" ``` -------------------------------- ### Example of using the slicer function Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/slicer.md An example demonstrating how to use the slicer function with the pbc2 dataset. ```R data(pbc2, package = "JMbayes2") data(pbc2.id, package = "JMbayes2") pbc2_slc <- slicer(n_slices = 2, id_var = "id", data_long = pbc2, data_surv = pbc2.id, seed = 123L) length(pbc2_slc$long) # 2 length(pbc2_slc$surv) # 2 ``` -------------------------------- ### Installation Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/index.html Instructions for installing the JMbayes2 package from CRAN and the development version from GitHub. ```R install.packages("JMbayes2") ``` ```R # install.packages("remotes") remotes::install_github("drizopoulos/jmbayes2") ``` -------------------------------- ### Example 2: Federated or sequential posterior sampling Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/articles/Consensus.html This example demonstrates a workflow for fitting joint models on partitioned datasets (simulating multi-center data or handling memory constraints) and then obtaining consensus posterior estimates. ```R # 1. Prepare the data pbc2.id$status2 <- as.numeric(pbc2.id$status != 'alive') # composite event indicator (1 = death or transplantation, 0 = otherwise) data_slc <- slicer(n_slices = 2, # target number of data subsamples id_var = "id", # clustering variable data_long = pbc2, # longitudinal dataset data_surv = pbc2.id, # survival dataset seed = 123L) # seed for reproducibility data.A <- list(long = data_slc$long[[1]], # data from Center A surv = data_slc$surv[[1]]) data.B <- list(long = data_slc$long[[2]], # data from Center B surv = data_slc$surv[[2]]) # 2. Fit the longitudinal and survival submodels on each subsample ## On Center A ## 2.1.A Fit the linear mixed-effects model lme_fit.A <- lme(data = data.A$long, # list with partitioned longitudinal dataset (class `sliced_data`) fixed = log(serBilir) ~ year * sex, # fixed-effects formula random = ~ year | id) # random-effects formula describing ## 2.2.A Fit the proportional-hazards model ph_fit.A <- coxph(data = data.A$surv, # list with partitioned survival dataset (class `sliced_data`) formula = Surv(years, status2) ~ sex) # model formula # 2.3.A Fit the joint model jm_fit.A <- jm(Surv_object = ph_fit.A, # list with survival submodels (class `sliced_coxph`) Mixed_objects = lme_fit.A, # list with longitudinal submodels (class `sliced_lme` or `MixMod`) time_var = "year") # time variable in the longitudinal submodel jm_fit.A$model_data <- NULL # remove the datasets ## On Center B ## 2.1.B Fit the linear mixed-effects model lme_fit.B <- lme(data = data.B$long, # longitudinal dataset (class `data.frame`) fixed = log(serBilir) ~ year * sex, # fixed-effects formula random = ~ year | id) # random-effects formula describing ## 2.2.B Fit the proportional-hazards model ph_fit.B <- coxph(data = data.B$surv, # survival dataset (class `data.frame`) formula = Surv(years, status2) ~ sex) # model formula # 2.3.B Fit the joint model jm_fit.B <- jm(Surv_object = ph_fit.B, # survival submodels (class `coxph`) Mixed_objects = lme_fit.B, # longitudinal submodels (class `lme` or `MixMod`) time_var = "year") # time variable in the longitudinal submodel jm_fit.B$model_data <- NULL # remove the datasets # 3. Obtain consensus posterior estimates fits <- list(jm_fit.A, jm_fit.B) class(fits) <- c("sliced_jm", class(fits)) # consensus() expects the class `sliced_jm` cons_fit2 <- consensus(fits, # list with joint models (class `sliced_jm`) parm = c("alphas", "betas1"), # parameters of interest method = "var_weight", # consensus algorithm seed = 123L) # seed for reproducibility cons_fit2 ``` -------------------------------- ### Example of Joint Model Fitting and Coefficient Extraction Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/methods.md This example demonstrates fitting a linear mixed model, a Cox model, and then a joint model using jm. It concludes by extracting the fixed effects for the event process using coef(). ```r # linear mixed model fits fit_lme1 <- lme(log(serBilir) ~ year:sex + age, random = ~ year | id, data = pbc2) fit_lme2 <- lme(prothrombin ~ sex, random = ~ year | id, data = pbc2) # cox model fit fit_cox <- coxph(Surv(years, status2) ~ age, data = pbc2.id) # joint model fit fit_jm <- jm(fit_cox, list(fit_lme1, fit_lme2), time_var = "year", n_chains = 1L, n_iter = 11000L, n_burnin = 1000L) # coef(): fixed effects for the event process coef(fit_jm) ``` -------------------------------- ### Dynamic Posterior-Posterior Predictive Checks Setup Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/articles/Posterior_Predictive_Checks.html Preparing the dataset for dynamic posterior-posterior predictive checks by setting a landmark time and adjusting variables. ```R t0 <- 3 prothro_t0 <- prothro[prothro$Time > t0 & prothro$time <= t0, ] prothro_t0$Time <- t0 prothro_t0$death <- 0 ``` -------------------------------- ### Installation from GitHub (Development Version) Source: https://github.com/drizopoulos/jmbayes2/blob/master/README.md Installs the development version of the JMbayes2 package from GitHub using the 'remotes' package. ```r # install.packages("remotes") remotes::install_github("drizopoulos/jmbayes2") ``` -------------------------------- ### Brier Weights Calculation Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/articles/Super_Learning.md Example of calculating Brier weights using the tvBrier function. ```r Brier_weights #> #> Cross-Validated Prediction Error using the Library of Joint Models 'Models_folds' #> #> Super Learning Estimated Integrated Brier score: 0.0548 #> In the time interval: [6, 8) #> For the 166 subjects at risk at time 6 #> Number of subjects with an event in [6, 8): 18 #> Number of subjects with a censored time in [6, 8): 44 #> Accounting for censoring using model-based weights #> #> Integrated Brier score per model: 0.0596 0.0588 0.0613 0.0523 0.0611 #> Weights per model: 0.199 0.1976 0.1925 0.214 0.1969 #> Number of folds: 5 ``` ``` -------------------------------- ### Statistical Model Output Example Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/jm.html Example output from a statistical model, likely from the jmbayes package, showing coefficients, standard deviations, confidence intervals, and p-values for longitudinal outcomes. ```R #> year 0.2605 0.0300 0.2020 0.3203 0.0000 #> sexfemale -0.2451 0.1481 -0.5341 0.0420 0.0988 #> year:sexfemale -0.0716 0.0305 -0.1324 -0.0125 0.0180 #> sigma 0.3483 0.0066 0.3357 0.3618 0.0000 #> #> Longitudinal Outcome: hepatomegaly (family = binomial, link = logit) #> Mean StDev 2.5% 97.5% P #> (Intercept) 0.2790 1.0209 -1.7124 2.2872 0.7926 #> sexfemale -0.9275 0.5312 -1.9833 0.0903 0.0756 #> age 0.0138 0.0165 -0.0190 0.0458 0.3974 #> year 0.2594 0.0730 0.1222 0.4067 0.0004 #> #> Longitudinal Outcome: ascites (family = binomial, link = logit) #> Mean StDev 2.5% 97.5% P #> (Intercept) -8.0466 0.8839 -9.8390 -6.3740 0 #> year 0.4513 0.0598 0.3147 0.5598 0 #> age 0.0733 0.0153 0.0443 0.1052 0 #> #> MCMC summary: #> chains: 1 #> iterations per chain: 11000 #> burn-in per chain: 1000 #> thinning: 1 #> time: 1.5 min ``` -------------------------------- ### Example of crisk_setup function usage Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/articles/Competing_Risks.html Demonstrates how to prepare data for competing risks analysis using the crisk_setup function, showing the transformation of a standard survival dataset into the competing risks long format. ```r pbc2.idCR <- crisk_setup(pbc2.id, statusVar = "status", censLevel = "alive", nameStrata = "CR") ``` -------------------------------- ### Function Signature Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/rc_setup.html The usage signature for the rc_setup function. ```R rc_setup(rc_data, trm_data, idVar = "id", statusVar = "status", startVar = "start", stopVar = "stop", trm_censLevel, nameStrata = "strata", nameStatus = "status") ``` -------------------------------- ### Usage Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/cr_setup.md The function signature for `crisk_setup`. ```R crisk_setup(data, statusVar, censLevel, nameStrata = "strata", nameStatus = "status2") ``` -------------------------------- ### Example 1: Speeding posterior sampling (parallel fitting) Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/articles/Consensus.md This code block demonstrates how to speed up posterior sampling by splitting the data into subsamples and fitting in parallel on a single machine. ```r library(JMbayes2) # Load the pbc2 and pbc2.id datasets data("pbc2") data("pbc2.id") # Define the model formula for the longitudinal and survival components long.formula <- Surv2(time, event, Ltime, Levent, data = pbc2, id = id, time_vary = NULL, nTstop = 1, max_Ltime = max(pbc2$time, na.rm = TRUE)) ~ sex + age + log(bili) + albumin + stage surv.formula <- Surv2(time, event, Ltime, Levent, data = pbc2, id = id, time_vary = NULL, nTstop = 1, max_Ltime = max(pbc2$time, na.rm = TRUE)) ~ sex + age + albumin + stage # Split the data into 4 subsamples # Note: For reproducibility, we set a seed. set.seed(123) subsamples <- split_data(pbc2, nsplit = 4, id = pbc2$id) # Fit the joint model on each subsample in parallel # We use the 'parallel' argument to specify the number of cores to use. # For this example, we use 4 cores. fit_parallel <- jointModelBayes2(long.formula, surv.formula, data = pbc2, id = pbc2$id, subsamples = subsamples, parallel = 4) # The 'fit_parallel' object now contains the consensus posterior draws. # We can summarize the results as usual. summary(fit_parallel) # We can also compare the consensus posterior with the full-data posterior # (if feasible to compute). # For this example, we fit the full-data model for comparison. fit_full <- jointModelBayes2(long.formula, surv.formula, data = pbc2, id = pbc2$id) # Compare the summary statistics of the consensus and full-data models. summary(fit_parallel)$statistics[1:5, ] summary(fit_full)$statistics[1:5, ] ``` -------------------------------- ### Data Setup for Joint Model Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/articles/Recurring_Events.html R code demonstrating the use of the `rc_setup` function to prepare data for fitting a joint model, combining terminal and recurrent event data. ```R cox_data <- rc_setup(rc_data = recu_data, trm_data = term_data, idVar = "id", statusVar = "status", startVar = "tstart", stopVar = "tstop", trm_censLevel = 0, nameStrata = "strata", nameStatus = "status") ``` -------------------------------- ### Installation from CRAN Source: https://github.com/drizopoulos/jmbayes2/blob/master/README.md Installs the JMbayes2 package from the Comprehensive R Archive Network (CRAN). ```r install.packages("JMbayes2") ``` -------------------------------- ### Data Setup for Cox Model Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/articles/Recurring_Events.md Using the rc_setup function to combine terminal and recurrent event data into a format suitable for a Cox model, including strata for event types. ```R cox_data <- rc_setup(rc_data = recu_data, trm_data = term_data, idVar = "id", statusVar = "status", startVar = "tstart", stopVar = "tstop", trm_censLevel = 0, nameStrata = "strata", nameStatus = "status") ``` -------------------------------- ### Example Usage of slicer function Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/slicer.html This example demonstrates how to use the slicer function to split longitudinal and survival data into two slices. It then checks the number of slices created for both longitudinal and survival data. ```r data(pbc2, package = "JMbayes2") data(pbc2.id, package = "JMbayes2") pbc2_slc <- slicer(n_slices = 2, id_var = "id", data_long = pbc2, data_surv = pbc2.id, seed = 123L) length(pbc2_slc$long) # 2 #> [1] 2 length(pbc2_slc$surv) # 2 #> [1] 2 ``` -------------------------------- ### Predicting with jmbayes2 Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/predict.html Examples of using the plot function for predictions with jmbayes2. ```R plot(predLong2, predSurv, outcomes = 1) # the first two longitudinal outcomes plot(predLong1, predSurv, outcomes = 1:2) # all three longitudinal outcomes, we display survival probabilities instead # of cumulative risk, and we transform serum bilirubin to the original scale plot(predLong2, predSurv, outcomes = 1:3, fun_event = function (x) 1 - x, fun_long = list(exp, identity, identity), ylab_event = "Survival Probabilities", ylab_long = c("Serum Bilirubin", "Prothrombin", "Ascites"), pos_ylab_long = c(1.9, 1.9, 0.08)) # } ``` -------------------------------- ### Data Filtering and Preparation Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/articles/Multi_State_Processes.md This code snippet demonstrates how to filter longitudinal data to the observed period and prepare the survival data for multi-state modeling. ```r df_long2 <- merge(df_long, df_surv[, c("id", "Tstop")], by = c("id")) df_long2 <- df_long2[!duplicated(df_long2, by = c('id', 'time')), ] df_long2 <- df_long2[df_long2$time <= df_long2$Tstop, ] df_surv$transition <- factor(df_surv$transition) ``` -------------------------------- ### Usage Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/coda_methods.html Usage of the coda methods for 'jm' objects. ```R traceplot(object, ...) # S3 method for class 'jm' traceplot(object, parm = c("all", "betas", "sigmas", "D", "bs_gammas", "tau_bs_gammas", "gammas", "alphas"), ...) ggtraceplot(object, ...) # S3 method for class 'jm' ggtraceplot(object, parm = c("all", "betas", "sigmas", "D", "bs_gammas", "tau_bs_gammas", "gammas", "alphas"), linewidth = 1, alpha = 0.8, theme = c('standard', 'catalog', 'metro', 'pastel', 'beach', 'moonlight', 'goo', 'sunset', 'custom'), grid = FALSE, gridrows = 3, gridcols = 1, custom_theme = NULL, ...) gelman_diag(object, ...) # S3 method for class 'jm' gelman_diag(object, parm = c("all", "betas", "sigmas", "D", "bs_gammas", "tau_bs_gammas", "gammas", "alphas"), ...) densplot(object, ...) # S3 method for class 'jm' densplot(object, parm = c("all", "betas", "sigmas", "D", "bs_gammas", "tau_bs_gammas", "gammas", "alphas"), ...) ggdensityplot(object, ...) # S3 method for class 'jm' ggdensityplot(object, parm = c("all", "betas", "sigmas", "D", "bs_gammas", "tau_bs_gammas", "gammas", "alphas"), linewidth = 1, alpha = 0.6, theme = c('standard', 'catalog', 'metro', 'pastel', 'beach', 'moonlight', 'goo', 'sunset', 'custom'), grid = FALSE, gridrows = 3, gridcols = 1, custom_theme = NULL, ...) cumuplot(object, ...) # S3 method for class 'jm' cumuplot(object, parm = c("all", "betas", "sigmas", "D", "bs_gammas", "tau_bs_gammas", "gammas", "alphas"), ...) ``` -------------------------------- ### Summarizing the Joint Model Fit Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/articles/Recurring_Events.html Example of summarizing the fitted joint model to view results. ```R summary(jm_fit) ``` -------------------------------- ### Fitting the Stratified Cox Model Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/articles/Recurring_Events.html Example of fitting a stratified Cox model using coxph() with transformed data. ```R cox_fit <- coxph(Surv(tstart, tstop, status) ~ (group + age):strata(strata), data = cox_data) ``` -------------------------------- ### Usage Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/reference/methods.html Usage of various methods for objects of class 'jm'. ```R coef(object, ...) # S3 method for class 'jm' coef(object, ...) fixef(object, ...) # S3 method for class 'jm' fixef(object, outcome = Inf, ...) ranef(object, ...) # S3 method for class 'jm' ranef(object, outcome = Inf, post_vars = FALSE, ...) terms(x, ...) # S3 method for class 'jm' terms(x, process = c("longitudinal", "event"), type = c("fixed", "random"), ...) model.frame(formula, ...) # S3 method for class 'jm' model.frame(formula, process = c("longitudinal", "event"), type = c("fixed", "random"), ...) model.matrix(object, ...) # S3 method for class 'jm' model.matrix(object, ...) family(object, ...) # S3 method for class 'jm' family(object, ...) compare_jm(..., type = c("marginal", "conditional"), order = c("WAIC", "DIC", "LPML", "none")) ``` -------------------------------- ### Data Preparation for Competing Risks Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/articles/Competing_Risks.md Example of preparing survival data for competing risks analysis using the `crisk_setup` function in JMbays2. It shows the transformation of standard survival data into a long format suitable for joint models, where each subject has multiple rows representing different competing event types. ```r 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 ``` -------------------------------- ### Fitting the Linear Mixed Model Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/articles/Recurring_Events.html Example of fitting a linear mixed model for the longitudinal outcome using lme(). ```R lme_fit <- lme(y ~ ns(time, k = c(1, 3), B = c(0, 7)), random = list(id = pdDiag(form = ~ ns(time, k = c(1, 3), B = c(0, 7)))), data = lme_data, control = lmeControl(opt = "optim", niterEM = 45)) ``` -------------------------------- ### Fitting the Joint Model Source: https://github.com/drizopoulos/jmbayes2/blob/master/docs/articles/Recurring_Events.html Example of fitting the joint model using jm() with specified functional forms and recurrent timescale. ```R jm_fit <- jm(cox_fit, lme_fit, time_var = "time", recurrent = "gap", functional_forms = ~ value(y):strata) ```