The Multivariate Linear Mixed Model for Markers

Why Model Several Markers Together?

A participant in a longitudinal study may contribute several outcomes, such as a laboratory marker and a symptom score, measured at different times. The marker panel is not necessarily observed on a common schedule: one marker may be measured at every visit while another is recorded only at selected visits. The earlier material on representing unequally spaced marker panels motivates keeping the marker identity and its actual measurement time explicit.

A multivariate linear mixed model gives each marker its own population trajectory and allows a participant’s departures from those trajectories to be represented jointly. The model is useful for describing the longitudinal panel and its dependence structure. Here the focus is its specification and interpretation, not a survival component or a particular parameterization of the covariance matrices.

Model Setup and Assumptions

Suppose there are \(n\) participants and \(K\) markers. Let \(Y_{ikj}\) be the measurement of marker \(k\) for participant \(i\) at that marker’s \(j\)th observed time, \(t_{ikj}\), where \(j=1,\ldots,n_{ik}\). The number and timing of observations may differ across markers and participants. Let \(x_i\) be a participant-level covariate, such as an indicator for treatment group.

For a linear time trend and a binary covariate, one useful specification is

The \(\beta\) coefficients are fixed effects: they describe the population-average trajectory conditional on the covariate. The \(b_{0ik}\) and \(b_{1ik}\) terms are participant-specific random intercept and slope departures for marker \(k\). The residual \(\varepsilon_{ikj}\) is the remaining measurement-level deviation. The time origin and units must be stated because the intercept is the expected outcome at \(t=0\), and the slope is change per unit of time.

Here \(b_i=(b_{0i1},b_{1i1},\ldots,b_{0iK},b_{1iK})^\mathsf{T}\) is a \(2K\)-dimensional random-effects vector. Its covariance matrix \(G\) describes variation and dependence among participant-specific marker trajectories. The vector \(\varepsilon_i\) stacks the residuals for participant \(i\), and \(R_i\) is its covariance matrix. Both covariance matrices must be symmetric and positive semidefinite; a nonsingular covariance is needed for an ordinary nondegenerate Gaussian density. The specific structures for these matrices should be selected and justified for the study rather than inferred from the mean model.

The Gaussian assumption applies to the random effects and residual vector conditional on the design. It does not require the observed outcomes themselves to be independent: shared random effects induce dependence. The time and covariate design is treated as given, and the model assumes the stated mean and covariance adequately describe the outcomes. In a simple working specification, residuals are independent with marker-specific variances, but that assumption is not automatic and should be checked against the measurement process.

Population Trajectories and Participant Departures

For marker \(k\), the population mean at covariate value \(x_i\) is \(\beta_{0k}+\beta_{1k}t+\beta_{2k}x_i+\beta_{3k}x_it\). Its time slope is \(\beta_{1k}+\beta_{3k}x_i\). A participant’s conditional trajectory adds \(b_{0ik}+b_{1ik}t\), so that participant’s slope is \(\beta_{1k}+\beta_{3k}x_i+b_{1ik}\). The random intercept is the participant’s departure at the chosen time origin; it is not necessarily a meaningful baseline value if \(t=0\) is outside the observed range.

The joint distribution of \(b_i\) permits participant departures for different markers to be related. For instance, a participant’s unusually high intercept for one marker may tend to accompany a high or low intercept for another. Such dependence is a property of the model’s random-effects covariance, not a claim that one marker causes another. Marker-specific fixed effects do not imply that the markers are statistically independent.

Stacked Form and Marginal Mean and Covariance

For compactness, stack all observed outcomes for participant \(i\) into a vector \(Y_i\), in a stated order, such as marker by marker and then time within marker. The fixed-effects design matrix \(X_i\) contains the covariate and time columns for each marker; the random-effects design matrix \(Z_i\) places the appropriate \(1\) and \(t_{ikj}\) entries in the columns for that marker’s random intercept and slope. The model is

This representation accommodates unequal schedules: each participant’s matrices have rows only for the outcomes actually observed, with their actual times. It also makes clear that marker-specific coefficients are placed in separate parts of \(\beta\); using a shared coefficient instead is a substantive restriction, not merely a different data layout.

To justify the result, take expectations and covariances conditional on the design. Since \(E(b_i)=0\) and \(E(\varepsilon_i)=0\), linearity gives \(E(Y_i)=X_i\beta\). For the covariance, subtract the mean and expand:

The two cross-covariance terms are zero because the random effects and residuals are independent. A linear transformation and sum of jointly Gaussian vectors is Gaussian, giving the stated marginal distribution. This derivation also shows why a conditional model with independent residuals can still imply correlated repeated outcomes through \(Z_iGZ_i^\mathsf{T}\).

Design Checks, Alternatives, and Limitations

Before interpreting or estimating a model, check whether its design can support its parameters. Each marker needs enough observations across participants to estimate its fixed time trend; estimating participant-specific slopes also requires useful within-participant time variation. The fixed-effects design must have full column rank: for example, a treatment indicator cannot be separated from an intercept if everyone is in the same group. Random-effects covariance parameters can be weakly identified when there are few participants, sparse repeated measurements, or little variation in times. Centering time at a clinically meaningful point can make intercepts interpretable and may reduce numerical dependence between intercept and slope parameters.

Residual plots and participant-level trajectory plots can reveal curvature, changing spread, outlying trajectories, or systematic departures by marker. A linear time trend is not suitable merely because it is convenient: inspect whether it captures the observed trajectory over the time range. Check whether residual dependence assumptions are plausible for the measurement schedule, and whether estimated covariance matrices are near singular or yield implausible variability. Such warning signs call for revisiting the mean or covariance specification, not simply reporting the fit unchanged.

Separate univariate mixed models are simpler to specify and may be appropriate when cross-marker dependence is not of interest. They do not estimate the joint random-effects distribution across markers. A pooled model that imposes common fixed effects can be more parsimonious, but it asserts that the markers share those effects on their modeled scales. The multivariate model here instead allows marker-specific means and a joint random-effects distribution, at the cost of more parameters and greater data demands. It does not, by itself, correct informative observation times, missingness, nonlinear trajectories, or measurement error beyond the specified model.