---
title: "Introduction to FastJM"
author: "FastJM Team"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to FastJM}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, echo=FALSE, results="hide", warning=FALSE}
suppressPackageStartupMessages({
  library(FastJM)
})
knitr::opts_chunk$set(
  fig.width = 12,
  fig.height = 12,
  out.width = "100%",
  dpi = 150,
  collapse = TRUE, 
  comment = "#>", 
  warning=FALSE, 
  message=FALSE, 
  eval=TRUE)
```

# FastJM

<!-- badges: start -->
[![R-CMD-check](https://github.com/shanpengli/FastJM/workflows/R-CMD-check/badge.svg)](https://github.com/shanpengli/FastJM/actions)
[![metacran downloads](https://cranlogs.r-pkg.org/badges/FastJM)](https://cran.r-project.org/package=FastJM)
[![](https://cranlogs.r-pkg.org/badges/grand-total/FastJM)](https://cran.r-project.org/package=FastJM)
[![CRAN_time_from_release](https://www.r-pkg.org/badges/ago/FastJM)](https://cran.r-project.org/package=FastJM)
[![CRAN_Status_Badge_version_last_release](https://www.r-pkg.org/badges/version-last-release/FastJM)](https://cran.r-project.org/package=FastJM)
<!-- badges: end -->

The `FastJM` package implements efficient computation of semi-parametric joint model of longitudinal and competing risks data. To view a brief guide on the purpose and use of this package, please refer to our [introductory video](https://youtu.be/sspYjUATICM?si=idTbVgT5DswN-yhe).

# Installation

To download the CRAN version of the FastJM package, please refer to the code chunk below. 

```{r, eval=FALSE}
# install.packages("FastJM")
```

Use the following in order to install the development version of the package. 

```{r, eval=FALSE}
# install.packages("remotes")
remotes::install_github("shanpengli/FastJM")
```

# Examples

## Single-biomarker joint model (`jmcs`)

The `FastJM` package comes with several simulated datasets. To fit a joint model, we use `jmcs` function. In the example below, we are using the following built-in data sets:

- ydata: longitudinal data for a **single** biomarker per patient
- cdata: competing risks time-to-event data per patient

```{r}
require(FastJM)
require(survival)
data(ydata)
data(cdata)
fit <- jmcs(ydata = ydata, cdata = cdata, 
            long.formula = response ~ time + gender + x1 + race, 
            surv.formula = Surv(surv, failure_type) ~ x1 + gender + x2 + race, 
            random =  ~ time| ID)
fit
```

The fitted `jmcs` object also provides a standard diagnostic plotting method through `plot()`. This function displays four panels: subject-specific residuals for the longitudinal process versus their corresponding fitted values; a normal Q-Q plot of the residuals of the standardized subject-specific residuals for the longitudinal process, an estimate of the marginal survival function for the event process, and an estimate of the marginal cumulative risk function for the event process. These plots provide a quick graphical summary of the longitudinal residual behavior and the estimated event-time process from the fitted joint model.

```{r}
plot(fit)
```
We can further examine the fitted joint model using diagnostic plots. The `timeplot()` function displays the longitudinal biomarker trajectories, the empirical log residual variance over follow-up time, and the event process. For competing-risk models, the event plot is shown as cumulative incidence curves for the specified event types.

The log residual variance plot is intended as an exploratory diagnostic. Apparent changes over time may reflect departures from constant residual variance, but they may also be affected by finite-sample variability, model fit, or changing subject composition over follow-up.
```{r}
crplot <- timeplot(
  object = fit,
  biomarker = response,
  id_col = ID,
  time_col = time,
  time_bin_width = 0.5,
  fail_code = 1,
  cr_code = 2,
  censor_code = 0,
  primary_event_label = "Event type 1",
  competing_event_label = "Event type 2",
  x_lab = "Visit time",
  event_x_lab = "Survival time",
  n.obs = 200
)
```

The `FastJM` package can make dynamic prediction given the longitudinal history information. Below is a toy example for competing risks data. Conditional cumulative incidence probabilities for each failure will be presented.

```{r}
ND <- ydata[ydata$ID %in% c(419, 218), ]
ID <- unique(ND$ID)
NDc <- cdata[cdata$ID  %in% ID, ]
survfit <- survfitJM(fit, 
                     ynewdata = ND, 
                     cnewdata = NDc, 
                     u = seq(3, 4.8, by = 0.2), 
                     method = "GH",
                     obs.time = "time")
survfit
```

To assess the prediction accuracy of the fitted joint model, we may run `DynPredAcc` to assess the prediction accuracy by calculating all available evaluation metrics. 
```{r}
res <- DynPredAcc(
  object = fit,
  landmark.time = 3,
  horizon.time = c(3.6, 4, 4.4),
  obs.time = "time",
  method = "GH",
  maxiter = 1000,
  n.cv = 3,
  quantile.width = 0.25,
  metrics = c("AUC", "Cindex", "Brier Score", "MAE", "MAEQ")
)

summary(res, metric = "Brier Score")
summary(res, metric = "MAE")
summary(res, metric = "MAEQ")
summary(res, metric = "AUC")
summary(res, metric = "Cindex")
```

Or we can calculate the overall, time-independent Cindex over the entire time period, evaluated by the linear predictor of the (cause-specific) Cox model.
```{r}
Concord <- Concordance(seed = 100, fit, n.cv = 3)
summary(Concord)
```

## Multi-biomarker Joint Model (`mvjmcs`)

To fit a joint model with multiple longitudinal outcomes and competing risks, we can use the `mvjmcs` function. In the example below, we are using the following built-in data sets:

- mvydata: longitudinal data for **multiple** biomarkers per patient
- mvcdata: competing risks time-to-event data per patient

```{r, eval=TRUE}
data(mvydata)
data(mvcdata)
mvfit <- mvjmcs(ydata = mvydata, cdata = mvcdata,
              long.formula = list(Y1 ~ X11 + X12 + time,
                                  Y2 ~ X11 + X12 + time),
              random = list(~ time | ID,
                            ~ 1 | ID),
              surv.formula = Surv(survtime, cmprsk) ~ X21 + X22)
mvfit
```

We can extract the components of the model as follows:

```{r, eval=TRUE}
# Longitudinal fixed effects
fixef(mvfit, process = "Longitudinal")
summary(mvfit, process = "Longitudinal")

# Survival fixed effects
fixef(mvfit, process = "Event")
summary(mvfit, process = "Event")

# Random effects for first few subjects
head(ranef(mvfit))
```

The `FastJM` package can now make dynamic prediction in the presence of multiple longitudinal outcomes. Below is a toy example for competing risks data. Conditional cumulative incidence probabilities for each failure will be presented.

```{r}
require(dplyr)
set.seed(08252025)
sampleID <- sample(mvcdata$ID, 5, replace = FALSE)

subcdata <- mvcdata %>%
  dplyr::filter(ID %in% sampleID)

subydata <- mvydata %>%
  dplyr::filter(ID %in% sampleID)

### Set up a landmark time of 4.75 and make predictions at time u
survmvfit <- survfitJM(mvfit, seed = 100, ynewdata = subydata, cnewdata = subcdata,
                       u = c(7, 8, 9), Last.time = 4.75, obs.time = "time")

survmvfit

```

To assess the prediction accuracy of the fitted joint model, we may run `DynPredAcc` to assess the prediction accuracy by calculating all available evaluation metrics. 
```{r}
res <- DynPredAcc(
  object = mvfit,
  landmark.time = 3,
  horizon.time = c(3.6, 4, 4.4),
  obs.time = "time",
  maxiter = 1000,
  n.cv = 3,
  quantile.width = 0.25,
  metrics = c("AUC", "Cindex", "Brier Score", "MAE", "MAEQ")
)

summary(res, metric = "Brier Score")
summary(res, metric = "MAE")
summary(res, metric = "MAEQ")
summary(res, metric = "AUC")
summary(res, metric = "Cindex")
```

### Landmark Multivariate Joint Model

An alternative approach to characterize flexible latent associations is a landmark multivariate joint model, which specifies a landmark time so that only subjects who remain event-free beyond the landmark time contribute to the model fitting. Here, we consider the current value of the latent process as the association structure in the survival sub-model.
```{r landmark-fit, eval=FALSE}
data(mvydata)
data(mvcdata)
fit.mvlm <- mvjmcs(ydata = mvydata, cdata = mvcdata,
                   long.formula = list(Y1 ~ X11 + X12 + time,
                                       Y2 ~ X11 + X12 + time),
                   random = list(~ time | ID,
                            ~ 1 | ID),
                   surv.formula = Surv(survtime, cmprsk) ~ X21 + X22,
                   control = mvjmcs_control(opt = "optim"),
                   latAsso = "presentlp",
                   landmark = TRUE,
                   s = 4,
                   ytime = "time")
fit.mvlm
```

## Single-biomarker joint model in the presence of heterogeneous within-subject variability (`JMMLSM`)

- ydatah: longitudinal data for a **single** biomarker per patient
- cdatah: competing risks time-to-event data per patient

```{r}
data(ydatah)
data(cdatah)
## fit a joint model
fit <- JMMLSM(cdata = cdatah, ydata = ydatah, 
              long.formula = Y ~ Z1 + Z2 + Z3 + time,
              surv.formula = Surv(survtime, cmprsk) ~ var1 + var2 + var3,
              variance.formula = ~ Z1 + Z2 + Z3 + time, 
              random = ~ 1|ID)
fit
```

```{r}
cnewdata <- cdatah[cdatah$ID %in% c(122, 152), ]
ynewdata <- ydatah[ydatah$ID %in% c(122, 152), ]
survfit <- survfitJM(fit, seed = 100, ynewdata = ynewdata, cnewdata = cnewdata, 
                     u = seq(5.2, 7.2, by = 0.5), Last.time = "survtime",
                     obs.time = "time", method = "GH")
survfit
oldpar <- par(mfrow = c(2, 2), mar = c(5, 4, 4, 4))
plot(survfit, include.y = TRUE)
par(oldpar)
```

To assess the prediction accuracy of the fitted joint model, we may run `DynPredAcc` to assess the prediction accuracy by calculating all available evaluation metrics. 
```{r}
res <- DynPredAcc(
  object = fit,
  landmark.time = 3, horizon.time = c(4:6),
  obs.time = "time",
  method = "GH",
  maxiter = 1000,
  n.cv = 3,
  quantile.width = 0.25,
  metrics = c("AUC", "Cindex", "Brier Score", "MAE", "MAEQ")
)

summary(res, metric = "Brier Score")
summary(res, metric = "MAE")
summary(res, metric = "MAEQ")
summary(res, metric = "AUC")
summary(res, metric = "Cindex")
```
### Simulate Data (Optional)

In order to create simulated data for `mvjmcs`, we can use the `simmvJMdata` function, which creates longitudinal and survival data as a nested list (which are unpacked the this example). When first calling the function, it provides censoring and risk rates. 

```{r}
# Simulate data
  sim <- simmvJMdata(seed = 100, N = 50) # returns list of cdata and ydata for a sample size of 50
  c_data <- sim$mvcdata # survival-side data, one row per ID
  y_data <- sim$mvydata # longitudinal measurements (multiple rows per ID)
```
Below is the simulated longitudinal data for **multiple** biomarkers, wherein Y1 and Y2 represent our biomarkers and X11 and X12 represent measurement-level predictors for the longitudinal submodel. 

```{r}
head(y_data)
```

Below is the simulated survival data wherein X21 and X22 represent patient-level predictors for the survival model.  

```{r}
head(c_data)
```