Run a parameter-recovery grid — recovery

The loop the validation scripts in bmm wrote by hand, with the requirements their overnight runs taught: one durable file per cell written as soon as the cell finishes, a resume that reads those files and needs no temporary state, a smoke mode in its own directory, a preflight fit that catches a compile or init error before any cell runs, and cells ordered so that the first completed block spans the design rather than exhausting one level.

Usage

recovery_grid(
  model,
  grid,
  pars,
  dir,
  reps = 1,
  sds = NULL,
  cors = NULL,
  covariates = NULL,
  tasks = NULL,
  task_col = "task",
  coding = c("cell", "contrast"),
  contrasts = NULL,
  formula = NULL,
  prior = NULL,
  prior_sd = NULL,
  trial_design = NULL,
  generator = NULL,
  seed = NULL,
  subjects = c("redraw", "fixed"),
  re_cor = c("none", "within", "all"),
  scale = c("natural", "link"),
  levels = c("population", "subject"),
  correlations = NULL,
  cor_scale = "link",
  convergence = NULL,
  ml = FALSE,
  smoke = FALSE,
  preflight = TRUE,
  cells = NULL,
  ...,
  .fitter = NULL
)

Arguments

model

A bmmodel, or a function of one grid row (a one-row data frame) returning one, for models whose constructor depends on the design. Or a list of recovery_component()s, or a function of the row returning such a list, to simulate several models for the same people and fit each on its own; see the section "Components".

grid

A data frame with the columns n_subjects and n_trials (with components, n_subjects only; see "Components" for their columns). A column named after a parameter gives that cell's population value on the link scale, overriding pars; a column sd_<parameter> overrides sds; a column cor_<a>__<b> sets the correlation of two parameters or covariates (the names in either order). With tasks, these columns may also use full task terms (kappa_task2, sd_kappa_task2, cor_kappa_task1__kappa_task2); a bare parameter column sets every task, and a full-term column in the same row overrides it for its task. A SimDesign design is a data frame and works as is.

pars, sds, cors

Defaults for every cell, as in simulate_recovery(). Each may also be a function(row) of the one-row grid data frame, evaluated under the cell's seed, which draws new hyperparameters for every data set; the grid columns above are applied to its result.

dir

Directory for the per-cell files; created if missing.

reps

Replications per cell.

covariates

As in simulate_recovery(), the same for every cell.

tasks, task_col

As in simulate_recovery(), the same for every cell. With tasks, the default formula is recovery_formula(model, re_cor = re_cor, task_col = task_col).

coding, contrasts

As in simulate_recovery(), the same for every cell. coding = "contrast" needs tasks, makes the default formula 1 + task rather than 0 + task, and scores the intercept and the contrasts in place of the cell means. What is simulated does not change, so a grid can be run under both codings from one set of generating values.

formula

A bmmformula; NULL means recovery_formula() of the cell's model; or a function(row) of the one-row grid data frame returning a bmmformula, called once per row, for designs whose formula depends on the row. For a custom m3 with its built-in generator the formula also holds the activations, so it reaches the generator too (see simulate_recovery()) and is written to the grid record.

For a model with alternative parameterizations (from bmm 1.4.1, criterion and bias of sdt_yn and sdt_rating), pars chooses the member that is simulated and recovered, and the default formula of each cell writes that member, as recovery_formula(pars = ) does. A given formula, or a formula or start in ml, must estimate the same member; one that names the other, or names neither while pars chose the non-default member, is refused before any cell is fitted. A function pars is checked as each row's first cell is simulated, so a mismatch in a later row stops the grid after the earlier rows' fits. With ml, each cell's ML fit estimates the cell's member.

prior

Passed to fit_cached(); with components, NULL or a list with a prior per component.

prior_sd

The prior SD each posterior is contracted against (#16): NULL (the default) for no contraction and no extra fit; "fit" to fit the prior alone (sample_prior = "only") once per distinct formula, model and prior among the grid's rows, cached in dir through fit_cached() with the backend named, and reused by every cell that shares it; "analytic" for the closed forms of prior_sd_table() applied to bmm::default_prior() with prior, which needs no fit and leaves NA where there is no closed form; or a table as recover() takes it, which may carry a condition column ("row-1", ...) to differ between rows. Subject rows have no contraction. The table used is written to <dir>/prior-sd.rds, which collect_grid() reads.

trial_design

As in simulate_recovery(), the same for every cell. A data frame without id needs every row of the grid to have its number of rows as n_trials; a function of n_trials serves a grid whose n_trials varies. Checked against every row before any cell runs, and recorded in <dir>/grid.rds.

generator

As in simulate_recovery().

seed

Master seed. Each cell derives its own from it and the cell's row and replication, so a cell is reproducible on its own: the cell's seed fixes its simulated data, the initial values of its chains and the sampler, so a rerun gives the same draws (see "Seed" in fit_cached() for what it does not fix: dynamic within-chain threading, the caller's RNGkind() and another platform's arithmetic). NULL leaves everything unseeded and records NA.

subjects

"redraw" draws new subject values in every replication; "fixed" draws them once per row and reuses them, so replications become a simulated retest of the same people. The realised population values, SDs, correlations and covariate values of replication 1 are reused too.

re_cor

Passed to recovery_formula() when formula is NULL.

scale

The scale recovery is scored on, as in recover(); natural by default.

levels

The estimate levels extracted from every fit and scored: one or more of "population", "subject" and "sd" (the between-subject standard deviations, always scored on the link scale).

correlations

NULL (the default) for no correlations, or one or more of the estimators of extract_correlations(), "model", "draws" and "point", to extract and score the between-subject correlations of every cell, pairs with covariates included.

cor_scale

The scale the correlations are extracted and scored on, "link" or "natural". The "model" estimator exists on the link scale only.

convergence

NULL, the default, gates every cell with check_convergence()'s defaults. A named list gives its arguments (rhat_max, ess_bulk_min, ess_tail_min, divergent_max, treedepth_max, variables), checked before any cell runs. The gate is computed once per cell and its whole diagnostics row is stored in the cell's sidecar, so a study can report why a cell failed without reading the fit.

ml

FALSE, the default, fits every cell hierarchically only. TRUE also fits every cell subject by subject with fit_ml() on the same simulated data, and scores both estimators against one truth at the subject level; the ML rows carry estimator = "ml". A named list gives arguments to fit_ml() (prior, draws, max_abs_link, ci_level, formula, .fitter, and through its ... the fitter's own, such as refresh = 0, silent = 2 to quiet the optimiser); model, data, file, refit, seed and by are the grid's. Needs "subject" in levels. See the section "Subject-wise ML".

smoke

TRUE runs the first two rows with two replications into <dir>/smoke, so a smoke run never overwrites a full one.

preflight

Run the first cell once with one chain and 200 iterations into <dir>/preflight before the loop. Skipped when the first cell already has a cached fit or a sidecar that can be used.

cells

Which cells this call runs, to split one grid over several processes. NULL, the default, runs every cell. integer(0) is the setup call: it runs no cell. Cell numbers (one element of grid_tasks()) run those cells of the full grid. See "Running a grid in several processes".

...

Passed to fit_cached() and so to the fitter: chains, iter, warmup, backend, cores, ...

.fitter

As in fit_cached().

Value

A bmmtools_recovery object over every cell that produced a fit, at the requested levels, with a condition column ("row-<i>") naming the grid row. summary() groups by condition. Its attributes:

cells: a tibble with one row per cell (condition, replication, n_subjects, n_trials, seed, file, status, elapsed in seconds, converged), followed by the cell's convergence record (max_rhat, min_ess_bulk, min_ess_tail, n_divergent, n_max_treedepth, n_variables, failed) and what it was run with (fit_seconds, chains, iter, threads). All of these come from the cell's sidecar, so a resumed grid reports them without reading a fit. elapsed is how long the cell took in this run, read time included; fit_seconds is how long the fit itself took when it was run.
grid: the grid as run. The same record, with everything else the grid was called with, is written to <dir>/grid.rds, which is what collect_grid() reads.
correlations: with correlations requested, a bmmtools_cor_recovery from recover_correlations() over every cell, with condition and replication; otherwise NULL.
subject_means: a tibble with one row per cell, subject and term: condition, replication, id, term, covariate (whether the term is a covariate), the posterior mean and median on the link scale, and true_value, the generating value on the link scale (a covariate's value is its own mean, median and true value). Its attribute links is the model's link table. subject_table() turns it into one row per subject.
ml_cells: with ml, a tibble with one row per cell (condition, replication, status, elapsed, n_subjects, n_converged, converged, file) for the ML fit; see "Subject-wise ML". cells stays one row per cell and describes the hierarchical fit.
run_info: one row per call that reached its cells in dir (this one included), in the order they started, read from <dir>/run-info/: call ("all", "setup" or "slice"), cells (as text, "1-3,7"), n_cells, n_reused (cells whose hierarchical results were read from their sidecar, not fitted), n_failed, started, ended, pid, the scheduler's scheduler, job_id, array_job_id, array_task_id and cores_granted (NA without one, as on a cloud VM), then every run_info() column. ended, n_reused and n_failed are NA while the call runs, and stay NA for a call that never finished its cells: a process that was stopped, or an R error outside a cell's fit once the cells had started (data that cannot be generated, say). A cell whose fit errors is counted in n_failed instead.

With cell numbers in cells, the object holds those cells only, and collect_grid() reads the whole grid. The setup call returns, invisibly, a tibble of every cell (condition, replication, n_subjects, n_trials, seed, file).

Details

A cell whose fit errors is recorded with status = "error" and the grid goes on; a warning at the end names the failed cells. A cell whose data cannot be generated stops the grid, because that is a design error rather than a sampler accident. Cell files are cell-<row>-rep-<rep>-sim.rds (data and truth) and cell-<row>-rep-<rep>.rds with its .key (the fit, through fit_cached()); a resume reads an existing simulation file rather than regenerating, so a cached fit and its truth always match.

The sidecar. While a fit is in memory, everything scored from it is written to cell-<row>-rep-<rep>-est.rds: the estimates at levels, the correlations, the subject means, and what they were extracted with (the cell's cache key, the bmmtools version, levels, correlations, cor_scale). A resume uses the sidecar without reading the fit when the key and the version match and it holds every requested level and estimator, so the fit files may be deleted once a grid has run. Otherwise the cell goes through fit_cached() again, which reuses a cached fit, and the sidecar is rewritten. A sidecar written before bmmtools stored the central 50 % interval counts as out of date for the same reason: it is re-extracted from the cached fit, not refitted.

Backend

Use cmdstanr (backend = "cmdstanr", or no backend and no brms.backend option with cmdstanr installed, which is what bmm then picks). It compiles a model once and every cell with the same Stan code reuses the executable; rstan compiles again for every fit. Measured 2026-10-04 on a mixture2p grid (bmm 1.3.2.9000, brms 2.23.0, rstan 2.32.7, cmdstanr 0.9.0, CmdStan 2.40.0, Apple M3 Max): under rstan each fit spent 19 to 23 s compiling, 77 to 85 % of a small cell (10 to 20 subjects, 50 trials, 2 chains of 500 iterations) and 22 to 23 % of a cell of 40 subjects, 100 trials and 4 chains of 2000 iterations; under cmdstanr the one compile took about 6 s and each later cell 0.03 to 0.05 s to find the executable. cmdstanr keeps executables in tempdir() unless told otherwise, so a grid resumed in a new R session compiles once more; options(cmdstanr_write_stan_file_dir = "<a directory that persists>") avoids that. The same holds for sbc().

Running a grid in several processes

A grid can be run as one setup call and then as slices, each in its own R process (one per core group of a large machine, or one per job), all writing into one dir:

The setup call, cells = integer(0), runs once. It writes grid.rds, runs the preflight, fits the priors of prior_sd and writes prior-sd.rds, and simulates every cell. It fits no cell.
Each slice is the same call with cells = grid_tasks(grid, reps, n)[[k]]. It checks its call against grid.rds, including the prior, model, formula, fitter arguments and prior_sd the setup call recorded, and stops when they differ or when no setup call wrote grid.rds. It fits its cells only and writes neither record.
collect_grid() reads the cells of all slices as one grid.

A full run (cells = NULL) in a dir that a setup call wrote, with the setup call's arguments, leaves grid.rds as it is, and slices can still run. One with another prior, model, formula, fitter arguments (other than cores) or prior_sd rewrites grid.rds without what the slices are checked against, since it may also have rewritten prior-sd.rds; every slice after it stops until the setup call is run again. One with another seed or reps rewrites grid.rds with its own seed and reps but keeps what the slices are checked against, so slices with its seed and reps are accepted. With prior_sd = "fit", prior-sd.rds holds the prior SDs of the old seed until that run has fitted its priors: start slices only after it has returned, and if it was stopped, run it or the setup call with its seed again first.

A cell has the seed and the files the full grid gives it, whichever slice runs it and in whatever order, so the slices together write the cells one call with cells = NULL writes. That needs a single number as seed (a call with cells and any other seed is refused), and each slice in its own R process under R's default RNGkind(): a worker of a parallel backend may run another generator (L'Ecuyer-CMRG under future.seed = TRUE), and its fits would differ. Slices given by grid_tasks() share no cell; do not run the same cell in two processes at once. Start the slices after the setup call has returned: it simulates the cells last. A slice that is stopped and run again resumes, like any grid. A slice may name fewer cores than the setup call. Its sidecars then record those cores among the sampler arguments where a sequential run's sidecars record its own; no result reads that field, and a cell's cache key leaves cores out, so the fits and estimates are the same. Each process compiles the Stan model once (see "Backend"). Do not point cmdstanr_write_stan_file_dir at one directory for processes that start together: they would compile the same executable at the same time.

The component form (a list of models) does not take cells, and neither does smoke = TRUE.

Subject-wise ML

With ml, every cell is fitted twice on the same simulated data: the hierarchical fit, and one fit_ml() call with no pooling, cached as cell-<row>-rep-<rep>-ml.rds under the cell's seed. The ML rows are bound to the hierarchical subject rows before scoring, so the result has both estimators at the subject level, "bayes" and "ml", against one truth; summary() reports them as separate rows and plot_recovery(color_by = "estimator") colours them. An ML fit has no population or sd level, so those levels stay hierarchical only, and its subject means are not added to subject_means. The ... of the grid (chains, iter, backend, refresh, ...) belong to the hierarchical fit and do not reach fit_ml(); ml = list(draws = , refresh = 0) and the other fit_ml() arguments do.

The sidecar records the ML request (its prior, draws, interval level, link-scale range and formula) with the ML rows, and a resume reuses them only when the request is the same; a change to max_abs_link or ci_level re-extracts from the cached ML fit without refitting, and a change to prior, draws or formula refits through fit_cached(). A cell whose ML fit fails keeps its hierarchical rows, is recorded with status = "error" in ml_cells, and is named in a warning at the end; its sidecar is written without ML rows, so a resume tries the ML fit again. There is no ML preflight: a compile error of the no-pooling model shows up as every cell's ML status.

Within a cell, a subject whose ML estimate leaves the link-scale range keeps its row with estimate = NA and converged = FALSE, as fit_ml() documents; the recovery rows carry that verdict per subject and term, and ml_cells counts the subjects per cell. Report n next to any comparison of the two estimators.

Components

With model a list of recovery_component()s (or a function(row) returning one, the same components by name in every row), each cell simulates one set with simulate_components() under the cell's seed and fits every component on its own through fit_cached(). pars, sds, generator, trial_design, tasks, task_col, re_cor and formula belong to the components and are errors here; cors, covariates, seed, reps, scale, levels, correlations, cor_scale, smoke, preflight and ... keep their meaning, and prior is NULL or a list by component, as in fit_components(). subjects = "fixed" is not supported for components yet. A component cannot be named sim, est, cor or sd.

Grid columns: n_subjects; n_trials_<comp> for a component's number of trials (a plain n_trials column is an error); <comp>_<par> or <comp>_<par>_<task_col><level> for population values, a bare parameter setting every task and a full term overriding it; sd_<comp>_<par> for SDs; cor_<a>__<b> with prefixed terms or covariates, in either order. They apply to numeric pars and sds of a component and to the result of a function-valued one. A missing value (NA) leaves the component's value, so a column may name a task that only some rows have. A column whose component or term is unknown is an error; so is a column whose prefix is no component but whose remainder is a term of one, such as a misspelt component name. Other columns are left for model and cors functions.

Files per cell: cell-<row>-rep-<rep>-sim.rds holds the set; cell-<row>-rep-<rep>-<comp>.rds (with its .key) each fit; a sidecar cell-<row>-rep-<rep>-<comp>-est.rds per component holds its estimates and, with "model" requested, its model correlations; and cell-<row>-rep-<rep>-cor.rds holds what needs every fit, the correlations of extract_correlations() on the set of fits (covariates included) and the subject means, with the components' keys. A resume in which every sidecar matches reads no fit; otherwise the fits are obtained through fit_cached(), which reuses cached ones, and the stale sidecars are rewritten. The preflight fits every component of the first cell into <dir>/preflight-<comp>, and is skipped when each of them has a cached fit or a usable sidecar.

The result is one bmmtools_recovery with prefixed terms (a_kappa, b_c_task1); converged is that component fit's verdict, and scale = "natural" uses every component's links. The cells attribute has one row per cell and component, with a component column; the correlations and subject_means attributes use prefixed terms, so subject_table() works on the result. A component whose fit fails is recorded as "error" for that cell, the cell's other components are still scored, its correlations and subject means are skipped, and the grid goes on.

Correlations between parameters of separately fitted components are attenuated by the reliabilities of both estimates; see "Separate fits" in extract_correlations(). The "model" estimator works within one fit only and needs a component formula with correlated terms.

Examples

if (FALSE) { # \dontrun{
out <- recovery_grid(
  bmm::mixture2p(resp_error = "y"),
  grid = expand.grid(n_subjects = c(20, 50), n_trials = c(30, 100)),
  pars = c(kappa = log(8), thetat = qlogis(0.75)),
  sds = c(kappa = 0.3, thetat = 0.5),
  dir = "fits/mixture2p", reps = 10, seed = 1,
  chains = 4, iter = 2000, backend = "cmdstanr"
)
summary(out)
} # }