Fitting a Basic Shared-Random-Effects Model in R

From Prepared Data to a Fitted Joint Model

Once the longitudinal measurements and event records are organized on a common time scale, the next task is to fit a model that links them. This tutorial uses a common starting point: a continuous Gaussian marker measured repeatedly, subject-specific random intercepts and slopes, and a time-to-event outcome. In R, the JM package provides one route to fitting this model by combining a longitudinal mixed-effects fit with a survival fit.

The key distinction is between fitting the two submodels separately and fitting them jointly. Separate fits can be useful for checking the inputs, but they do not estimate the association between the marker trajectory and event risk. A joint fit estimates the model components together, accounting for the fact that the same subject-specific random effects contribute to both the repeated measurements and the event process.

Model and Assumptions

Definition: Let $Y_{ij}$ be the Gaussian marker for subject $i$ at measurement time $t_{ij}$, and let $T_i$ denote the event or censoring time with event indicator $\delta_i$. A random-intercept-and-slope longitudinal model is $$ Y_{ij} = \beta_0 + \beta_1 t_{ij} + \boldsymbol{\beta}_x^\top \boldsymbol{x}_i + b_{0i} + b_{1i}t_{ij} + \epsilon_{ij}, $$ where $\boldsymbol{x}_i$ contains fixed covariates, $(b_{0i},b_{1i})^\top$ are subject-specific random effects, and $\epsilon_{ij}$ is measurement error. In a current-value shared-random-effects model, the event hazard is $$ h_i(t\mid b_i) = h_0(t)\exp\left\{\boldsymbol{\gamma}^\top \boldsymbol{w}_i + \alpha m_i(t)\right\}, \qquad m_i(t)=\beta_0+\beta_1t+\boldsymbol{\beta}_x^\top\boldsymbol{x}_i+b_{0i}+b_{1i}t. $$ The shared random effects enter the hazard through the subject’s latent marker trajectory $m_i(t)$. The parameter $\alpha$ describes the association between current latent marker level and hazard, conditional on the included covariates and random effects.

This formulation relies on substantive and distributional assumptions. The random effects are usually taken to follow a multivariate normal distribution with mean zero and an estimated covariance matrix. Conditional on those effects and the fixed effects, the longitudinal errors are Gaussian, have mean zero, and are independent across measurements with a common variance in this basic model. The event submodel assumes proportional hazards conditional on the random effects and covariates, with a chosen baseline hazard. The longitudinal and event processes are conditionally independent given the shared random effects and modeled covariates. The event indicator must correctly distinguish observed events from censoring, and each subject’s measurement and event times must use compatible time origins and units.

These assumptions are not guaranteed by choosing software options. They are modeling choices to assess using subject-level data checks, plots, scientific knowledge, and fit diagnostics. A random slope, for example, is most meaningful when subjects have measurements at multiple times and the time range supports estimating between-subject slope variation.

Prepare and Match the Two Data Sets

The longitudinal data ordinarily have one row per measurement, with an identifier, measurement time, marker value, and any longitudinal covariates. The event data ordinarily have one row per subject, with the same identifier, follow-up time, event indicator, and baseline covariates. The identifier is the link between the tables; rows should not be matched merely by their order.

Worked Example: Checking Subject-Level Alignment

Consider a hypothetical three-subject cohort. The longitudinal table contains these records: subject A has marker values 4.2 at time 0 and 5.1 at time 1; subject B has 3.8 at time 0 and 4.0 at time 1; subject C has 5.0 at time 0 and 5.4 at time 1. The event table records A with follow-up time 1.6 and event indicator 1, B with time 2.0 and indicator 0, and C with time 1.2 and indicator 1.

There are six measurement rows but three event rows, as expected for these table layouts. A’s event indicator of 1 means the event was observed at 1.6 time units; B’s 0 means follow-up ended at 2.0 without an observed event. These event records do not change the measurement times: for example, A’s measurement at time 1 occurred before the event at time 1.6. Before fitting, verify that each identifier in the measurement table has exactly one corresponding event record, that there are no duplicated subject rows in the event table, and that event times and measurement times refer to the same origin and units.

This tiny table illustrates structure, not a suitable sample for estimating a model with several variance and association parameters. In an actual analysis, the number and timing of subjects’ measurements and events must support the model being fit.

The JM package accepts a fitted longitudinal mixed model and a fitted Cox model as inputs. The Cox fit must retain its model matrix, which is why the example below uses x = TRUE. For this basic specification, the survival data provide one record per subject and the longitudinal data provide repeated records.

Fit the Basic Model in R

Worked Example: Fitting with JM

Suppose a hypothetical study has data frames named long_dat and event_dat. The former has columns id, time, marker, and trt; the latter has id, event_time, event, and trt. Here event is 1 for an observed event and 0 for right censoring, and trt is a categorical treatment indicator. Assuming the IDs and time variables have passed the checks above, the following code fits the model:

R code:
library(nlme)
library(survival)
library(JM)

lme_fit <- lme(fixed = marker ~ time + trt,
    random = ~ time | id,
    data = long_dat,
    method = "ML")

cox_fit <- coxph(Surv(event_time, event) ~ trt,
    data = event_dat,
    x = TRUE)

joint_fit <- jointModel(lme_fit, cox_fit,
    timeVar = "time",
    method = "Cox-PH-aGH",
    parameterization = "value")

summary(joint_fit)

The mixed model specifies fixed effects for time and treatment and subject-specific random intercepts and slopes for time. The Cox model supplies the survival outcome and its treatment covariate. In the joint fit, timeVar = "time" identifies the longitudinal time variable used to connect the marker trajectory to event follow-up; parameterization = "value" specifies the current-value association. The selected Cox baseline-hazard method and the package version should be recorded with the analysis.

The two component fits are inputs to the joint fitting procedure, not the final joint estimates. The joint model estimates the longitudinal and survival parameters in the linked model. In particular, the association estimate is not obtained by simply inserting predictions from the separate mixed model into a Cox regression.

Inspect Estimates in Context

Start by checking that the fitted object corresponds to the intended data and specification. Inspect the summary for the fixed longitudinal effects, random-effects and residual variation, survival covariate effects, and marker–hazard association. Confirm that the time variable, event coding, covariates, and association structure are what the analysis plan calls for. Check that estimates and standard errors are finite, and consider whether the estimated random-effects variability is plausible given the observed trajectories.

Worked Example: Interpreting an Illustrative Association Estimate

Suppose, purely for illustration, a joint-model summary reports a current-value association estimate of $\hat{\alpha}=0.35$ with standard error $0.14$. If the marker is measured in units of 10 points, then, conditional on the model’s covariates and random effects, a one-unit higher latent marker value corresponds to a hazard ratio of

\exp(\hat{\alpha})=\exp(0.35)=1.419.

Thus the fitted model estimates a hazard about 1.42 times as high for a one-unit higher latent marker value at the same time, holding the included covariates fixed. This is an estimated conditional association, not evidence that increasing the marker causes an increase in hazard.

A rough Wald calculation gives $z=0.35/0.14=2.50$, with a two-sided normal-approximation p-value of about 0.0124. A corresponding approximate 95% interval for $\alpha$ is

0.35 \mathbin{\pm} 1.96(0.14) = [0.0756,\ 0.6244].

Exponentiating the endpoints gives an approximate hazard-ratio interval of $[\exp(0.0756),\exp(0.6244)]=[1.079,\ 1.867]$ per 10-point marker unit. Under the model and the Wald approximation, this interval is compatible with hazard ratios from about 1.08 to 1.87. It is not a probability statement about where the fixed parameter lies. The reported estimate, standard error, and interval here are hypothetical values used to demonstrate interpretation; they are not claimed as output from the code above.

The fixed treatment effect in the survival submodel has a different interpretation from $\alpha$. If its estimate is $\hat{\gamma}_{\mathrm{trt}}=-0.40$, the fitted conditional hazard ratio comparing the coded treatment group with the reference group is $\exp(-0.40)=0.670$, at the same latent marker value. The reference category and covariate coding must be checked before describing that comparison.

Justification for the Shared-Effects Interpretation

The software’s current-value association links the hazard to the latent trajectory, rather than directly to a noisy observed measurement. The shared-effects interpretation follows by substituting the trajectory into the hazard model:

\begin{aligned} h_i(t\mid b_i) &=h_0(t)\exp\{\boldsymbol{\gamma}^\top\boldsymbol{w}_i+\alpha m_i(t)\}\\ &=h_0(t)\exp\{\boldsymbol{\gamma}^\top\boldsymbol{w}_i+ \alpha(\beta_0+\beta_1t+\boldsymbol{\beta}_x^\top\boldsymbol{x}_i +b_{0i}+b_{1i}t)\}. \end{aligned}

Both random effects, $b_{0i}$ and $b_{1i}$, now appear in the hazard’s linear predictor. A subject’s latent intercept and slope therefore influence both the distribution of repeated marker values and the event hazard. The association is through the current latent value, so this specification should not be described as a separate direct effect of the slope.

Checks, Alternatives, and Limitations

Basic checks can catch problems before estimates are given a scientific interpretation. Confirm that the measurement table has one row per actual measurement, the event table has one row per subject, IDs match, and event coding agrees with the survival definition. Plot individual trajectories and look for patterns inconsistent with a roughly linear mean trajectory or common Gaussian residual variation. Inspect whether the estimated random-intercept-and-slope covariance is plausible and whether the association estimate has an interpretable scale. A warning, extreme estimate, or unusually large uncertainty is a reason to investigate, not a reason to report the output uncritically.

The main alternative at this stage is to fit a longitudinal model and survival model separately. Those fits may be useful for descriptive checks, but they do not account for the shared latent effects when estimating event association. Landmarking is another approach to prediction, with a different conditioning strategy; it does not fit this same shared-effects likelihood. The joint model is attractive when the longitudinal trajectory and event process are linked through latent subject heterogeneity, but it relies on its assumptions and requires enough information to estimate that link.

A linear time trend, Gaussian errors, a random intercept and slope, and a current-value association can all be inadequate. Nonlinear trajectories, non-Gaussian markers, different association structures, and alternative baseline hazards require deliberate model choices; the basic code is not a general solution to them. Subjects with few measurements, limited follow-up, or few events may provide weak information about some parameters. Model fit, sensitivity to assumptions, and the intended prediction or association target should guide whether this basic model is appropriate.

Common mistake: Treating the two separate component fits as if they were already a joint model. The correction is to fit the linking model and interpret its joint estimates, not to infer a marker–hazard association from unrelated summaries.
Common mistake: Interpreting $\exp(\hat{\alpha})$ as the effect of a one-unit increase in an observed noisy reading. It is the conditional hazard ratio for a one-unit increase in the latent current marker value, under the specified model.
Common mistake: Assuming matching row order proves the data sets align. Match by subject ID and verify one event record per subject.

Key takeaway: A basic R joint-model fit combines repeated Gaussian measurements with a subject-level survival record and links the hazard to the latent marker trajectory. Check the data alignment and model specification before interpreting the fitted estimates.

Check Your Understanding

Answer each question, then compare with the concise solutions below.

Why does the longitudinal data set usually have more rows than the event data set?
In the model shown, what does $\alpha$ describe?
If $\hat{\alpha}=0.20$ per marker unit, what is the estimated hazard ratio for a one-unit higher latent marker value?
Why should IDs be matched explicitly rather than assuming the tables have the same row order?
Name one assumption of the basic Gaussian random-intercept-and-slope model that should be assessed.

Solutions. (1) Each subject can contribute multiple measurement rows but typically contributes one event or censoring record. (2) It is the conditional association between current latent marker value and event hazard on the log-hazard scale. (3) $\exp(0.20)=1.221$, an estimated hazard about 1.22 times as high per one-unit increase, conditional on the model. (4) Row order does not establish subject identity; incorrect matches link measurements to the wrong event records. (5) Examples include Gaussian random effects, Gaussian conditionally independent measurement errors, or a linear mean trajectory; each is a modeling assumption, not guaranteed by the fitting function.