Mixed Models For Repeated Omics Measurements

Why Use a Mixed Model for Repeated Omics Measurements?

Suppose an assay measures the same molecular feature in each participant at several visits. Measurements from one participant are often more alike than measurements from different participants: they share biology, sample handling, and other subject-specific influences. Treating every measurement as independent ignores that structure and can understate uncertainty. A linear mixed model represents both the average feature trajectory and the clustering of observations within subjects.

The focus here is one approximately continuous omics feature, such as a normalized expression value or log-intensity. The model can describe how its mean changes with time and how that change differs by a group. It is not a model for several features jointly, and its use does not by itself establish a causal effect of group or time.

Model, Parameters, and Assumptions

Let \(y_{ij}\) denote the measured value of one feature for subject \(i\) at visit \(j\), and let \(t_{ij}\) be time since a chosen reference visit. Let \(g_i\) equal 0 for the reference group and 1 for the other group. A linear random-intercept and random-slope model is

y_{ij}=\beta_0+\beta_1t_{ij}+\beta_2g_i+\beta_3t_{ij}g_i+b_{0i}+b_{1i}t_{ij}+\epsilon_{ij}.

The \(\beta\) parameters are fixed effects: they describe the average trajectory across the population. The subject-specific random intercept \(b_{0i}\) shifts subject \(i\)'s trajectory, while the random slope \(b_{1i}\) allows that subject's rate of change to differ from the population slope. The residual \(\epsilon_{ij}\) is the remaining visit-specific deviation.

Definition: In a linear mixed model, the response is written as a fixed-effects mean plus random effects and residual error. Here the fixed-effects mean is \(\beta_0+\beta_1t_{ij}+\beta_2g_i+\beta_3t_{ij}g_i\). The random effects \((b_{0i},b_{1i})\) describe subject-level deviations from that mean, and \(\epsilon_{ij}\) describes the remaining measurement-level deviation.

A standard version of this model assumes that, across subjects, the random-effects vectors are independent and identically distributed with mean zero and covariance matrix \(G\):

\begin{pmatrix}b_{0i}\\b_{1i}\end{pmatrix} \sim N\left( \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}\tau_0^2&\tau_{01}\\\tau_{01}&\tau_1^2\end{pmatrix} \right).

It also assumes that residuals have mean zero, variance \(\sigma^2\), and are independent across visits conditional on the random effects; the residuals are independent of the random effects and of other subjects. For the usual Gaussian likelihood, random effects and residuals are normally distributed. The covariance matrix \(G\) must be positive semidefinite. The fixed-effects mean should be appropriately specified for the time range and design; the straight-line time terms above assume a linear mean trajectory. These assumptions concern the model, not a guarantee that every observed trajectory is linear or normally shaped.

Subjects are the independent clusters in this setup. The model presumes that different subjects' random effects and residual errors are independent. It also presumes that the measurement process and missing visits do not create unmodelled dependence; in particular, likelihood-based analysis with incomplete outcomes needs an appropriate missingness assumption, commonly missing at random conditional on observed data and included covariates. The model does not automatically correct informative dropout.

Interpreting Fixed Effects and Predictions

At \(t=0\), the reference group's population mean is \(\beta_0\), and its change per unit time is \(\beta_1\). The other group's population mean at \(t=0\) is \(\beta_0+\beta_2\), and its time slope is \(\beta_1+\beta_3\). Thus \(\beta_2\) is the between-group mean difference at the reference time, while \(\beta_3\) is the difference in mean time slopes. These interpretations depend on the coding of group and the chosen time origin. Since the random effects have mean zero, averaging over subjects gives the fixed-effects mean.

Worked Example: Comparing Group Trajectories

Hypothetical data and model. Suppose a normalized feature is measured in two groups over time, with \(t\) in years and \(y\) in log2 assay units. A fitted model gives \(\beta_0=5.00\), \(\beta_1=0.20\), \(\beta_2=0.30\), and \(\beta_3=-0.10\). Interpret the coefficients and calculate the population-average mean for group 1 at \(t=2\).

For group 0, the mean at \(t=0\) is \(5.00\) log2 units and the slope is \(0.20\) log2 units per year. For group 1, the mean at \(t=0\) is \(5.00+0.30=5.30\) log2 units, and its slope is \(0.20+(-0.10)=0.10\) log2 units per year. The groups' mean slopes differ by \(-0.10\) log2 units per year.

At \(t=2\), the group 1 population-average mean is

5.00+(0.20)(2)+0.30+(-0.10)(2) =5.00+0.40+0.30-0.20 =5.50\text{ log2 units}.

The corresponding group 0 mean is \(5.00+(0.20)(2)=5.40\) log2 units, so the group 1 minus group 0 contrast at two years is \(5.50-5.40=0.10\) log2 units. This is an estimated difference in mean feature values, not evidence that group membership caused the difference.

How the Model Represents Within-Subject Correlation

Write \(z(t)=(1,t)^\mathsf{T}\) for the random-effects design vector. For two different visits \(j\) and \(k\) on the same subject, the shared random effects induce covariance. Under the stated residual independence assumption,

\operatorname{Cov}(y_{ij},y_{ik}) =z(t_{ij})^\mathsf{T}Gz(t_{ik}) =\tau_0^2+\tau_{01}(t_{ij}+t_{ik})+\tau_1^2t_{ij}t_{ik}, \qquad j\ne k.

At one visit, the residual variance also contributes, so \(\operatorname{Var}(y_{ij})=z(t_{ij})^\mathsf{T}Gz(t_{ij})+\sigma^2\). This structure allows repeated measurements to be correlated without treating them as independent replicates. It also shows that correlation need not be constant across pairs of visits when random slopes are included.

Justification: Write \(y_{ij}=x_{ij}^\mathsf{T}\beta+z(t_{ij})^\mathsf{T}b_i+\epsilon_{ij}\), where \(x_{ij}\) contains the fixed-effects terms. For \(j\ne k\), the fixed part is nonrandom, the residuals are independent, and residuals are independent of \(b_i\). Therefore, \(\operatorname{Cov}(y_{ij},y_{ik})=\operatorname{Cov}(z(t_{ij})^\mathsf{T}b_i,z(t_{ik})^\mathsf{T}b_i)=z(t_{ij})^\mathsf{T}Gz(t_{ik})\). Expanding this matrix product gives the covariance formula above. For \(j=k\), add \(\operatorname{Var}(\epsilon_{ij})=\sigma^2\).

Worked Example: Calculating Within-Subject Covariance

Hypothetical fitted variance components. Suppose \(\tau_0^2=0.36\), \(\tau_1^2=0.04\), \(\tau_{01}=-0.03\), and \(\sigma^2=0.25\), where the variances are in squared log2 units and time is in years. Calculate the covariance between measurements at years 0 and 2, and their correlation.

First, the random-effects covariance matrix is admissible: its determinant is \((0.36)(0.04)-(-0.03)^2=0.0135>0\), and its diagonal entries are positive. The covariance between distinct visits at 0 and 2 years is

0.36+(-0.03)(0+2)+(0.04)(0)(2) =0.36-0.06+0 =0.30\text{ squared log2 units}.

At year 0, the total variance is \(0.36+0.25=0.61\) squared log2 units. At year 2, it is \(0.36+2(-0.03)(2)+(0.04)(2^2)+0.25=0.36-0.12+0.16+0.25=0.65\) squared log2 units. Hence the correlation is

\frac{0.30}{\sqrt{(0.61)(0.65)}}= \frac{0.30}{0.6297}\approx0.476.

The two measurements are positively correlated under this fitted model. The residual variance contributes to each measurement's variance, but not to covariance across distinct visits under the independent-residual assumption.

Subject-Specific Predictions and Estimation

A population-average prediction sets the random effects to their mean, zero. A subject-specific fitted trajectory instead adds estimates of that subject's random effects. These estimates are often called empirical Bayes estimates or best linear unbiased predictions (BLUPs); they are shrunk toward zero because the subject's data are noisy and the random-effects distribution provides population-level information.

Worked Example: Shrinkage of a Subject Intercept

Hypothetical residuals. Consider a random-intercept-only model with intercept variance \(\tau_0^2=0.36\) and residual variance \(\sigma^2=0.25\), both in squared log2 units. For one subject, suppose two measurements have residuals \(+0.40\) and \(-0.20\) log2 units after subtracting the fixed-effects means. Calculate the estimated random intercept.

For this balanced random-intercept case, the conditional estimate is \(\hat b_{0i}=\frac{\tau_0^2}{\sigma^2+m\tau_0^2}\sum_j r_{ij}\), where \(m=2\) and \(r_{ij}\) are the fixed-mean residuals. Substitution gives

\hat b_{0i} =\frac{0.36}{0.25+2(0.36)}(0.40-0.20) =\frac{0.36}{0.97}(0.20) \approx0.0742\text{ log2 units}.

The subject's average residual is \(0.10\) log2 units, but the estimated random intercept is smaller because it is shrunk toward the population value zero. This conditional estimate is not a new fixed group effect.

Fixed effects and variance components are commonly estimated by maximum likelihood (ML) or restricted maximum likelihood (REML). REML accounts for the degrees of freedom used to estimate fixed effects when estimating variance components; ML is often preferred when comparing models with different fixed-effects structures using likelihood-based criteria. In either case, estimation relies on the specified mean and covariance structure. A more elaborate random-effects structure is not automatically better, especially with few subjects or few visits.

Diagnostics, Alternatives, and Limitations

Diagnostics should assess whether the fitted structure is plausible, not merely whether the software converged. Plot residuals against fitted values and time to look for changing variance, systematic curvature, or patterns left in the residuals. Inspect subject-specific trajectories and estimated random effects for unusual influential subjects. Check whether the estimated covariance matrix is valid and whether variance estimates are near zero or correlations near their boundaries. Such boundary estimates can signal that the data do not support the chosen random-effects complexity.

A linear time effect may be inadequate if trajectories bend; systematic residual patterns over time are a warning. Random slopes may be weakly identified with few visits per subject, and a complicated random-effects covariance may be unstable with a small number of subjects. Normality of residuals and random effects matters most for small-sample likelihood-based inference and prediction, although severe outliers or heavy tails can also distort estimates. Use scientific knowledge and plots alongside formal checks; do not select a structure solely because a test is nonsignificant.

Approach	Useful when	Trade-off
Linear mixed model	Subject-specific trajectories and unbalanced visit times are of interest.	Requires a credible random-effects and residual covariance model.
Repeated-measures ANOVA	Visits are common, scheduled, and the response structure is relatively simple.	Can be restrictive for irregular timing, incomplete visits, or individual slopes.
Generalized estimating equations	The target is a population-average association and robust cluster-level inference is desired.	Does not provide the same subject-specific random-effect predictions; robust variance estimation needs enough independent subjects.

A mixed model's fixed effects describe the mean conditional on random effects and, because those effects have mean zero, also give the model's population-average mean. That does not make mixed models and generalized estimating equations interchangeable: their assumptions, estimation, and handling of covariance differ. The model also does not solve confounding, assay batch effects, or informative visit and dropout processes unless these are appropriately represented. In particular, a single-feature mixed model is not a substitute for modelling several omics layers jointly.

Common Mistakes

Counting visits as independent subjects. Repeated measurements from one person share random effects. The number of independent clusters is the number of subjects, not the total number of rows.
Calling a random effect a fixed group effect. A random intercept is a subject-level deviation modeled as a distribution; a group coefficient is a systematic mean contrast defined by the fixed-effects design.
Interpreting \(\beta_2\) without its time origin. With a time-by-group interaction, \(\beta_2\) is the group contrast at \(t=0\), not at every time.
Assuming a random slope fixes any trajectory problem. Random slopes represent between-subject variation in linear rates; they do not make a curved population mean linear model adequate.
Assuming missing visits are harmless. Likelihood-based analysis does not automatically account for dropout related to unobserved feature values.

Key takeaway: Specify the population mean and subject-level covariance separately. Fixed effects describe the average trajectory; random effects induce within-subject dependence and represent individual deviations. Both components need to be scientifically plausible and checked against the observed data.

Check Your Understanding

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

In the stated model, what does \(\beta_3\) measure, and in what units?
If two distinct visits have residuals independent conditional on random effects, why can their outcomes still be correlated?
For a random-intercept-only model with \(\tau_0^2=0.25\), \(\sigma^2=0.25\), and two visits, what is the covariance between the visits?
What is the population-average mean when \(t=1\), \(g=1\), and the coefficients are \(5.00, 0.20, 0.30,-0.10\), respectively?
Name one diagnostic sign that a random-slope model may be inadequate or too complex.

Solutions. (1) \(\beta_3\) is the group difference in mean time slope, in log2 feature units per unit time. (2) Both outcomes contain the same subject's random effects, so their covariance includes \(z(t_{ij})^\mathsf{T}Gz(t_{ik})\). (3) Their covariance is the shared random-intercept variance, \(0.25\) squared response units. (4) \(5.00+(0.20)(1)+0.30+(-0.10)(1)=5.40\) log2 units. (5) A near-zero slope variance or a boundary random-effects correlation can indicate that the data do not support the added complexity; systematic residual curvature instead suggests the linear mean may be inadequate.