Repeated-measures mixed models

What does the input look like?

The easiest route is a stacked long table with one row per sample, subject, and repeated-measure level. Here we build a small ROI example with three visits:

subjects <- paste0("sub-", sprintf("%02d", 1:6))
visits <- c("pre", "mid", "post")
time <- c(-1, 0, 1)

rm_df <- expand.grid(
  sample = c("ROI_1", "ROI_2"),
  subject = subjects,
  contrast = visits,
  KEEP.OUT.ATTRS = FALSE,
  stringsAsFactors = FALSE
)

subject_intercept <- c(-0.25, -0.10, 0.00, 0.10, 0.20, 0.30)
names(subject_intercept) <- subjects
subject_slope <- c(-0.10, -0.04, 0.00, 0.02, 0.06, 0.12)
names(subject_slope) <- subjects
sample_shift <- c(ROI_1 = 0.00, ROI_2 = 0.35)
sample_slope <- c(ROI_1 = 0.45, ROI_2 = 0.70)
sample_slope_scale <- c(ROI_1 = 0.5, ROI_2 = 2.0)

rm_df$time <- time[match(rm_df$contrast, visits)]
rm_signal <- 0.9 +
  sample_shift[rm_df$sample] +
  subject_intercept[rm_df$subject] +
  (sample_slope[rm_df$sample] +
     sample_slope_scale[rm_df$sample] * subject_slope[rm_df$subject]) * rm_df$time
rm_df$beta <- rm_signal + rnorm(nrow(rm_df), sd = 0.06)

knitr::kable(head(rm_df), row.names = FALSE)

The time column is constant within each visit, so it can be collapsed into contrast-level metadata during ingestion. The two simulated ROIs deliberately have different amounts of subject-to-subject slope variation; that makes the pooled-versus-voxelwise assumption visible later.

The response is beta. Despite the name, it need not be a raw regression coefficient: it can be any approximately Gaussian subject-level effect estimate with a consistent definition and scale across subjects and visits. Choose the reducer from the information that is actually available:

A t or z statistic is therefore not interchangeable with an effect estimate: its scale also depends on uncertainty and, for t, degrees of freedom. When an unthresholded t map was computed exactly as beta / se for the same contrast, the effect-scale variance can be recovered as (beta / t)^2; verify that provenance before using it. A z map may instead be a signed transformation of a p value, in which case beta / z is not the first-level standard error. If a trustworthy se assay is available, the known-variance reducers derive var from it automatically.

If the scientific target is evidence combination rather than an effect-scale repeated-measures model, a method such as combine:stouffer may be appropriate after accounting for dependence. Repeated observations from the same subject are not independent inputs to an ordinary Stouffer combination.

What model is fitted?

For sample \(v\), subject \(i\), and repeated-measure level \(r\), every random intercept plus slope reducer uses the same conditional mean:

\[ y_{vir} = x_{ir}^{\mathsf T}\beta_v + b_{0iv} + t_r b_{1iv} + \epsilon_{vir}. \]

Here, \(x_{ir}\) is the row of the fixed-effects design built from formula, and \(t_r\) is the numeric variable named by options$slope. The simpler random-intercept reducers drop \(b_{1iv}\).

The two reducer families differ in their marginal covariance. The unweighted methods assume

\[ \begin{aligned} (b_{0iv}, b_{1iv})^{\mathsf T} &\sim N(0,\;\sigma_v^2 \Lambda_v),\\ \epsilon_{vir} &\sim N(0,\;\sigma_v^2). \end{aligned} \]

This parameterization separates the sample-specific residual scale — sigma2 or vc_resid — from relative random-effect variances in λ (lambda_intercept, lambda_slope, and, for a full covariance, lambda_cov_intercept_slope). The absolute random-effect variances reported as vc_intercept and vc_slope are λ multiplied by sigma2. The simpler lmm:ri reducer drops \(b_{1iv}\) and estimates one intercept-to-residual variance ratio, lambda.

The *_knownvar methods instead use the supplied sampling-variance vector \(d_v\):

\[ V_v = \operatorname{diag}(d_v) + \sigma_{\text{resid},v}^2 I + ZG_vZ^{\mathsf T}. \]

Here, \(G_v\) contains absolute random-effect variances and \(\sigma_{\text{resid},v}^2\) is additional observation-level heterogeneity. The known sampling variances establish an absolute scale, so these reducers report vc_intercept, vc_slope, and vc_resid but not relative lambda_* assays. They treat var as fixed input rather than estimating a common multiplier for it.

covariance = "diag" fixes the intercept–slope covariance to zero. covariance = "full" estimates that covariance through a Cholesky parameterization and also reports corr_intercept_slope. The full model is more flexible, but it spends an additional variance parameter and can be weakly identified when there are few subjects or few repeated-measure levels. Centering the slope variable before ingestion makes the random intercept refer to a meaningful time point and usually makes the intercept–slope covariance easier to interpret.

What is the fastest end-to-end path?

Use gds(..., contrast_data_cols = "time") to put the repeated-measure metadata on the contrast axis, then call reduce() with one of the restricted LMM reducers. The first example has effect estimates but no sampling variances, so it uses the unweighted method:

g_rm <- gds(rm_df, contrast_data_cols = "time") |>
  reduce(
    method = "lmm:ri_slope1",
    formula = ~ time,
    options = list(
      slope = "time",
      covariance = "diag",
      fit = "REML",
      theta_mode = "voxelwise"
    )
  ) |>
  compute()

The result is still a GDS, but the subject and contrast axes have collapsed to the model level. A compact summary is usually easier to read than printing individual assay arrays:

knitr::kable(lmm_summary_display, row.names = FALSE)

Fixed-effect outputs still follow the same naming convention as the regression reducers (coef:<term>, se_coef:<term>, t_coef:<term>, p_coef:<term>), while variance components are exposed as ordinary assays such as vc_intercept, vc_slope, and lambda_slope.

How are the random effects estimated?

All four reducers stack the subject-by-visit observations for every sample, build one shared fixed-effects matrix and one shared random-effects design, and integrate the subject-specific random effects out of the Gaussian marginal likelihood. They estimate the covariance of the random effects; they do not return subject-level conditional modes or BLUPs.

The unweighted reducers then:

1. optimize the relative covariance parameters λ by profile maximum likelihood (fit = "ML") or restricted maximum likelihood (fit = "REML");
2. conditional on those fitted parameters, solve the generalized least-squares problem for each sample’s fixed effects; and
3. estimate a residual scale for each sample, then convert the relative covariance parameters into the absolute vc_* components.

For a known-variance reducer, diag(var) remains fixed while the optimizer estimates the absolute random-effect covariance \(G\) and the additional residual variance vc_resid. At every trial covariance, the reducer profiles out the fixed effects with generalized least squares. There is no single residual scale that can be factored out of the whole covariance, which is why these methods do not report relative lambda_* components.

REML is the default because it accounts for the fixed-effects degrees of freedom while estimating variance components. ML is mainly useful when models with different fixed-effects structures must be compared on the same data. For lmm:ri, the optimizer explicitly considers the zero random-intercept boundary. With covariance = "full", the slope methods obtain positive-definite covariance matrices from a Cholesky factor. Any of the fitted variance components can be near zero at the numerical boundary, especially with few subjects or little replication.

The reported fixed-effect standard errors are plug-in GLS standard errors: they condition on the fitted variance parameters. The reported t tests use n_observations - n_fixed_effects residual degrees of freedom. They do not use Satterthwaite or Kenward–Roger corrections and do not propagate uncertainty in the estimated variance components. With few subjects or weakly identified random slopes, treat p_coef:* as an approximate screening summary rather than small-sample-calibrated inference. The restricted implementation is designed for fast, common-design imaging fits; it is not a drop-in replacement for every inferential feature of lmer(). The p_coef:* assays are also uncorrected for the number of ROIs or voxels tested.

How do first-level variance estimates enter the model?

The original lmm:ri and lmm:ri_slope1 reducers are unweighted: they require only beta. If the input also contains var or se, it is carried by the source GDS but is not used in their likelihood, estimates, or weights. Their sigma2 and vc_* outputs are estimated entirely from the supplied effects.

The following check deliberately adds implausibly heterogeneous first-level variances. The LMM result is unchanged because var is outside this reducer’s input contract:

rm_df_with_var <- rm_df
rm_df_with_var$var <- rep(c(1e-6, 100), length.out = nrow(rm_df_with_var))

g_var_probe <- gds(rm_df_with_var, contrast_data_cols = "time") |>
  reduce(
    method = "lmm:ri_slope1",
    formula = ~ time,
    options = list(
      slope = "time",
      covariance = "diag",
      fit = "REML",
      theta_mode = "voxelwise"
    )
  ) |>
  compute()

variance_probe <- data.frame(
  output = c("time coefficient", "time standard error", "slope variance ratio"),
  largest_absolute_change = c(
    max(abs(assay(g_var_probe, "coef:time") - assay(g_rm, "coef:time"))),
    max(abs(assay(g_var_probe, "se_coef:time") - assay(g_rm, "se_coef:time"))),
    max(abs(assay(g_var_probe, "lambda_slope") - assay(g_rm, "lambda_slope")))
  )
)

knitr::kable(
  transform(variance_probe, largest_absolute_change = signif(largest_absolute_change, 3)),
  row.names = FALSE
)

For an unweighted fit, this has two distinct implications:

• A biased first-level variance estimate does not directly create an excessive inverse-variance weight.
• Differences in first-level precision are nevertheless unmodeled. Variation due to scan duration, motion, censoring, design efficiency, or first-level noise modeling is absorbed into the LMM residual and is not separated from other within-subject variation. When precision differs materially across observations, the homoscedastic residual assumption can be poor and the nominal fixed-effect standard errors can be misleading.

Use lmm:ri_knownvar or lmm:ri_slope1_knownvar when the input variances should enter the model. For each sample, these reducers fit

\[ V = \operatorname{diag}(v) + \sigma_{\text{resid}}^2 I + ZGZ^{\mathsf T}, \]

where \(v\) is the supplied var assay, \(G\) is the subject random-effect covariance, and vc_resid is additional observation-level heterogeneity not explained by the reported sampling variance. The supplied variances are treated as known: they are not multiplied by a fitted calibration factor.

This is the same covariance contract accepted by metafor::rma.mv() when V contains the known sampling variances and random effects represent the subject and observation-level heterogeneity. It is not equivalent to passing weights = 1 / var to lme4::lmer(): in lmer(), the inverse weights are still multiplied by an estimated common residual scale.

The next example makes the variance contract part of the data-generating process. Subjects have different first-level standard errors, visits differ slightly in precision, and every observation also contains residual variation with variance 0.0064 that is not included in var.

rm_df_knownvar <- rm_df
subject_se <- c(0.02, 0.03, 0.05, 0.07, 0.09, 0.12)
names(subject_se) <- subjects
visit_precision_factor <- c(pre = 1.15, mid = 0.85, post = 1.00)

rm_df_knownvar$var <- (
  subject_se[rm_df_knownvar$subject] *
    visit_precision_factor[rm_df_knownvar$contrast]
)^2

true_extra_residual_var <- 0.0064
set.seed(11)
rm_df_knownvar$beta <- rm_signal +
  rnorm(nrow(rm_df_knownvar), sd = sqrt(rm_df_knownvar$var)) +
  rnorm(nrow(rm_df_knownvar), sd = sqrt(true_extra_residual_var))

g_knownvar_unweighted <- gds(rm_df_knownvar, contrast_data_cols = "time") |>
  reduce(
    method = "lmm:ri_slope1",
    formula = ~ time,
    options = list(
      slope = "time",
      covariance = "diag",
      fit = "REML",
      theta_mode = "voxelwise"
    )
  ) |>
  compute()

g_knownvar <- gds(rm_df_knownvar, contrast_data_cols = "time") |>
  reduce(
    method = "lmm:ri_slope1_knownvar",
    formula = ~ time,
    options = list(
      slope = "time",
      covariance = "diag",
      fit = "REML"
    )
  ) |>
  compute()

knownvar_compare <- data.frame(
  sample = sample_labels(g_knownvar),
  unweighted_time = as.numeric(assay(g_knownvar_unweighted, "coef:time")),
  knownvar_time = as.numeric(assay(g_knownvar, "coef:time")),
  unweighted_se = as.numeric(assay(g_knownvar_unweighted, "se_coef:time")),
  knownvar_se = as.numeric(assay(g_knownvar, "se_coef:time")),
  mean_sampling_var = as.numeric(assay(g_knownvar, "sampling_var_mean")),
  fitted_extra_residual_var = as.numeric(assay(g_knownvar, "vc_resid"))
)

knitr::kable(knownvar_compare, digits = 4, row.names = FALSE)

Both fits use exactly the same effect estimates. Their coefficients and standard errors differ because the known-variance fit recognizes that some observations are measured less precisely. fitted_extra_residual_var estimates the simulated 0.0064 component; with a data set this small it need not recover that value exactly. A near-zero estimate in real data can mean either that the reported sampling variances explain the remaining observation-level variation or that the component is weakly identified.

There are two important limits. First, var must contain strictly positive sampling variances on the squared scale of beta; if the input supplies a standard error, square it during ingestion or map it to se so that fmrigds can derive var. Second, this version assumes the sampling-error covariance is diagonal. Contrasts estimated from the same first-level GLM can have correlated errors. Marginal variance maps alone cannot recover that covariance, so do not interpret the diagonal model as a full solution for dependent first-level contrasts.

Because the variances are treated as known, their accuracy now matters directly. A variance biased downward gives that observation too much precision; vc_resid can absorb general excess heterogeneity but cannot repair an observation-specific miscalibration. Inspect the distribution and provenance of var, and use sensitivity analyses when first-level noise modeling is doubtful.

That distinction matters in fMRI. Under a correctly specified mean model with mean-zero errors unrelated to the regressors, ordinary least squares coefficients need not be biased merely because residuals are positively autocorrelated, but an IID variance calculation generally overstates the effective information and can make standard errors too small. Prewhitening choices can materially change first-level uncertainty and, in some settings, the resulting activation estimates. See Woolrich et al. (2001) and the multi-package evaluation by Olszowy et al. (2019).

Before treating a beta map as a comparable response, record how it was made: the first-level design and contrast, effect scaling, temporal autocorrelation model and prewhitening, censoring or nuisance regression, residual degrees of freedom, and any transformations. If first-level sampling variances need to enter the group model, use a *_knownvar method and audit the variance model; the unweighted methods are not measurement-error or inverse-variance meta-analytic LMMs. For the distinction between ordinary and precision-weighted two-stage group analysis, see Mumford and Nichols (2009) and Chen et al. (2012).

Should variance components be pooled or fitted voxelwise?

This choice currently applies to the unweighted reducers. The known-variance reducers are voxelwise only: sharing an absolute random-effect covariance across voxels would usually be scientifically indefensible, while sharing a relative covariance requires an explicit cross-voxel scaling rule that is not yet part of the API.

The choice is a bias–variance tradeoff, not just a speed switch. In pooled mode, all complete samples contribute to one profile-likelihood estimate of the relative covariance λ. Each sample still gets its own residual variance sigma2, so its absolute vc_intercept and vc_slope may differ. In voxelwise mode, both the relative covariance and the residual scale are estimated separately for every sample.

Pooled mode is not spatial smoothing and does not shrink voxelwise estimates toward a mean; it fits one common covariance shape by summing profile-likelihood contributions across samples. It also does not model spatial correlation among samples. A large set of highly correlated voxels therefore supplies less independent information than its voxel count suggests, and a common λ can be a poor description when tissue classes or regions have genuinely different variance ratios.

g_pooled <- gds(rm_df, contrast_data_cols = "time") |>
  reduce(
    method = "lmm:ri_slope1",
    formula = ~ time,
    options = list(
      slope = "time",
      covariance = "diag",
      fit = "REML",
      theta_mode = "pooled"
    )
  ) |>
  compute()

theta_compare <- data.frame(
  sample = sample_labels(g_rm),
  pooled_lambda_slope = as.numeric(assay(g_pooled, "lambda_slope")),
  voxelwise_lambda_slope = as.numeric(assay(g_rm, "lambda_slope")),
  pooled_vc_slope = as.numeric(assay(g_pooled, "vc_slope")),
  voxelwise_vc_slope = as.numeric(assay(g_rm, "vc_slope")),
  pooled_se_time = as.numeric(assay(g_pooled, "se_coef:time")),
  voxelwise_se_time = as.numeric(assay(g_rm, "se_coef:time")),
  stringsAsFactors = FALSE
)

knitr::kable(
  theta_compare,
  digits = 3,
  row.names = FALSE
)

The pooled λ is identical across rows by construction. The absolute pooled slope variance need not be identical because it is the shared λ multiplied by each ROI’s residual variance. Differences between the final two columns show the inferential consequence of the choice: the fitted covariance changes the GLS standard error, not merely a diagnostic variance map. With only six subjects, the voxelwise covariance estimates in this toy example should be read as a sensitivity demonstration, not as trustworthy ROI-specific biology.

knitr::kable(transform(bench_times, elapsed = round(elapsed, 3)), row.names = FALSE)

Runtime is a practical reason to pool, but not a statistical justification. For a real analysis, start from the covariance assumption you can defend and then inspect sensitivity: compare coefficient and standard-error maps, the distribution of lambda_*, and the frequency of near-boundary components. The current converged assay distinguishes an available numerical fit from an unavailable one; it is not a rich optimizer-diagnostic code. A large change in scientific conclusions is evidence that the random-effects structure is consequential; it is not evidence that the more favorable mode is correct.

What if your repeated-measure metadata is already separate?

If your time, visit order, or other repeated-measure annotations already live outside the long table, attach them explicitly with with_contrast_data().

contrast_df <- data.frame(time = time, row.names = visits)

g_attached <- gds(rm_df[c("sample", "subject", "contrast", "beta")]) |>
  with_contrast_data(contrast_df) |>
  reduce(
    method = "lmm:ri_slope1",
    formula = ~ time,
    options = list(
      slope = "time",
      covariance = "diag",
      fit = "REML",
      theta_mode = "pooled"
    )
  ) |>
  compute()

knitr::kable(
  data.frame(
    sample = sample_labels(g_attached),
    time = round(as.numeric(assay(g_attached, "coef:time")), 3),
    stringsAsFactors = FALSE
  ),
  row.names = FALSE
)

If you only need a random intercept, use method = "lmm:ri" instead:

g_ri <- gds(rm_df, contrast_data_cols = "time") |>
  reduce(
    method = "lmm:ri",
    formula = ~ time,
    options = list(fit = "REML", theta_mode = "pooled")
  ) |>
  compute()

knitr::kable(
  data.frame(
    sample = sample_labels(g_ri),
    time = round(as.numeric(assay(g_ri, "coef:time")), 3),
    lambda = round(as.numeric(assay(g_ri, "lambda")), 3),
    stringsAsFactors = FALSE
  ),
  row.names = FALSE
)

What is supported, and what is intentionally out of scope?

The repeated-measures LMM reducers are deliberately narrower than lmer():

• one grouping factor only
• Gaussian outcomes only
• shared observation layout across all samples
• random intercept only or random intercept plus one slope, each available as an unweighted method or a *_knownvar method
• unweighted methods estimate one residual variance per sample and ignore observation-level precision
• known-variance methods add diag(var) to the marginal covariance and estimate additional residual heterogeneity; they currently assume diagonal sampling covariance and voxelwise variance components
• complete observations within a sample; a non-finite subject-by-repeat value makes that sample’s model outputs unavailable
• plug-in GLS inference with residual degrees of freedom, without Satterthwaite or Kenward–Roger correction
• no crossed random effects
• no general random-effects formula grammar
• no spatial covariance model or spatial regularization of voxelwise variance components

That restriction is the point: it lets fmrigds reuse one design across many responses and produce a whole-brain or whole-ROI result with one call to reduce(). It also defines the boundary of the scientific claim: the reducer models dependence among repeated effect estimates from the same subject. The known-variance variants additionally model reported marginal sampling uncertainty, but none of these reducers models dependence between neighboring samples or sampling-error covariance between contrasts.

Where do you go next?

• See ?reduce for the reducer contract and options.
• See ?with_contrast_data if your repeated-measure metadata is already separate from the long table.
• See inst/examples/fmrigds_lmm_benchmark.R for a benchmark harness that compares pooled and voxelwise theta modes.
• See inst/examples/fmrigds_lmm_knownvar_benchmark.R for a benchmark of the variance-aware random-slope reducer.