Package 'multilevelcoda'

Description

Internal function only.

Usage

.meanvar.complr(x, weight = c("equal", "proportional"), parts = 1, ...)

Arguments

Description

This is a constructor function for a complr object. It checks that the input list has the necessary structure and components to be considered a complr object, and if so, it assigns the class "complr" to it. This allows for method dispatch on complr objects.

Usage

as.complr(x)

Arguments

Value

A complr object.

Description

Extract amounts and compositions in conventional formats as data.frames, matrices, or arrays.

Usage

## S3 method for class 'complr'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

## S3 method for class 'complr'
as.matrix(x, ...)

Arguments

Description

This is a constructor function for a diagnostics object. It checks that the input list has the necessary structure and components to be considered a diagnostics object, and if so, it assigns the class "diagnostics" to it. This allows for method dispatch on diagnostics objects.

Usage

as.diagnostics(x)

Arguments

Value

A diagnostics object.

Description

Compute Bayes factors from marginal likelihoods

Usage

## S3 method for class 'brmcoda'
bayes_factor(x1, x2, ...)

Arguments

See Also

bayes_factor.brmsfit

Description

Fit a brm model with multilevel ILR coordinates

Usage

brmcoda(complr, formula, ...)

Arguments

Value

A brmcoda with two elements

Examples

if(requireNamespace("cmdstanr")){
  x1 <- complr(data = mcompd, sbp = sbp,
                 parts = c("TST", "WAKE", "MVPA", "LPA", "SB"), idvar = "ID")

  # inspect variables before passing to brmcoda
  get_variables(x1)

  ## model with compositional predictor at between and within-person levels
  m1 <- brmcoda(complr = x1,
                formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 +
                                   wz1_1 + wz2_1 + wz3_1 + wz4_1 + (1 | ID),
                chain = 1, iter = 500,
                backend = "cmdstanr")

  ## model with compositional outcome
  m2 <- brmcoda(complr = x1,
                formula = mvbind(z1_1, z2_1, z3_1, z4_1) ~ Stress + Female + (1 | ID),
                chain = 1, iter = 500,
                backend = "cmdstanr")

  ## model with compositional predictor and outcome
  x2 <- complr(data = mcompd,
                parts = list(c("TST", "WAKE"), c("MVPA", "LPA", "SB")),
                total = list(c(480), c(960)),
                idvar = "ID",
                transform = "ilr")

  m3 <- brmcoda(complr = x2,
                formula = mvbind(z1_2, z2_2) ~ z1_1 + Female + (1 | ID),
                chain = 1, iter = 500,
                backend = "cmdstanr")
  }

Description

This function is an alias of substitution to estimates the difference in an outcome when compositional parts are substituted for specific unit(s) at between level using a single reference composition (e.g., compositional mean at sample level). It is recommended that users run substitution model using the substitution function.

Usage

bsub(
  object,
  delta,
  ref = "grandmean",
  level = "between",
  summary = TRUE,
  aorg = TRUE,
  at = NULL,
  parts = 1,
  base,
  type = "one-to-one",
  weight = "equal",
  scale = c("response", "linear"),
  cores = NULL,
  ...
)

Arguments

Value

A list containing the results of multilevel compositional substitution model. The first six lists contain the results of the substitution estimation for a compositional part.

See Also

substitution

Examples

if(requireNamespace("cmdstanr")){
cilr <- complr(data = mcompd, sbp = sbp, 
               parts = c("TST", "WAKE", "MVPA", "LPA", "SB"), idvar = "ID", total = 1440)

# model with compositional predictor at between and between-person levels
m <- brmcoda(complr = cilr, 
             formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 + 
                                wz1_1 + wz2_1 + wz3_1 + wz4_1 + Female + (1 | ID), 
             chain = 1, iter = 500,
             backend = "cmdstanr")
subm <- bsub(object = m, base = psub, delta = 5)
}

Description

This function is an alias of substitution to estimates the difference in an outcome when compositional parts are substituted for specific unit(s) at between level using cluster mean (e.g., compositional mean at individual level) as reference composition. It is recommended that users run substitution model using the substitution function.

Usage

bsubmargin(
  object,
  delta,
  ref = "clustermean",
  level = "between",
  summary = TRUE,
  at = NULL,
  parts = 1,
  base,
  type = "one-to-one",
  weight = "proportional",
  scale = c("response", "linear"),
  cores = NULL,
  ...
)

Arguments

Value

A list containing the results of multilevel compositional substitution model. The first six lists contain the results of the substitution estimation for a compositional part.

See Also

substitution

Examples

if(requireNamespace("cmdstanr")){
cilr <- complr(data = mcompd[ID %in% 1:10, .SD[1:3], by = ID], sbp = sbp, 
               parts = c("TST", "WAKE", "MVPA", "LPA", "SB"), idvar = "ID", total = 1440)

m <- brmcoda(complr = cilr, 
             formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 + 
                                wz1_1 + wz2_1 + wz3_1 + wz4_1 + 
                                Female + (1 | ID), 
             chains = 1, iter = 500,
             backend = "cmdstanr")
             
subm <- bsubmargin(object = m, base = psub, delta = 5)
}

Description

Make a data set of all possible pairwise substitution of a composition which can be used as the base for substitution models.

Usage

build.base(parts, type = NULL)

Arguments

Value

A data table of all possible pairwise substitution.

Examples

ps1 <- build.base(parts = c("TST", "WAKE", "MVPA", "LPA", "SB"))
print(ps1)

ps2 <- build.base(c("WAKE", "MVPA", "LPA", "SB"), type = "one-to-all")
print(ps2)

Description

Build a dataset for fitted.brmcoda used in substitution model

Usage

build.rg(object, ref, at, parts, level, weight, fill = FALSE)

Arguments

Value

A reference grid consisting of a combination of covariates in brmcoda

Description

Build a default sequential binary partition for complr object. The default sequential binary partition is a pivot balance that allows the effect of this first balance coordinate to be interpreted as the change in the prediction for the dependent variable when that given part increases while all remaining parts decrease by a common proportion.

Usage

build.sbp(parts)

Arguments

Value

A matrix sequential binary partition.

Examples

sbp1 <- build.sbp(parts = c("TST", "WAKE", "MVPA", "LPA", "SB"))
print(sbp1)

sbp2 <- build.sbp(c("WAKE", "MVPA", "LPA", "SB"))
print(sbp2)

Description

Extract model coefficients, which are the sum of population-level effects and corresponding group-level effects of the brmsfit object in a brmcoda object.

Usage

## S3 method for class 'brmcoda'
coef(object, ...)

Arguments

Value

A list of 3D arrays (one per grouping factor). If summary is TRUE, the 1st dimension contains the factor levels, the 2nd dimension contains the summary statistics (see posterior_summary), and the 3rd dimension contains the group-level effects. If summary is FALSE, the 1st dimension contains the posterior draws, the 2nd dimension contains the factor levels, and the 3rd dimension contains the group-level effects.

See Also

coef.brmsfit

Examples

## fit a model
if(requireNamespace("cmdstanr")){
  m <- brmcoda(complr = complr(data = mcompd, sbp = sbp,
                               parts = c("TST", "WAKE", "MVPA", "LPA", "SB"),
                               idvar = "ID", total = 1440),
               formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 +
                                  wz1_1 + wz2_1 + wz3_1 + wz4_1 + (1 | ID),
                                  chain = 1, iter = 500,
                                  backend = "cmdstanr")

  ## extract population and group-level coefficients separately
  fixef(m)
  ranef(m)

  ## extract combined coefficients
  coef(m)
}

Description

Indices from a (dataset of) Multilevel Composition(s) (deprecated.)

Usage

compilr(...)

Arguments

Value

A complr object with at least the following elements.

See Also

complr

Description

Compute sets of compositions and log ratio transformation for multilevel compositional data

Usage

complr(data, parts, sbp = NULL, total = 1, idvar = NULL, transform = "ilr")

Arguments

Details

The ilr-transform maps the D-part compositional data from the simplex into non-overlapping subgroups in the (D-1)-dimension Euclidean space isometrically by using an orthonormal basis, thereby preserving the compositional properties and yielding a full-rank covariance matrix. ilr transformation should be preferred. However, the alr and clr are alternatives. The alr-transform maps a D-part composition in the Aitchison-simplex non-isometrically to a (D-1)-dimension Euclidian vectors, commonly treating the last part as the common denominator of the others. alr transformation does not rely on distance which breaks the constraint of compositional data. clr-transform maps a D-part composition in the Aitchison-simplex isometrically to a D-dimensional Euclidian vector subspace. clr transformation is not injetive, resulting in singular covariance matrices.

Value

A complr object with at least the following elements.

Examples

x1 <- complr(data = mcompd,
                parts = c("TST", "WAKE", "MVPA", "LPA", "SB"),
                idvar = "ID", total = 1440)
str(x1)

x2 <- complr(data = mcompd,
                parts = list(c("TST", "WAKE"), c("MVPA", "LPA", "SB")),
                total = list(c(480), c(960)),
                idvar = "ID",
                transform = "ilr")
str(x2)

x3 <- complr(data = mcompd, sbp = sbp,
                parts = c("TST", "WAKE", "MVPA", "LPA", "SB"),
                idvar = "ID",
                 transform = "ilr")
str(x3)

x_wide <- complr(data = mcompd[!duplicated(ID)], sbp = sbp,
                parts = c("TST", "WAKE", "MVPA", "LPA", "SB"))
str(x_wide)

Description

Create generator specifications for gamma and beta variables.

Usage

gen_gamma(
  vars,
  level = c("single", "level2", "multilevel"),
  fixed_intercept = NULL,
  ...,
  scale_fixed_intercept = NULL,
  random_cov = NULL
)

gen_beta(
  vars,
  level = c("single", "level2", "multilevel"),
  fixed_intercept = NULL,
  ...,
  scale_fixed_intercept = NULL,
  random_cov = NULL
)

Arguments

Value

An mlsim_generator_spec object for use in simulate_data().

See Also

Other predictor generators: count-generators, gen_categorical(), gen_custom(), gen_mvn(), gen_outcome(), gen_template()

Examples

sim <- simulate_data(
  n = 5,
  seed = 5,
  generators = list(
    positive = gen_gamma(
      "positive",
      fixed_intercept = log(1.5),
      scale_fixed_intercept = log(2)
    ),
    proportion = gen_beta(
      "proportion",
      fixed_intercept = stats::qlogis(0.4),
      scale_fixed_intercept = log(10)
    )
  )
)
sim$data

Description

Create generator specifications for binomial, Poisson, and negative binomial count variables.

Usage

gen_binomial(
  vars,
  level = c("single", "level2", "multilevel"),
  size = NULL,
  fixed_intercept = NULL,
  random_cov = NULL
)

gen_poisson(
  vars,
  level = c("single", "level2", "multilevel"),
  fixed_intercept = NULL,
  random_cov = NULL
)

gen_negbin(
  vars,
  level = c("single", "level2", "multilevel"),
  fixed_intercept = NULL,
  ...,
  scale_fixed_intercept = NULL,
  random_cov = NULL
)

Arguments

Details

Count generators use link-scale intercepts at every level. Negative-binomial size is exp(scale_fixed_intercept).

Value

An mlsim_generator_spec object for use in simulate_data().

See Also

Other predictor generators: continuous-generators, gen_categorical(), gen_custom(), gen_mvn(), gen_outcome(), gen_template()

Examples

sim <- simulate_data(
  n = 5,
  seed = 4,
  generators = list(
    events = gen_poisson("events", fixed_intercept = log(2)),
    successes = gen_binomial(
      "successes",
      size = 4,
      fixed_intercept = stats::qlogis(0.5)
    ),
    overdispersed = gen_negbin(
      "overdispersed",
      fixed_intercept = log(2),
      scale_fixed_intercept = log(3)
    )
  )
)
sim$data

Description

Extract Diagnostic Quantities from brmsfit Models in brmcoda

Usage

## S3 method for class 'brmcoda'
log_posterior(object, ...)

## S3 method for class 'brmcoda'
nuts_params(object, ...)

## S3 method for class 'brmcoda'
rhat(x, ...)

## S3 method for class 'brmcoda'
neff_ratio(object, ...)

Arguments

Value

The exact form of the output depends on the method.

See Also

log_posterior.brmsfit

nuts_params.brmsfit

rhat.brmsfit

neff_ratio.brmsfit

Description

This is a generic function for diagnostics methods in multilevelcoda. See diagnostics.complr for the method for complr objects.

Usage

diagnostics(x, ...)

Arguments

Value

diagnostics class object

Description

Uses robust minimum covariance determinant estimates to calculate squared Mahalanobis distances for the total, between, or within log-ratio coordinates stored in a complr object. This can be used to identify extreme values (outliers).

Usage

## S3 method for class 'complr'
diagnostics(
  x,
  level = c("total", "between", "within"),
  parts = 1,
  ev.perc = 0.001,
  extremevalues = c("theoretical", "empirical"),
  ...
)

Arguments

Value

A diagnostics object with the raw composition matrix in x, a data.table of robust distances in distance, the extreme-value cutoff in cutoff, the requested upper-tail proportion in ev.perc, the cutoff type in extremevalues, and the diagnosed coordinate level in levels. The distance element includes an is_extremevalue column indicating whether each fitted observation is flagged as an extreme value. For level = "between", diagnostics are fitted using one observation per idvar.

Examples

data(mcompd)
data(sbp)

ids <- unique(mcompd$ID)[1:20]
cilr <- complr(
  data = mcompd[ID %in% ids, .SD[1:3], by = ID],
  sbp = sbp,
  parts = c("TST", "WAKE", "MVPA", "LPA", "SB"),
  idvar = "ID",
  total = 1440
)

# One diagnostic object per coordinate level
total_diag <- diagnostics(cilr, level = "total")
between_diag <- diagnostics(cilr, level = "between")
within_diag <- diagnostics(cilr, level = "within")

# Use an empirical cutoff and a larger tail proportion
empirical_diag <- diagnostics(
  cilr,
  level = "between",
  ev.perc = .05,
  extremevalues = "empirical"
)

is.diagnostics(between_diag)
head(between_diag$distance)

Description

Index brmcoda objects

Usage

## S3 method for class 'brmcoda'
variables(x, ...)

## S3 method for class 'brmcoda'
nvariables(x, ...)

## S3 method for class 'brmcoda'
niterations(x)

## S3 method for class 'brmcoda'
nchains(x)

## S3 method for class 'brmcoda'
ndraws(x)

Arguments

See Also

variables.brmsfit

nvariables.brmsfit

niterations.brmsfit

nchains.brmsfit

ndraws.brmsfit

Description

Compute posterior draws of the expected value of the posterior predictive distribution of a brmsfit model in the brmcoda object. Can be performed for the data used to fit the model (posterior predictive checks) or for new data. By definition, these predictions have smaller variance than the posterior predictions performed by the predict.brmcoda method. This is because only the uncertainty in the expected value of the posterior predictive distribution is incorporated in the draws computed by fitted while the residual error is ignored there. However, the estimated means of both methods averaged across draws should be very similar.

Usage

## S3 method for class 'brmcoda'
fitted(object, scale = c("linear", "response"), parts = 1, summary = TRUE, ...)

Arguments

Value

An array of predicted mean response values. If summary = FALSE the output resembles those of posterior_epred.brmsfit.

If summary = TRUE the output depends on the family: For categorical and ordinal families, the output is an N x E x C array, where N is the number of observations, E is the number of summary statistics, and C is the number of categories. For all other families, the output is an N x E matrix. The number of summary statistics E is equal to 2 + length(probs): The Estimate column contains point estimates (either mean or median depending on argument robust), while the Est.Error column contains uncertainty estimates (either standard deviation or median absolute deviation depending on argument robust). The remaining columns starting with Q contain quantile estimates as specified via argument probs.

In multivariate models, an additional dimension is added to the output which indexes along the different response variables.

See Also

fitted.brmsfit

Examples

## fit a model
if(requireNamespace("cmdstanr")){
  ## compute composition and ilr coordinates
  x <- complr(data = mcompd, sbp = sbp,
                 parts = c("TST", "WAKE", "MVPA", "LPA", "SB"),
                 idvar = "ID", total = 1440)

  ## fit a model
  m1 <- brmcoda(complr = x,
                formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 +
                                   wz1_1 + wz2_1 + wz3_1 + wz4_1 + (1 | ID),
                chain = 1, iter = 500,
                backend = "cmdstanr")

  ## compute expected predictions
  epred <- fitted(m1)
  head(epred)

  ## fit a model with compositional outcome
  m2 <- brmcoda(complr = x,
                formula = mvbind(z1_1, z2_1, z3_1, z4_1) ~ Stress + Female + (1 | ID),
                chain = 1, iter = 500,
                backend = "cmdstanr")

  ## expected predictions on compositional scale
  epredcomp <- fitted(m2, scale = "response")
  head(epredcomp)
}

Description

Extract the population-level ('fixed') effects from the brmsfit object in a brmcoda object.

Usage

## S3 method for class 'brmcoda'
fixef(object, ...)

Arguments

Value

If summary is TRUE, a matrix returned by posterior_summary for the population-level effects. If summary is FALSE, a matrix with one row per posterior draw and one column per population-level effect.

See Also

fixef.brmsfit

Examples

## fit a model
if(requireNamespace("cmdstanr")){
  ## fit a model
  m <- brmcoda(complr = complr(data = mcompd, sbp = sbp,
                               parts = c("TST", "WAKE", "MVPA", "LPA", "SB"),
                               idvar = "ID", total = 1440),
               formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 +
                                  wz1_1 + wz2_1 + wz3_1 + wz4_1 + (1 | ID),
                                  chain = 1, iter = 500,
                                  backend = "cmdstanr")

  ## extract population-Level coefficients
  fixef(m)
}

Description

Create a generator specification for binary or multicategory categorical variables.

Usage

gen_categorical(
  vars,
  level = c("single", "level2", "multilevel"),
  categories = NULL,
  fixed_intercept = NULL,
  random_cov = NULL,
  reference = NULL,
  output = c("factor", "character", "integer"),
  ordered = FALSE
)

Arguments

Details

Multilevel categorical generators record their group-level random-effect draws as per-row truth columns named ⁠.mlsim_<var>_random_intercept_<category>⁠, one per non-reference category, with ⁠<var>⁠ and ⁠<category>⁠ sanitized by make.names(). See the Value section of simulate_data() for the full naming contract.

Value

An mlsim_generator_spec object for use in simulate_data().

See Also

Other predictor generators: continuous-generators, count-generators, gen_custom(), gen_mvn(), gen_outcome(), gen_template()

Examples

sim <- simulate_data(
  n = 6,
  seed = 3,
  generators = list(
    arm = gen_categorical(
      "arm",
      categories = c("control", "treatment"),
      fixed_intercept = stats::qlogis(0.5)
    )
  )
)
sim$data

Description

Wrap a custom generator function so it can participate in simulate_data().

Usage

gen_custom(
  vars,
  generator,
  level = c("single", "level2", "multilevel"),
  ...,
  metadata = list()
)

Arguments

Value

An mlsim_generator_spec object for use in simulate_data().

See Also

Other predictor generators: continuous-generators, count-generators, gen_categorical(), gen_mvn(), gen_outcome(), gen_template()

Examples

constant_generator <- function(context, vars, level, value = 1) {
  list(data = data.frame(value = rep(value, context$n_rows)), names = vars)
}

sim <- simulate_data(
  n = 3,
  generators = list(x = gen_custom("x", constant_generator, value = 7))
)
sim$data

Description

Create a generator specification for univariate or multivariate Gaussian variables, optionally transformed from ILR coordinates to compositional parts.

Usage

gen_mvn(
  vars,
  level = c("single", "level2", "multilevel"),
  fixed_intercept = NULL,
  residual_cov = NULL,
  random_cov = NULL,
  ...,
  scale_fixed_intercept = NULL,
  residual_cor = NULL,
  compositional = FALSE,
  parts = NULL,
  total = 1,
  keep_ilr = TRUE,
  sbp = NULL
)

Arguments

Details

Compositional MVN generators simulate ILR coordinates first, then use the SBP basis to transform the coordinates into positive composition parts that sum to total.

Multilevel generators additionally emit the latent between- and within-group components of each generated variable as visible columns. For non-compositional variables these are ⁠<var>_between⁠ (the group-level mean, fixed plus random intercept) and ⁠<var>_within⁠ (the row-level residual), with ⁠<var> = <var>_between + <var>_within⁠ exactly. For compositional generators the same decomposition is emitted on the ILR scale (⁠<ilr>_between⁠, ⁠<ilr>_within⁠), together with the true between composition ⁠<part>_between⁠, the closed back-transform of the group-level ILR means. All generated variables are labelled with column roles so that between() and within() terms in gen_outcome() formulas resolve to the latent components; roles (and hence between()/within()) apply to the ILR coordinates of a compositional generator, not to its parts. When keep_ilr = FALSE the observed ILR columns are dropped, but the between/within columns and their roles are still emitted for multilevel generators.

Multilevel generators also record their group-level random-effect draws as per-row truth columns: ⁠.mlsim_<var>_random_intercept⁠ for the location intercept and, when random_cov is a joint location-scale covariance, ⁠.mlsim_<var>_scale_random_intercept⁠ for the scale intercept. The other multilevel scalar generators (gen_negbin(), gen_gamma(), gen_beta(), gen_poisson(), gen_binomial()) share this machinery and emit the same columns. See the Value section of simulate_data() for the full naming contract.

Value

An mlsim_generator_spec object for use in simulate_data().

See Also

Other predictor generators: continuous-generators, count-generators, gen_categorical(), gen_custom(), gen_outcome(), gen_template()

Examples

sim <- simulate_data(
  n = 4,
  seed = 2,
  generators = list(
    x = gen_mvn("x", fixed_intercept = 0, residual_cov = 1),
    z = gen_mvn(
      c("z1", "z2"),
      fixed_intercept = c(0, 0),
      residual_cov = diag(2),
      compositional = TRUE,
      parts = c("sleep", "activity", "sedentary")
    )
  )
)
sim$data
rowSums(as.matrix(sim$data[, c("sleep", "activity", "sedentary"), with = FALSE]))

Description

Create an outcome generator for simulate_data(). With the default family = "gaussian", outcomes are simulated on a Gaussian scale, optionally as ILR coordinates that are back-transformed to strictly positive compositions, and ar1() in the location formula defines a residual VAR(1) process, not an observed lagged-outcome predictor. Non-Gaussian families ("poisson", "binomial", "negbin", "gamma", "beta") simulate a single outcome from the family distribution with the location linear predictor mapped through the family link.

Usage

gen_outcome(
  formula,
  scale = NULL,
  params,
  burnin = NULL,
  family = "gaussian",
  trials = NULL,
  composition = list(parts = NULL, total = 24, sbp = NULL, keep_ilr = TRUE),
  ar_stability = c("resample", "shrink", "error"),
  max_stability_attempts = 1000,
  shrink_target_radius = 0.98
)

Arguments

Value

An mlsim_generator_spec object for use in simulate_data().

Non-Gaussian families

Non-Gaussian outcomes are drawn independently per row from the family distribution after mapping the location linear predictor through the inverse link, with the auxiliary parameter (size, shape, or precision) taken as exp() of the scale linear predictor. ar1() is currently only supported for family = "gaussian". Group-level effects on the location and scale linear predictors are drawn jointly from params$random[[group]]$covariance exactly as for Gaussian outcomes. Because the links are nonlinear, fixed effects are conditional (subject-specific) effects: the marginal mean of the outcome is not the inverse link of the fixed-effect predictor alone. Extreme linear-predictor values that produce non-finite means or auxiliary parameters are an error; draws that land on a boundary the family likelihood cannot support (beta outcomes numerically at 0 or 1, gamma outcomes at 0) produce a warning.

Time spacing

When ar1() is present, time must be complete and equally spaced within each participant (or within the single series). Participants may have different start times, end times, and numbers of observations, and the simulator does not check or enforce equal spacing between participants: different participants may also use different step sizes. AR and VAR coefficients are defined per observation step, not per unit of real time, so dynamic parameters are only comparable across participants in real-time units when all participants share the same step size.

AR stability and realized moderator draws

Stability is enforced through the row-wise spectral radius of the assembled AR coefficient matrices for every observed row. When ar1() interacts with predictors (for example within(stress):ar1()), the row-specific AR matrix depends on the moderator values realized earlier in the generator pipeline. Stability is therefore a property of the AR parameters jointly with the realized predictor data, not of the parameters alone: the same AR parameter values may be accepted in one simulated data set and rejected in another with more extreme moderator draws, and the chance of an unstable row grows with the number of rows. ar_stability = "resample" redraws only group-level effects, so it cannot repair instability caused by the population-level part of a moderated AR term; that case errors instead.

Two further diagnostics guard the validity of accepted draws. First, when the AR matrices vary over time within a series (moderated ar1() with more than one outcome), a spectral radius below 1 for every row is necessary but not sufficient for stability of the time-varying process: products of individually stable matrices can still diverge. The simulator therefore also computes the spectral norm (largest singular value) of each row-specific AR matrix; a largest norm below 1 guarantees stability. When the matrices vary within a series and the largest spectral norm is 1 or more, stability cannot be guaranteed and a warning is issued; the result is recorded in the stability metadata as time_varying_stability_guaranteed. Second, ar_stability = "resample" truncates the group-level effect distribution: each group's entire random-effect vector (location, scale, and AR effects jointly) is redrawn until the AR part is stable, so when the acceptance rate is below 1 the realized random-effect covariance is smaller than params$random[[group]]$covariance requested. The realized covariance is always recorded in the generator metadata as realised_group_level_effect_covariance, and a message is emitted when the acceptance rate falls below 0.95.

Performance

The simulator is written to scale to large designs without changing the data-generating model. Innovations are drawn in one vectorized step from the fixed conditional correlation matrix and scaled by the row-wise conditional SDs (an exact draw from the row-specific innovation distribution, since its covariance is the SD-scaled correlation matrix). Row-specific AR matrices are assembled with vectorized array operations, spectral radii are computed once per unique AR matrix so row-constant AR designs cost one eigendecomposition per participant, and stability resampling or shrinkage re-evaluates only the affected participant's rows. Because the order in which random numbers are consumed is part of the implementation, a given seed maps to a particular realization only within a package version; the distribution of simulated data is unaffected.

See Also

Other predictor generators: continuous-generators, count-generators, gen_categorical(), gen_custom(), gen_mvn(), gen_template()

Examples

beta_location <- matrix(
  c(0, 0, 0.2, -0.1),
  nrow = 2,
  dimnames = list(c("(Intercept)", "treatmenttreatment"), c("ilr1", "ilr2"))
)
beta_scale <- matrix(
  log(c(0.4, 0.35)),
  nrow = 1,
  dimnames = list("(Intercept)", c("ilr1", "ilr2"))
)
beta_ar <- array(
  c(0.25, 0.02, -0.01, 0.2),
  dim = c(1, 2, 2),
  dimnames = list("ar1()", c("ilr1", "ilr2"), c("ilr1", "ilr2"))
)
corr <- diag(2)
dimnames(corr) <- list(c("ilr1", "ilr2"), c("ilr1", "ilr2"))

sim <- simulate_data(
  n_groups = 4,
  n_per_group = 4,
  group_id = "ID",
  time_id = "day",
  seed = 1,
  generators = list(
    treatment = gen_categorical(
      "treatment",
      level = "level2",
      categories = c("control", "treatment"),
      fixed_intercept = stats::qlogis(0.5)
    ),
    outcome = gen_outcome(
      mvbind(ilr1, ilr2) ~ treatment + ar1(),
      scale = sigma ~ 1,
      params = list(
        location = list(beta = beta_location),
        scale = list(beta = beta_scale, correlation = corr),
        ar = list(beta = beta_ar)
      ),
      burnin = 10,
      composition = list(parts = c("sleep", "active", "sedentary"), total = 24)
    )
  )
)
head(sim$data)

Description

Build a generator that parses a gen_outcome() specification and records a complete parameter template without simulating outcome values. Use this inside simulate_data() after any generators that define predictors used by the outcome formula, especially predictors referenced by between() or within().

Usage

gen_template(
  formula,
  scale = NULL,
  burnin = NULL,
  family = "gaussian",
  composition = list(parts = NULL, total = 24, sbp = NULL, keep_ilr = TRUE),
  ar_stability = c("resample", "shrink", "error"),
  max_stability_attempts = 1000,
  shrink_target_radius = 0.98
)

Arguments

Details

The generated template is stored at sim$generator_metadata[[name]]$params, where name is the name given to this generator in the generators list. The template must be created with the same simulation design, previous generators, factor levels, formula, scale formula, and composition settings that will be used for the later gen_outcome() call.

Template values are initialized to zero for regression, AR, and group-level covariance parameters, and to identity for residual ILR correlations. The object is ready to edit and pass as params to gen_outcome().

Templates are family-aware: params$scale is only included for families with an auxiliary scale model ("gaussian", "negbin", "gamma", "beta"), params$scale$correlation only for "gaussian", and params$ar only when ar1() appears (which requires family = "gaussian"). trials is not needed to build a binomial template because no parameter block depends on it.

Value

An mlsim_generator_spec object for use in simulate_data(). It emits no data columns and stores the parameter template in generator metadata.

See Also

Other predictor generators: continuous-generators, count-generators, gen_categorical(), gen_custom(), gen_mvn(), gen_outcome()

Examples

template_sim <- simulate_data(
  n_groups = 3,
  n_per_group = 2,
  group_id = "ID",
  seed = 1,
  generators = list(
    treatment = gen_categorical(
      "treatment",
      level = "level2",
      categories = c("control", "treatment"),
      fixed_intercept = stats::qlogis(0.5)
    ),
    outcome_template = gen_template(
      mvbind(ilr1, ilr2) ~ treatment,
      scale = sigma ~ 1,
      composition = list(parts = c("sleep", "active", "sedentary"), total = 24)
    )
  )
)

params <- template_sim$generator_metadata$outcome_template$params
params$location$beta["treatmenttreatment", "ilr1"] <- 0.2

Description

Extract Sequential Binary Partition from a complr object.

Usage

get_sbp(object)

Arguments

Description

Generic function to extract variable names from a supported object.

Usage

get_variables(object)

## S3 method for class 'complr'
get_variables(object)

## S3 method for class 'brmcoda'
get_variables(object)

Arguments

Value

A list of variable names.

Examples

# For a complr object:
# get_variables(complr_object)

# For a brmcoda object:
# get_variables(brmcoda_object)

Description

Functions used only internally to estimate substitution model

Description

Checks if argument is a brmcoda object

Usage

is.brmcoda(x)

Arguments

Description

Checks if argument is a complr object

Usage

is.complr(x)

Arguments

Description

Checks if argument is a diagnostics object

Usage

is.diagnostics(x)

Arguments

Description

Checks if argument is a substitution object

Usage

is.substitution(x)

Arguments

Description

Create lagged versions of one or more columns using the actual time values rather than row position. Each row's lag is taken from the row (within the same group, when group is supplied) whose time value equals time - n * time_step. When no observation exists at that time – for example, a skipped day in a daily diary – the lag is NA instead of the value from before the gap. This makes the lags robust to missing time points, unlike positional shifts such as data.table::shift().

Usage

lag_by_time(
  data,
  cols,
  time,
  group = NULL,
  time_step = NULL,
  n = 1L,
  suffix = "_lag"
)

Arguments

Details

Rows are matched on a per-group integer time grid: within each group, time values are converted to step counts from the earliest observation, (time - min(time)) / time_step, and a row's lag source is the row n grid points earlier. A value is treated as on the grid when it is within 1e-8 * max(1, time_step) (on the time scale) of a whole number of steps. Off-grid rows – irregular spacing that is not a multiple of the step – silently receive NA lags and are counted in n_irregular; if this affects many rows, supply the intended time_step explicitly.

Gaps are silent by design: no message or warning is emitted when a lag is NA because the previous time point is missing. Inspect the "lag_by_time" attribute to see how many rows were affected.

When time_step is inferred, at least one group must contain two or more time points; otherwise an error asks for an explicit time_step. Duplicate time values within a group are an error.

Value

data with one new column per entry of cols, named paste0(cols, suffix), appended in the original row order. The input is never modified in place; a data frame input returns a data frame and a data.table input returns a new data.table. The result carries an attribute "lag_by_time", a list with:

See Also

prep_sim_analysis(), which uses lag_by_time() to build its ar1() lag columns.

Examples

# day 3 is missing, so the lag on day 4 is NA rather than the day-2 value
d <- data.frame(
  id = c(1, 1, 1, 2, 2, 2),
  day = c(1, 2, 4, 1, 2, 3),
  y = c(10, 11, 12, 20, 21, 22)
)
lag_by_time(d, cols = "y", time = "day", group = "id")

# Date time with an explicit one-day step
d2 <- data.frame(date = as.Date("2026-01-01") + c(0, 1, 3), y = 1:3)
lag_by_time(d2, cols = "y", time = "date", time_step = 1)

Description

Provide an interface to shinystan for models fitted with brms

Usage

## S3 method for class 'brmcoda'
launch_shinystan(object, ...)

Arguments

Value

An S4 shinystan object

See Also

launch_shinystan

Description

Perform approximate leave-one-out cross-validation based on the posterior likelihood using the loo package. For more details see loo.

Usage

## S3 method for class 'brmcoda'
loo(x, ...)

Arguments

Value

If just one object is provided, an object of class loo. If multiple objects are provided, an object of class loolist.

See Also

loo.brmsfit

Description

Call MCMC plotting functions implemented in the bayesplot package.

Usage

## S3 method for class 'brmcoda'
mcmc_plot(object, ...)

Arguments

Value

A ggplot object that can be further customized using the ggplot2 package.

See Also

mcmc_plot.brmsfit

Examples

## Not run: 
cilr <- complr(data = mcompd, sbp = sbp,
        parts = c("TST", "WAKE", "MVPA", "LPA", "SB"), idvar = "ID")

# model with compositional predictor at between and within-person levels
fit <- brmcoda(complr = cilr,
                formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 +
                                   wz1_1 + wz2_1 + wz3_1 + wz4_1 + (1 | ID),
               chain = 1, iter = 500)
mcmc_plot(fit)

## End(Not run)

Description

A simulated dataset containing multiple days of compositional data.

Usage

mcompd

Format

A data table containing 10 variables.

ID

A unique identifier for each individual

Time

Recurrence time of repeated measures by individual

Stress

Self report stress measures on a 0 to 10 scale — repeated measure

TST

Total Sleep Time (minutes) — repeated measure

WAKE

Wake time while in bed, trying to sleep (minutes) — repeated measure

MVPA

Moderate to Vigorous Physical Activity (minutes) — repeated measure

LPA

Light Physical Activity (minutes) — repeated measure

SB

Sedentary Behavior (minutes) — repeated measure

Age

Age in years — baseline measure only

Female

Binary: whether participants identified as female (1) or not (0) — baseline measure only

Description

Mean amounts and mean compositions presented in a complr object.

Usage

## S3 method for class 'complr'
mean(x, weight = c("equal", "proportional"), parts = 1, ...)

Arguments

Examples

x <- complr(data = mcompd, sbp = sbp,
                parts = c("TST", "WAKE", "MVPA", "LPA", "SB"),
                idvar = "ID")
mean(x)

Description

Extracting the Model Frame from a Formula or Fit from brmcoda object

Usage

## S3 method for class 'brmcoda'
model.frame(formula, ...)

Arguments

Description

Provide the full results for a simulation study testing the performance of multilevelcoda

Usage

multilevelcoda_sim()

Value

An S4 shiny object

Description

Extract Number of Observations from brmcoda object

Usage

## S3 method for class 'brmcoda'
nobs(object, ...)

Arguments

Create a matrix of output plots from a brmcoda's brmsfit object

Description

A pairs method that is customized for MCMC output.

Usage

## S3 method for class 'brmcoda'
pairs(x, ...)

Arguments

See Also

pairs.brmsfit

Examples

## Not run: 
cilr <- complr(data = mcompd, sbp = sbp,
        parts = c("TST", "WAKE", "MVPA", "LPA", "SB"), idvar = "ID")

# model with compositional predictor at between and within-person levels
fit <- brmcoda(complr = cilr,
                formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 +
                                   wz1_1 + wz2_1 + wz3_1 + wz4_1 + (1 | ID),
               chain = 1, iter = 500)
pairs(fit)

## End(Not run)

Description

This function estimates pivot balance coordinates for each compositional part by either "rotate" the sequential binary partition using the same brmcoda object or "refit" the brmcoda object.

Usage

pivot_coord(
  object,
  summary = TRUE,
  method = c("rotate", "refit"),
  parts = 1,
  ...
)

Arguments

Value

Estimated pivot balance coordinates representing the effect of increasing one compositional part relative to the remaining compositional parts.

Examples

if(requireNamespace("cmdstanr")){
  x <- complr(data = mcompd, sbp = sbp,
                 parts = c("TST", "WAKE", "MVPA", "LPA", "SB"), idvar = "ID",
                 total = 1440)
  
  # inspects ILRs before passing to brmcoda
  names(x$between_logratio)
  names(x$within_logratio)
  names(x$logratio)
  
  # model with compositional predictor at between and within-person levels
  m <- brmcoda(complr = x,
                formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 +
                                   wz1_1 + wz2_1 + wz3_1 + wz4_1 + (1 | ID),
                chain = 1, iter = 500,
                backend = "cmdstanr")
  
  m_pivot_coord <- pivot_coord(m)
  summary(m_pivot_coord)
  }

Description

This function is an alias of pivot_coord to estimates the pivot balance coordinates by "refit" the brmcoda object.

Usage

pivot_coord_refit(object, summary = TRUE, parts = 1, ...)

Arguments

Value

Estimated pivot balance coordinates representing the effect of increasing one compositional part relative to the remaining compositional parts.

See Also

pivot_coord

Examples

if(requireNamespace("cmdstanr", "brms")){
  x <- complr(data = mcompd, sbp = sbp,
                 parts = c("TST", "WAKE", "MVPA", "LPA", "SB"), idvar = "ID",
                 total = 1440)
  
  m <- brmcoda(complr = x,
                formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 +
                                   wz1_1 + wz2_1 + wz3_1 + wz4_1 + (1 | ID),
                chain = 1, iter = 500,
                backend = "cmdstanr")
  
  m_pivot_coord_refit <- pivot_coord_refit(m)
  summary(m_pivot_coord_refit)
  
  m_pivot_coord_raw <-  pivot_coord_refit(m, summary = FALSE)
  lapply(m_pivot_coord_raw$output, brms::posterior_summary)
  
  }

Description

This function is an alias of pivot_coord to estimates the pivot balance coordinates by "rotate" the sequential binary partition on the same brmcoda object.

Usage

pivot_coord_rotate(object, summary = TRUE, parts = 1, ...)

Arguments

Value

Estimated pivot balance coordinates representing the effect of increasing one compositional part relative to the remaining compositional parts.

See Also

pivot_coord

Examples

if(requireNamespace("cmdstanr")){
  x <- complr(data = mcompd, sbp = sbp,
                 parts = c("TST", "WAKE", "MVPA", "LPA", "SB"), idvar = "ID",
                 total = 1440)
  
  m <- brmcoda(complr = x,
                formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 +
                                   wz1_1 + wz2_1 + wz3_1 + wz4_1 + (1 | ID),
                chain = 1, iter = 500,
                backend = "cmdstanr")
  
  m_pivot_coord_rotate <- pivot_coord_rotate(m)
  summary(m_pivot_coord_rotate)
  
  m_pivot_coord_raw <-  pivot_coord_rotate(m, summary = FALSE)
  lapply(m_pivot_coord_raw$output, brms::posterior_summary)
  
  }

Description

Make a plot of brmcoda model results.

Usage

## S3 method for class 'brmcoda'
plot(x, ...)

Arguments

Value

An invisible list of gtable objects.

See Also

plot.brmsfit

Examples

## Not run: 
cilr <- complr(data = mcompd, sbp = sbp,
        parts = c("TST", "WAKE", "MVPA", "LPA", "SB"), idvar = "ID")

# model with compositional predictor at between and within-person levels
fit <- brmcoda(complr = cilr,
                formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 +
                                   wz1_1 + wz2_1 + wz3_1 + wz4_1 + (1 | ID),
               chain = 1, iter = 500)
plot(fit)

## End(Not run)

Description

Make a scatterplot matrix of diagnostics results from a complr object, with density and rug plots on the diagonal.

Usage

## S3 method for class 'diagnostics'
plot(x, transform = c("clr", "raw"), ...)

Arguments

Value

A ggplot graph object showing pairwise scatterplots for each composition part, with point colour indicating extreme value status. Diagonal panels show marginal densities and rug marks for observations flagged as extreme values.

Examples

data(mcompd)
data(sbp)

ids <- unique(mcompd$ID)[1:20]
cilr <- complr(
  data = mcompd[ID %in% ids, .SD[1:3], by = ID],
  sbp = sbp,
  parts = c("TST", "WAKE", "MVPA", "LPA", "SB"),
  idvar = "ID",
  total = 1440
)

between_diag <- diagnostics(cilr, level = "between", ev.perc = .05)

# Plot centered logratio transformed composition values
plot(between_diag)

# Plot on the raw composition scale
plot(between_diag, transform = "raw")

Description

Make a plot of substitution model results.

Usage

## S3 method for class 'substitution'
plot(x, to, ref, level, ...)

Arguments

Value

A ggplot graph object showing the estimated difference in outcome when each pair of compositional variables are substituted for a specific time.

Description

Perform posterior predictive checks with the help of the bayesplot package.

Usage

## S3 method for class 'brmcoda'
pp_check(
  object,
  type = "dens_overlay",
  ndraws = NULL,
  prefix = c("ppc", "ppd"),
  newdata = NULL,
  draw_ids = NULL,
  nsamples = NULL,
  subset = NULL,
  scale = c("linear", "response"),
  parts = NULL,
  ...
)

Arguments

See Also

pp_check.brmsfit

Examples

if(requireNamespace("cmdstanr")){
  ## fit a model
  x <- complr(data = mcompd, sbp = sbp,
                 parts = c("TST", "WAKE", "MVPA", "LPA", "SB"),
                 idvar = "ID", total = 1440)

  m1 <- brmcoda(complr = x,
                formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 +
                                   wz1_1 + wz2_1 + wz3_1 + wz4_1 + (1 | ID),
                chain = 1, iter = 500,
                backend = "cmdstanr")

  ## posterior predictive checks
  bayesplot::pp_check(m1, ndraws = 5)

  ## fit a model with compositional outcome
  m2 <- brmcoda(complr = x,
                formula = mvbind(z1_1, z2_1, z3_1, z4_1) ~ 
                          bz1_1 + bz2_1 + bz3_1 + bz4_1 + Female + (1 | ID),
                chain = 1, iter = 500,
                backend = "cmdstanr")

  ## posterior predictive checks for compositional outcome -- linear scale
  bayesplot::pp_check(m2, resp = "z11", ndraws = 10)
  ## posterior predictive checks for compositional outcome -- original response scale
  bayesplot::pp_check(m2, parts = "WAKE", scale = "response", ndraws = 10)
}

Description

Compute posterior draws of the posterior predictive distribution of a brmsfit model in the brmcoda object. Can be performed for the data used to fit the model (posterior predictive checks) or for new data. By definition, these draws have higher variance than draws of the expected value of the posterior predictive distribution computed by fitted.brmcoda. This is because the residual error is incorporated in posterior_predict. However, the estimated means of both methods averaged across draws should be very similar.

Usage

## S3 method for class 'brmcoda'
predict(object, scale = c("linear", "response"), parts = 1, ...)

Arguments

Value

An array of predicted response values. If summary = FALSE the output resembles those of posterior_predict.brmsfit.

If summary = TRUE the output depends on the family: For categorical and ordinal families, the output is an N x C matrix, where N is the number of observations, C is the number of categories, and the values are predicted category probabilities. For all other families, the output is a N x E matrix where E = 2 + length(probs) is the number of summary statistics: The Estimate column contains point estimates (either mean or median depending on argument robust), while the Est.Error column contains uncertainty estimates (either standard deviation or median absolute deviation depending on argument robust). The remaining columns starting with Q contain quantile estimates as specified via argument probs.

See Also

predict.brmsfit

Examples

if(requireNamespace("cmdstanr")){
  ## fit a model
  x <- complr(data = mcompd, sbp = sbp,
                 parts = c("TST", "WAKE", "MVPA", "LPA", "SB"),
                 idvar = "ID", total = 1440)

  m1 <- brmcoda(complr = x,
                formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 +
                                   wz1_1 + wz2_1 + wz3_1 + wz4_1 + (1 | ID),
                chain = 1, iter = 500,
                backend = "cmdstanr")

  ## predicted responses
  pred <- predict(m1)
  head(pred)

  ## fit a model with compositional outcome
  m2 <- brmcoda(complr = x,
                formula = mvbind(z1_1, z2_1, z3_1, z4_1) ~ 
                          bz1_1 + bz2_1 + bz3_1 + bz4_1 + Female + (1 | ID),
                chain = 1, iter = 500,
                backend = "cmdstanr")

  ## predicted responses on ilr scale
  predilr <- predict(m2, scale = "linear")
  head(predilr)

  ## predicted responses on compositional scale
  predcomp <- predict(m2, scale = "response")
  head(predcomp)
}

Description

Convert an simulate_data() result with a gen_outcome() generator into an analysis-ready data set and inferred brms formula. The helper recomputes observed between- and within-person predictors from the raw simulated columns, creates lagged outcome predictor columns for ar1() terms (person-mean-centered by default; see lag_center), and creates a complr() object when the outcome is compositional or when between()/within() terms reference ILR coordinates of a compositional predictor generator.

Usage

prep_sim_analysis(
  sim,
  outcome = NULL,
  drop_lag_na = FALSE,
  time_step = NULL,
  lag_center = c("within", "none"),
  centering = c("manifest", "latent"),
  link_random = TRUE
)

Arguments

Details

The inferred analysis model is a deliberately pragmatic default estimator built from observed data, not the matched model for the simulated data-generating process. See the "Pragmatic default estimator" section before interpreting parameter recovery results.

The analysis formula is inferred from the stored gen_outcome() formula. Simulation terms between(x) and within(x) become observed-data columns named x_between and x_within, computed from x by the simulation group identifier. These columns are recomputed even when columns with the same names already exist in sim$data.

When between()/within() reference ILR coordinates of a compositional gen_mvn() predictor, the helper instead builds a complr() object for that predictor composition (using the simulator's parts, SBP, and total) and maps the terms to the complr between/within coordinates (⁠bz<k>_<j>⁠ and ⁠wz<k>_<j>⁠, where k is the coordinate and j the composition index), exactly as in the standard multilevelcoda workflow. When the outcome is also compositional, one multi-composition complr object covers both the predictor and outcome compositions.

The returned brms formula carries the family recorded by the outcome generator (gaussian(), poisson(), binomial(), brms::negbinomial(), Gamma(link = "log"), or brms::Beta()). When the generator has a scale model, it becomes the matching distributional formula (sigma, shape, or phi). Binomial outcomes use ⁠y | trials(<outcome>_trials)⁠ with the trials column generated by gen_outcome(). When the same grouping factor appears in both the mean and the scale formulas, the random effects are emitted with brms ID-linked syntax (for example (1 | p1 | ID)) by default, because the simulator draws these group-level effects from one joint covariance and the linked syntax lets the analysis model estimate their correlation; see link_random to opt out for large models fitted with approximate algorithms.

For dynamic formulas, ar1() is translated to lagged response columns whose construction and name follow the two centering arguments:

• lag_center = "within", centering = "manifest" (both defaults): the lagged observed response centered at each person's observed mean of all their values (not the mean of the lagged values), ⁠lag_<response>_within⁠.

• lag_center = "within", centering = "latent": the lagged response centered at the latent conditional mean, which is the recorded latent residual state, ⁠lag_<response>_latent⁠.

• lag_center = "none", either centering: the raw lagged observed response, ⁠lag_<response>⁠.

For compositional outcomes, the helper rebuilds the ILR coordinates through complr() using the simulator's parts and SBP metadata, then lags the generated z coordinates used by brmcoda().

The lags are built by lag_by_time() as time-based lags on the simulation time_id: a row's lag comes from the observation at time - time_step, not from the previous row. When a time point is missing – as can happen when analysing real data with skipped observations – the following row's lag is NA rather than the value from before the gap. The autoregressive coefficient therefore refers to one time_step. For data simulated with ar1() the time grid is complete and equally spaced, so the result is identical to a positional shift. The number of gap-affected rows is recorded in metadata$lag_gaps.

Because a single time_step is applied to all groups, groups whose own spacing is wider than the inferred step end up with only NA lags, and drop_lag_na = TRUE would then remove them from the analysis data entirely. prep_sim_analysis() warns in both situations and records the diagnostics in the metadata: metadata$time_step_by_group (each group's smallest positive spacing, NA for single-observation groups), metadata$time_step_heterogeneous (whether those spacings differ), metadata$lag_na_groups (groups without a single usable lagged row), and metadata$dropped_groups (groups removed by drop_lag_na).

Value

An mlsim_analysis object, a list with:

Pragmatic default estimator

gen_outcome() simulates latent residual AR/VAR dynamics around the model-implied mean and resolves between()/within() from latent generating components supplied by upstream predictor generators. By default (centering = "manifest") prep_sim_analysis() instead constructs the model applied analysts commonly fit to observed data:

• between(x) and within(x) become x_between and x_within, recomputed from realised person means of the observed x (manifest centering), not from the latent generating components.

• ar1() becomes, by default, person-mean-centered lagged observed response predictors (⁠lag_<response>_within⁠), not the latent residual state.

• For ILR coordinates of compositional predictors, between(ilr) and within(ilr) become the complr() coordinates ⁠bz*⁠/⁠wz*⁠, where the between composition is the ILR of each person's closed arithmetic-mean composition of the observed parts. The simulator's between(ilr), in contrast, is the latent group-level ILR mean, whose back-transform is the closed geometric (Aitchison) center. These are different estimands, and the gap has two parts: a structural Jensen gap between the arithmetic and geometric centers, which does not shrink as the number of observations per group grows, plus ordinary sampling error, which does. Manifest centering is therefore biased for latent between-person effects even in long series, unlike the scalar case where only sampling error is involved.

These observed-data constructions target different estimands from the simulation truth. Manifest person-mean centering and observed-score lagged regression are known to yield biased estimates of between-person effects and of inertia and cross-lag parameters relative to the latent generating values, especially with short series (Ludtke et al. 2008; Hamaker & Grasman 2014). This mismatch is intentional: it lets simulation studies quantify the bias of the pragmatic estimator. Do not interpret systematic discrepancies between estimates from this default analysis model and the gen_outcome() truth parameters as errors in the simulator.

Latent centering

centering = "latent" replaces the manifest construction with the generating components themselves, so for the affected between()/within() terms the fitted model targets exactly the quantity the simulator recorded:

• Scalar predictors take x_between from the column the generator labelled with component "between" – for level = "level2" predictors that is the variable itself – and set x_within <- x - x_between.

• Compositional predictors keep the returned complr() object, but its between and within blocks (bX, wX, bZ, wZ, and the matching ⁠b<part>⁠/⁠w<part>⁠/⁠bz*⁠/⁠wz*⁠ columns) are rebuilt from the latent between composition emitted by the generator. Column names are unchanged, so brmcoda() and substitution() work as usual.

This is an oracle: the latent components exist only in simulated data, so "latent" has no counterpart in a real analysis. Use it to separate estimator bias from simulator error – if a parameter is recovered under "latent" but not under "manifest", the gap is the estimator's, not the simulator's. It requires a grouped design and a generator that emitted latent components; gen_custom() generators that label roles without emitting them can only be analysed with "manifest".

Latent centering and ar1()

When the outcome formula contains ar1() and the lag is centered (lag_center = "within", the default), centering = "latent" centers the lag with latent values, which replaces the lag predictor entirely: instead of a lagged observed response centered at the person's observed mean, the analysis uses the latent residual state the simulator recorded, in columns named ⁠lag_<response>_latent⁠. The two constructions are the same operation applied with different values, because the recorded residual is y_t - mu_t, the response centered at its latent conditional mean. The latent state needs no further centering: person-mean centering it on top would reintroduce exactly the attenuation this construction avoids.

With lag_center = "none" there is no lag centering for centering to act on, so the raw lagged observed response is used under "latent" just as it is under "manifest". That combination is the point of keeping the two arguments separate: a model can carry between()/within() predictors that have nothing to do with the lag, and those can be built from the latent generating components while the AR part keeps the non-centered parametrization discussed below.

A latent lag makes the AR, scale and rescor parameters estimate exactly what the simulator recorded only when the analysis model has no group-level effects. The reason is not the proxy problem that lag_center addresses. Writing the generating model as ⁠z_t = X_t beta + Phi e_{t-1} + eps_t⁠, the regressor ⁠e_{t-1}⁠ is predetermined but not strictly exogenous: Cov(eps_t, e_t) equals the innovation variance, and e_t is the next row's regressor. Any estimator that mixes rows within a group – which is what the random-effects GLS transformation implied by a (1 | group) term does – picks up that cross-row covariance and incurs an O(1/T) dynamic-panel (Nickell) bias. Substituting the true latent state removes the measurement problem but not this one. So with any group-level effect, no parameter of a dynamic model is guaranteed to keep its generating estimand, under either centering.

For a matched-model recovery of a multilevel dynamic model, the appropriate analysis is not lag regression at all but a residual autocorrelation structure, for example brms::bf(y ~ 1 + (1 | ID)) + brms::ar(time = day, gr = ID, p = 1, cov = TRUE). That form is outside what this helper builds – brms ar() supports only a scalar unmoderated AR per response and parametrizes sigma as the marginal rather than the innovation standard deviation – so construct it by hand.

Note that under "manifest" the latent columns of the same name are overwritten in analysis$data; sim$data always keeps the originals.

Whether the ar1() lag columns are centered at all is controlled by lag_center, independently of centering. Hamaker & Grasman (2014) show that in multilevel autoregressive models, cluster-mean centering the lagged outcome (the "within" default) attenuates the average autoregressive coefficient, whereas the non-centered lagged outcome ("none") recovers it nearly unbiased. The trade-off is that under "none" the whole location mean structure is reparametrized: the intercept estimates ⁠(I - Phi_i) mu_i⁠ rather than the person mean, and any location covariate effects absorb omitted ⁠-Phi x_{t-1} beta⁠ terms. Centering does not fully resolve this either: the centered observed lag still carries the lagged mean structure, so when a predictor varies within series its current coefficient can absorb omitted lagged-covariate terms under "within" as well. Whenever the formula contains ar1() and the lag is an observed response, then, no parameter is guaranteed to keep its generating estimand, under either lag_center setting. Use "none" when the average autoregressive (inertia and cross-lag) coefficients are the estimands of interest, and the default "within" when person means and their predictors must remain interpretable; the data-simulation vignette compares parameter recovery under both.

References

Ludtke, O., et al. (2008). The multilevel latent covariate model. Psychological Methods, 13(3), 203-229.

Hamaker, E. L., & Grasman, R. P. P. P. (2014). To center or not to center? Investigating inertia with a multilevel autoregressive model. Frontiers in Psychology, 5, 1492.

Examples

params <- list(
  location = list(beta = matrix(0, nrow = 1, dimnames = list("(Intercept)", "y"))),
  scale = list(beta = matrix(log(0.2), nrow = 1, dimnames = list("(Intercept)", "y")))
)
sim <- simulate_data(
  n = 5,
  seed = 1,
  generators = list(
    outcome = gen_outcome(
      y ~ 1,
      scale = sigma ~ 1,
      params = params
    )
  )
)
analysis <- prep_sim_analysis(sim)
analysis$formula

Description

Print a Summary for a fitted brmsfit model in a brmcoda object

Usage

## S3 method for class 'brmcoda'
print(x, ...)

Arguments

See Also

summary.brmcoda

Examples

if(requireNamespace("cmdstanr")){
  m <- brmcoda(complr = complr(data = mcompd, sbp = sbp,
                               parts = c("TST", "WAKE", "MVPA", "LPA", "SB"),
                               idvar = "ID", total = 1440),
  formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 +
                     wz1_1 + wz2_1 + wz3_1 + wz4_1 + (1 | ID),
  chain = 1, iter = 500,
  backend = "cmdstanr")

 print(m)
}

Description

Print a Summary for a complr object

Usage

## S3 method for class 'complr'
print(x, ...)

Arguments

See Also

summary.complr

Examples

cilr <- complr(data = mcompd, sbp = sbp,
                parts = c("TST", "WAKE", "MVPA", "LPA", "SB"),
                idvar = "ID")
print(cilr)

Description

Print a compact overview of an mlsim_data object.

Usage

## S3 method for class 'mlsim_data'
print(x, ...)

Arguments

Value

x, invisibly.

Examples

sim <- simulate_data(
  n = 3,
  generators = list(
    x = gen_mvn("x", fixed_intercept = 0, residual_cov = 1)
  )
)
print(sim)

Description

Print a Summary for a substitution object

Usage

## S3 method for class 'substitution'
print(x, ...)

Arguments

See Also

summary.substitution

Examples

if(requireNamespace("cmdstanr")){
  ## fit a model with compositional predictor at between and between-person levels
  m <- brmcoda(complr = complr(data = mcompd, sbp = sbp,
                               parts = c("TST", "WAKE", "MVPA", "LPA", "SB"),
                               idvar = "ID", total = 1440),
  formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 +
                     wz1_1 + wz2_1 + wz3_1 + wz4_1 + (1 | ID),
  chain = 1, iter = 500,
  backend = "cmdstanr")

  subm <- substitution(object = m, delta = 5)
  print(subm)
}

Description

Print a Simulated Data Summary

Usage

## S3 method for class 'summary.mlsim_data'
print(x, ...)

Arguments

Value

x, invisibly.

Description

Compute Bayes factors from marginal likelihoods

Usage

## S3 method for class 'brmcoda'
prior_summary(object, ...)

Arguments

See Also

prior_summary.brmsfit

Description

A dataset containing possible pairwise subsitutions.

Usage

psub

Format

A data table containing 5 variables.

TST

first compositional variable

WAKE

second compositional variable

MVPA

third compositional variable

LPA

fourth compositional variable

SB

fifth compositional variable

Description

Extract the group-level ('random') effects of each level of the brmsfit object in a brmcoda object.

Usage

## S3 method for class 'brmcoda'
ranef(object, ...)

Arguments

Value

A list of 3D arrays (one per grouping factor). If summary is TRUE, the 1st dimension contains the factor levels, the 2nd dimension contains the summary statistics (see posterior_summary), and the 3rd dimension contains the group-level effects. If summary is FALSE, the 1st dimension contains the posterior draws, the 2nd dimension contains the factor levels, and the 3rd dimension contains the group-level effects.

See Also

ranef.brmsfit

Examples

## fit a model
if(requireNamespace("cmdstanr")){
  m <- brmcoda(complr = complr(data = mcompd, sbp = sbp,
                               parts = c("TST", "WAKE", "MVPA", "LPA", "SB"),
                               idvar = "ID", total = 1440),
               formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 +
                                  wz1_1 + wz2_1 + wz3_1 + wz4_1 + (1 | ID),
                                  chain = 1, iter = 500,
                                  backend = "cmdstanr")

  ## extract group-level coefficients
  ranef(m)
}

Description

Compute posterior draws of residuals/predictive errors

Usage

## S3 method for class 'brmcoda'
residuals(object, ...)

Arguments

Value

An array of predictive error/residual draws. If summary = FALSE the output resembles those of predictive_error.brmsfit. If summary = TRUE the output is an N x E matrix, where N is the number of observations and E denotes the summary statistics computed from the draws.

See Also

residuals.brmsfit

Examples

## fit a model
if(requireNamespace("cmdstanr")){
  m <- brmcoda(complr = complr(data = mcompd, sbp = sbp,
                               parts = c("TST", "WAKE", "MVPA", "LPA", "SB"),
                               idvar = "ID", total = 1440),
               formula = Stress ~ bz1_1 + bz2_1 + bz3_1 + bz4_1 +
                                  wz1_1 + wz2_1 + wz3_1 + wz4_1 + (1 | ID),
                                  chain = 1, iter = 500,
                                  backend = "cmdstanr")

  ## extract residuals
  res <- residuals(m)
  head(res)
}

Description

A matrix of sequential binary partition.

Usage

sbp

Format

A matrix with 5 columns and 4 rows.

TST

first compositional variable

WAKE

second compositional variable

MVPA

third compositional variable

LPA

fourth compositional variable

SB

fifth compositional variable

Description

A list of 4 components

Usage

sim

Format

A list with 5 columns and 4 rows.

brmcoda_tab

Simulation results for brmcoda() for tables

sub_tab

Simulation results for substitution() for tables

brmcoda_plot

Simulation results for brmcoda() for graphs

sub_plot

Simulation results for substitution() for graphs

Description

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

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.

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 by brms::VarCorr().

• group: grouping factor for random-effect rows, NA otherwise.

• parameter: the brms-style parameter label.

• truth: the true generating value from the simulator (sqrt() of covariance diagonals for random_sd rows, correlations for random_cor rows).

• estimate, est_error, lower, upper: posterior summary of the fitted parameter at probs.

• bias: estimate - truth.

• covered: whether truth lies inside ⁠[lower, upper]⁠.

• simulator_name: the simulator-side parameter identifier the row was derived from, for provenance.

See Also

prep_sim_analysis() for the analysis object and its ⁠$truth⁠ table, gen_template() for authoring the truth values.

Examples

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)
}

Description

Build a single-level, grouped, or longitudinal simulation design and evaluate generator specifications into one shared data table.

Usage

simulate_data(
  generators,
  n = NULL,
  n_groups = NULL,
  n_per_group = NULL,
  group_id = "group_id",
  obs_id = "obs_id",
  time_id = NULL,
  time_values = NULL,
  time_truncate = TRUE,
  seed = NULL
)

Arguments

Details

Generators are evaluated in list order. Each generator can depend on columns produced by earlier generators through the simulation context. Generated column names must be unique and must not collide after make.names().

Value

An object of class mlsim_data, a list with:

data

A data.table::data.table() containing design columns and generated variables.

metadata

Simulation design metadata, seed, index metadata, and generator order.

generator_specs

Generator specifications with simulator closures removed.

generator_metadata

Per-generator metadata recorded during simulation.

For multilevel predictor generators, data also contains the group-level random-effect draws expanded to one value per row, under a stable .mlsim_ naming scheme:

⁠.mlsim_<var>_random_intercept⁠

The location random intercept of a continuous variable (e.g. from gen_mvn()).

⁠.mlsim_<var>_scale_random_intercept⁠

The scale random intercept of a continuous variable, present only when the generator was given a joint location-scale random_cov.

⁠.mlsim_<var>_random_intercept_<category>⁠

The location random intercept of a gen_categorical() variable, one column per non-reference category.

Both ⁠<var>⁠ and ⁠<category>⁠ are sanitized with make.names(). These are truth columns, not observed data: they record the generating group-level effects so the data-generating process can be validated, and they are never used as predictors by prep_sim_analysis().

Examples

sim <- simulate_data(
  n_groups = 3,
  n_per_group = 2,
  seed = 10,
  generators = list(
    group_x = gen_mvn("group_x", level = "level2", fixed_intercept = 0, residual_cov = 1),
    y = gen_poisson("y", fixed_intercept = log(2))
  )
)
sim$data
summary(sim)

longitudinal <- simulate_data(
  n_groups = 2,
  n_per_group = 3,
  time_id = "visit",
  generators = list(x = gen_mvn("x", fixed_intercept = 0, residual_cov = 1))
)
longitudinal$metadata$index

Description

This function is an alias of substitution to estimates the difference in an outcome when compositional parts are substituted for specific unit(s) using a aggregate reference composition (e.g., compositional mean at sample level, not seperated by between- and within effects). It is recommended that users run substitution model using the substitution function.

Usage

sub(
  object,
  delta,
  ref = "grandmean",
  level = "aggregate",
  summary = TRUE,
  aorg = TRUE,
  at = NULL,
  parts = 1,
  base,
  type = "one-to-one",
  weight = "equal",
  scale = c("response", "linear"),
  cores = NULL,
  ...
)