Worked Repeated-Measures Time Analysis

Why Model Repeated Visits?

When a donor contributes pseudobulk samples at several visits, those samples are not independent biological replicates. Measurements from the same donor can resemble one another because of stable donor characteristics, even after accounting for time. Treating every visit as independent can therefore understate uncertainty about a time trend. A donor random-intercept model represents that within-donor dependence while estimating the average change over time.

The worked analysis below uses a gene’s normalized log2 pseudobulk expression at three equally spaced visits. The outcome is expression for a donor-by-visit sample, not an individual cell count. The model is intended for a continuous, approximately linear time trend on the modeled expression scale.

Definition: For donor $i$ at visit $j$, a linear-time random-intercept model is $$ Y_{ij}=\beta_0+\beta_1 t_{ij}+b_i+\varepsilon_{ij}, $$ where $Y_{ij}$ is the transformed pseudobulk expression, $t_{ij}$ is time in the chosen unit, $\beta_0$ is the population-average expression at time zero, $\beta_1$ is the population-average change in expression per time unit, $b_i$ is a donor-specific random intercept, and $\varepsilon_{ij}$ is the residual. The standard model assumes $b_i\sim N(0,\tau^2)$, $\varepsilon_{ij}\sim N(0,\sigma^2)$, and that random intercepts and residuals are mutually independent. Residuals are independent across visits conditional on the donor effect and the fixed effects.

The random intercept induces positive dependence between visits from the same donor. Under these assumptions, two distinct visits from one donor have covariance $\tau^2$, while a visit has variance $\tau^2+\sigma^2$. The intraclass correlation is therefore $\tau^2/(\tau^2+\sigma^2)$. Different donors are assumed independent. The model also assumes that the mean trajectory is linear in time and that the time slope is common across donors; it does not say that every donor follows an identical observed trajectory.

For a gene-level analysis, this model is usually applied to an appropriately normalized and transformed pseudobulk outcome whose residual distribution is reasonably described by the Gaussian model. It should not be applied mechanically to untransformed counts as though they were normally distributed. The interpretation of $\beta_1$ is on the outcome scale: for log2 expression, it is a change in log2 units per time unit.

Worked Analysis With Three Visits

Consider a clearly hypothetical complete dataset with four donors, each measured at months 0, 1, and 2. The following summary describes the twelve donor-visit expression values. Donor means are means across the three visits; the pooled residual sum of squares is calculated after fitting the common intercept and linear time slope, including donor intercepts.

Summary from hypothetical data	Value
Number of donors and visits per donor	4 donors; 3 visits each
Times	0, 1, and 2 months
Donor means (log2-expression units)	5.3000, 5.0333, 5.4333, 5.2667
Estimated common slope	0.3125 log2-expression units per month
Residual sum of squares after fitting the slope	0.05875 squared log2-expression units

Worked Example: Estimating the Variance Components

For these balanced data, estimate the residual variance from the within-donor residual contrasts, then estimate the donor-intercept variance from the variation among donor means.

There are twelve observations. Fitting four donor intercepts and one common time slope leaves $12-5=7$ within-donor residual degrees of freedom. The restricted maximum likelihood (REML) residual variance estimate is

\hat{\sigma}^2=\frac{0.05875}{7}=0.0083929.

The sample variance of the four donor means is $0.0276852$ squared log2-expression units. Each donor mean averages three residuals, so its residual contribution to the variance is $\sigma^2/3$. Thus, for these balanced data, the REML donor-intercept variance estimate is the sample variance of the donor means minus the estimated variance of a three-visit residual mean:

\hat{\tau}^2=0.0276852-\frac{0.0083929}{3} =0.0276852-0.0027976 =0.0248876.

Both estimated variance components are positive. The estimated within-donor correlation is $\hat{\tau}^2/(\hat{\tau}^2+\hat{\sigma}^2)=0.0248876/(0.0248876+0.0083929)\approx0.7479$. This indicates substantial estimated similarity among measurements from the same donor in this hypothetical dataset; it is not evidence that time caused the similarity.

Estimating and Testing the Time Effect

In a complete balanced design with times 0, 1, and 2, the centered visit times are $-1,0,1$. A donor’s random intercept contributes equally to all three observations and cancels from the within-donor time contrast. The slope estimate is the total expression contrast between the last and first visits, divided by the corresponding time contrast, with the donor contributions averaged across the donors. In this example the fitted common slope is $0.3125$ log2-expression units per month.

Result: With complete, balanced repeated measurements at fixed equally spaced times, a random-intercept model’s common linear slope is informed by within-donor time contrasts. Under the normal-error assumptions, the estimated slope’s standard error in this design is $\sqrt{\hat{\sigma}^2/\sum_i\sum_j(t_{ij}-\bar{t}_i)^2}$. For four donors measured at 0, 1, and 2, the denominator is $4(1+0+1)=8$. The residual contrasts used to estimate $\sigma^2$ have $12-4-1=7$ degrees of freedom; the slope test therefore uses $df=7$ in this balanced example.

Worked Example: Testing the Average Monthly Change

Test $H_0:\beta_1=0$ against the two-sided alternative $H_a:\beta_1\ne0$, using $\alpha=0.05$. The assumptions are the stated normal random-intercept model, independent donors, a common linear time trend, and a complete balanced schedule. The residual variance estimate is positive, and the example has no missing visits.

The sum of squared centered visit times is $4(1^2+0^2+1^2)=8$. Therefore,

SE(\hat{\beta}_1)=\sqrt{\frac{0.0083929}{8}} =\sqrt{0.0010491} \approx0.03239 \quad\text{log2-expression units per month}.

The test statistic is

t=\frac{0.3125-0}{0.03239}\approx9.648,\qquad df=7.

The two-sided $p$-value is approximately $0.0000274$, which is below $0.05$. The data provide evidence against a zero average linear time slope under this model. A 95% confidence interval uses $t_{0.975,7}\approx2.3646$:

0.3125\pm2.3646(0.03239) \approx0.3125\pm0.0766 =[0.2359,\,0.3891] \quad\text{log2-expression units per month}.

Under the model, this interval estimates the population-average monthly change in log2 expression. Its repeated-sampling interpretation is that intervals constructed by this method would cover the true slope in 95% of repeated samples meeting the model assumptions. It does not assign a 95% probability to the fixed slope being in this particular interval.

The degrees of freedom matter. In this balanced design, the slope is based on within-donor contrasts, and the estimated residual variance comes from seven residual contrasts orthogonal to the fitted donor intercepts and common slope. Using $df=7$ reflects that variance estimation. Substituting the number of visits or the total number of observations as the denominator degrees of freedom would not represent the uncertainty in this slope test.

Interpreting the Coefficient and Checking the Fit

Worked Example: Translating the Slope to a Two-Month Change

The grand mean of the donor means is $(5.3000+5.0333+5.4333+5.2667)/4=5.2583$ log2-expression units. Because the visits are centered around month 1, this is the fitted population mean at month 1. The fitted intercept at month 0 is therefore $5.2583-0.3125=4.9458$, and the fitted mean at month 2 is $4.9458+2(0.3125)=5.5708$.

The estimated change from month 0 to month 2 is $2(0.3125)=0.6250$ log2-expression units. On the corresponding ratio scale, this is $2^{0.6250}\approx1.542$: the model’s fitted expression level at month 2 is about 1.54 times its fitted level at month 0. This is an estimated average trend in the modeled population, not a claim that every donor’s expression increased by that amount.

Before interpreting the coefficient, assess whether the assumptions are plausible. Plot each donor’s measurements against time to see whether a straight-line mean trajectory is reasonable and whether donor trajectories show pronounced departures from a common slope. Examine residuals against fitted values and time for changing spread or systematic patterns, and use a residual quantile plot to look for severe departures from normality. Check for influential donor samples and confirm that donor and visit labels were matched correctly. These diagnostics can reveal problems, but they cannot establish the model assumptions as true.

A complete design is especially helpful for separating the within-donor time contrast from donor-level variation. With missing visits or markedly unequal schedules, the simple balanced-design standard-error calculation and $df=7$ justification above do not automatically apply. Use the fitted model’s appropriate degrees-of-freedom procedure for the actual design, and inspect whether missingness or visit timing could be informative. A near-zero estimated variance component, convergence warnings, or a singular fit also calls for caution; the numerical output may not support the assumed dependence structure.

Alternatives and Limitations

For two visits, a paired analysis of within-donor changes is a direct alternative: it focuses on the mean paired difference and also respects the donor as the unit linking observations. With three or more visits, a linear-time mixed model summarizes change using one slope, making the trend and its units easy to report. It can be misleading if the mean trajectory bends or changes abruptly. Treating visit as a categorical factor avoids imposing a straight line, but estimates separate visit means and generally requires a different test for an overall time effect; it answers a less compact question than a single linear slope.

A donor fixed-effect regression can also account for stable donor differences by giving each donor a separate intercept. It targets within-donor change but does not estimate a donor-intercept variance in the same way as the random-intercept model. The random-intercept approach instead treats donors as draws from a population with a shared variance component, an assumption relevant to inference beyond the observed donors. Neither approach repairs confounding of time with batch, systematic changes in sample processing, or selective dropout. Those design and measurement problems require investigation rather than a different standard error alone.

A common mistake is to count twelve donor-visit samples as twelve independent donors and fit ordinary regression. That ignores the shared donor component: observations from the same donor are correlated, so the independence-based standard error can be too small. The correction is to represent the donor grouping in the model and base the test’s uncertainty on the relevant residual information. Another mistake is to describe a positive slope as a causal effect of time. The model estimates an association between visit time and expression under its design and assumptions; causal interpretation needs an appropriate study design and additional assumptions.

Key takeaway: A repeated-measures time analysis needs both a clear time estimand and a dependence model. In the balanced example, the random intercept accounts for stable donor-level differences, while the common slope estimates average monthly change. The result remains conditional on a plausible linear trend, an adequate outcome model, and valid donor-level independence.

Check Your Understanding

Answer each question, then compare with the solutions below.

In a random-intercept model, what does $\beta_1$ measure when time is recorded in months and the outcome is log2 expression?
For six donors with three complete visits at times 0, 1, and 2, what is the sum of squared centered visit times?
Why is the balanced example’s slope test based on seven residual degrees of freedom rather than eleven?
Give one diagnostic sign that the common linear time trend may not fit well, and name one suitable alternative representation of time.

Solutions. (1) It is the population-average change in log2-expression units per month. (2) Each donor contributes $1+0+1=2$, so the total is $6(2)=12$. (3) Twelve observations minus four donor intercepts and one common slope leave $12-4-1=7$ residual contrasts for estimating the residual variance. (4) A systematic curved pattern in donor trajectories or residuals against time suggests nonlinearity; treating visit as a categorical factor allows separate mean expression at each visit rather than imposing one straight-line slope.