sim_recovery() joins the posterior summaries of a model fitted to
simulated data with the true generating parameter values recorded by
simulate_data(), using the brms parameter names on both sides.
It replaces the error-prone manual name matching otherwise needed to
compare gen_outcome() truth with brms::fixef() and
brms::VarCorr() output.
Usage
sim_recovery(fit, analysis, probs = c(0.025, 0.975))Arguments
- fit
A
brmcoda()model object or brms::brmsfit fitted toanalysis$data(oranalysis$complr) usinganalysis$formula.- analysis
An
mlsim_analysisobject fromprep_sim_analysis()with a non-NULL$truthelement.- probs
A length-two numeric vector of lower and upper posterior interval probabilities used for the coverage indicator. Defaults to
c(0.025, 0.975).
Value
A data.table::data.table of class mlsim_recovery with one
row per model parameter:
type:"fixed","random_sd","random_cor", or"rescor". Random-effect rows are reported on the standard deviation / correlation scale used bybrms::VarCorr().group: grouping factor for random-effect rows,NAotherwise.parameter: the brms-style parameter label.truth: the true generating value from the simulator (sqrt()of covariance diagonals forrandom_sdrows, correlations forrandom_corrows).estimate,est_error,lower,upper: posterior summary of the fitted parameter atprobs.bias:estimate - truth.covered: whethertruthlies inside[lower, upper].simulator_name: the simulator-side parameter identifier the row was derived from, for provenance.
Details
The truth table is constructed by prep_sim_analysis() and stored in
its $truth element; sim_recovery() matches every truth row to
exactly one fitted parameter and vice versa. Any parameter that cannot
be matched in either direction is an error, never a silently wrong
row: truth values are attached purely by name, and duplicate or
ambiguous names (for example response names that collide after brms
removes _ and . characters) abort with an explanatory message.
All gen_outcome() families are supported, not only Gaussian models.
Rows for a modeled distributional parameter are labeled with the
family's brms parameter name (sigma_* for "gaussian", shape_*
for "negbin" and "gamma", phi_* for "beta"), and their truth
values are on the log link scale that both gen_outcome() and brms
use for these parameters; location rows are on the family's location
link (identity, log, or logit), again matching on both sides.
"poisson" and "binomial" families have no distributional
parameter, so their tables contain no scale rows, and a binomial
trials() term contributes no parameter. Because ar1(),
multivariate mvbind() responses, and compositional outcomes are
Gaussian-only simulator features, recovery tables for the other
families never contain rescor rows.
Two situations yield no recovery row at all: cross-block random-effect
correlations when the analysis was prepared with
prep_sim_analysis()(link_random = FALSE) (the fitted model never
samples them), and analysis models whose structure cannot be aligned
with the simulator's single joint covariance (recorded in
metadata$truth_unavailable_reason).
Interpreting the bias and covered columns requires knowing what the
fitted model estimates. prep_sim_analysis() deliberately builds the
pragmatic observed-data model that applied analysts commonly fit, not
the matched model for the generating process, so some parameters
estimate a related quantity rather than the one the simulator recorded
– most notably whenever the formula contains ar1(), and for
between()/within() terms under the default
centering = "manifest". Systematic bias in those parameters is an
expected property of the estimator, not an error in the simulator. See
the "Pragmatic default estimator" section of prep_sim_analysis()
before drawing conclusions from a recovery table.
See also
prep_sim_analysis() for the analysis object and its
$truth table, gen_template() for authoring the truth values.
Examples
# \donttest{
if (requireNamespace("cmdstanr", quietly = TRUE)) {
params <- list(
location = list(beta = matrix(
c(0.5, 0.2),
nrow = 2,
dimnames = list(c("(Intercept)", "x"), "y")
)),
scale = list(beta = matrix(
log(0.5),
nrow = 1,
dimnames = list("(Intercept)", "y")
))
)
sim <- simulate_data(
n = 200,
seed = 1,
generators = list(
x = gen_mvn("x", fixed_intercept = 0, residual_cov = 1),
y = gen_outcome(y ~ x, scale = sigma ~ 1, params = params)
)
)
analysis <- prep_sim_analysis(sim)
analysis$truth
fit <- brms::brm(
analysis$formula,
data = analysis$data,
backend = "cmdstanr",
refresh = 0
)
sim_recovery(fit, analysis)
}
#> Error: CmdStan path has not been set yet. See ?set_cmdstan_path.
# }