#knitr::opts_chunk$set(collapse = TRUE, comment = "#>") library(Matrix) library(nlme) library(survival) library(joineRML) if (requireNamespace('joineR', quietly = TRUE)) { library('joineR') } else { message("'joineR' not available") }
The joineRML
package implements methods for analyzing data from multiple longitudinal studies in which the responses from each subject consists of time-sequences of repeated measurements and a possibly censored time-to-event outcome. The modelling framework for the repeated measurements is the multivariate linear mixed effects model. The model for the time-to-event outcome is a Cox proportional hazards model with log-Gaussian frailty. Stochastic dependence is captured by allowing the Gaussian random effects of the linear model to be correlated with the frailty term of the Cox proportional hazards model. For full details of the model, please consult the technical vignette by running
vignette("technical", package = "joineRML")
The simplest way to explain the concepts of the package is through an example. joineRML
comes with the data set heart.valve
. Details of this data can be found in the help file by running the command
help("heart.valve", package = "joineRML")
This data is in so-called long or unbalanced format:
library("joineRML") data("heart.valve") head(heart.valve)
The data refer to r length(unique(heart.valve$num))
patients and are stored in the unbalanced format, which is convenient here because measurement times were unique to each subject. The data are stored as a single R object, heart.valve
, which is a data frame of dimension r nrow(heart.valve)
by r ncol(heart.valve)
. The average number of repeated measurements per subject is therefore r nrow(heart.valve)
/r length(unique(heart.valve$num))
= r round(nrow(heart.valve) / length(unique(heart.valve$num)), 2)
. As with any unbalanced data set, values of time-constant variables are repeated over all rows that refer to the same subject. The dimensionality of the data set can be confirmed by a call to the dim()
function, whilst the names of the 25 variables can be listed by a call to the names()
function:
dim(heart.valve) names(heart.valve)
We will only analyse a subset of this data, namely records with case-complete data for heart valve gradient (grad
) and left ventricular mass index (lvmi
):
hvd <- heart.valve[!is.na(heart.valve$grad) & !is.na(heart.valve$lvmi), ]
Strictly speaking, this is not necessary because joineRML
can handle the situation of different measurement schedules within subjects That is, a subject does not need to have all multiple longitudinal outcomes recorded at each visit. It is conceivable that some biomarkers will be measured more or less frequently than others. For example, invasive measurements may only be recorded annually, whereas a simple biomarker measurement might be recorded more frequently. joineRML
can handle this situation by specifying each longitudinal outcome its own data frame.
The main function in the joineRML
package is the mjoint()
function. Its main (required) arguments are:
formLongFixed
: a list (of length equal to the number of longitudinal outcome types considered) of two-sided formulae specifying the response on the left-hand side and the mean linear predictor terms for the fixed effects in the linear mixed models on the right-hand side.
formLongRandom
: a list (of same length as formLongFixed
) of one-sided formulae specifying the model for random effects in the linear mixed models.
formSurv
: a formula specifying the proportional hazards regression model for the time-to-event outcome in the same structure as for survival::coxph
.
data
: a list (of same length as formLongFixed
) of data.frames; one for each longitudinal outcome. It is assumed that the event time data is in the first data.frame (i.e. data[[1]]
), unless the argument survData
(which defaults to NULL
) is specified. If $K>1$ and the data are balanced within patients (i.e. multiple markers measured at common measurement times), then one can specify data
as a data frame rather than as a list.
timeVar
: the column name indicating the time variable in the linear mixed effects model. If $K>1$ and the data frames have different column names for time, then timeVar
can alternatively be specified as a vector of strings of length $K$.
We can fit a bivariate joint model to the log-transformed valve gradient and LVMI indices in the hvd
subset using
set.seed(12345) fit <- mjoint( formLongFixed = list("grad" = log.grad ~ time + sex + hs, "lvmi" = log.lvmi ~ time + sex), formLongRandom = list("grad" = ~ 1 | num, "lvmi" = ~ time | num), formSurv = Surv(fuyrs, status) ~ age, data = list(hvd, hvd), inits = list("gamma" = c(0.11, 1.51, 0.80)), timeVar = "time")
Details on the model estimation algorithm are provided in the technical details vignette. We note here that this is not necessarily the most appropriate model for the data, and is included only for the purposes of demonstration. There are a number of other useful arguments in the mjoint
function; for example, inits
for specifying (partial) initial values, control
for controlling the optimization algorithm, and verbose
for monitoring the convergence output in real-time. A full list of all arguments with explanation are given in the help documentation, accessed by running help("mjoint")
.
Once we have a fitted mjoint
object, we can begin to extract relevant information from it. Most summary statistics are available from the summary
function:
summary(fit)
One can also extract the coefficients, fixed effects, and random effects using standard generic functions:
coef(fit) fixef(fit, process = "Longitudinal") fixef(fit, process = "Event") head(ranef(fit))
Although a model fit may indicate convergence, it is generally a good idea to examine the convergence plots. These can be viewed using the plot
function for each group of model parameters.
plot(fit, params = "gamma") plot(fit, params = "beta")
Once an mjoint
model has converged, and assuming the pfs
argument is TRUE
(default), then approximated standard errors are calculated based on the empirical information matrix of the profile likelihood at the maximizer. Theoretically, these standard errors will be underestimated (see the technical vignette). In principle, residual Monte Carlo error will oppose this through an increase in uncertainty.
fit.se <- bootSE(fit, nboot = 100)
Bootstrapping is a computationally intensive method, possibly taking many hours to fit. For this reason, one can relax the control parameter constraints on the optimization algorithm for each bootstrap model; however, this will be at the possible expense of inflated standard errors due to Monte Carlo error.
We can call the bootSE
object to interrogate it
fit.se
or alternatively re-run the summary
command, passing the additional argument of bootSE = fit.se
summary(fit, bootSE = fit.se)
joineRML
versus joineR
There are a growing number of software options for fitting joint models of a single longitudinal outcome and a single time-to-event outcome; what we call here univariate joint models. joineR
(version 1.2.7) is one package available in R for fitting such models, however joineRML
can fit these models too, since the univariate model is simply a special case of the multivariate model. It is useful to contrast these two packages. There are theoretical and practical implementation differences between the packages beyond just univariate versus multivariate capability:
joineR
uses Gauss-Hermite quadrature (with 3 nodes) for numerical integration, whereas joineRML
uses Monte Carlo integration with an automated selection of sample size.
joineR
only allows for random-intercept models, random-intercept and random-slope models, or a quadratic model. joineRML
, on the other hand, allows for any random effects structure.
joineR
only allows for specification of convergence based on an absolute difference criterion.
joineR
does not calculate approximate standard errors, and instead requires a bootstrap approach be used after the model fit.
The current version of joineR
requires a data pre-processing step in order to generate a joint.data
object, whereas joineRML
can work straight from the data frame.
To fit a univariate model in joineR
we run the following code for the hvd
data
joineR_available <- require(joineR)
library(joineR, quietly = TRUE) hvd.surv <- UniqueVariables(hvd, var.col = c("fuyrs", "status"), id.col = "num") hvd.cov <- UniqueVariables(hvd, "age", id.col = "num") hvd.long <- hvd[, c("num", "time", "log.lvmi")] hvd.jd <- jointdata(longitudinal = hvd.long, baseline = hvd.cov, survival = hvd.surv, id.col = "num", time.col = "time") fit.joiner <- joint(data = hvd.jd, long.formula = log.lvmi ~ time + age, surv.formula = Surv(fuyrs, status) ~ age, model = "intslope") summary(fit.joiner)
To fit a univariate model in joineRML
we run the following code for the hvd
data
set.seed(123) fit.joinerml <- mjoint(formLongFixed = log.lvmi ~ time + age, formLongRandom = ~ time | num, formSurv = Surv(fuyrs, status) ~ age, data = hvd, timeVar = "time") summary(fit.joinerml)
In addition to just comparing model parameter estimates, we can also extract the predicted (or posterior) random effects from each model and plot them.
id <- as.numeric(row.names(fit.joiner$coefficients$random)) id.ord <- order(id) # joineR rearranges patient ordering during EM fit par(mfrow = c(1, 2)) plot(fit.joiner$coefficients$random[id.ord, 1], ranef(fit.joinerml)[, 1], main = "Predicted random intercepts", xlab = "joineR", ylab = "joineRML") grid() abline(a = 0, b = 1, col = 2, lty = "dashed") plot(fit.joiner$coefficients$random[id.ord, 2], ranef(fit.joinerml)[, 2], main = "Predicted random slopes", xlab = "joineR", ylab = "joineRML") grid() abline(a = 0, b = 1, col = 2, lty = "dashed")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.