---
title: Repeated-measures mixed models
output:
  rmarkdown::html_vignette:
    toc: yes
    toc_depth: 2.0
    css: albers.css
    includes:
      in_header: albers-header.html
params:
  family: red
  preset: interaction
resource_files:
- albers.css
- albers.js
- albers-header.html
- fonts

vignette: |
  %\VignetteIndexEntry{Repeated-measures mixed models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
if (requireNamespace("ragg", quietly = TRUE)) knitr::opts_chunk$set(dev = "ragg_png")
if (requireNamespace("albersdown", quietly = TRUE)) {
  if (requireNamespace("systemfonts", quietly = TRUE)) try(albersdown::albers_register_fonts(), silent = TRUE)
  if (requireNamespace("ggplot2", quietly = TRUE)) try(ggplot2::theme_set(albersdown::theme_albers(family = params$family, preset = params$preset)), silent = TRUE)
}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE
)
library(fmrigds)
set.seed(1)
```

```{r albers-classes, echo=FALSE, results='asis'}
cat(sprintf(
  paste0(
    '<script>document.addEventListener("DOMContentLoaded",function(){',
    'document.body.classList.remove("palette-red","palette-lapis","palette-ochre","palette-teal","palette-green","palette-violet","preset-homage","preset-interaction","preset-study","preset-structural","preset-adobe","preset-midnight");',
    'document.body.classList.add("palette-%s","preset-%s");',
    '});</script>'
  ),
  params$family,
  params$preset
))
```

Repeated-measures fMRI analyses are often statistically simple but operationally
awkward. You have one effect estimate per ROI or voxel, per subject, per visit
or condition, and you want a group model that respects the repeated structure
without hand-writing a loop over samples. `fmrigds` gives you a narrow, fast
path for that case.

This vignette shows the supported contract: one grouping factor, Gaussian
responses, shared observation layout across samples, and either a random
intercept or a random intercept plus one within-subject slope. By the end, you
will have one model-level GDS containing fixed-effect estimates and
variance-component assays that you can inspect or write out like any other
group result. You will also know when to use an unweighted model, when to treat
first-level variances as known sampling variances, what `pooled` actually
pools, and which claims the available summary statistics can support.

## 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:

```{r make-long-table}
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:

| Available input | Recommended reducer | What the model can claim |
|---|---|---|
| Effect estimate only | `lmm:ri` or `lmm:ri_slope1` | Models the empirical second-level covariance, without separating first-level sampling error |
| Effect estimate plus variance or standard error | Matching `*_knownvar` reducer | Treats the supplied marginal sampling variance as known and estimates additional heterogeneity |
| Effect estimate plus *t* or *z* | Derive or supply `var` only when the effect-scale standard error can be recovered | Same interpretation as `*_knownvar`; the test statistic alone is not the variance |
| *t* or *z* only | No effect-scale repeated-measures LMM | The effect magnitude and its sampling variance cannot generally be recovered from a standardized statistic alone |

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:

```{r quick-win}
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()
```

```{r quick-win-checks, include = FALSE}
lmm_summary <- data.frame(
  sample = sample_labels(g_rm),
  intercept = as.numeric(assay(g_rm, "coef:(Intercept)")),
  time = as.numeric(assay(g_rm, "coef:time")),
  vc_intercept = as.numeric(assay(g_rm, "vc_intercept")),
  vc_slope = as.numeric(assay(g_rm, "vc_slope")),
  lambda_slope = as.numeric(assay(g_rm, "lambda_slope")),
  stringsAsFactors = FALSE
)

lmm_summary_display <- transform(
  lmm_summary,
  intercept = round(intercept, 3),
  time = round(time, 3),
  vc_intercept = round(vc_intercept, 3),
  vc_slope = round(vc_slope, 3),
  lambda_slope = round(lambda_slope, 3)
)

stopifnot(
  identical(subjects(g_rm), "meta"),
  identical(contrasts(g_rm), "model"),
  all(is.finite(lmm_summary$time)),
  all(lmm_summary$vc_slope > 0),
  all(lmm_summary$lambda_slope > 0)
)
```

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:

```{r inspect-summary}
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:

```{r first-level-variance-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
)
```

```{r first-level-variance-contract-check, include = FALSE}
stopifnot(all(variance_probe$largest_absolute_change < 1e-10))
```

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()`](https://wviechtb.github.io/metafor/reference/rma.mv.html)
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`.

```{r known-variance-fit}
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)
```

```{r known-variance-fit-check, include = FALSE}
stopifnot(
  all(is.finite(as.matrix(knownvar_compare[-1]))),
  any(abs(knownvar_compare$knownvar_se - knownvar_compare$unweighted_se) > 1e-6),
  all(knownvar_compare$fitted_extra_residual_var > 0)
)
```

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)](https://ora.ox.ac.uk/objects/uuid%3Ae26cdb91-02a9-4252-b2f6-cc57a8e25844)
and the multi-package evaluation by
[Olszowy et al. (2019)](https://www.nature.com/articles/s41467-019-09230-w).

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)](https://pmc.ncbi.nlm.nih.gov/articles/PMC2719771/)
and
[Chen et al. (2012)](https://pmc.ncbi.nlm.nih.gov/articles/PMC3404516/).

## 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.

| Question | `theta_mode = "pooled"` | `theta_mode = "voxelwise"` |
|---|---|---|
| What is shared? | Random-effect-to-residual variance ratios, and the intercept–slope correlation for `covariance = "full"` | Nothing beyond the fixed/random design |
| Main advantage | More stable covariance shape and much less optimization | Allows genuine spatial or ROI heterogeneity in covariance shape |
| Main risk | A common shape can misrepresent regions with different longitudinal variability | Noisy or boundary variance estimates with small samples; spatially speckled component maps |
| Effect on fixed effects | Shared covariance determines the within-subject GLS weighting | Sample-specific covariance gives sample-specific GLS weighting |
| Sensible starting point | Whole-brain data with few repeats, or when covariance heterogeneity is not the scientific target | A modest number of ROIs, stronger replication, or a planned sensitivity analysis |

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.

```{r pooled-vs-voxelwise}
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.

```{r pooled-vs-voxelwise-checks, include = FALSE}
bench_subjects <- paste0("sub-", sprintf("%02d", 1:10))
bench_visits <- c("pre", "mid", "post")
bench_time <- c(-1, 0, 1)
bench_samples <- paste0("ROI_", sprintf("%02d", 1:24))

bench_df <- expand.grid(
  sample = bench_samples,
  subject = bench_subjects,
  contrast = bench_visits,
  KEEP.OUT.ATTRS = FALSE,
  stringsAsFactors = FALSE
)

bench_subject_intercept <- seq(-0.25, 0.25, length.out = length(bench_subjects))
names(bench_subject_intercept) <- bench_subjects
bench_subject_slope <- seq(-0.04, 0.04, length.out = length(bench_subjects))
names(bench_subject_slope) <- bench_subjects
bench_sample_shift <- seq(0.00, 0.45, length.out = length(bench_samples))
names(bench_sample_shift) <- bench_samples
bench_sample_slope <- seq(0.35, 0.85, length.out = length(bench_samples))
names(bench_sample_slope) <- bench_samples

bench_df$time <- bench_time[match(bench_df$contrast, bench_visits)]
bench_df$beta <- 0.8 +
  bench_sample_shift[bench_df$sample] +
  bench_subject_intercept[bench_df$subject] +
  (bench_sample_slope[bench_df$sample] + bench_subject_slope[bench_df$subject]) *
    bench_df$time +
  rnorm(nrow(bench_df), sd = 0.05)
bench_pooled_plan <- gds(bench_df, contrast_data_cols = "time") |>
  reduce(
    method = "lmm:ri_slope1",
    formula = ~ time,
    options = list(
      slope = "time",
      covariance = "diag",
      fit = "REML",
      theta_mode = "pooled"
    )
  )

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

bench_times <- data.frame(
  theta_mode = c("pooled", "voxelwise"),
  elapsed = c(
    median(replicate(3, unname(system.time(compute(bench_pooled_plan))[["elapsed"]]))),
    median(replicate(3, unname(system.time(compute(bench_voxel_plan))[["elapsed"]])))
  ),
  stringsAsFactors = FALSE
)

stopifnot(
  length(unique(as.numeric(assay(g_pooled, "lambda_slope")))) == 1L,
  all(is.finite(as.numeric(assay(g_rm, "lambda_slope")))),
  bench_times$elapsed[1] < bench_times$elapsed[2]
)
```

```{r pooled-vs-voxelwise-timing}
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()`.

```{r attach-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
)
```

```{r attach-contrast-data-checks, include = FALSE}
stopifnot(
  max(abs(as.numeric(assay(g_attached, "coef:time")) -
            as.numeric(assay(g_pooled, "coef:time")))) < 1e-10,
  max(abs(as.numeric(assay(g_attached, "lambda_slope")) -
            as.numeric(assay(g_pooled, "lambda_slope")))) < 1e-10
)
```

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

```{r random-intercept}
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.