JMbayes2

https://drizopoulos.github.io/JMbayes2/
Admin
JMbayes2 is an R package for fitting joint models to analyze longitudinal and time-to-event data...
Tokens:64,491
Snippets:590
Trust Score:9.7
Update:4 months ago
Show doc for...
Tokens:
Raw
### Example Usage

Source: https://drizopoulos.github.io/JMbayes2/reference/cr_setup.html

An example demonstrating the usage of the `crisk_setup` function with the `pbc2.id` dataset.

```R
head(crisk_setup(pbc2.id, "status", "alive"))
```

--------------------------------

### Examples

Source: https://drizopoulos.github.io/JMbayes2/reference/slicer.html

Example of 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
#> [1] 2
length(pbc2_slc$surv) # 2
#> [1] 2
```

--------------------------------

### Installation from GitHub (Development Version)

Source: https://drizopoulos.github.io/JMbayes2/index.html

To install the development version of JMbayes2 from GitHub, use the following commands.

```r
# install.packages("remotes")
remotes::install_github("drizopoulos/jmbayes2")
```

--------------------------------

### Example of fitting a joint model and using coda methods

Source: https://drizopoulos.github.io/JMbayes2/reference/coda_methods.html

This example demonstrates how to fit a linear mixed model, a Cox model, and then combine them into a joint model using jm(). It also shows how to use the coda methods on the fitted joint model object.

```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)


```

--------------------------------

### Example Usage

Source: https://drizopoulos.github.io/JMbayes2/reference/sliced_model_generics.html

Example demonstrating the use of lme, coxph, and mixed_model with a 'sliced_data' object.

```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 Methods

Source: https://drizopoulos.github.io/JMbayes2/reference/methods.html

Illustrates how to use JMBayes2 methods by fitting example linear mixed models and a Cox model, which can then be used with joint modeling functions.

```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)
```

--------------------------------

### Installation from CRAN

Source: https://drizopoulos.github.io/JMbayes2/index.html

To install JMbayes2 from CRAN, use the following command.

```bash
install.packages("JMbayes2")
```

--------------------------------

### Examples

Source: https://drizopoulos.github.io/JMbayes2/reference/ppcheck.html

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)


# }

```

--------------------------------

### Example 1: Speeding posterior sampling (parallel fitting)

Source: https://drizopoulos.github.io/JMbayes2/articles/Consensus.html

This code demonstrates how to split a dataset, fit longitudinal and survival submodels in parallel, and then obtain consensus posterior estimates for a joint model.

```R
# 1. Split the dataset into `n_slices` subsamples
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

# 2. Specify the number of available cores
n_cores <- max(1L, parallelly::availableCores(omit = 1L))

# 3. Fit the longitudinal and survival submodels on each subsample
## 3.1 Fit a linear mixed-effects model 
lme_fit <- lme(data   = data_slc$long, # list with partitioned longitudinal dataset (class `sliced_data`)
               fixed  = log(serBilir) ~ year * sex, # fixed-effects formula  
               random = ~ year | id, # random-effects formula describing 
               cores  = n_cores) # number of available cores

## 3.2 Fit a proportional-hazards model
ph_fit <- coxph(data    = data_slc$surv, # list with partitioned survival dataset (class `sliced_data`)
                formula = Surv(years, status2) ~ sex, # model formula
                cores   = n_cores) # number of available cores

# 4. Fit joint model(s) on each subsample
jm_fit <- jm(Surv_object   = ph_fit, # list with survival submodels (class `sliced_coxph`)
             Mixed_objects = lme_fit, # list with longitudinal submodels (class `sliced_lme`)
             time_var      = "year", # time variable in the longitudinal submodel
             cores         = n_cores) # number of available cores

# 5. Obtain consensus posterior estimates
cons_fit <- consensus(jm_fit, # 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_fit
#> 
#> Consensus summary (2 splits)
#> Method (var_weight): iteration-wise weighted average across splits using inverse-variance weights.
#> 
#> Parameter block: alphas (9000 draws)
#>                        Mean  StDev   2.5%  97.5% P   w_S1   w_S2
#> value(log(serBilir)) 1.2503 0.0885 1.0817 1.4297 0 0.4261 0.5739
#> 
#> Parameter block: betas1 (9000 draws)
#>                   Mean  StDev    2.5%   97.5%      P   w_S1   w_S2
#> (Intercept)     0.7235 0.1693  0.3870  1.0565 0.0000 0.6853 0.3147
#> year            0.2648 0.0381  0.1911  0.3407 0.0000 0.5906 0.4094
#> sexfemale      -0.2410 0.1844 -0.5961  0.1181 0.1880 0.6668 0.3332
#> year:sexfemale -0.0875 0.0406 -0.1678 -0.0072 0.0329 0.5706 0.4294
```
```

--------------------------------

### Data Setup for Joint Model

Source: https://drizopoulos.github.io/JMbayes2/articles/Recurring_Events.html

This R code uses the `rc_setup` function to combine terminal and recurrent event data into a format suitable for fitting a joint model, using the counting process formulation.

```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")
```

--------------------------------

### Minimal Example

Source: https://drizopoulos.github.io/JMbayes2/index.html

This example demonstrates fitting a joint model with three longitudinal outcomes using the PBC dataset. It involves fitting a Cox model for the event process and linear mixed-effects models for the longitudinal outcomes, then combining them using the jm() function.

```r
# Cox model for the composite event death or transplantation
pbc2.id$status2 <- as.numeric(pbc2.id$status != 'alive')
CoxFit <- coxph(Surv(years, status2) ~ sex, data = pbc2.id)

# a linear mixed model for log serum bilirubin
fm1 <- lme(log(serBilir) ~ year * sex, data = pbc2, random = ~ year | id)

# a linear mixed model for the prothrombin time
fm2 <- lme(prothrombin ~ year * sex, data = pbc2, random = ~ year | id)

# a mixed effects logistic regression for ascites
fm3 <- mixed_model(ascites ~ year + sex, data = pbc2,
 random = ~ year | id, family = binomial())

# the joint model that links all sub-models
jointFit <- jm(CoxFit, list(fm1, fm2, fm3), time_var = "year",
 n_iter = 12000L, n_burnin = 2000L, n_thin = 5L)
summary(jointFit)
```

--------------------------------

### Variogram function usage

Source: https://drizopoulos.github.io/JMbayes2/reference/variogram.html

Demonstrates the basic usage of the variogram function with example data.

```R
ind <- aids$patient == 2
yy <- aids$CD4[ind]
tt <- aids$obstime[ind]
olds <- aids$patient[ind]
variogram(yy, tt, ids)
#&gt; $svar
#&gt;      time_lag     diffs2
#&gt; [1,]        6 1.61906969
#&gt; [2,]       12 1.51724651
#&gt; [3,]       18 0.87722340
#&gt; [4,]        6 6.27097906
#&gt; [5,]       12 4.87980798
#&gt; [6,]        6 0.08712153
#&gt; 
#&gt; $sigma2
#&gt; [1] NA
#&gt; 
#&gt; attr(,"class")
#&gt; [1] "vrgm"
```

--------------------------------

### Data Preprocessing and Simulation Setup

Source: https://drizopoulos.github.io/JMbayes2/articles/Non_Gaussian_Mixed_Models.html

This code block preprocesses the data by handling potential outliers and setting up variables for simulating event times based on a survival model. It includes defining parameters, creating a design matrix, and simulating event times using inverse transform sampling.

```R
# we assume that values below -4 are not observed, and set equal to -4
DF$ind <- as.numeric(DF$y < -4)
DF$y <- pmax(DF$y, -4)

upp_Cens <- 15 # fixed Type I censoring time
shape_wb <- 5 # shape Weibull
alpha <- 0.8 # association coefficients
gammas <- c("(Intercept)" = -9, "sex" = 0.5)
W <- model.matrix(~ sex, data = DF[!duplicated(DF$id), ])
# linear predictor for the survival model
eta_t <- as.vector(W %*% gammas)
# to simulate event times we use inverse transform sampling
# (https://en.wikipedia.org/wiki/Inverse_transform_sampling). Namely, we want 
# to find t, such that S(t) = u, where S(.) is the survival function, and u a 
# number from the Unif(0, 1) distribution. The function below calculates 
# log(u) - log(S(t)), and for a given u, we want to find t for which it equals
# zero. We do that below using the uniroot() function
invS <- function (t, i) {
  # i denotes the subject
  sex_i <- W[i, 2L]
  # h() is the hazard function and we assume a Weibull baseline hazard
  h <- function (s) {
    X_at_s <- cbind(1, sex_i, s, sex_i * s)
    Z_at_s <- cbind(1, s)
    # the linear predictor of the mixed model evaluated at time s
    f <- as.vector(X_at_s %*% betas + 
                     rowSums(Z_at_s * b[rep(i, nrow(Z_at_s)), ]))
    exp(log(shape_wb) + (shape_wb - 1) * log(s) + eta_t[i] + f * alpha)
  }
  # -log(S(t)) = H(t), where H(t) is the cumulative hazard function
  integrate(h, lower = 0, upper = t)$value + log(u[i])
}
# we simulate the event times
u <- runif(n)
trueTimes <- numeric(n)
for (i in seq_len(n)) {
    Up <- 100
    Root <- try(uniroot(invS, interval = c(1e-05, Up), i = i)$root, TRUE)
    trueTimes[i] <- if (!inherits(Root, "try-error")) Root else 150
}

# we use fixed Type I right censoring denoting the end of the trial.
Ctimes <- upp_Cens
Time <- pmin(trueTimes, Ctimes)
event <- as.numeric(trueTimes <= Ctimes) # event indicator

# we keep the longitudinal measurements before the event times
DF$Time <- Time[DF$id]
DF$event <- event[DF$id]
DF <- DF[DF$time <= DF$Time, ]

```

--------------------------------

### Example of eta calculation

Source: https://drizopoulos.github.io/JMbayes2/articles/Posterior_Predictive_Checks.html

This snippet shows the calculation of eta, which is likely a linear predictor in a statistical model, and its return as a matrix.

```R
eta <- c(X %*% betas[[1]]) + rowSums(Z * bi)
    # we return as a matrix
    cbind(eta)
}
```

--------------------------------

### Example of fitting a multivariate joint model and preparing for predictions

Source: https://drizopoulos.github.io/JMbayes2/reference/predict.html

This code snippet demonstrates how to fit a multivariate joint model using the jm function, which is a prerequisite for making predictions.

```R
# We fit a multivariate joint model
pbc2.id$status2 <- as.numeric(pbc2.id$status != 'alive')
CoxFit <- coxph(Surv(years, status2) ~ sex, data = pbc2.id)
fm1 <- lme(log(serBilir) ~ ns(year, 3) * sex, data = pbc2,
           random = ~ ns(year, 3) | id, control = lmeControl(opt = 'optim'))
fm2 <- lme(prothrombin ~ ns(year, 2) * sex, data = pbc2,
           random = ~ ns(year, 2) | id, control = lmeControl(opt = 'optim'))
fm3 <- mixed_model(ascites ~ year * sex, data = pbc2,
                   random = ~ year | id, family = binomial())

jointFit <- jm(CoxFit, list(fm1, fm2, fm3), time_var = "year", n_chains = 1L)

# we select the subject for whom we want to calculate predictions

```

--------------------------------

### Example Data Structure for Subject 1

Source: https://drizopoulos.github.io/JMbayes2/articles/Recurring_Events.html

This R code displays the structure of the prepared data for subject 1, illustrating the counting process formulation with recurrent ('R') and terminal ('T1') events.

```R
cox_data[cox_data$id == 1, c("id", "tstart", "tstop", "status", "strata")]
```

--------------------------------

### rc_setup Function Usage

Source: https://drizopoulos.github.io/JMbayes2/reference/rc_setup.html

The usage signature for the rc_setup function, detailing its arguments.

```R
rc_setup(rc_data, trm_data,
    idVar = "id", statusVar = "status",
    startVar = "start", stopVar = "stop",
    trm_censLevel,
    nameStrata = "strata", nameStatus = "status")
```

--------------------------------

### Coda Package Methods

Source: https://drizopoulos.github.io/JMbayes2/reference/index.html

Functions for visualizing and diagnosing MCMC convergence, adapted from the 'coda' package.

```R
traceplot()
ggtraceplot()
gelman_diag()
densplot()
ggdensityplot()
cumuplot()
```

--------------------------------

### Calculate Cross-Validated Brier Score Average

Source: https://drizopoulos.github.io/JMbayes2/articles/Dynamic_Predictions.html

Computes the average of the Brier scores across all folds to get the cross-validated estimate.

```R
average_Brier <- mean(sapply(Brier_per_fold, "||", "Brier"))
average_Brier
#> [1] 0.0862818
```

--------------------------------

### Function Signature

Source: https://drizopoulos.github.io/JMbayes2/reference/cr_setup.html

The signature of the `crisk_setup` function.

```R
crisk_setup(data, statusVar,
    censLevel,
    nameStrata = "strata", nameStatus = "status2")
```

--------------------------------

### Prepare data for competing risks

Source: https://drizopoulos.github.io/JMbayes2/articles/Competing_Risks.html

Illustrates the use of the crisk_setup() function to transform standard survival data into the competing risks long format, using the PBC dataset where death and transplantation are treated as competing risks.

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

--------------------------------

### Plotting predictions

Source: https://drizopoulos.github.io/JMbayes2/reference/predict.html

Various examples of plotting longitudinal and event outcome predictions, including combined plots and customized labels.

```r
plot(predLong1)

# for subject 25, outcomes in reverse order
plot(predLong2, outcomes = 3:1, subject = 25)




# prediction for the event outcome
plot(predSurv)


# combined into one plot, the first longitudinal outcome and cumulative risk
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))
```

--------------------------------

### Prepare the data

Source: https://drizopoulos.github.io/JMbayes2/articles/Consensus.html

This code snippet prepares the data by creating a composite event indicator and splitting the dataset into two subsamples for longitudinal and survival analysis.

```R
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]]) 
```

--------------------------------

### Expected Predictive Cross-Entropy (EPCE) Calculation

Source: https://drizopoulos.github.io/JMbayes2/articles/Super_Learning.html

Demonstrates the calculation of EPCE weights using the tvEPCE function.

```R
EPCE_weights <- tvEPCE(Models_folds, newdata = CVdats$testing, 
                       Tstart = tstr, Thoriz = thor)
EPCE_weights
```

--------------------------------

### Multivariate Joint Model Fitting and ROC Curve Calculation

Source: https://drizopoulos.github.io/JMbayes2/reference/accuracy.html

This example demonstrates fitting a multivariate joint model using the `jm` function and then calculating time-dependent ROC curves using `tvROC` and `tvAUC`.

```R
# \donttest{
# We fit a multivariate joint model
pbc2.id$status2 <- as.numeric(pbc2.id$status != 'alive')
CoxFit <- coxph(Surv(years, status2) ~ sex, data = pbc2.id)
fm1 <- lme(log(serBilir) ~ ns(year, 3) * sex, data = pbc2,
           random = ~ ns(year, 3) | id, control = lmeControl(opt = 'optim'))
fm2 <- lme(prothrombin ~ ns(year, 2) * sex, data = pbc2,
           random = ~ ns(year, 2) | id, control = lmeControl(opt = 'optim'))
fm3 <- mixed_model(ascites ~ year * sex, data = pbc2,
                   random = ~ year | id, family = binomial())

jointFit <- jm(CoxFit, list(fm1, fm2, fm3), time_var = "year", n_chains = 1L)

roc <- tvROC(jointFit, newdata = pbc2, Tstart = 4, Dt = 3, cores = 1L)
roc
#>
#>
#>	Time-dependent Sensitivity and Specificity for the Joint Model jointFit
#>
#>
#>
#>	At time: 7
#>
#>	Using information up to time: 4 (225 subjects still at risk)
#>
#>	Accounting for censoring using model-based weights
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
#>
A very detailed explanation of the `tvROC` output is provided. The output is a data frame with columns `cut-off`, `SN` (Sensitivity), and `SP` (Specificity). The `cut-off` column represents the threshold for the predicted risk score. The `SN` column shows the sensitivity at each cut-off, which is the proportion of actual positive cases that are correctly identified by the model. The `SP` column shows the specificity at each cut-off, which is the proportion of actual negative cases that are correctly identified by the model. The `tvAUC` function then calculates the time-dependent Area Under the ROC Curve (AUC), which is a measure of the model's overall ability to discriminate between positive and negative cases over time. An AUC of 1 indicates perfect discrimination, while an AUC of 0.5 indicates discrimination no better than chance. The `plot` function with `optimal_cutoff = "Youden"` visualizes the ROC curve and highlights the optimal cut-off point based on the Youden index, which maximizes the sum of sensitivity and specificity.

```

--------------------------------

### Usage of Standard Generic Functions

Source: https://drizopoulos.github.io/JMbayes2/reference/methods.html

Demonstrates the usage of standard generic functions like coef, fixef, ranef, terms, model.frame, model.matrix, family, and compare_jm for 'jm' objects.

```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"))
```

--------------------------------

### Fit joint model with custom functional forms

Source: https://drizopoulos.github.io/JMbayes2/reference/jm.html

This example shows how to extend a joint model by specifying custom functional forms for the association between longitudinal outcomes and the survival outcome using the 'functional_forms' argument in jm(). It includes value, slope, and area terms.

```R
# [1] Fit the mixed-effects models using lme() for continuous outcomes
# and mixed_model() for categorical outcomes.
fm1 <- lme(fixed = log(serBilir) ~ year * sex, random = ~ year | id, data = pbc2)

fm2 <- mixed_model(hepatomegaly ~ sex + age + year, data = pbc2,
                   random = ~ year | id, family = binomial())

fm3 <- mixed_model(ascites ~ year + age, data = pbc2,
                   random = ~ year | id, family = binomial())

# [2] Save all the fitted mixed-effects models in a list.
Mixed <- list(fm1, fm2, fm3)

# [3] Fit a Cox model, specifying the baseline covariates to be included in the
# joint model.
fCox1 <- coxph(Surv(years, status2) ~ drug + age, data = pbc2.id)
```

--------------------------------

### Data Filtering and Preparation

Source: https://drizopoulos.github.io/JMbayes2/articles/Multi_State_Processes.html

This code snippet filters longitudinal data to the observed period and prepares the survival data by converting the transition factor.

```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)
```

--------------------------------

### Univariate joint model for serum bilirubin

Source: https://drizopoulos.github.io/JMbayes2/reference/jm.html

Example of fitting a univariate joint model with one continuous outcome (serum bilirubin). It includes fitting a mixed-effects model using lme(), a Cox model using coxph(), and then the joint model using jm().

```R
# [1] Fit the mixed model using lme().
fm1 <- lme(fixed = log(serBilir) ~ year * sex + I(year^2) +
           age + prothrombin, random =  ~ year | id, data = pbc2)

# [2] Fit a Cox model, specifying the baseline covariates to be included in the
# joint model.
fCox1 <- coxph(Surv(years, status2) ~ drug + age, data = pbc2.id)

# [3] The basic joint model is fitted using a call to jm() i.e.,
joint_model_fit_1 <- jm(fCox1, fm1, time_var = "year",
        n_chains = 1L, n_iter = 11000L, n_burnin = 1000L)
summary(joint_model_fit_1)
#> 
#> Call:
#> JMbayes2::jm(Surv_object = fCox1, Mixed_objects = fm1, time_var = "year", 
#>     n_chains = 1L, n_iter = 11000L, n_burnin = 1000L)
#> 
#> Data Descriptives:
#> Number of Groups: 312		Number of events: 140 (44.9%)
#> Number of Observations:
#>   log(serBilir): 1945
#> 
#>                  DIC     WAIC      LPML
#> marginal    4204.643 5075.364 -2931.823
#> conditional 3334.409 3165.829 -1817.475
#> 
#> Random-effects covariance matrix:
#>                     
#>        StdDev   Corr
#> (Intr) 0.9752 (Intr)
#> year   0.1772 0.3429
#> 
#> Survival Outcome:
#>                         Mean  StDev    2.5%  97.5%     P
#> drugD-penicil        -0.0294 0.2359 -0.4760 0.4417 0.913
#> age                   0.0639 0.0092  0.0465 0.0824 0.000
#> value(log(serBilir))  1.4254 0.0976  1.2296 1.6264 0.000
#> 
#> Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
#>                   Mean  StDev    2.5%   97.5%      P
#> (Intercept)     0.2444 0.3605 -0.4688  0.9355 0.4880
#> year            0.2281 0.0370  0.1568  0.3005 0.0000
#> sexfemale      -0.2421 0.1816 -0.6013  0.1134 0.1850
#> I(year^2)       0.0026 0.0010  0.0007  0.0045 0.0066
#> age            -0.0017 0.0054 -0.0122  0.0092 0.7452
#> prothrombin     0.0529 0.0085  0.0364  0.0695 0.0000
#> year:sexfemale -0.0881 0.0385 -0.1644 -0.0137 0.0206
#> sigma           0.3452 0.0068  0.3322  0.3594 0.0000
#> 
#> MCMC summary:
#> chains: 1 
#> iterations per chain: 11000 
#> burn-in per chain: 1000 
#> thinning: 1 
#> time: 33 sec
traceplot(joint_model_fit_1)

```

--------------------------------

### ggdensityplot with grid

Source: https://drizopoulos.github.io/JMbayes2/reference/coda_methods.html

Generates a density plot with a grid for the fixed effects (betas) using ggdensityplot.

```R
ggdensityplot(fit_jm, parm = "betas", grid = TRUE)
```

--------------------------------

### ggdensityplot with custom theme

Source: https://drizopoulos.github.io/JMbayes2/reference/coda_methods.html

Generates a density plot with a custom theme for the fixed effects (betas) using ggdensityplot.

```R
ggdensityplot(fit_jm, parm = "betas", custom_theme = c('1' = 'black'))
```

--------------------------------

### Prepare Data for Landmark Time

Source: https://drizopoulos.github.io/JMbayes2/articles/Posterior_Predictive_Checks.html

Preparing the dataset for dynamic posterior-posterior predictive checks at a landmark time tL=3 years.

```R
t0 <- 3
prothro_t0 <- prothro[prothro$Time > t0 & prothro$time <= t0, ]
prothro_t0$Time <- t0
prothro_t0$death <- 0
```

--------------------------------

### Constructing Dataset for Prediction

Source: https://drizopoulos.github.io/JMbayes2/articles/Super_Learning.html

Constructs a dataset with subjects at risk at year six, considering longitudinal measurements before that time, and setting event time and status.

```R
ND <- pbc2[pbc2$years > tstr & pbc2$year <= tstr, ]
ND$id <- ND$id[, drop = TRUE]
ND$years <- tstr
ND$status2 <- 0
```

--------------------------------

### JMbayes2 Package

Source: https://drizopoulos.github.io/JMbayes2/reference/index.html

The main package object and its associated documentation.

```R
JMbayes2-package
JMbayes2
```
Tokens:
Raw