Why Repeated Measurements Need Their Own Analysis
During treatment, a molecular feature can change within a patient, and patients can also differ in their typical molecular levels. Those are different sources of variation. For example, patients with greater average clinical improvement might have different average protein abundance, even if a patient's protein abundance does not track that patient's clinical changes over time. Treating all visits as independent observations blurs these patterns and understates uncertainty because measurements from the same patient are generally correlated.
A useful analysis therefore identifies the patient as the unit that links repeated observations, represents treatment time, and separates patient-level differences from time-varying associations. The methods below concern one molecular feature at a time; applying them across many features requires attention to multiplicity and validation.
A Model That Separates Within- and Between-Patient Associations
For a continuous molecular abundance outcome, a linear mixed model can include both response components and treatment time:
Here \(\beta_W\) is the time-adjusted within-patient association: the expected difference in molecular abundance for a one-unit increase in response above a patient's own average, holding the modeled time trend and other terms fixed. \(\beta_B\) is the between-patient association: the expected abundance difference associated with a one-unit difference in patients' average response. The patient random intercept \(b_i\) captures residual patient-level differences not explained by the fixed effects. These coefficients are associations, not causal effects.
A standard version assumes \(b_i\sim N(0,\tau^2)\), independent across patients and independent of the included predictors; conditional on \(b_i\) and the predictors, \(\epsilon_{ij}\) has mean zero, variance \(\sigma^2\), and is independent across visits and patients. Normality is commonly used for likelihood-based inference in small samples; residual plots and sensitivity analyses help assess whether it is reasonable. The linear specification also assumes that the mean time trend and response associations are adequately represented by the chosen terms. If the trajectory is curved, time can instead be represented by a suitable nonlinear function or prespecified time categories.
Including \(\bar{R}_i\) alongside the patient-centered response deviation is called a within-between decomposition. It prevents the within-patient coefficient from being silently interpreted as a comparison between different patients. When repeated values are irregularly spaced, use actual visit times rather than visit number if elapsed time is scientifically relevant.
Worked Example: Separating Patient-Level and Within-Patient Patterns
Suppose two hypothetical patients have response scores and molecular abundances at three visits. The response score is higher when clinical status is better. For Patient A, responses are 2, 4, and 6, and abundances are 5, 7, and 9 units. For Patient B, responses are 6, 8, and 10, and abundances are 12, 14, and 16 units.
For A, \(\bar{R}_A=(2+4+6)/3=4\) and \(\bar{Y}_A=(5+7+9)/3=7\). The response deviations are \(-2,0,2\), and abundance deviations are \(-2,0,2\). The within-patient slope for these centered observations is
Patient B has the same centered values, so the pooled within-patient slope is also 1 abundance unit per response-score unit. Between patients, the mean response differs by \(8-4=4\) score units and the mean abundance differs by \(14-7=7\) units, giving a two-point between-patient slope of \(7/4=1.75\) abundance units per response-score unit. The two slopes differ: the within-patient pattern is not the same as the between-patient comparison. This tiny illustration has no uncertainty estimate and no time adjustment; moreover, its response values rise in lockstep with visit order, so it cannot distinguish response association from a time trend.
Time Adjustment and Trajectory Comparisons
A response measure often changes during treatment, as do many molecular features. If both rise or fall over time, an apparent association can result from their shared time trend. Adjusting for time helps separate that trend from the response association, but it does not guarantee that all time-varying confounding has been removed. The time function should reflect the study schedule and scientific question; with few visits, prespecified visit indicators may be safer than an elaborate curve.
A related question is whether a molecular feature follows different trajectories in patients classified by a baseline clinical response group. For a binary baseline group \(G_i\), a model may include a time-by-group interaction:
With a linear time term, the reference group's expected slope is \(\beta_t\), and the other group's slope is \(\beta_t+\beta_{tG}\). Thus \(\beta_{tG}\) is the difference in slopes. A baseline group is patient-level and does not vary from visit to visit. It must not be confused with a time-varying response score, which can be decomposed into within- and between-patient parts as above.
Worked Example: A Time-Adjusted Within-Patient Association
Consider hypothetical data from two patients measured at days 0, 1, and 2. Let \(x_{ij}=R_{ij}-\bar{R}_i\). Patient A has \(x=(-1,0,1)\) and molecular abundance \(Y=(7,12,17)\) units. Patient B has \(x=(1,0,-1)\) and \(Y=(23,22,21)\) units. Fit a model with patient intercepts, a common linear time effect, and a within-patient response effect. In this constructed example, the patient and time terms are nuisance terms to be removed before estimating the response slope.
The patient means are 12 units for A and 22 units for B. The pooled time means are 15, 17, and 19 units at days 0, 1, and 2. Subtracting the relevant patient mean and time deviation from each observation leaves residualized abundances \(Y^* = (-3,0,3)\) for A and \((3,0,-3)\) for B. For example, A's day-0 residualized value is \(7-12-(15-17)=-3\). The centered response values are orthogonal to patient indicators and to the common time trend: each patient's \(x\) values sum to zero, and the cross-patient sums at each visit are also zero.
The time-adjusted within-patient slope is therefore
In these hypothetical data, a one-point increase above a patient's mean response is associated with 3 additional abundance units, after removing the modeled patient intercept and linear time trend. The exact fit is intentionally simple and supplies no useful estimate of population uncertainty. In an actual study, the coefficient and its interval or test would be obtained from the mixed-model fit, with uncertainty reflecting the number of patients and the repeated-measure structure.
Worked Analysis Choices
Worked Example: Comparing Baseline Response-Group Slopes
Suppose a hypothetical study measures one protein in two baseline groups at days 0 and 14. The mean abundance is 10 and 7 units for the responder group, and 10 and 9 units for the non-responder group. Use the two-visit linear change to compare descriptive slopes.
The responder-group slope is \((7-10)/(14-0)=-3/14=-0.2143\) abundance units per day. The non-responder-group slope is \((9-10)/14=-1/14=-0.07143\) abundance units per day. Their slope difference, responder minus non-responder, is
The group means show a faster decline in this protein in the responder group over the observed 14 days. This is a descriptive calculation, not a significance test: the group sizes, patient-level variability, and within-patient correlation are not provided. A mixed model with patient-level repeated observations and a time-by-group interaction would estimate the slope difference and its uncertainty. Baseline response-group membership is not randomized, so the contrast does not establish that response status caused the different trajectory.
The repeated-measure model is not the only possible approach. The choice depends on the estimand, number of patients, visit schedule, and whether the goal is estimation, prediction, or description.
| Approach | Useful when | Main caution |
|---|---|---|
| Mixed model with patient random effects | Repeated observations are unbalanced or visit times differ; patient-level dependence must be represented. | Results depend on model structure and assumptions, including the random-effects specification. |
| Patient-level summaries or change scores | A small number of prespecified visits makes a simple patient-level contrast scientifically meaningful. | Summarizing can discard timing information and may be inefficient or sensitive to the chosen visits. |
| Population-average model with a working correlation | The target is an average association across the population rather than a patient-specific random effect. | With few independent patients, robust standard errors can be unreliable; correlation structure still matters for efficiency. |
Diagnostics, Assumptions, and Failure Modes
A fitted model is not validated merely because software returns coefficients. Check the data and model at both the patient and visit levels. Plot each patient's molecular measurements against actual time, with clinical response marked or shown in a companion panel. Look for nonlinear trajectories, outlying patients, abrupt changes at assay batches, and patterns that are hidden by a cohort-average curve. Examine residuals against fitted values and time for nonconstant variance, systematic curvature, or remaining trends. Check whether observed residual patterns suggest dependence across visits beyond the random intercept; if so, a richer correlation structure may be needed, provided the number of patients supports it.
Check that response and time are not so strongly aligned that their separate effects are poorly identified. If every patient's response improves monotonically at exactly the same rate as time, the data may not contain enough independent variation to estimate distinct time and response effects. Also report how many patients and visits contribute to each analysis, the amount and pattern of missingness, and whether missing visits are plausibly explained by observed information. Likelihood-based mixed models can use incomplete repeated outcomes under a missing-at-random assumption conditional on included observed data and the specified model; they do not automatically correct informative dropout related to unobserved outcomes.
Assay batch and treatment time can be confounded. If samples from early visits were assayed in one batch and late visits in another, an apparent trajectory may be a technical shift. Include known batch terms when estimable, inspect quality-control information, and be explicit when design confounding prevents separation. A random intercept does not fix batch confounding, an incorrectly specified time trend, or unmeasured time-varying clinical factors.
- Wrong: “There are 60 observations, so the analysis has 60 independent patients.” Correction: repeated visits from one patient are correlated; the number of independent patient units is 20 if 20 patients supplied those observations.
- Wrong: “The between-patient coefficient shows that a patient's molecular abundance rises when that patient's response improves.” Correction: \(\beta_B\) compares patients' average response levels. The within-patient coefficient \(\beta_W\) addresses deviations from an individual's average.
- Wrong: “The protein changed during treatment, so treatment caused the change.” Correction: time trends in an observational treatment cohort can reflect treatment, disease course, concomitant care, batch, or other factors. The model describes associations unless the design and assumptions support a causal interpretation.
When many genes, proteins, or metabolites are tested, single-feature p-values do not account for the number of analyses. Prespecify features when possible and use an appropriate multiplicity strategy for exploratory screens. Report effect sizes in interpretable abundance units (or clearly state any transformation or standardization), uncertainty intervals, and the scale of the clinical response. Avoid selecting and reporting only the most striking features without disclosing the screening process.
Interpretation and Scope
Repeated multi-omics measurements can reveal whether molecular variation tracks clinical variation within patients, whether patients with different average response also differ molecularly, and whether baseline groups have different molecular trajectories. These are distinct questions, and the corresponding coefficients should be reported separately. A molecular signal may also be an assay artifact, a correlate of another clinical process, or a consequence rather than a driver of response.
The linear mixed model is a useful framework, not a universal solution. It may be unstable with very few patients, sparse visits, severe outliers, highly nonlinear trajectories, or a random-effects structure too complex for the data. A continuous response is assumed in the examples; binary or count outcomes require an appropriate generalized mixed model and a suitable link function. More complicated models do not compensate for weak design, missing time points, or poor assay quality.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- What does \(R_{ij}-\bar{R}_i\) represent, and what does \(\bar{R}_i\) represent?
- A fitted model estimates \(\beta_W=0.8\) abundance units per response-score unit. State the interpretation, assuming time and other terms are held fixed.
- Why can an unadjusted association between response and molecular abundance be misleading if both change steadily during treatment?
- A patient has abundances of 4, 7, and 10 units at days 0, 3, and 6. Calculate the linear change per day between the first and last visit.
- What does a time-by-baseline-group interaction estimate in a linear trajectory model?
Solutions. (1) \(R_{ij}-\bar{R}_i\) is the visit-specific response deviation from that patient's observed mean; \(\bar{R}_i\) is the patient's average response across included visits. (2) At a visit when a patient's response is one score unit above that patient's own mean, expected molecular abundance is 0.8 units higher, conditional on the modeled time trend and other terms; this is an association, not a causal effect. (3) Shared time trends can induce an association even if there is no response-specific relationship, so time must be represented and the remaining within-patient variation assessed. (4) The change is \((10-4)/(6-0)=6/6=1\) abundance unit per day. (5) It estimates the difference in expected molecular-abundance slopes per unit time between the baseline groups, under the specified trajectory model.