Skip to contents

Data Simulation

Source: vignettes/F-data-simulation.Rmd
F-data-simulation.Rmd

NOTE THIS IS EXPERIMENTAL FOR NOW

library(multilevelcoda)
library(knitr)
library(data.table)
#> data.table 1.18.4 using 9 threads (see ?getDTthreads).  Latest news: r-datatable.com
#> 
#> Attaching package: 'data.table'
#> 
#> The following object is masked from 'package:base':
#> 
#>     %notin%
library(cmdstanr)
#> This is cmdstanr version 0.9.0
#> - CmdStanR documentation and vignettes: mc-stan.org/cmdstanr
#> - CmdStan path: /Users/wileyj/.cmdstan/cmdstan-2.39.0
#> - CmdStan version: 2.39.0
library(brms)
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.23.0). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
#> 
#> Attaching package: 'brms'
#> 
#> The following object is masked from 'package:stats':
#> 
#>     ar
library(JWileymisc)
#> Registered S3 method overwritten by 'lme4':
#>   method           from
#>   na.action.merMod car
library(ggplot2)

options(digits = 3, pillar.sigfig = 3)

This vignette shows how to use the data simulation features that are part of the multilevelcoda package, particularly the simulate_data() function. Data simulation is useful for multiple purposes, such as:

  • Power analysis and sample size planning.
  • Method development and testing.
  • Generating synthetic data to test models or write analysis code before data collection finishes.

The simulate_data() function is designed to support these efforts. Broadly, simulate_data() builds an indexed design, runs named generators in order, and returns an mlsim_data object with major pieces:

  • data: the generated data.table.
  • metadata: design metadata such as group sizes, time indexing, and generator order.
  • generator_specs: the generator settings.
  • generator_metadata: e.g., parameters, random effects, covariance matrices, helper columns.

The first part of this vignette covers study designs and the predictor generators. The second part turns to model-based outcomes: we simulate outcomes from known parameters, fit the matching models, and check how well the true values are recovered.

Note that while part of the aspirations of data simulation here are to support things like Monte Carlo power analyses and simulation studies, these remain difficult for two reasons. Firstly, specifying all necessary parameters is not trivial, and it is our experience as researchers that such information is unlikely to be available prior to data collection. Secondly, Bayesian multilevel models are rather computationally demanding. Many random effects are possible, even sensible, and so at least in 2026, the computational demands of simulating and running hundreds of models exceeds most local devices and would likely require access to a high performance server.

Features at a Glance

The following tables summarise current features.

Design options

Design options control the highest level structure of the simulated data.

Option Use
n Total rows for a single-level design, or total rows to divide across groups if multilevel.
n_groups Number of groups in a grouped or multilevel design.
n_per_group Group sizes: scalar, vector, function, or count-distribution list (the latter support unequal or variable sample sizes per group as is common in real data).
group_id Name of the grouping column.
obs_id Name of the within-design observation index.
time_id Optional time column name.
time_values Optional time vector or function, including Date/POSIXct values.
time_truncate Whether shared time values are truncated for shorter groups.
seed Reproducible simulation seed that restores the caller RNG state.

Generator families

Generators add simulated data columns. Each generator draws its own values. However, gen_outcome() can generate outcomes that depend on columns created by earlier generators through the simulation context. The current generators are listed below.

Generator Role
gen_mvn() Univariate and multivariate Gaussian predictors, including ILR compositions and back-transforms to original compositional scale.
gen_categorical() Binary, unordered categorical, and ordered categorical variables.
gen_binomial() Binomial counts with logit-scale multilevel support.
gen_poisson() Poisson counts with log-scale multilevel support.
gen_negbin() Negative-binomial counts with optional location-scale multilevel support.
gen_gamma() Positive continuous variables with optional location-scale multilevel support.
gen_beta() Bounded (0, 1) continuous variables with optional location-scale support.
gen_custom() User-supplied generators that receive the active simulation context.
gen_outcome() Gaussian, dynamic VAR(1), ILR compositional, and GLM-family (Poisson, binomial, negative binomial, gamma, beta) outcomes generated from prior predictors.

Advanced simulation options

These are some of the features supported by generators. Not all features are necessarily supported by all generators.

Feature Use
level Choose row-level, group-level, or multilevel generation.
fixed_intercept Set the fixed location or link-scale intercept in any predictor generator.
random_cov Add group-specific random intercepts; a joint location-scale matrix also adds scale random effects.
residual_cov Set row-level residual variation for Gaussian and MVN generators.
scale_fixed_intercept Enable a scale model and set its fixed log-scale intercept (log residual SD, size, shape, or precision).
compositional, parts, sbp, total, keep_ilr Generate compositions from ILR coordinates and store SBP metadata.
between(), within(), ar1() Build outcome models from labelled predictor components (including ILR coordinates of compositional predictors) and residual VAR(1) dynamics.

Study Designs: Groups, Sizes, and Time

To start, let’s look at a simple, single level example. Here we generate one normally distributed variable, one categorical variable with three levels, and one binary categorical variable. Each variable in this first example is generated independently of the others.

single <- simulate_data(
  n = 20,
  seed = 2026,
  generators = list(
    x = gen_mvn("x", fixed_intercept = 0, residual_cov = 1),
    edu = gen_categorical("edu", categories = c("< Bachelor", "Bachelor", "> Bachelor"),
      fixed_intercept = c("Bachelor" = 0, "> Bachelor" = log(4 / 3))),
    tx = gen_categorical("tx", categories = c("control", "treatment"), fixed_intercept = stats::qlogis(.5))
  )
)

kable(single$data)
obs_id x edu tx
1 0.521 Bachelor control
2 -1.080 > Bachelor treatment
3 0.139 Bachelor treatment
4 -0.085 Bachelor control
5 -0.667 < Bachelor control
6 -2.516 < Bachelor control
7 -0.735 > Bachelor treatment
8 -1.020 < Bachelor control
9 0.114 Bachelor treatment
10 -0.474 < Bachelor treatment
11 -0.408 > Bachelor treatment
12 -0.730 > Bachelor control
13 -0.221 > Bachelor treatment
14 -0.226 > Bachelor treatment
15 -2.547 Bachelor control
16 1.347 > Bachelor treatment
17 0.616 > Bachelor treatment
18 0.218 < Bachelor treatment
19 -0.805 > Bachelor control
20 0.690 < Bachelor treatment 
single
#> <mlsim_data>
#>   rows: 20
#>   columns: 4 (generated: 3)
#>   seed: 2026
#>   grouping: none
#>   generators: 3
#> 
#> Generators:
#>  generator distribution  level vars
#>          x          mvn single    x
#>        edu  categorical single  edu
#>         tx  categorical single   tx

Printing an mlsim_data object gives an overview of the design. If we want the design and generator metadata as data.table objects, we can use summary().

summary(single)$design
#>        n n_cols n_generated_cols n_groups group_id group_size_min group_size_median
#>    <int>  <int>            <int>    <int>   <char>          <int>             <num>
#> 1:    20      4                3       NA     <NA>             NA                NA
#>    group_size_max obs_id time_id  seed n_generators
#>             <int> <char>  <char> <int>        <int>
#> 1:             NA obs_id    <NA>  2026            3
summary(single)$generators
#>    generator distribution  level   vars n_vars parameter_level parameter_count
#>       <char>       <char> <char> <char>  <int>          <char>           <int>
#> 1:         x          mvn single      x      1             row              20
#> 2:       edu  categorical single    edu      1             row              20
#> 3:        tx  categorical single     tx      1             row              20
#>    has_row_parameters has_group_parameters has_fixed_parameters has_random_cov
#>                <lgcl>               <lgcl>               <lgcl>         <lgcl>
#> 1:               TRUE                FALSE                 TRUE          FALSE
#> 2:               TRUE                FALSE                 TRUE          FALSE
#> 3:               TRUE                FALSE                 TRUE          FALSE
#>    has_random_effects has_residuals has_scale_model has_composition has_custom_output
#>                <lgcl>        <lgcl>          <lgcl>          <lgcl>            <lgcl>
#> 1:              FALSE          TRUE           FALSE           FALSE             FALSE
#> 2:              FALSE         FALSE           FALSE           FALSE             FALSE
#> 3:              FALSE         FALSE           FALSE           FALSE             FALSE

We can generate multilevel data as well. A single scalar n_per_group creates a balanced grouped design, with exactly the same number of observations per group. To show it more easily, we make just one variable.

balanced <- simulate_data(
  n_groups = 3,
  n_per_group = 4,
  seed = 2026,
  generators = list(
    x = gen_mvn("x", fixed_intercept = 0, residual_cov = 1)
  )
)

kable(balanced$data)
group_id obs_id x
1 1 0.521
1 2 -1.080
1 3 0.139
1 4 -0.085
2 1 -0.667
2 2 -2.516
2 3 -0.735
2 4 -1.020
3 1 0.114
3 2 -0.474
3 3 -0.408
3 4 -0.730

However, while we generate the labels for multilevel data, the variable x was still single level, Gaussian data. To address this, we speicfy the parameters of a random intercept only multilevel model, which is used to generate matching data. The fixed_intercept is the overall mean, random_cov is the group-level variance, and residual_cov is the residual variance. The level argument controls where a variable varies: level = "multilevel" gives row-level values with group-specific random effects, and level = "level2" gives one value per group, constant within groups but varying between them. The same approach broadly applies to other generator families.

Note that in the data output are the simulated data as well as the “true” between and within values. These are useful because in a multilevel structure, with a small, finite number of samples for a participant, the sample mean can differ from that person’s true mean. These values instead show the “true” values. This can be important for subsequent simulations (e.g., outcome generation where we would want it generated from the truth, not from the sample).

balanced2 <- simulate_data(
  n_groups = 3,
  n_per_group = 4,
  seed = 2026,
  generators = list(
    x = gen_mvn("x", level = "multilevel",
      fixed_intercept = 0, random_cov = 1, residual_cov = .01),
    y = gen_mvn("y", level = "level2", fixed_intercept = 0, residual_cov = 1)
  )
)

kable(balanced2$data)
group_id obs_id x x_between x_within .mlsim_x_random_intercept y
1 1 0.512 0.521 -0.008 0.521 1.347
1 2 0.454 0.521 -0.067 0.521 1.347
1 3 0.269 0.521 -0.252 0.521 1.347
1 4 0.447 0.521 -0.074 0.521 1.347
2 1 -1.182 -1.080 -0.102 -1.080 0.616
2 2 -1.068 -1.080 0.011 -1.080 0.616
2 3 -1.127 -1.080 -0.047 -1.080 0.616
2 4 -1.121 -1.080 -0.041 -1.080 0.616
3 1 0.066 0.139 -0.073 0.139 0.218
3 2 0.117 0.139 -0.022 0.139 0.218
3 3 0.117 0.139 -0.023 0.139 0.218
3 4 -0.115 0.139 -0.255 0.139 0.218

If we specify n_per_group as a vector, we can get an unbalanced design.

unbalanced <- simulate_data(
  n_groups = 4,
  n_per_group = c(3, 5, 4, 2),
  seed = 2026,
  generators = list(
    x = gen_mvn("x", level = "multilevel",
      fixed_intercept = 0, random_cov = 1, residual_cov = .01),
    y = gen_mvn("y", level = "level2", fixed_intercept = 0, residual_cov = 1)
  )
)

kable(unbalanced$data)
group_id obs_id x x_between x_within .mlsim_x_random_intercept y
1 1 0.454 0.521 -0.067 0.521 -0.805
1 2 0.269 0.521 -0.252 0.521 -0.805
1 3 0.447 0.521 -0.074 0.521 -0.805
2 1 -1.182 -1.080 -0.102 -1.080 0.690
2 2 -1.068 -1.080 0.011 -1.080 0.690
2 3 -1.127 -1.080 -0.047 -1.080 0.690
2 4 -1.121 -1.080 -0.041 -1.080 0.690
2 5 -1.153 -1.080 -0.073 -1.080 0.690
3 1 0.117 0.139 -0.022 0.139 -0.329
3 2 0.117 0.139 -0.023 0.139 -0.329
3 3 -0.115 0.139 -0.255 0.139 -0.329
3 4 0.274 0.139 0.135 0.139 -0.329
4 1 -0.023 -0.085 0.062 -0.085 -0.165
4 2 -0.063 -0.085 0.022 -0.085 -0.165

Group sizes can also be drawn from a count distribution. We give a list with the distribution name and its parameters, along with optional minimum and maximum sizes. This is an easier way to specify that there may be differences without having to manually specify the number for each participant.

drawn_sizes <- simulate_data(
  n_groups = 5,
  n_per_group = list(
    distribution = "poisson",
    lambda = 4,
    minimum = 2,
    maximum = 6
  ),
  seed = 2026,
  generators = list(
    x = gen_mvn("x", level = "multilevel",
      fixed_intercept = 0, random_cov = 1, residual_cov = .01),
    y = gen_mvn("y", level = "level2", fixed_intercept = 0, residual_cov = 1)
  )
)

kable(drawn_sizes$data)
group_id obs_id x x_between x_within .mlsim_x_random_intercept y
1 1 -1.994 -1.958 -0.037 -1.958 0.339
1 2 -2.261 -1.958 -0.304 -1.958 0.339
1 3 -2.169 -1.958 -0.211 -1.958 0.339
1 4 -1.998 -1.958 -0.040 -1.958 0.339
1 5 -2.112 -1.958 -0.154 -1.958 0.339
2 1 0.995 1.085 -0.089 1.085 1.207
2 2 0.977 1.085 -0.108 1.085 1.207
2 3 1.111 1.085 0.026 1.085 1.207
2 4 0.992 1.085 -0.093 1.085 1.207
3 1 -0.010 0.204 -0.214 0.204 0.571
3 2 0.185 0.204 -0.019 0.204 0.571
4 1 0.495 0.501 -0.006 0.501 -1.686
4 2 0.488 0.501 -0.013 0.501 -1.686
4 3 0.616 0.501 0.115 0.501 -1.686
5 1 1.015 1.030 -0.015 1.030 -0.838
5 2 0.848 1.030 -0.182 1.030 -0.838
5 3 0.976 1.030 -0.054 1.030 -0.838
5 4 0.958 1.030 -0.072 1.030 -0.838

Longitudinal designs can add a time index. When groups have different lengths, time_truncate = TRUE uses the first values for shorter groups. Adding a time index may not matter for some designs, but it is useful for creating lags and similar derived variables later.

dated <- simulate_data(
  n_groups = 3,
  n_per_group = c(5, 3, 4),
  time_id = "date",
  time_values = as.Date("2026-01-01") + 0:4,
  seed = 2026,
  generators = list(
    x = gen_mvn("x", level = "multilevel",
      fixed_intercept = 0, random_cov = 1, residual_cov = .01),
    y = gen_mvn("y", level = "level2", fixed_intercept = 0, residual_cov = 1)
  )
)

kable(dated$data)
group_id obs_id date x x_between x_within .mlsim_x_random_intercept y
1 1 2026-01-01 0.512 0.521 -0.008 0.521 1.347
1 2 2026-01-02 0.454 0.521 -0.067 0.521 1.347
1 3 2026-01-03 0.269 0.521 -0.252 0.521 1.347
1 4 2026-01-04 0.447 0.521 -0.074 0.521 1.347
1 5 2026-01-05 0.419 0.521 -0.102 0.521 1.347
2 1 2026-01-01 -1.068 -1.080 0.011 -1.080 0.616
2 2 2026-01-02 -1.127 -1.080 -0.047 -1.080 0.616
2 3 2026-01-03 -1.121 -1.080 -0.041 -1.080 0.616
3 1 2026-01-01 0.066 0.139 -0.073 0.139 0.218
3 2 2026-01-02 0.117 0.139 -0.022 0.139 0.218
3 3 2026-01-03 0.117 0.139 -0.023 0.139 0.218
3 4 2026-01-04 -0.115 0.139 -0.255 0.139 0.218

Distribution Families

All predictor generators use fixed location or link-scale intercepts, and the multilevel forms add random covariance terms. Here is a quick demonstration of different distribution families that we support with single level data.

families <- simulate_data(
  n = 10,
  seed = 2026,
  generators = list(
    normal = gen_mvn("normal", fixed_intercept = 5, residual_cov = 1),
    successes = gen_binomial("successes", size = 5, fixed_intercept = stats::qlogis(0.40)),
    visits = gen_poisson("visits", fixed_intercept = log(2)),
    events = gen_negbin("events", fixed_intercept = log(2), scale_fixed_intercept = log(3)),
    cost = gen_gamma("cost", fixed_intercept = log(20), scale_fixed_intercept = log(4)),
    adherence = gen_beta("adherence", fixed_intercept = stats::qlogis(0.70), scale_fixed_intercept = log(12))
  )
)

kable(families$data)
obs_id normal successes visits events cost adherence
1 5.52 2 4 1 10.63 0.875
2 3.92 2 1 0 36.47 0.822
3 5.14 1 3 5 38.06 0.657
4 4.92 0 0 4 41.29 0.664
5 4.33 2 2 2 33.78 0.500
6 2.48 1 2 4 10.91 0.766
7 4.26 2 1 0 24.79 0.897
8 3.98 1 2 2 13.94 0.841
9 5.11 0 3 1 14.41 0.568
10 4.53 2 2 2 8.27 0.846

Here is an example of a multilevel count generator.

link_scale_counts <- simulate_data(
  n_groups = 3,
  n_per_group = 5,
  seed = 2026,
  generators = list(
    count = gen_poisson(
      "count",
      level = "multilevel",
      fixed_intercept = log(2.5),
      random_cov = 0.20
    )
  )
)

link_scale_counts
#> <mlsim_data>
#>   rows: 15
#>   columns: 4 (generated: 2)
#>   seed: 2026
#>   grouping: group_id (3 groups; size min/median/max: 5/5/5)
#>   generators: 1
#> 
#> Generators:
#>  generator distribution      level  vars
#>      count      poisson multilevel count

Categorical Variables

gen_categorical() supports binary variables, labeled categories, ordered factors, character output, and integer codes.

categorical <- simulate_data(
  n = 12,
  seed = 2026,
  generators = list(
    binary_factor = gen_categorical("binary_factor", fixed_intercept = stats::qlogis(0.35)),
    education = gen_categorical(
      "education",
      categories = c("< bachelor", "bachelor", "> bachelor"),
      fixed_intercept = c("bachelor" = log(0.35 / 0.45), "> bachelor" = log(0.20 / 0.45)),
      ordered = TRUE
    ),
    education_code = gen_categorical(
      "education_code",
      categories = c("arts", "science", "psychology"),
      fixed_intercept = c("science" = log(0.35 / 0.45), "psychology" = log(0.20 / 0.45)),
      output = "integer"
    )
  )
)

kable(categorical$data)
obs_id binary_factor education education_code
1 1 < bachelor 0
2 0 > bachelor 0
3 0 < bachelor 0
4 0 < bachelor 0
5 0 bachelor 0
6 0 < bachelor 1
7 0 < bachelor 2
8 1 < bachelor 0
9 0 < bachelor 1
10 0 < bachelor 0
11 0 < bachelor 1
12 1 < bachelor 0

Multilevel categorical generators use baseline-category logits and optional group random effects. With two categories, this is effectively a multilevel Bernoulli model. With more than two categories, it is a multilevel multinomial model with the first category as the baseline. The random intercepts can have their own variances and covariances, and a full random effect covariance matrix must be specified. Use 0s on the off diagonal for uncorrelated random intercepts.

categorical_ml <- simulate_data(
  n_groups = 4,
  n_per_group = 4,
  seed = 2026,
  generators = list(
    education = gen_categorical(
      "education",
      level = "multilevel",
      categories = c("< bachelor", "bachelor", "> bachelor"),
      fixed_intercept = c("bachelor" = 0.15, "> bachelor" = -0.55),
      random_cov = matrix(c(0.20, 0.04, 0.04, 0.12), nrow = 2),
      ordered = TRUE
    )
  )
)

kable(categorical_ml$data)
group_id obs_id education .mlsim_education_random_intercept_bachelor .mlsim_education_random_intercept_X..bachelor
1 1 bachelor -0.306 0.105
1 2 < bachelor -0.306 0.105
1 3 < bachelor -0.306 0.105
1 4 < bachelor -0.306 0.105
2 1 bachelor 0.155 0.940
2 2 bachelor 0.155 0.940
2 3 < bachelor 0.155 0.940
2 4 < bachelor 0.155 0.940
3 1 bachelor -0.150 0.194
3 2 < bachelor -0.150 0.194
3 3 bachelor -0.150 0.194
3 4 < bachelor -0.150 0.194
4 1 < bachelor -0.089 0.318
4 2 bachelor -0.089 0.318
4 3 > bachelor -0.089 0.318
4 4 < bachelor -0.089 0.318

Correlated and Compositional Predictors

gen_mvn() simulates correlated blocks of variables. With compositional = TRUE, the MVN columns are interpreted as ILR coordinates and back-transformed into compositional parts. Here we are back to single level data for the first illustration. With MVN data, we need to specify the residual covariance matrix, which is used to control any correlation.

correlated_comp <- simulate_data(
  n = 8,
  seed = 2026,
  generators = list(
    affect_block = gen_mvn(
      c("affect", "energy"),
      fixed_intercept = c(0, 1),
      residual_cov = matrix(c(1.0, 0.5, 0.5, 1.2), nrow = 2)
    ),
    time_use = gen_mvn(
      c("ilr_1", "ilr_2"),
      fixed_intercept = c(0.4, -1.5),
      residual_cov = diag(c(0.05, 0.05)),
      compositional = TRUE,
      parts = c("sleep", "active", "sedentary"),
      total = 24,
      keep_ilr = TRUE
    )
  )
)

kable(correlated_comp$data)
obs_id affect energy ilr_1 ilr_2 sleep active sedentary
1 0.351 1.566 0.389 -1.36 8.38 1.99 13.6
2 -0.587 -0.290 -0.027 -1.45 5.64 2.09 16.3
3 0.355 0.938 0.013 -1.68 5.30 1.59 17.1
4 0.366 0.561 0.387 -1.35 8.41 2.02 13.6
5 -0.405 0.238 0.256 -1.57 6.93 1.67 15.4
6 -1.890 -1.579 0.014 -1.54 5.65 1.87 16.5
7 0.922 -0.962 0.518 -1.81 7.86 1.16 15.0
8 -1.621 0.655 0.363 -1.17 8.73 2.44 12.8 
correlated_comp$generator_metadata$time_use$sbp
#>           sleep active sedentary
#> balance_1     1     -1        -1
#> balance_2     0      1        -1
correlated_comp$generator_metadata$time_use$ilr_coordinate_map
#>       ilr sbp_row positive_parts   negative_parts
#>    <char>   <int>         <AsIs>           <AsIs>
#> 1:  ilr_1       1          sleep active,sedentary
#> 2:  ilr_2       2         active        sedentary

If we only want the parts in the output, we can set keep_ilr = FALSE.

parts_only <- simulate_data(
  n = 5,
  seed = 2026,
  generators = list(
    time_use = gen_mvn(
      c("z1", "z2"),
      fixed_intercept = c(0.4, -1.5),
      residual_cov = diag(c(0.05, 0.05)),
      compositional = TRUE,
      parts = c("sleep", "active", "sedentary"),
      total = 24,
      keep_ilr = FALSE
    )
  )
)

kable(parts_only$data)
obs_id sleep active sedentary
1 12.41 1.44 10.2
2 8.38 1.23 14.4
3 9.70 1.59 12.7
4 7.83 1.69 14.5
5 8.29 1.39 14.3

Location-Scale Predictors

scale_fixed_intercept enables a scale model and sets its fixed log-scale intercept. For a univariate gen_mvn() generator, the scale is the log residual standard deviation. For the negative binomial, gamma, and beta generators, it is the size, shape, and precision, respectively.

On its own, scale_fixed_intercept gives every group the same scale. To make groups differ not only in their average level (location) but also in their variability (scale), we supply a joint location-scale random_cov with dimension twice the number of variables, ordered with all location intercepts first and all scale intercepts second. A location-only random_cov leaves the scale constant between groups even when scale_fixed_intercept is set. The generators below are univariate, so their 2x2 random_cov matrices are the joint form: one location row and one scale row each, and both location and scale vary by group.

There are two places to inspect the group-varying scale that was actually generated: the realized scale parameters in generator_metadata (shown below), and the group-level scale draws themselves, which appear in sim$data as .mlsim_<variable>_scale_random_intercept columns.

location_scale <- simulate_data(
  n_groups = 4,
  n_per_group = 5,
  seed = 2026,
  generators = list(
    symptom = gen_mvn(
      "symptom",
      level = "multilevel",
      fixed_intercept = 10,
      scale_fixed_intercept = log(1.2),
      random_cov = matrix(
        c(1.00, 0.15,
          0.15, 0.08),
        nrow = 2
      )
    ),
    event_count = gen_negbin(
      "event_count",
      level = "multilevel",
      fixed_intercept = log(2),
      scale_fixed_intercept = log(6),
      random_cov = matrix(
        c(0.20, 0.03,
          0.03, 0.06),
        nrow = 2
      )
    )
  )
)

kable(data.table(
  group_id = location_scale$data$group_id,
  symptom_residual_sd = location_scale$generator_metadata$symptom$residual_sd[, "symptom"],
  event_count_size = location_scale$generator_metadata$event_count$size
))
group_id symptom_residual_sd event_count_size
1 1.29 5.99
1 1.29 5.99
1 1.29 5.99
1 1.29 5.99
1 1.29 5.99
2 2.57 4.75
2 2.57 4.75
2 2.57 4.75
2 2.57 4.75
2 2.57 4.75
3 1.39 6.04
3 1.39 6.04
3 1.39 6.04
3 1.39 6.04
3 1.39 6.04
4 1.54 4.54
4 1.54 4.54
4 1.54 4.54
4 1.54 4.54
4 1.54 4.54

Custom Generators

Sometimes the built-in generators are not enough. Custom generators receive the simulation context and return generated data, column names, and optional metadata. They can use the active design columns, previously generated variables, and any extra arguments supplied to gen_custom(). Their metadata is captured in the same way as for the built-in generators.

Below is an examplethat can be useful when trying to simulate synthetic data for use in building or testing data analysis code. For example, this could support pre-registration, registered reports where data analysis code is written in advance, or power analyses where there is pilot data on variables that follow distributions that do not conveniently match any of the currently supported distributions. We use density_inversion() function which first calculates the density of a set of values passed in, and then inverts to generate synthetic data with approximately the same density. We apply it to the iris data built into R for which the petal lengths are famously bimodal. While not exact, we can see that the simulated data broadly captures and reflects the petal length data present in the iris dataset. You may want more metadata saved, but here we just save the source (density inversion) and the number of pilot samples passed in.

custom_sampler <- function(context, vars, level, pilot_values) {
  values <- density_inversion(pilot_values, n = context$n_rows,
    KDEn = 50)

  list(
    data = data.frame(values = values),
    names = vars,
    metadata = list(source = "density inversion", n_pilot = length(pilot_values))
  )
}

custom <- simulate_data(
  n = 500,
  seed = 2026,
  generators = list(
    custom = gen_custom(
      "custom",
      generator = custom_sampler,
      pilot_values = iris$Petal.Length
    )
  )
)

dcompare <- data.table(
  value = c(iris$Petal.Length, custom$data$custom),
  group = rep(c("original", "simulated"), times = c(
    length(iris$Petal.Length),
    length(custom$data$custom)
  )))

ggplot(dcompare, aes(value, colour = group, linetype = group)) +
  geom_density(linewidth = 1) + theme_classic()
plot of chunk custom-generators

plot of chunk custom-generators

A First Outcome Model and Parameter Recovery

So far, we have only simulated predictors. gen_outcome() adds a model-based outcome generator to the same simulate_data() workflow.

The rest of this vignette focuses on the workflow that builds on it: simulate data with known parameters, prepare an analysis data set, fit the matching model, and check that the true values are recovered.

Before the more complex models in the later sections, here is the complete loop on a simple example: a Gaussian outcome with one predictor and a random intercept. We need to specify all of the parameters and parameter matrices. The location matrix has one row per term and one column per outcome, the scale matrix holds the log residual standard deviation, and the random-effect covariance uses structured names to label which parameter each random effect belongs to. Even in this simple example, the number of parameters and the specific structure to hold them is rather complicated. The next section will show how gen_template() can create these structures for us, which is probably easier for most use cases.

first_random_name <- "location|outcome=y|term=(Intercept)"
first_params <- list(
  location = list(beta = matrix(
    c(2, 0.5),
    nrow = 2,
    dimnames = list(c("(Intercept)", "x"), "y")
  )),
  scale = list(beta = matrix(
    log(0.8),
    nrow = 1,
    dimnames = list("(Intercept)", "y")
  )),
  random = list(ID = list(covariance = matrix(
    0.3^2,
    dimnames = list(first_random_name, first_random_name)
  )))
)

first_sim <- simulate_data(
  n_groups = 200,
  n_per_group = 10,
  group_id = "ID",
  seed = 2026,
  generators = list(
    x = gen_mvn("x", fixed_intercept = 0, residual_cov = 1),
    y = gen_outcome(y ~ x + (1 | ID), scale = sigma ~ 1, params = first_params)
  )
)

kable(head(first_sim$data))
ID obs_id x y
1 1 0.521 1.76
1 2 -1.080 2.40
1 3 0.139 1.11
1 4 -0.085 1.40
1 5 -0.667 1.59
1 6 -2.516 1.57

prep_sim_analysis() turns the simulated data into an analysis-ready data set and infers the matching brms formula. It also returns the true parameter values in $truth, labeled with the parameter names the fitted model will use. The goal of this is to support cases where we want to test if the true value was recovered (e.g., simulation studies on bias, coverage rates).

first_analysis <- prep_sim_analysis(first_sim)
first_analysis$formula
#> y ~ x + (1 | ID) 
#> sigma ~ 1

There is no composition in this example, so we fit the model directly with brms.

first_fit <- brm(
  formula = first_analysis$formula,
  data = first_analysis$data,
  backend = "cmdstanr",
  seed = 2026,
  chains = 2,
  cores = 2,
  iter = 1000,
  refresh = 0,
  silent = 2
)

This object is a regular brms object so we could summary, predict, etc However, there are some more helper functions that may be useful here. sim_recovery() joins the fitted estimates with the simulation truth by parameter name, and adds bias and interval-coverage columns.

first_recovery <- sim_recovery(first_fit, first_analysis)
kable(first_recovery[, !"simulator_name"], digits = 3)
type group parameter truth estimate est_error lower upper bias covered
fixed NA Intercept 2.000 1.977 0.029 1.918 2.035 -0.023 TRUE
fixed NA x 0.500 0.495 0.020 0.455 0.533 -0.005 TRUE
fixed NA sigma_Intercept -0.223 -0.218 0.017 -0.251 -0.185 0.005 TRUE
random_sd ID Intercept 0.300 0.304 0.026 0.259 0.356 0.004 TRUE

Each row is one model parameter. The truth column holds the generating value, bias is the difference between the posterior estimate and the truth, and covered indicates whether the truth falls inside the credible interval. In this simple model the analysis matches the data-generating process, so each fitted parameter estimates exactly the quantity the simulator recorded. That will not always hold: as we get into time series and multilevel predictors below, the analysis we fit no longer matches the simulation model exactly, and some parameters then estimate a related quantity instead.

Dynamic Outcome Simulation

A major use case we are interested in for multilevelcoda are intensive longitudinal studies, such as wearables or diet or time use measured in people across consecutive days, for example for 1-2 weeks. In these cases, the data are multilevel, but importantly, they also are time series. Even if not all analyses are fundamentally focused on or even include a time series component, the underlying data would or could. To support that and simulate data that ideally match real data as much as feasible, we support dynamic outcome simulation, or outcomes that are assumed to be time series.

gen_outcome() supports residual autoregressive dynamics and multivariate compositional outcomes, which we combine in this section. As before, we focus on the default family = "gaussian".

For Gaussian outcomes, the scale model is required. Use scale = sigma ~ 1 for a constant conditional standard deviation. Dynamics are supported through an autoregressive component, implemented using ar1(). When ar1() is present, the scale model controls the variability of the innovations. When ar1() is absent, it controls ordinary residual variability. Note that the AR process is applied to the residual ILR states, not to lagged observed outcomes. ar1() is currently only supported for family = "gaussian". This is because building models with residual autoregression and innovations becomes substantially more complicated for categorical data and other not Gaussian distributions.

When ar1() is present, burnin is required. Each series starts with zero residuals and is iterated burnin steps before the first observed row. Do not use burnin = 0 or other small values with ar1(). If the burn-in is too short, the first observations will be under-dispersed relative to the stationary distribution of the AR process. Choose a burn-in long enough for the process to forget its start. A simple decision rule can be at least log(0.01) / log(rho) steps, where rho is the largest spectral radius of the AR matrices, which works out to about 90 steps at rho = 0.95. Models without ar1() have no burn-in phase, so burnin can simply be omitted.

In the first outcome example, we wrote the parameter matrices by hand. However, that is quite difficult with many names. To help with this, gen_template() can be used. Basically, it is a stand-in for gen_outcome() which does not require parameters to be specified but will generate the full structure so we can extract the necessary parameter matrices. We run the template with the same formulas, grouping structure, and composition settings that we plan to use in the final simulation. The design size does not need to match: a small template design is enough to obtain the parameter structure, which we then fill in and reuse for the full-size simulation.

Here we see our first dynamic outcome formula. The setup here is a 3 component compositional dyanmic outcome. Much of the formula is a “standard” multivariate, multilevel location scale formula specification for brms. What is new is the use of ar1(). This gets parsed by multilevelcoda and automatically expanded to mean the lag1 outcome variable, and this can be done as fixed and random effects, as done below. Now once we have the template, we can extract the relevant parameter matrices and then specify the parameter values we want to use for the simulation. This saves us manually writing / figuring out the structure and the very specific names needed for simulate_data().

dynamic_outcome_formula <- mvbind(ilr1, ilr2) ~ ar1() + (1 + ar1() | ID)
dynamic_scale_formula <- sigma ~ 1 + (1 | ID)
dynamic_composition <- list(parts = c("sleep", "activity", "sedentary"), total = 24)

dynamic_template <- simulate_data(
  n_groups = 4,
  n_per_group = 5,
  group_id = "ID",
  time_id = "day",
  seed = 2026,
  generators = list(
    outcome_template = gen_template(
      dynamic_outcome_formula,
      scale = dynamic_scale_formula,
      composition = dynamic_composition,
      burnin = 20
    )
  )
)

dynamic_params <- dynamic_template$generator_metadata$outcome_template$params
dynamic_params$location$beta["(Intercept)", ] <- c(0.2, -1.5)
dynamic_params$scale$beta["(Intercept)", ] <- log(c(0.2, 0.15))
dynamic_params$ar$beta["ar1()", , ] <- matrix(c(0.50, 0.10, 0.10, 0.30), 2, 2, byrow = TRUE)

dynamic_random_names <- rownames(dynamic_params$random$ID$covariance)
dynamic_random_sd_values <- c(
  "location|outcome=ilr1|term=(Intercept)" = 0.20,
  "location|outcome=ilr2|term=(Intercept)" = 0.18,
  "ar|term=ar1()|to=ilr1|from=ilr1" = 0.02,
  "ar|term=ar1()|to=ilr1|from=ilr2" = 0.01,
  "ar|term=ar1()|to=ilr2|from=ilr1" = 0.01,
  "ar|term=ar1()|to=ilr2|from=ilr2" = 0.02,
  "scale|outcome=ilr1|term=(Intercept)" = 0.16,
  "scale|outcome=ilr2|term=(Intercept)" = 0.14
)
dynamic_random_sd <- dynamic_random_sd_values[dynamic_random_names]

## we initialise all the random effect correlations at 0
dynamic_random_cor <- diag(length(dynamic_random_names))
dimnames(dynamic_random_cor) <- list(dynamic_random_names, dynamic_random_names)
set_dynamic_random_cor <- function(x, y, value) {
  dynamic_random_cor[x, y] <<- value
  dynamic_random_cor[y, x] <<- value
}
## we set a few select correlations to non zero values
set_dynamic_random_cor(
  "location|outcome=ilr1|term=(Intercept)",
  "ar|term=ar1()|to=ilr1|from=ilr1",
  0.20
)
set_dynamic_random_cor(
  "ar|term=ar1()|to=ilr1|from=ilr1",
  "ar|term=ar1()|to=ilr1|from=ilr2",
  0.15
)
set_dynamic_random_cor(
  "location|outcome=ilr2|term=(Intercept)",
  "ar|term=ar1()|to=ilr2|from=ilr1",
  0.20
)
set_dynamic_random_cor(
  "ar|term=ar1()|to=ilr2|from=ilr1",
  "ar|term=ar1()|to=ilr2|from=ilr2",
  0.15
)

if (anyNA(dynamic_random_sd) ||
    min(eigen(dynamic_random_cor, symmetric = TRUE, only.values = TRUE)$values) <= 0) {
  stop("The dynamic outcome random-effect covariance is not valid.")
}

dynamic_params$random$ID$covariance <- diag(dynamic_random_sd) %*%
  dynamic_random_cor %*%
  diag(dynamic_random_sd)
dimnames(dynamic_params$random$ID$covariance) <- list(dynamic_random_names, dynamic_random_names)

We use clearly visible dynamics for autoregression to make the centering comparison later in this section clearer. The fixed autoregressive (diagonal) coefficients are 0.50 and 0.30, and the cross-lagged (off-diagonal) coefficients are 0.10. The off diagonal coefficients could be zero, but keeping them probably makes some sense as it is plausible, depending on the sequential binary partition used, that one ILR would have some predictive utility to the subsequent value of other ILRs, especially as two ILRs commonly will include values from at least one common component. In any case, these values give the fixed AR matrix a spectral radius of about 0.54. A burn-in of 150 steps is far above the log(0.01) / log(rho) approximate rule. For simulating individual datasets, we do not see a reason not to put a very generous burnin because the computational cost, especially vis a vis the cost of fitting any analytic model to such complex data is negligible in our view. That said, it does add some time just to the data simulation, and we can see that if a very large number of people or hundreds or thousdands of datasets were being simulated, there may be some utility in striving for the “minimal” burnin needed. Another reason for our larger choice is that with random slopes, while the fixed matrix may have a low spectral radius, unless we had distinct burnins for different people, some participants, due to random effects, could have a larger spectral radius.

We keep the group-level AR standard deviations small: 0.02 for the autoregressive terms and 0.01 for the cross-lagged terms. Larger random slopes can push individual AR matrices toward the stationarity boundary.

If we get towards a stationary point, that will trigger resampling and truncate the random-effect distribution (because internally, the code will ensure that stationary matrices for any given participant is not used). Finally, we simulate 50 observations per person. Series length matters here and probably quite a bit. We simulate data based on a residual process. However, for analyses it is quite common to just use a lagged outcome value possibly person mean centred. In these cases, it becomes much more precise as the time series length increases, because otherwise, we are centering on something rather noisy.

dynamic_comp <- simulate_data(
  n_groups = 200,
  n_per_group = 50,
  group_id = "ID",
  time_id = "day",
  seed = 2026,
  generators = list(
    outcome = gen_outcome(
      dynamic_outcome_formula,
      scale = dynamic_scale_formula,
      params = dynamic_params,
      composition = dynamic_composition,
      burnin = 150
    )
  )
)

kable(head(dynamic_comp$data[, c("ID", "day", "ilr1", "ilr2", "sleep", "sedentary", "activity"), with = FALSE]))
ID day ilr1 ilr2 sleep sedentary activity
1 1 0.600 -1.43 9.62 12.7 1.68
1 2 0.461 -1.14 9.48 12.1 2.40
1 3 0.292 -1.32 7.86 14.0 2.16
1 4 0.230 -1.14 7.93 13.4 2.67
1 5 0.454 -1.36 8.82 13.2 1.93
1 6 0.521 -1.12 9.96 11.7 2.38

Before fitting anything, we can run a few quick checks on the simulation. For example, we can look at the realized maximum spectral radius across the individual AR matrices, which should be below 1. Here it is around 0.6, only a little above the fixed-matrix value of 0.54, because we kept the AR random effects small.

dynamic_comp$generator_metadata$outcome$ar$stability$max_spectral_radius_overall
#> [1] 0.596

Next, we prepare the generated data into an analysis data set that can be passed to brmcoda(). prep_sim_analysis() rebuilds a complr object using the same SBP basis as the simulator and translates ar1() into within-person centered, lag 1 ILR predictors. The fitted model includes random level terms, random inertia terms through the lagged ILR slopes, and random conditional variability through response-specific sigma models.

The lag columns are built by the exported lag_by_time() helper. These are time-based lags on time_id, so a row’s lag comes from the observation at time - time_step. This matters when working with real data that has skipped time points, such as a missed diary day. The row after a gap gets an NA lag. The step is inferred as the smallest positive within-person time difference, or it can be set explicitly with prep_sim_analysis(sim, time_step = ). The number of gap-affected rows is recorded in metadata$lag_gaps. lag_by_time() can also be used directly on observed data sets to construct lags for hand-built analysis models.

By default, prep_sim_analysis() emits brms ID-linked random effects, such as (1 | p1 | ID), whenever the same grouping factor appears in both the mean and the scale formulas. The simulator draws those group-level effects from one joint covariance, and the linked syntax lets the analysis model estimate their correlations. Here, that joins all eight random effects (mean levels, lagged slopes, and sigma intercepts for both responses) into a single correlated block, which matches the generative model. One practical caveat is that large linked random-effect blocks are demanding. For very large designs, we can opt out with prep_sim_analysis(sim, link_random = FALSE), which emits separate, uncorrelated random-effect blocks instead. Of course, if the simulation included correlations (any non zero true correlations), this will result in some degree of bias. Of course, even this may be interesting to study. The possible number of random effects grows very quickly. For example, with a 5 part composition, there are 4 ILR terms, resulting in: 4 location intercepts, 4 scale intercepts, 4 autoregressive slopes, and 12 cross-lagged slopes for a total of 24 parameters that could be all random effects and if so a 24 x 24 random effect covariance matrix would have 300 parameters (24 variances and 276 unique covariances). Such a large matrix would be extremely unstable with small sample sizes. Even with large sample sizes, it would likely be very computationally demanding. So an interesting simulation question could be: how much bias in the fixed effect autoregression estimates is introduced by forcing the random effect covariances to be zero?

In any case, back to our actual example, we can prepare the analysis and see the columns and formulae. All of the data and appropriate names for the data and formulae are generated automatically.

dynamic_analysis <- prep_sim_analysis(dynamic_comp, centering = "latent")
dynamic_analysis$metadata$lag_columns
#> [1] "lag_z1_1_latent" "lag_z2_1_latent"
dynamic_analysis$formula
#> z1_1 ~ lag_z1_1_latent + lag_z2_1_latent + (1 + lag_z1_1_latent + lag_z2_1_latent | p1 | ID) 
#> sigma ~ 1 + (1 | p1 | ID)
#> z2_1 ~ lag_z1_1_latent + lag_z2_1_latent + (1 + lag_z1_1_latent + lag_z2_1_latent | p1 | ID) 
#> sigma ~ 1 + (1 | p1 | ID)

These are large models, with 10,000 rows, a bivariate response, and a linked random-effect block. To speed things up, the fits below use brms within-chain threading on top of the two parallel chains, via threads = brms::threading(4, static = TRUE). Using static = TRUE keeps the results reproducible for a given number of threads.

dynamic_fit <- brmcoda(
  complr = dynamic_analysis$complr,
  formula = dynamic_analysis$formula,
  backend = "cmdstanr",
  seed = 2026,
  chains = 2,
  cores = 2,
  threads = brms::threading(4, static = TRUE),
  iter = 1000,
  refresh = 0
)
summary(dynamic_fit)
#> Warning: Parts of the model have not converged (some Rhats are > 1.05). Be careful when
#> analysing the results! We recommend running more iterations and/or setting stronger
#> priors.
#>  Family: MV(gaussian, gaussian) 
#>   Links: mu = identity; sigma = log
#>          mu = identity; sigma = log 
#> Formula: z1_1 ~ lag_z1_1_latent + lag_z2_1_latent + (1 + lag_z1_1_latent + lag_z2_1_latent | p1 | ID) 
#>          sigma ~ 1 + (1 | p1 | ID)
#>          z2_1 ~ lag_z1_1_latent + lag_z2_1_latent + (1 + lag_z1_1_latent + lag_z2_1_latent | p1 | ID) 
#>          sigma ~ 1 + (1 | p1 | ID)
#>    Data: complr$dataout (Number of observations: 9800) 
#>   Draws: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
#>          total post-warmup draws = 1000
#> 
#> Multilevel Hyperparameters:
#> ~ID (Number of levels: 200) 
#>                                              Estimate Est.Error l-95% CI u-95% CI Rhat
#> sd(z11_Intercept)                                0.20      0.01     0.18     0.21 1.01
#> sd(z11_lag_z1_1_latent)                          0.03      0.02     0.00     0.06 1.03
#> sd(z11_lag_z2_1_latent)                          0.04      0.02     0.00     0.08 1.05
#> sd(sigma_z11_Intercept)                          0.15      0.01     0.13     0.18 1.00
#> sd(z21_Intercept)                                0.19      0.01     0.17     0.21 1.00
#> sd(z21_lag_z1_1_latent)                          0.01      0.01     0.00     0.04 1.01
#> sd(z21_lag_z2_1_latent)                          0.02      0.02     0.00     0.06 1.00
#> sd(sigma_z21_Intercept)                          0.14      0.01     0.12     0.16 1.00
#> cor(z11_Intercept,z11_lag_z1_1_latent)           0.22      0.26    -0.36     0.67 1.00
#> cor(z11_Intercept,z11_lag_z2_1_latent)           0.35      0.26    -0.23     0.76 1.00
#> cor(z11_lag_z1_1_latent,z11_lag_z2_1_latent)     0.11      0.33    -0.55     0.70 1.01
#> cor(z11_Intercept,sigma_z11_Intercept)          -0.16      0.08    -0.31    -0.00 1.00
#> cor(z11_lag_z1_1_latent,sigma_z11_Intercept)    -0.03      0.28    -0.53     0.54 1.03
#> cor(z11_lag_z2_1_latent,sigma_z11_Intercept)     0.11      0.28    -0.46     0.62 1.09
#> cor(z11_Intercept,z21_Intercept)                -0.07      0.07    -0.21     0.07 1.01
#> cor(z11_lag_z1_1_latent,z21_Intercept)          -0.08      0.21    -0.66     0.30 1.02
#> cor(z11_lag_z2_1_latent,z21_Intercept)          -0.04      0.25    -0.44     0.57 1.10
#> cor(sigma_z11_Intercept,z21_Intercept)           0.07      0.08    -0.08     0.23 1.01
#> cor(z11_Intercept,z21_lag_z1_1_latent)          -0.12      0.29    -0.64     0.46 1.01
#> cor(z11_lag_z1_1_latent,z21_lag_z1_1_latent)    -0.06      0.32    -0.66     0.60 1.00
#> cor(z11_lag_z2_1_latent,z21_lag_z1_1_latent)    -0.03      0.33    -0.66     0.60 1.00
#> cor(sigma_z11_Intercept,z21_lag_z1_1_latent)    -0.06      0.31    -0.62     0.56 1.00
#> cor(z21_Intercept,z21_lag_z1_1_latent)          -0.03      0.28    -0.55     0.55 1.00
#> cor(z11_Intercept,z21_lag_z2_1_latent)           0.16      0.30    -0.48     0.68 1.00
#> cor(z11_lag_z1_1_latent,z21_lag_z2_1_latent)     0.08      0.33    -0.57     0.65 1.01
#> cor(z11_lag_z2_1_latent,z21_lag_z2_1_latent)     0.09      0.35    -0.58     0.75 1.01
#> cor(sigma_z11_Intercept,z21_lag_z2_1_latent)    -0.04      0.28    -0.57     0.52 1.00
#> cor(z21_Intercept,z21_lag_z2_1_latent)          -0.03      0.30    -0.61     0.53 1.00
#> cor(z21_lag_z1_1_latent,z21_lag_z2_1_latent)    -0.04      0.33    -0.64     0.62 1.00
#> cor(z11_Intercept,sigma_z21_Intercept)           0.17      0.09     0.01     0.34 1.00
#> cor(z11_lag_z1_1_latent,sigma_z21_Intercept)     0.16      0.29    -0.43     0.66 1.08
#> cor(z11_lag_z2_1_latent,sigma_z21_Intercept)    -0.02      0.28    -0.54     0.58 1.05
#> cor(sigma_z11_Intercept,sigma_z21_Intercept)     0.19      0.10    -0.02     0.38 1.00
#> cor(z21_Intercept,sigma_z21_Intercept)          -0.09      0.09    -0.26     0.08 1.01
#> cor(z21_lag_z1_1_latent,sigma_z21_Intercept)    -0.17      0.28    -0.66     0.41 1.02
#> cor(z21_lag_z2_1_latent,sigma_z21_Intercept)    -0.04      0.29    -0.63     0.52 1.02
#>                                              Bulk_ESS Tail_ESS
#> sd(z11_Intercept)                                  70      245
#> sd(z11_lag_z1_1_latent)                            71      348
#> sd(z11_lag_z2_1_latent)                            43      360
#> sd(sigma_z11_Intercept)                           449      683
#> sd(z21_Intercept)                                 288      394
#> sd(z21_lag_z1_1_latent)                           244      557
#> sd(z21_lag_z2_1_latent)                           251      445
#> sd(sigma_z21_Intercept)                           559      816
#> cor(z11_Intercept,z11_lag_z1_1_latent)           1162      718
#> cor(z11_Intercept,z11_lag_z2_1_latent)            849      449
#> cor(z11_lag_z1_1_latent,z11_lag_z2_1_latent)      360      455
#> cor(z11_Intercept,sigma_z11_Intercept)            564      683
#> cor(z11_lag_z1_1_latent,sigma_z11_Intercept)       31       77
#> cor(z11_lag_z2_1_latent,sigma_z11_Intercept)       18       46
#> cor(z11_Intercept,z21_Intercept)                   90      245
#> cor(z11_lag_z1_1_latent,z21_Intercept)             24       48
#> cor(z11_lag_z2_1_latent,z21_Intercept)             23       37
#> cor(sigma_z11_Intercept,z21_Intercept)            126      234
#> cor(z11_Intercept,z21_lag_z1_1_latent)           1845      468
#> cor(z11_lag_z1_1_latent,z21_lag_z1_1_latent)      560      590
#> cor(z11_lag_z2_1_latent,z21_lag_z1_1_latent)      898      793
#> cor(sigma_z11_Intercept,z21_lag_z1_1_latent)     1380      826
#> cor(z21_Intercept,z21_lag_z1_1_latent)           1559      744
#> cor(z11_Intercept,z21_lag_z2_1_latent)           1243      719
#> cor(z11_lag_z1_1_latent,z21_lag_z2_1_latent)      759      705
#> cor(z11_lag_z2_1_latent,z21_lag_z2_1_latent)      303      538
#> cor(sigma_z11_Intercept,z21_lag_z2_1_latent)     1100      912
#> cor(z21_Intercept,z21_lag_z2_1_latent)           1486      685
#> cor(z21_lag_z1_1_latent,z21_lag_z2_1_latent)      890      768
#> cor(z11_Intercept,sigma_z21_Intercept)            780      652
#> cor(z11_lag_z1_1_latent,sigma_z21_Intercept)       26       61
#> cor(z11_lag_z2_1_latent,sigma_z21_Intercept)       35       38
#> cor(sigma_z11_Intercept,sigma_z21_Intercept)      466      657
#> cor(z21_Intercept,sigma_z21_Intercept)           1011      841
#> cor(z21_lag_z1_1_latent,sigma_z21_Intercept)       58       87
#> cor(z21_lag_z2_1_latent,sigma_z21_Intercept)       69      310
#> 
#> Regression Coefficients:
#>                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> z11_Intercept           0.23      0.01     0.20     0.25 1.01       43      112
#> sigma_z11_Intercept    -1.61      0.01    -1.64    -1.59 1.00      421      643
#> z21_Intercept          -1.52      0.01    -1.54    -1.49 1.01      213      481
#> sigma_z21_Intercept    -1.90      0.01    -1.93    -1.88 1.01      856      796
#> z11_lag_z1_1_latent     0.47      0.01     0.46     0.49 1.00     1474      584
#> z11_lag_z2_1_latent     0.10      0.01     0.08     0.13 1.00      970      805
#> z21_lag_z1_1_latent     0.10      0.01     0.08     0.11 1.00     1941      657
#> z21_lag_z2_1_latent     0.27      0.01     0.25     0.29 1.00     2349      787
#> 
#> Residual Correlations: 
#>                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> rescor(z11,z21)    -0.00      0.01    -0.02     0.02 1.00     1624      702
#> 
#> Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

As in our first outcome example, the last step of a recovery study is aligning the fitted estimates with the simulation truth. That alignment does more work here. The simulator labels its parameters with structured names, for example ar|term=ar1()|to=ilr1|from=ilr2. However, brms reports estimates under its own coefficient names, here z11_lag_z2_1_within, and matching the two by hand is error-prone. To avoid this, prep_sim_analysis() returns the truth table pre-aligned to the analysis model’s parameter names in $truth. sim_recovery() then joins it with the fitted model’s posterior summaries (fixef(), VarCorr(), and the residual correlations) by name, and adds bias and interval-coverage columns. Random effects are reported on the standard deviation and correlation scale used by VarCorr().

dynamic_recovery <- sim_recovery(dynamic_fit, dynamic_analysis)
kable(dynamic_recovery[, !"simulator_name"], digits = 3)
type group parameter truth estimate est_error lower upper bias covered
fixed NA z11_Intercept 0.20 0.225 0.013 0.201 0.248 0.025 FALSE
fixed NA z21_Intercept -1.50 -1.515 0.013 -1.540 -1.489 -0.015 TRUE
fixed NA sigma_z11_Intercept -1.61 -1.614 0.013 -1.639 -1.588 -0.005 TRUE
fixed NA sigma_z21_Intercept -1.90 -1.905 0.012 -1.927 -1.882 -0.008 TRUE
fixed NA z11_lag_z1_1_latent 0.50 0.475 0.009 0.456 0.492 -0.025 FALSE
fixed NA z11_lag_z2_1_latent 0.10 0.104 0.013 0.079 0.131 0.004 TRUE
fixed NA z21_lag_z1_1_latent 0.10 0.097 0.007 0.083 0.111 -0.003 TRUE
fixed NA z21_lag_z2_1_latent 0.30 0.273 0.010 0.252 0.293 -0.027 FALSE
random_sd ID z11_Intercept 0.20 0.196 0.010 0.179 0.214 -0.004 TRUE
random_sd ID z21_Intercept 0.18 0.185 0.010 0.167 0.207 0.005 TRUE
random_sd ID z11_lag_z1_1_latent 0.02 0.025 0.016 0.002 0.059 0.005 TRUE
random_sd ID z11_lag_z2_1_latent 0.01 0.036 0.021 0.002 0.079 0.026 TRUE
random_sd ID z21_lag_z1_1_latent 0.01 0.014 0.011 0.001 0.040 0.004 TRUE
random_sd ID z21_lag_z2_1_latent 0.02 0.021 0.015 0.001 0.058 0.001 TRUE
random_sd ID sigma_z11_Intercept 0.16 0.152 0.012 0.131 0.176 -0.008 TRUE
random_sd ID sigma_z21_Intercept 0.14 0.140 0.011 0.119 0.161 0.000 TRUE
random_cor ID cor(z11_Intercept,z21_Intercept) 0.00 -0.072 0.071 -0.214 0.073 -0.072 TRUE
random_cor ID cor(z11_Intercept,z11_lag_z1_1_latent) 0.20 0.218 0.258 -0.362 0.672 0.018 TRUE
random_cor ID cor(z11_Intercept,z11_lag_z2_1_latent) 0.00 0.345 0.257 -0.228 0.757 0.345 TRUE
random_cor ID cor(z11_Intercept,z21_lag_z1_1_latent) 0.00 -0.125 0.289 -0.636 0.458 -0.125 TRUE
random_cor ID cor(z11_Intercept,z21_lag_z2_1_latent) 0.00 0.162 0.300 -0.482 0.682 0.162 TRUE
random_cor ID cor(z11_Intercept,sigma_z11_Intercept) 0.00 -0.156 0.080 -0.306 -0.001 -0.156 FALSE
random_cor ID cor(z11_Intercept,sigma_z21_Intercept) 0.00 0.171 0.085 0.011 0.344 0.171 FALSE
random_cor ID cor(z21_Intercept,z11_lag_z1_1_latent) 0.00 -0.084 0.213 -0.658 0.300 -0.084 TRUE
random_cor ID cor(z21_Intercept,z11_lag_z2_1_latent) 0.00 -0.042 0.251 -0.439 0.574 -0.042 TRUE
random_cor ID cor(z21_Intercept,z21_lag_z1_1_latent) 0.20 -0.026 0.279 -0.554 0.549 -0.226 TRUE
random_cor ID cor(z21_Intercept,z21_lag_z2_1_latent) 0.00 -0.029 0.296 -0.610 0.527 -0.029 TRUE
random_cor ID cor(z21_Intercept,sigma_z11_Intercept) 0.00 0.067 0.081 -0.081 0.226 0.067 TRUE
random_cor ID cor(z21_Intercept,sigma_z21_Intercept) 0.00 -0.093 0.089 -0.258 0.085 -0.093 TRUE
random_cor ID cor(z11_lag_z1_1_latent,z11_lag_z2_1_latent) 0.15 0.113 0.326 -0.546 0.702 -0.037 TRUE
random_cor ID cor(z11_lag_z1_1_latent,z21_lag_z1_1_latent) 0.00 -0.055 0.324 -0.662 0.597 -0.055 TRUE
random_cor ID cor(z11_lag_z1_1_latent,z21_lag_z2_1_latent) 0.00 0.081 0.326 -0.574 0.646 0.081 TRUE
random_cor ID cor(z11_lag_z1_1_latent,sigma_z11_Intercept) 0.00 -0.035 0.278 -0.526 0.538 -0.035 TRUE
random_cor ID cor(z11_lag_z1_1_latent,sigma_z21_Intercept) 0.00 0.158 0.294 -0.431 0.663 0.158 TRUE
random_cor ID cor(z11_lag_z2_1_latent,z21_lag_z1_1_latent) 0.00 -0.028 0.331 -0.658 0.597 -0.028 TRUE
random_cor ID cor(z11_lag_z2_1_latent,z21_lag_z2_1_latent) 0.00 0.085 0.351 -0.583 0.749 0.085 TRUE
random_cor ID cor(z11_lag_z2_1_latent,sigma_z11_Intercept) 0.00 0.114 0.283 -0.462 0.622 0.114 TRUE
random_cor ID cor(z11_lag_z2_1_latent,sigma_z21_Intercept) 0.00 -0.016 0.279 -0.542 0.575 -0.016 TRUE
random_cor ID cor(z21_lag_z1_1_latent,z21_lag_z2_1_latent) 0.15 -0.040 0.328 -0.636 0.617 -0.190 TRUE
random_cor ID cor(z21_lag_z1_1_latent,sigma_z11_Intercept) 0.00 -0.064 0.305 -0.621 0.558 -0.064 TRUE
random_cor ID cor(z21_lag_z1_1_latent,sigma_z21_Intercept) 0.00 -0.167 0.281 -0.663 0.409 -0.167 TRUE
random_cor ID cor(z21_lag_z2_1_latent,sigma_z11_Intercept) 0.00 -0.043 0.283 -0.570 0.522 -0.043 TRUE
random_cor ID cor(z21_lag_z2_1_latent,sigma_z21_Intercept) 0.00 -0.035 0.290 -0.626 0.517 -0.035 TRUE
random_cor ID cor(sigma_z11_Intercept,sigma_z21_Intercept) 0.00 0.186 0.102 -0.016 0.379 0.186 TRUE
rescor NA rescor(z11,z21) 0.00 -0.001 0.010 -0.019 0.018 -0.001 TRUE

Most parameters are recovered well, but two patterns are worth noting. First, the AR “inertia” coefficients are estimated somewhat below their true values, with credible intervals that can exclude the truth. Second, the location random-intercept standard deviations are estimated above their true values. Under person-mean centering (the default), the fitted intercept is each person’s observed mean. The variance of those observed means is the trait variance plus the variance of each person’s sample mean of the autoregressive residual, which inflates the estimated standard deviation.

Both patterns are expected and are not bugs in the simulator or the model. prep_sim_analysis() deliberately builds the pragmatic analysis model that researchers typically fit to observed data: ar1() becomes person-mean-centered lagged observed responses, whereas gen_outcome() simulates a latent residual VAR(1). Lagged regression on observed scores with person-mean centering is known to attenuate inertia and cross-lag estimates relative to the latent generating values, especially with short series (Nickell bias; see Hamaker & Grasman, 2014, and the “Pragmatic default estimator” section of ?prep_sim_analysis). Note that in prep_sim_analysis(), we used centering = "latent" to build the lagged predictors, this means that the person mean centering is done on the true latent person means, rather than observed person means. However, these are still not the same as the latent residuals.

This mismatch is intentional as latent residual VAR(1) models, at least that can be easily fit, have various limitations (e.g., in our simulation we support the autoregressive effects being moderated by other variables). Because of this, once ar1() is in the model, to our knowledge no parameter is guaranteed to keep its generating estimand, including the location effects: the observed lag carries the lagged mean structure as well as the lagged residual, and person-mean centering does not generally remove it. Static models are the more forgiving case. There most parameters do estimate the quantity the simulator recorded, and those are the ones where systematic bias would indicate a problem.

The next section fits the non-centered alternative to the same simulated data, which shows how much of the AR attenuation comes from the centering itself.

To Center or Not to Center: Non-Centered Lagged Predictors

The person-mean-centered lagged predictor used above is only one of two common ways to parametrize a multilevel autoregressive model. Hamaker & Grasman (2014) compare them directly. With cluster-mean centering (CMC, the prep_sim_analysis() default also called person-mean centred), the average autoregressive coefficient is attenuated toward zero. Regressing on the raw, non-centred lagged outcome (NC) instead recovers the average autoregressive coefficient nearly unbiased.

The trade-off is the interpretability of the level parameters. In the NC parametrization, the model intercept corresponds to (1−ϕi)μi(1 - \phi_i)\,\mu_i (multivariate: (𝐈−𝚽i)𝛍i(\mathbf{I} - \boldsymbol{\Phi}_i)\,\boldsymbol{\mu}_i) rather than the person mean μi\mu_i. As a result, the location intercepts, and their random-effect standard deviations and correlations, no longer estimate the simulated data’s mean-structure truth. Hamaker & Grasman’s practical recommendation is therefore to fit both parametrizations: NC when the average inertia and cross-lag effects are the parameters of interest, and CMC when person means must remain meaningful, for example for interpretable intercepts or for cross-level effects of person-level predictors on ϕi\phi_i.

prep_sim_analysis() supports the NC parametrization through lag_center = "none", which builds raw lag_<response> columns instead of the person-mean-centered lag_<response>_within columns. The caveat above does not change: with ar1() in the model, no parameter is guaranteed to keep its generating estimand under either parametrization. What changes is which parameters we expect to be recovered nearly without bias. Under NC, the AR rows. Under CMC, the intercepts and means. We reuse the simulated data from above, so the two analysis models see exactly the same observations.

dynamic_analysis_nc <- prep_sim_analysis(dynamic_comp, lag_center = "none")
dynamic_analysis_nc$metadata$lag_columns
#> [1] "lag_z1_1" "lag_z2_1"
dynamic_analysis_nc$formula
#> z1_1 ~ lag_z1_1 + lag_z2_1 + (1 + lag_z1_1 + lag_z2_1 | p1 | ID) 
#> sigma ~ 1 + (1 | p1 | ID)
#> z2_1 ~ lag_z1_1 + lag_z2_1 + (1 + lag_z1_1 + lag_z2_1 | p1 | ID) 
#> sigma ~ 1 + (1 | p1 | ID)
dynamic_fit_nc <- brmcoda(
  complr = dynamic_analysis_nc$complr,
  formula = dynamic_analysis_nc$formula,
  backend = "cmdstanr",
  seed = 2026,
  chains = 2,
  cores = 2,
  threads = brms::threading(4, static = TRUE),
  iter = 1000,
  refresh = 0
)

The full recovery table for the non-centered model parallels the one shown above for the centered default. Note the lag_<response> parameter names, rather than lag_<response>_within:

dynamic_recovery_nc <- sim_recovery(dynamic_fit_nc, dynamic_analysis_nc)
kable(dynamic_recovery_nc[, !"simulator_name"], digits = 3)
type group parameter truth estimate est_error lower upper bias covered
fixed NA z11_Intercept 0.20 0.270 0.020 0.234 0.311 0.070 FALSE
fixed NA z21_Intercept -1.50 -1.082 0.018 -1.119 -1.048 0.418 FALSE
fixed NA sigma_z11_Intercept -1.61 -1.613 0.013 -1.639 -1.588 -0.003 TRUE
fixed NA sigma_z21_Intercept -1.90 -1.905 0.012 -1.929 -1.882 -0.008 TRUE
fixed NA z11_lag_z1_1 0.50 0.508 0.010 0.489 0.527 0.008 TRUE
fixed NA z11_lag_z2_1 0.10 0.108 0.012 0.085 0.132 0.008 TRUE
fixed NA z21_lag_z1_1 0.10 0.098 0.007 0.085 0.112 -0.002 TRUE
fixed NA z21_lag_z2_1 0.30 0.300 0.011 0.279 0.321 0.000 TRUE
random_sd ID z11_Intercept 0.20 0.105 0.012 0.086 0.135 -0.095 FALSE
random_sd ID z21_Intercept 0.18 0.132 0.011 0.110 0.155 -0.048 FALSE
random_sd ID z11_lag_z1_1 0.02 0.018 0.013 0.000 0.049 -0.002 TRUE
random_sd ID z11_lag_z2_1 0.01 0.013 0.010 0.001 0.036 0.003 TRUE
random_sd ID z21_lag_z1_1 0.01 0.012 0.010 0.000 0.036 0.002 TRUE
random_sd ID z21_lag_z2_1 0.02 0.016 0.011 0.001 0.039 -0.004 TRUE
random_sd ID sigma_z11_Intercept 0.16 0.151 0.011 0.129 0.174 -0.009 TRUE
random_sd ID sigma_z21_Intercept 0.14 0.139 0.011 0.119 0.160 -0.001 TRUE
random_cor ID cor(z11_Intercept,z21_Intercept) 0.00 -0.305 0.107 -0.490 -0.087 -0.305 FALSE
random_cor ID cor(z11_Intercept,z11_lag_z1_1) 0.20 0.137 0.302 -0.495 0.698 -0.063 TRUE
random_cor ID cor(z11_Intercept,z11_lag_z2_1) 0.00 0.258 0.357 -0.489 0.808 0.258 TRUE
random_cor ID cor(z11_Intercept,z21_lag_z1_1) 0.00 -0.099 0.284 -0.628 0.478 -0.099 TRUE
random_cor ID cor(z11_Intercept,z21_lag_z2_1) 0.00 0.192 0.312 -0.463 0.717 0.192 TRUE
random_cor ID cor(z11_Intercept,sigma_z11_Intercept) 0.00 -0.157 0.106 -0.350 0.062 -0.157 TRUE
random_cor ID cor(z11_Intercept,sigma_z21_Intercept) 0.00 0.175 0.112 -0.071 0.383 0.175 TRUE
random_cor ID cor(z21_Intercept,z11_lag_z1_1) 0.00 -0.129 0.268 -0.687 0.389 -0.129 TRUE
random_cor ID cor(z21_Intercept,z11_lag_z2_1) 0.00 0.167 0.286 -0.450 0.659 0.167 TRUE
random_cor ID cor(z21_Intercept,z21_lag_z1_1) 0.20 0.013 0.306 -0.592 0.599 -0.187 TRUE
random_cor ID cor(z21_Intercept,z21_lag_z2_1) 0.00 0.050 0.303 -0.576 0.606 0.050 TRUE
random_cor ID cor(z21_Intercept,sigma_z11_Intercept) 0.00 0.080 0.101 -0.140 0.269 0.080 TRUE
random_cor ID cor(z21_Intercept,sigma_z21_Intercept) 0.00 -0.121 0.103 -0.326 0.075 -0.121 TRUE
random_cor ID cor(z11_lag_z1_1,z11_lag_z2_1) 0.15 0.006 0.337 -0.610 0.639 -0.144 TRUE
random_cor ID cor(z11_lag_z1_1,z21_lag_z1_1) 0.00 -0.014 0.329 -0.602 0.623 -0.014 TRUE
random_cor ID cor(z11_lag_z1_1,z21_lag_z2_1) 0.00 0.076 0.338 -0.582 0.696 0.076 TRUE
random_cor ID cor(z11_lag_z1_1,sigma_z11_Intercept) 0.00 -0.214 0.302 -0.750 0.396 -0.214 TRUE
random_cor ID cor(z11_lag_z1_1,sigma_z21_Intercept) 0.00 0.048 0.266 -0.534 0.533 0.048 TRUE
random_cor ID cor(z11_lag_z2_1,z21_lag_z1_1) 0.00 -0.010 0.331 -0.620 0.631 -0.010 TRUE
random_cor ID cor(z11_lag_z2_1,z21_lag_z2_1) 0.00 0.028 0.325 -0.567 0.650 0.028 TRUE
random_cor ID cor(z11_lag_z2_1,sigma_z11_Intercept) 0.00 -0.052 0.322 -0.656 0.558 -0.052 TRUE
random_cor ID cor(z11_lag_z2_1,sigma_z21_Intercept) 0.00 0.046 0.329 -0.572 0.672 0.046 TRUE
random_cor ID cor(z21_lag_z1_1,z21_lag_z2_1) 0.15 -0.006 0.319 -0.631 0.588 -0.156 TRUE
random_cor ID cor(z21_lag_z1_1,sigma_z11_Intercept) 0.00 -0.059 0.314 -0.651 0.572 -0.059 TRUE
random_cor ID cor(z21_lag_z1_1,sigma_z21_Intercept) 0.00 -0.091 0.294 -0.633 0.506 -0.091 TRUE
random_cor ID cor(z21_lag_z2_1,sigma_z11_Intercept) 0.00 -0.104 0.306 -0.665 0.533 -0.104 TRUE
random_cor ID cor(z21_lag_z2_1,sigma_z21_Intercept) 0.00 -0.034 0.286 -0.578 0.522 -0.034 TRUE
random_cor ID cor(sigma_z11_Intercept,sigma_z21_Intercept) 0.00 0.179 0.103 -0.040 0.378 0.179 TRUE
rescor NA rescor(z11,z21) 0.00 0.000 0.011 -0.021 0.021 0.000 TRUE

The two recovery tables label the same simulator parameters with different analysis-model names, lag_z1_1_within versus lag_z1_1. However, the simulator_name column is identical in both, so we can use it as a stable key for a side-by-side comparison of the autoregressive fixed effects:

ar_cmc <- dynamic_recovery[type == "fixed" & startsWith(simulator_name, "ar:")]
ar_nc <- dynamic_recovery_nc[type == "fixed" & startsWith(simulator_name, "ar:")]
ar_comparison <- merge(
  ar_cmc[, .(simulator_name, truth, estimate_cmc = estimate,
             ci_cmc = sprintf("[%.3f, %.3f]", lower, upper), bias_cmc = bias)],
  ar_nc[, .(simulator_name, estimate_nc = estimate,
            ci_nc = sprintf("[%.3f, %.3f]", lower, upper), bias_nc = bias)],
  by = "simulator_name"
)
kable(ar_comparison, digits = 3)
simulator_name truth estimate_cmc ci_cmc bias_cmc estimate_nc ci_nc bias_nc
ar:ar1()@to=ilr1,from=ilr1 0.5 0.475 [0.456, 0.492] -0.025 0.508 [0.489, 0.527] 0.008
ar:ar1()@to=ilr1,from=ilr2 0.1 0.104 [0.079, 0.131] 0.004 0.108 [0.085, 0.132] 0.008
ar:ar1()@to=ilr2,from=ilr1 0.1 0.097 [0.083, 0.111] -0.003 0.098 [0.085, 0.112] -0.002
ar:ar1()@to=ilr2,from=ilr2 0.3 0.273 [0.252, 0.293] -0.027 0.300 [0.279, 0.321] 0.000

The inertia (diagonal) coefficients show the pattern Hamaker & Grasman (2014) predict. The CMC estimates sit below their true values, while the NC estimates land close to the truth, with credible intervals that cover it.

The level parameters tell a different story, and neither parametrization recovers them exactly. Longer series shrink the CMC inflation, because the residual mean averages away as TT grows, but they do not change the NC shrinkage. So, generally the CMC will be closer than the NC for the intercepts.

level_cmc <- dynamic_recovery[
  type %in% c("fixed", "random_sd") & grepl("^location", simulator_name)]
level_nc <- dynamic_recovery_nc[
  type %in% c("fixed", "random_sd") & grepl("^location", simulator_name)]
level_comparison <- merge(
  level_cmc[, .(simulator_name, type, truth, estimate_cmc = estimate,
                ci_cmc = sprintf("[%.3f, %.3f]", lower, upper), bias_cmc = bias)],
  level_nc[, .(simulator_name, estimate_nc = estimate,
               ci_nc = sprintf("[%.3f, %.3f]", lower, upper), bias_nc = bias)],
  by = "simulator_name"
)
kable(level_comparison, digits = 3)
simulator_name type truth estimate_cmc ci_cmc bias_cmc estimate_nc ci_nc bias_nc
location:(Intercept)@ilr1 fixed 0.20 0.225 [0.201, 0.248] 0.025 0.270 [0.234, 0.311] 0.070
location:(Intercept)@ilr2 fixed -1.50 -1.515 [-1.540, -1.489] -0.015 -1.082 [-1.119, -1.048] 0.418
location|outcome=ilr1|term=(Intercept) random_sd 0.20 0.196 [0.179, 0.214] -0.004 0.105 [0.086, 0.135] -0.095
location|outcome=ilr2|term=(Intercept) random_sd 0.18 0.185 [0.167, 0.207] 0.005 0.132 [0.110, 0.155] -0.048

Which one to fit, or whether to fit both as Hamaker & Grasman (2014) recommend, depends on which parameters matter to you. See the “Pragmatic default estimator” section of ?prep_sim_analysis for more on which parameters keep their generating estimand under each parametrization.

Simulating the Between-Within (brmcoda) Model

The between-within multilevelcoda model (see the earlier vignettes on multilevel models with compositional predictors and on substitution analyses) predicts a scalar outcome from the between- and within-person ILR coordinates of a composition. We can simulate this model directly. A multilevel compositional gen_mvn() generator emits the latent decomposition of each ILR coordinate as visible columns (ilr1_between, ilr1_within, and so on) and the true between composition on the parts scale (sleep_between, …), and it labels the ILR coordinates with column roles. between() and within() terms in a gen_outcome() formula then resolve to those latent components.

As before, gen_template() provides the exact parameter structure required. Note the role-based names, such as between(ilr1), in the location matrix.

bw_formula <- y ~ between(ilr1) + between(ilr2) +
  within(ilr1) + within(ilr2) + (1 | ID)

bw_generator <- function() {
  gen_mvn(
    c("ilr1", "ilr2"),
    level = "multilevel",
    fixed_intercept = c(0.2, -0.1),
    random_cov = diag(2) * 0.2,
    residual_cov = diag(2) * 0.4,
    compositional = TRUE,
    parts = c("sleep", "active", "sedentary"),
    total = 24
  )
}

bw_template <- simulate_data(
  n_groups = 150,
  n_per_group = 10,
  group_id = "ID",
  seed = 2026,
  generators = list(
    sleep_comp = bw_generator(),
    outcome_template = gen_template(bw_formula, scale = sigma ~ 1)
  )
)

bw_params <- bw_template$generator_metadata$outcome_template$params
bw_params$location$beta[, "y"] <- c(0, 0.5, -0.3, 0.2, 0.1)
bw_params$scale$beta["(Intercept)", "y"] <- log(0.3)
bw_params$random$ID$covariance[] <- 0.1
bw_params$location$beta
#>                  y
#> (Intercept)    0.0
#> between(ilr1)  0.5
#> between(ilr2) -0.3
#> within(ilr1)   0.2
#> within(ilr2)   0.1
bw_sim <- simulate_data(
  n_groups = 150,
  n_per_group = 10,
  group_id = "ID",
  seed = 2026,
  generators = list(
    sleep_comp = bw_generator(),
    y = gen_outcome(bw_formula, scale = sigma ~ 1, params = bw_params)
  )
)

kable(head(bw_sim$data[, c(
  "ID", "sleep", "active", "sedentary",
  "ilr1_between", "ilr1_within", "sleep_between", "y"
), with = FALSE]))
ID sleep active sedentary ilr1_between ilr1_within sleep_between y
1 16.30 3.96 3.75 0.189 0.989 9.25 -0.071
1 11.90 6.12 5.98 0.189 0.363 9.25 0.297
1 9.46 5.22 9.32 0.189 0.060 9.25 0.904
1 11.86 8.98 3.16 0.189 0.465 9.25 -0.123
1 4.20 7.19 12.62 0.189 -0.858 9.25 0.481
1 8.44 5.90 9.65 0.189 -0.098 9.25 0.435

prep_sim_analysis() recognises that between() and within() reference ILR coordinates of a compositional predictor. Instead of generic person-mean centered columns, it builds a complr() object for the predictor composition, using the simulator’s parts, SBP, and total, and maps the terms to the standard multilevelcoda between and within coordinates bz* and wz*. The returned complr object plugs straight into brmcoda(). Because the fitted model is a regular brmcoda model, substitution() works on it too.

bw_analysis <- prep_sim_analysis(bw_sim)
bw_analysis$metadata$special_term_map
#> between(ilr1) between(ilr2)  within(ilr1)  within(ilr2) 
#>       "bz1_1"       "bz2_1"       "wz1_1"       "wz2_1"
bw_analysis$formula
#> y ~ bz1_1 + bz2_1 + wz1_1 + wz2_1 + (1 | ID) 
#> sigma ~ 1
bw_fit <- brmcoda(
  complr = bw_analysis$complr,
  formula = bw_analysis$formula,
  backend = "cmdstanr",
  seed = 2026,
  chains = 2,
  cores = 2,
  iter = 1000,
  refresh = 0
)

sim_recovery() works unchanged here. The truth table maps the role-based simulator names (between(ilr1), …) to the fitted bz* and wz* coefficients, and it also covers the random-intercept standard deviation automatically.

bw_recovery <- sim_recovery(bw_fit, bw_analysis)
kable(bw_recovery[, !"simulator_name"], digits = 3)
type group parameter truth estimate est_error lower upper bias covered
fixed NA Intercept 0.000 0.006 0.032 -0.055 0.072 0.006 TRUE
fixed NA bz1_1 0.500 0.409 0.066 0.271 0.539 -0.091 TRUE
fixed NA bz2_1 -0.300 -0.292 0.057 -0.402 -0.187 0.008 TRUE
fixed NA wz1_1 0.200 0.209 0.014 0.182 0.236 0.009 TRUE
fixed NA wz2_1 0.100 0.092 0.013 0.066 0.115 -0.008 TRUE
fixed NA sigma_Intercept -1.204 -1.191 0.019 -1.227 -1.156 0.013 TRUE
random_sd ID Intercept 0.316 0.321 0.022 0.281 0.367 0.005 TRUE

The same caveat about what is being estimated applies here as for the dynamic model. The fitted bz* coefficients use the ILR of each person’s closed arithmetic-mean composition, computed from the observed parts. The simulation truth between(ilr) is instead the latent group-level ILR mean. These differ, and observed between-person predictors are biased for latent between-person effects when the series are short (Ludtke et al., 2008). Within-person effects are typically recovered well. See the “Pragmatic default estimator” section of ?prep_sim_analysis for details.

However, if desired we can use the simulated latent variables directly.

bw_analysis_latent <- prep_sim_analysis(bw_sim, centering = "latent")
bw_fit_latent <- brmcoda(
  complr = bw_analysis_latent$complr,
  formula = bw_analysis_latent$formula,
  backend = "cmdstanr",
  seed = 2026,
  chains = 2,
  cores = 2,
  iter = 1000,
  refresh = 0
)

And the recovery

bw_recovery_latent <- sim_recovery(bw_fit_latent, bw_analysis_latent)
kable(bw_recovery_latent[, !"simulator_name"], digits = 3)
type group parameter truth estimate est_error lower upper bias covered
fixed NA Intercept 0.000 -0.005 0.028 -0.061 0.053 -0.005 TRUE
fixed NA bz1_1 0.500 0.453 0.072 0.310 0.594 -0.047 TRUE
fixed NA bz2_1 -0.300 -0.338 0.064 -0.454 -0.211 -0.038 TRUE
fixed NA wz1_1 0.200 0.209 0.014 0.183 0.237 0.009 TRUE
fixed NA wz2_1 0.100 0.092 0.013 0.067 0.116 -0.008 TRUE
fixed NA sigma_Intercept -1.204 -1.190 0.019 -1.229 -1.153 0.014 TRUE
random_sd ID Intercept 0.316 0.309 0.020 0.270 0.348 -0.008 TRUE

The bz* rows are the ones to compare against the previous table. Here they estimate the latent between-person composition the simulator generated from, rather than an observed proxy for it, so any remaining bias in them is not attributable to the arithmetic-versus-geometric centering gap. The wz* rows were already recovered well and change little.

Finally, because the fit is a regular brmcoda model with a complr object, substitution analyses work directly on the simulated data:

bw_sub <- substitution(bw_fit, delta = 1, level = "within")
summary(bw_sub)
#>    Estimate Est.Error CI_low CI_high Delta      From        To  Level Reference   Resp
#>       <num>     <num>  <num>   <num> <num>    <char>    <char> <char>    <char> <char>
#> 1:     0.02         0   0.02    0.03     1    active     sleep within grandmean      y
#> 2:     0.04         0   0.03    0.04     1 sedentary     sleep within grandmean      y
#> 3:    -0.02         0  -0.03   -0.02     1     sleep    active within grandmean      y
#> 4:     0.02         0   0.01    0.02     1 sedentary    active within grandmean      y
#> 5:    -0.04         0  -0.04   -0.03     1     sleep sedentary within grandmean      y
#> 6:    -0.01         0  -0.02   -0.01     1    active sedentary within grandmean      y

Non-Gaussian Outcomes

Beyond the Gaussian default, gen_outcome() supports five univariate GLM-style families: "poisson", "binomial", "negbin", "gamma", and "beta". The location formula works exactly as before, with ordinary terms, between(), within(), and one grouping term. It defines the linear predictor on the family link scale. That is a log link for Poisson, negative binomial, and gamma, and a logit link for binomial and beta.

The scale formula uses the family’s brms distributional-parameter name on the left-hand side, so the simulation formula matches the analysis formula:

  • "negbin": scale = shape ~ ... models the log size.
  • "gamma": scale = shape ~ ... models the log shape.
  • "beta": scale = phi ~ ... models the log precision.
  • "poisson" and "binomial" have no auxiliary parameter, so scale must be omitted.

Binomial outcomes require trials (a scalar, vector, function of n, or count-distribution list). The resolved trial sizes are written into the simulated data as a <outcome>_trials column so that the analysis model can use y | trials(y_trials).

There are two restrictions to keep in mind. First, ar1() is currently only supported for family = "gaussian", and combining it with any other family is an error. Second, because the links are nonlinear, the fixed effects are conditional (subject-specific) effects when random effects are present. The marginal mean of the outcome is not the inverse link of the fixed-effect predictor alone. For example, the marginal probability of a binary outcome is not plogis(b0) under random intercepts. Do not validate simulated data against marginal summaries. The recovery checks below fit the matching conditional model instead.

The examples below simulate each family with known parameters, show the brms formula that prep_sim_analysis() infers, and briefly fit each model to check that the true values are recovered. A single simulated data set only identifies parameters up to sampling error, so the tables report 95% credible intervals around each estimate next to the true value, rather than expecting point equality.

glm_intercept_beta <- function(value, terms = "(Intercept)") {
  matrix(value, nrow = length(terms), dimnames = list(terms, "y"))
}

glm_fit <- function(analysis, seed) {
  brm(
    formula = analysis$formula,
    data = analysis$data,
    backend = "cmdstanr",
    chains = 2,
    iter = 1000,
    seed = seed,
    refresh = 0,
    silent = 2
  )
}

glm_show_recovery <- function(fit, analysis) {
  kable(sim_recovery(fit, analysis)[, !"simulator_name"], digits = 3)
}

Poisson

We start with a multilevel Poisson outcome with a row-level predictor and a random intercept. The random intercept enters the log-linear predictor, so counts are conditionally Poisson given the group effect.

poisson_random_name <- "location|outcome=y|term=(Intercept)"
poisson_random_cov <- matrix(
  0.4^2,
  dimnames = list(poisson_random_name, poisson_random_name)
)
poisson_params <- list(
  location = list(beta = glm_intercept_beta(c(1, 0.3), c("(Intercept)", "x"))),
  random = list(ID = list(covariance = poisson_random_cov))
)

poisson_sim <- simulate_data(
  n_groups = 100,
  n_per_group = 10,
  group_id = "ID",
  seed = 2027,
  generators = list(
    x = gen_mvn("x", fixed_intercept = 0, residual_cov = 1),
    y = gen_outcome(
      y ~ x + (1 | ID),
      params = poisson_params,
      family = "poisson"
    )
  )
)

poisson_analysis <- prep_sim_analysis(poisson_sim)
kable(head(poisson_sim$data))
ID obs_id x y
1 1 -0.935 0
1 2 -0.694 1
1 3 -0.548 3
1 4 0.766 3
1 5 1.852 4
1 6 0.754 1 
poisson_analysis$formula
#> y ~ x + (1 | ID)
poisson_fit <- glm_fit(poisson_analysis, seed = 2027)
glm_show_recovery(poisson_fit, poisson_analysis)
type group parameter truth estimate est_error lower upper bias covered
fixed NA Intercept 1.0 0.932 0.045 0.842 1.013 -0.068 TRUE
fixed NA x 0.3 0.311 0.019 0.272 0.347 0.011 TRUE
random_sd ID Intercept 0.4 0.435 0.038 0.366 0.516 0.035 TRUE

The random-intercept standard deviation (true value 0.4) is included in the recovery table automatically.

Binomial

Binomial outcomes need trials. Here every row has ten trials. The resolved sizes are stored in the generated y_trials column, which prep_sim_analysis() wires into the trials() term of the analysis formula.

binomial_params <- list(
  location = list(beta = glm_intercept_beta(c(-0.5, 0.4), c("(Intercept)", "x")))
)

binomial_sim <- simulate_data(
  n = 500,
  seed = 2028,
  generators = list(
    x = gen_mvn("x", fixed_intercept = 0, residual_cov = 1),
    y = gen_outcome(
      y ~ x,
      params = binomial_params,
      family = "binomial",
      trials = 10
    )
  )
)

binomial_analysis <- prep_sim_analysis(binomial_sim)
kable(head(binomial_sim$data))
obs_id x y y_trials
1 0.096 3 10
2 0.292 5 10
3 -2.573 2 10
4 1.898 5 10
5 2.186 6 10
6 0.002 5 10 
binomial_analysis$formula
#> y | trials(y_trials) ~ x
binomial_fit <- glm_fit(binomial_analysis, seed = 2028)
glm_show_recovery(binomial_fit, binomial_analysis)
type group parameter truth estimate est_error lower upper bias covered
fixed NA Intercept -0.5 -0.497 0.03 -0.555 -0.441 0.003 TRUE
fixed NA x 0.4 0.401 0.03 0.341 0.460 0.001 TRUE

Negative Binomial

The negative binomial family uses scale = shape ~ ... as the log size. The same shape name appears as the distributional formula in the inferred brms model, so the simulated auxiliary parameter is directly comparable to the fitted shape_Intercept.

negbin_params <- list(
  location = list(beta = glm_intercept_beta(c(1, 0.3), c("(Intercept)", "x"))),
  scale = list(beta = glm_intercept_beta(log(2)))
)

negbin_sim <- simulate_data(
  n = 500,
  seed = 2038,
  generators = list(
    x = gen_mvn("x", fixed_intercept = 0, residual_cov = 1),
    y = gen_outcome(
      y ~ x,
      scale = shape ~ 1,
      params = negbin_params,
      family = "negbin"
    )
  )
)

negbin_analysis <- prep_sim_analysis(negbin_sim)
kable(head(negbin_sim$data))
obs_id x y
1 0.177 0
2 1.145 0
3 1.197 7
4 1.097 10
5 1.902 7
6 -1.700 1 
negbin_analysis$formula
#> y ~ x 
#> shape ~ 1
negbin_fit <- glm_fit(negbin_analysis, seed = 2038)
glm_show_recovery(negbin_fit, negbin_analysis)
type group parameter truth estimate est_error lower upper bias covered
fixed NA Intercept 1.000 1.020 0.043 0.930 1.102 0.020 TRUE
fixed NA x 0.300 0.304 0.040 0.222 0.380 0.004 TRUE
fixed NA shape_Intercept 0.693 0.625 0.115 0.412 0.856 -0.068 TRUE

Gamma

Gamma outcomes use a log link for the mean and scale = shape ~ ... as the log shape. Draws use rate shape / mu, so the simulated mean is exp(eta).

gamma_params <- list(
  location = list(beta = glm_intercept_beta(c(0.5, 0.2), c("(Intercept)", "x"))),
  scale = list(beta = glm_intercept_beta(log(3)))
)

gamma_sim <- simulate_data(
  n = 500,
  seed = 2030,
  generators = list(
    x = gen_mvn("x", fixed_intercept = 0, residual_cov = 1),
    y = gen_outcome(
      y ~ x,
      scale = shape ~ 1,
      params = gamma_params,
      family = "gamma"
    )
  )
)

gamma_analysis <- prep_sim_analysis(gamma_sim)
kable(head(gamma_sim$data))
obs_id x y
1 1.226 1.183
2 -1.026 1.727
3 -1.339 2.081
4 0.579 1.347
5 1.018 4.489
6 0.194 0.931 
gamma_analysis$formula
#> y ~ x 
#> shape ~ 1
gamma_fit <- glm_fit(gamma_analysis, seed = 2030)
glm_show_recovery(gamma_fit, gamma_analysis)
type group parameter truth estimate est_error lower upper bias covered
fixed NA Intercept 0.5 0.516 0.025 0.466 0.565 0.016 TRUE
fixed NA x 0.2 0.201 0.023 0.159 0.242 0.001 TRUE
fixed NA shape_Intercept 1.1 1.145 0.061 1.027 1.265 0.046 TRUE

Beta

Beta outcomes use a logit link for the mean and scale = phi ~ ... as the log precision. This example puts correlated random intercepts on both the mean and the precision. The simulator draws all of a group’s random effects from one joint covariance. prep_sim_analysis() therefore emits brms ID-linked random effects, (1 | p1 | ID), whenever the same grouping factor appears in both formulas, which lets the fitted model estimate the cross-parameter correlation.

beta_random_names <- c(
  "location|outcome=y|term=(Intercept)",
  "phi|outcome=y|term=(Intercept)"
)
beta_random_sd <- c(0.3, 0.2)
beta_random_cor <- matrix(c(1, 0.3, 0.3, 1), 2, 2)
beta_random_cov <- diag(beta_random_sd) %*% beta_random_cor %*% diag(beta_random_sd)
dimnames(beta_random_cov) <- list(beta_random_names, beta_random_names)

beta_params <- list(
  location = list(beta = glm_intercept_beta(0.2)),
  scale = list(beta = glm_intercept_beta(log(20))),
  random = list(ID = list(covariance = beta_random_cov))
)

beta_sim <- simulate_data(
  n_groups = 150,
  n_per_group = 8,
  group_id = "ID",
  seed = 2031,
  generators = list(
    y = gen_outcome(
      y ~ 1 + (1 | ID),
      scale = phi ~ 1 + (1 | ID),
      params = beta_params,
      family = "beta"
    )
  )
)

beta_analysis <- prep_sim_analysis(beta_sim)
kable(head(beta_sim$data))
ID obs_id y
1 1 0.621
1 2 0.513
1 3 0.532
1 4 0.655
1 5 0.758
1 6 0.455 
beta_analysis$formula
#> y ~ 1 + (1 | p1 | ID) 
#> phi ~ 1 + (1 | p1 | ID)
beta_fit <- glm_fit(beta_analysis, seed = 2031)
glm_show_recovery(beta_fit, beta_analysis)
type group parameter truth estimate est_error lower upper bias covered
fixed NA Intercept 0.2 0.210 0.029 0.154 0.265 0.010 TRUE
fixed NA phi_Intercept 3.0 2.972 0.049 2.876 3.072 -0.023 TRUE
random_sd ID Intercept 0.3 0.336 0.024 0.291 0.385 0.036 TRUE
random_sd ID phi_Intercept 0.2 0.166 0.081 0.022 0.328 -0.034 TRUE
random_cor ID cor(Intercept,phi_Intercept) 0.3 0.484 0.297 -0.169 0.966 0.184 TRUE

The random-effect standard deviations (true values 0.3 for the mean intercept and 0.2 for the phi intercept) and their correlation (true value 0.3) all appear in the recovery table. This works because the ID-linked random effects make the cross-parameter correlation estimable.