Package 'bsitar'

An object of class bgmfit representing the model to which the fit criteria will be added.

Names of model fit criteria to compute. Currently supported are "loo", "waic", "kfold", "loo_subsample", "bayes_R2" (Bayesian R-squared), "loo_R2" (LOO-adjusted R-squared), and "marglik" (log marginal likelihood).

A positive integer indicating the number of posterior draws to use in estimation. If NULL (default), all draws are used.

An integer specifying the specific posterior draw(s) to use in estimation (default NULL).

A flag indicating if the information criteria of the models should be compared to each other via loo_compare.

A flag indicating whether to compute the full log-likelihood matrix at once or separately for each observation. The latter approach is usually considerably slower but requires much less working memory. Accordingly, if one runs into memory issues, pointwise = TRUE is the way to go.

Optional name of the model. If NULL (the default) the name is taken from the call to x.

Logical; Indicates if already stored fit indices should be overwritten. Defaults to FALSE. Setting it to TRUE is useful for example when changing additional arguments of an already stored criterion.

Logical; only relevant if file is specified and ignored otherwise. If TRUE, the fitted model object will be saved regardless of whether new criteria were added via add_criterion.

Either NULL or a character string. In the latter case, the fitted model object including the newly added criterion values is saved via saveRDS in a file named after the string supplied in file. The .rds extension is added automatically. If x was already stored in a file before, the file name will be reused automatically (with a message) unless overwritten by file. In any case, file only applies if new criteria were actually added via add_criterion or if force_save was set to TRUE.

A logical value indicating whether only the estimate should be computed (TRUE), or whether the estimate along with SE and CI should be returned (FALSE, default). Setting summary to FALSE will increase computation time. Note that summary = FALSE is required to obtain correct estimates when re_formula = NULL.

A logical value to specify the summary options. If FALSE (default), the mean is used as the measure of central tendency and the standard deviation as the measure of variability. If TRUE, the median and median absolute deviation (MAD) are applied instead. Ignored if summary is FALSE.

The percentiles to be computed by the quantile function. Only used if summary is TRUE.

An optional data frame for estimation. If NULL (default), newdata is retrieved from the model.

A character string (default NULL) to specify the response variable when processing posterior draws for univariate_by and multivariate models. See bsitar() for details on univariate_by and multivariate models.

An integer specifying the number of cores for parallel execution.

• If NULL (default), the number of cores is set to future::availableCores() - 1.

• On non-Windows systems, this can be controlled globally via the mc.cores option.

A logical value specifying whether to estimate the velocity curve from the derivative function or by differentiating the distance curve. Set model_deriv = TRUE for functions that require the velocity curve, such as growthparameters() and plot_curves(). Set it to NULL for functions that use the distance curve (i.e., fitted values), such as loo_validation() and plot_ppc().

A logical (default TRUE) to indicate whether to return the model object or just assign the criteria to the model. Mainly for internal use.

A logical (default FALSE) to indicate whether to return the criterion. Mainly for internal use.

A logical argument (default FALSE) to specify whether to print information collected during the setup of the object(s).

A logical value or NULL (default FALSE). Controls whether Stan functions are exposed for post-processing.

• TRUE: Explicitly exposes Stan functions saved from the model fit. This is required when calculating fit criteria or bayes_R2 during model optimization.

• FALSE (default): Stan functions are not exposed.

• NULL: The setting is inherited from the original bsitar() model fit configuration.

Note: In the optimize_model() function, the default is NULL (inheriting behavior), whereas other post-processing functions default to FALSE.

A logical value (default NULL) indicating whether to use already exposed and saved Stan functions. This is typically set automatically based on the expose_functions argument from the bsitar() call. Manual specification of usesavedfuns is rarely needed and is intended for internal testing, as improper use can lead to unreliable estimates.

A logical value indicating whether to clear the exposed Stan functions from the environment (TRUE) or not (FALSE). If NULL, clearenvfuns is set based on the value of usesavedfuns: TRUE if usesavedfuns = TRUE, or FALSE if usesavedfuns = FALSE.

The environment used for function evaluation. The default is NULL, which sets the environment to parent.frame(). Since most post-processing functions rely on brms, it is recommended to set envir = globalenv() or envir = .GlobalEnv, especially for derivatives like velocity curves.

Further arguments passed on to the functions from the brms

Predictor variable (typically age in years). For a univariate model, x is a single variable. For univariate_by and multivariate models, x can either be the same for all sub-models, or different for each sub-model. For example, when fitting a bivariate model, x = list(x1, x2) specifies that x1 is the predictor variable for the first sub-model, and x2 is the predictor for the second sub-model. To use x1 as a common predictor variable for both sub-models, user can specify x = list(x1) or simply x = x1.

Response variable (e.g., repeated height measurements). For univariate and univariate_by models, y is specified as a single variable. In the case of a univariate_by model, the response vector for each sub-model is created and named internally based on the factor levels of the variable used to define the sub-model. For example, specifying univariate_by = sex creates response vectors Female and Male when Female is the first level and Male is the second level of the sex variable. In a multivariate model, the response variables are provided as a list, such as y = list(y1, y2), where y1 is the response variable for the first sub-model, and y2 is the response for the second sub-model. Note that for the multivariate model, the data are not stacked; instead, response vectors are separate variables in the data and must be of equal length.

A factor variable uniquely identifying the groups (e.g., individuals) in the data frame. For univariate_by and multivariate models, the id can be the same (typically) for all sub-models, or different for each sub-model (see the x argument for details on setting different arguments for sub-models).

A data frame containing variables such as x, y, id, etc.

Degrees of freedom for the natural cubic spline design matrix (default 4). The df is used internally to construct knots (quantiles of the x distribution) for the spline design matrix. For univariate_by and multivariate models, the df can be the same across sub-models (e.g., df = 4) or different for each sub-model, such as df = list(4, 5), where df = 4 applies to the first sub-model and df = 5 applies to the second sub-model.

A numeric vector specifying the knots for the natural cubic spline design matrix (default NULL). Note that user cannot specify both df and knots at the same time, nor can both be NULL. In other words, either df or knots must be specified. Like df, the knots can be the same for all sub-models, or different for each sub-model when fitting univariate_by and multivariate models (see df for details).

A named list to control knot optimization during model fitting. This feature explores a larger candidate set (typically df + 4) and selects an optimal subset that minimizes the chosen information criterion. The following elements are available:

select

Selection strategy for knots and/or degrees of freedom. Options are "knots", "df", or "both". "knots" (fully integrated into the bsitar workflow) automatically optimizes knot locations for a given df and passes the result to internal functions. The "both" option first determines the optimal degree of freedom, then selects knot locations using this df, which may result in a different (optimized) df from the user's initial choice. Users should note this potential change. Alternatively, select = "df" can be combined with return = TRUE to recover the optimal degree of freedom, which may then be reused as an argument in bsitar().

all_scores

Logical (default FALSE). When TRUE and select = "df", returns the information criterion value for each evaluated df. Ignored for other selection types.

plot_all_scores

Logical (default FALSE). When TRUE and select = "df" with all_scores = TRUE, plots the criterion for each df. Ignored otherwise.

nsearch

Number of candidate knots in the search space. If NULL (default), set automatically to df + 4 to provide sufficient coverage for optimization.

criteria

Model selection criterion. Options include "AIC", "BIC", and cross-validation ("CV"). Lower values are preferred. Defaults to "AIC". Note: Cross-validation often requires considerably more computation than "AIC" or "BIC".

cvk

Number of folds for cross-validation (default 10). Only used if criteria = "CV".

cviter

Number of cross-validation iterations to average (default 100). Only used if criteria = "CV".

bkrange

Logical (default TRUE). If TRUE, uses the full predictor range for boundary knots. Otherwise, boundary knots are set via quantiles (see Hmisc::rcspline.eval()). This option is relevant only for stype = "rcs" and overrides any later bkrange value globally for knots_selection.

fix_bknots

Logical (default TRUE). If TRUE, fixes boundary knots during internal optimization. If FALSE, boundary knots are also candidates for optimization. Applies only when stype = "rcs"; see notes on bkrange.

kspace

Knot spacing approach:

• "un": Uniform quantile-based spacing (default)

• "nu": Non-uniform spacing (algorithm-selected). Note: the actual number of knots may differ from the specified df.

method

Algorithm for knot optimization:

• "bs": built-in systematic workflow for knot optimization (default)

• "rs": Recursive search based on user-defined function. Not implemented yet.

when

Stage of knot optimization relative to predictor centering:

• "bc": Before centering (raw x values; default)

• "ac": After centering (i.e., x - xoffset)

what

Diagnostics or plotting output type. Default "knots". Available options:

• "knots": Optimized internal and boundary knots

• "df": Selected degree of freedom

• "plot1": Data and curve (original knots)

• "plot2": Data and curve (optimized knots)

• "plot3": Data and curves for both knot sets

• "plot4": Data, curve, and confidence bands (original knots)

• "plot5": Data, curve, and confidence bands (optimized knots)

• "plot6": Data, curve, and confidence bands (both knots)

• "plot7": Residual plot (original knots)

• "plot8": Residual plot (optimized knots)

• "plot9": Residual plots (comparison)

return

Logical. If TRUE, returns the diagnostic plot object (as specified by what). Default is FALSE.

print

Logical. If TRUE, displays the diagnostic plot (as specified by what). Default is FALSE.

Note that when both return and print are FALSE, the optimized knots are passed to the bsitar() function for model fitting. If return = TRUE, the knots are returned immediately, and model fitting is not performed. When return = FALSE and print = TRUE, the selected knots are displayed via a plot, and the model fitting proceeds using bsitar().

Example usage: bsitar( x = age, y = height, id = id, data = berkeley_exdata, knots_selection = list( select = "knots", nsearch = NULL, criteria = "AIC", bkrange = TRUE, fix_bknots = TRUE, method = "bs", when = "bc", what = "plot3", return = FALSE, print = TRUE ), seed = 123)

A character string specifying the fixed effects structure (default 'a+b+c'). For univariate_by and multivariate models, user can specify different fixed effect structures for each sub-model. For example, fixed = list('a+b+c', 'a+b') implies that the fixed effects structure for the first sub-model is 'a+b+c', and for the second sub-model it is 'a+b'.

A character string specifying the random effects structure (default 'a+b+c'). The approach to setting the random is the same as for the fixed effects structure (see fixed).

An optional character string or numeric value to set the origin of the predictor variable, x (i.e., centering of x). Available options include:

• 'mean': The mean of x (i.e., mean(x)),

• 'max': The maximum value of x (i.e., max(x)),

• 'min': The minimum value of x (i.e., min(x)),

• 'apv': Age at peak velocity estimated from the velocity curve derived from a simple linear model fit to the data,

• Any real number (e.g., xoffset = 12). Note that the value should be on the original scale of the x variable because this number will be automatically transformed based on the xfunxoffset which is set same as xfun unless xfunxoffset is explicitly set as FALSE. See xfunxoffset for details. The default is xoffset = 'mean'.

For univariate_by and multivariate models, xoffset can be the same or different for each sub-model (see x for details on setting different arguments for sub-models). If xoffset is a numeric value, it will be transformed internally (e.g., log or sqrt) depending on the xfun argument. Similarly, when xoffset is 'mean', 'min', or 'max', these values are calculated after applying the log or sqrt transformation to x.

An optional character string or numeric value to set the origin of the fixed effect parameter b. The bstart argument is used to establish the location parameter for location-scale based priors (such as normal()) via the b_prior_beta argument, and/or the initial value via the b_init_beta argument. The available options for bstart are:

• 'mean': The mean of x (i.e., mean(x)),

• 'min': The minimum value of x (i.e., min(x)),

• 'max': The maximum value of x (i.e., max(x)),

• 'apv': Age at peak velocity estimated from the velocity curve derived from a simple linear model fit to the data,

• Any real number (e.g., bstart = 12).

The default is bstart = 'xoffset' (i.e., the same value as xoffset). For univariate_by and multivariate models, bstart can be the same for all sub-models (typically), or different for each sub-model (refer to x for details on setting different arguments for sub-models).

An optional character string or numeric value to set the origin of the fixed effect parameter c. The cstart argument is used to establish the location parameter for location-scale based priors (such as normal()) via the c_prior_beta argument, and/or the initial value via the c_init_beta argument. The available options for cstart are:

• 'pv': Peak velocity estimated from the velocity curve derived from the simple linear model fit to the data,

• Any real number (e.g., cstart = 1).

Note that since parameter c is estimated on the exponential scale, the cstart should be adjusted accordingly. The default cstart is '0' (i.e., cstart = '0'). For univariate_by and multivariate models, cstart can be the same for all sub-models (typically), or different for each sub-model (refer to x for details on setting different arguments for sub-models).

An optional character string specifying the transformation of the predictor variable x. The default value is NULL, indicating that no transformation is applied and the model is fit to the data with the original scale of x. The available transformation options are:

• 'log': Logarithmic transformation,

• 'sqrt': Square root transformation.

For univariate_by and multivariate models, the xfun can be the same for all sub-models (typically), or different for each sub-model (refer to x for details on setting different arguments for sub-models).

An optional character string specifying the transformation of the response variable y. The default value is NULL, indicating that no transformation is applied and the model is fit to the data with the original scale of y. The available transformation options are:

• 'log': Logarithmic transformation,

• 'sqrt': Square root transformation.

For univariate_by and multivariate models, the yfun can be the same for all sub-models (typically), or different for each sub-model (refer to x for details on setting different arguments for sub-models).

Transformation applied to xoffset for x variable (default TRUE). See xfun for details. The default xfunxoffset = TRUE sets its value to same as xfun. Users rarely need to specify xfunxoffset themselves. One potential application is when a user wants to turn it off by setting xfunxoffset = FALSE. Note that xfunxoffset is called only when xoffset is a numeric values and not data based such as xoffset = mean. This is because even when xfunxoffset is not TRUE, xoffset is still automatically adjusted by xfun, as it is inferred from the transformed x variable.

An optional real number specifying the value by which the span of the predictor variable x should be extended (default is 0.04). This extension can help in modeling edge cases. For more details, refer to the sitar package documentation. For univariate_by and multivariate models, the bound can be the same for all sub-models (typically), or different for each sub-model (refer to x for details on setting different arguments for sub-models).

A character string or a named list specifying the spline type to be used. The available options are:

• 'rcs' (default): Constructs the spline design matrix using the truncated power basis (Harrell's method), as implemented in Hmisc::rcspline.eval().

• 'nsk': Implements a B-spline based natural cubic spline method, similar to splines2::nsk().

• 'nsp': Implements a B-spline based natural cubic spline method, similar to splines2::nsp().

The 'rcs' method uses a truncated power basis, whereas 'nsk' and 'nsp' are B-spline-based methods. Unlike splines2::nsp() and splines2::nsk(), which normalize the spline basis by default, 'nsk' and 'nsp' return the non-normalized version of the spline. If normalization is desired, the user can specify normalize = TRUE in a list. For example, to use a normalized 'nsp', one can specify stype = list(type = 'nsp', normalize = TRUE).

For more details, see Hmisc::rcspline.eval(), splines2::nsk(), and splines2::nsp().

An optional character string (default NULL) specifying terms on the right-hand side of the response variable, but before the formula tilde sign ~. The terms_rhs is used when fitting a measurement error model.

For example, when fitting a model with measurement error in the response variable, the formula in brms::brmsformula() could be specified as brmsformula(y | mi(sdy) ~ ...). In this case, mi(sdy) is passed to the formula via terms_rhs = 'mi(sdy)'.

For a multivariate model, each outcome can have its own measurement error variable. For instance, the terms_rhs can be specified as a list: terms_rhs = list(mi(sdy1), mi(sdy2)).

Note that brms::brmsformula() does not allow combining mi() with the subset() formulation which used in fitting univariate_by models.

Formula for the fixed effect parameter, a (default ~ 1). Users can specify different formulas when fitting univariate_by and multivariate models.

For example, a_formula = list(~1, ~1 + cov) specifies that the a_formula for the first sub-model includes only an intercept, while the second sub-model includes both an intercept and a covariate cov. The covariate(s) can be either continuous or factor variables. For factor covariates, dummy variables are created internally using stats::model.matrix().

The formula can include a combination of continuous and factor variables, as well as their interactions.

Formula for the fixed effect parameter, b (default ~ 1). See a_formula for details on how to specify the formula. The behavior and structure of b_formula are similar to a_formula.

Formula for the fixed effect parameter, c (default ~ 1). See a_formula for details on how to specify the formula. The behavior and structure of c_formula are similar to a_formula.

Formula for the fixed effect parameter, d (default ~ 1). See a_formula for details on how to specify the formula. The behavior and structure of d_formula are similar to a_formula.

Formula for the fixed effect parameter, s (default ~ 1). The s_formula sets up the spline design matrix. Typically, covariates are not included in the s_formula to limit the population curve to a single curve for the entire data. In fact, the sitar package does not provide an option to include covariates in the s_formula. However, the bsitar package allows the inclusion of covariates. In such cases, the user must justify the modeling of separate curves for each category when the covariate is a factor variable.

Formula for the random effect parameter a (default ~ 1). As with a_formula, users may specify different formulas when fitting univariate_by and multivariate models. The formula can include continuous and/or factor covariates, as well as their interactions (see a_formula for details).

In addition to defining the design matrix for the random effect parameter a, this argument allows specification of the grouping structure and correlation among random effects using vertical bar (|) notation.

For example, to include only random intercepts for parameters a, b, and c, one can specify:

a_formula_gr = ~1, b_formula_gr = ~1, c_formula_gr = ~1.

To include a grouping variable (e.g., id) and specify an unstructured correlation among the random effects, use:

a_formula_gr = ~ (1|i|id)
b_formula_gr = ~ (1|i|id)
c_formula_gr = ~ (1|i|id)

Here, i is a shared identifier indicating that the random effects for a, b, and c are modeled as jointly correlated with an unstructured covariance matrix. See [brms::brm()] for further details on this notation.

Alternatively, the grouping structure and correlation can be specified using the group_by argument. The equivalent specification to the above can be written as:

a_formula_gr = ~1, b_formula_gr = ~1, c_formula_gr = ~1,

together with group_by = list(groupvar = id, cor = un), where id defines the grouping variable and un specifies an unstructured correlation. See group_by for additional details.

Formula for the random effect parameter b (default ~ 1). See a_formula_gr for details.

Formula for the random effect parameter c (default ~ 1). See a_formula_gr for details.

Formula for the random effect parameter d (default ~ 1). See a_formula_gr for details.

Formula for the random effect parameter a (default NULL), used for specifying hierarchical models with three or more levels. For example, this may arise in data with repeated measurements (level 1) nested within individuals (level 2), which are further nested within studies (level 3).

For a_formula_gr_str, only the vertical bar (|) notation (see a_formula_gr) can be used to define grouping structures and correlation patterns. An example specification for a three-level model with random effect parameters a, b, and c is:

a_formula_gr_str = ~ (1|i|id:study) + (1|i2|study)
b_formula_gr_str = ~ (1|i|id:study) + (1|i2|study)
c_formula_gr_str = ~ (1|i|id:study) + (1|i2|study)

In this setup, |i| and |i2| define unstructured correlation among random effects at the individual and study levels, respectively. These identifiers must be distinct, as random effects are not correlated across different hierarchical levels.

More generally, users can specify models with additional hierarchical levels and include covariates within the random effect formulas as needed.

Formula for the random effect parameter b (default NULL), used for specifying hierarchical models with three or more levels. See a_formula_gr_str for details.

Formula for the random effect parameter c (default NULL), used for specifying hierarchical models with three or more levels. See a_formula_gr_str for details.

Formula for the random effect parameter d (default NULL), used for specifying hierarchical models with three or more levels. See a_formula_gr_str for details.

Logical indicator for adjusting the scale of the predictor variable x when including the random effect parameter d. The coefficient for d is modeled as a linear function of x, i.e., d * x.

If FALSE (default), the original x is used. If TRUE, x is adjusted for the timing (b) and intensity (c) parameters as (x - b) * exp(c), so that the model includes d * ((x - b) * exp(c)).

This adjusted scale represents individual developmental age rather than chronological age, making the parameter d more sensitive to differences in growth timing (e.g., pubertal timing) across individuals. See [sitar::sitar()] for further details.

Formula for the distributional (fixed effect) parameter sigma. This argument defines the fixed-effect design matrix for sigma, allowing inclusion of continuous and/or factor covariates, as well as their interactions (see a_formula for details). The parameter sigma is modeled on the log scale.

By default, sigma_formula = NULL, in which case [brms::brm()] estimates sigma as a residual standard deviation (RSD) parameter on the link scale. When specified, sigma_formula, together with sigma_formula_gr and sigma_formula_gr_str, enables modeling sigma with both fixed and random effects, analogous to other model parameters such as a, b, and c.

An alternative approach is to use the dpar_formula argument to specify the design matrix for sigma. The advantage of dpar_formula is its flexibility in supporting both linear and nonlinear specifications via [brms::lf()] and [brms::nlf()]. This allows additional control, such as centering predictors or removing the intercept using a 0 + formulation. However, a limitation of dpar_formula is that it does not support random effects for sigma.

The arguments sigma_formula and dpar_formula are mutually exclusive. When either is specified, the default estimation of sigma as an RSD parameter in [brms::brm()] is disabled.

External functions (e.g., poly, splines::ns) can be used within sigma_formula. Two common approaches are:

• Direct use of standard functions, e.g., sigma_formula = ~ 1 + splines::ns(age, df = 3).

• Use of a custom wrapper function with a single argument (the predictor), e.g., sigma_formula = mypoly(age). Additional arguments can be set by modifying the function definition, for example: mypoly <- poly; formals(mypoly)[['degree']] <- 3; mypoly(age). This approach can produce more concise coefficient names.

Formula for the random effect parameter sigma (default NULL). See a_formula_gr for details on specifying grouping structures and correlation using vertical bar notation. As with sigma_formula, external functions (e.g., poly, splines::ns) may be included. See sigma_formula for further details.

Formula for the random effect parameter sigma, used when fitting hierarchical models with three or more levels. See a_formula_gr_str for details on specifying multi-level structures.

As with sigma_formula, external functions can be incorporated. For example: sigma_formula_gr_str = (1 + splines::ns(age, df = 3) | i | id:study).

See sigma_formula for additional details.

A custom formula for advanced modeling of the distributional parameter sigma using the brms::nlf() and brms::lf() functions. This is an advanced and experimental feature for specifying complex variance structures. The sigma_formula_manual can be can be a character string or a list of formulas.

Use Cases: The primary use cases for sigma_formula_manual are:

1. Basic variance modelling: A simple formula for modelling heteroscedasticity (sigma).

2. Advanced variance modelling: Calculating variance weights where the variance of residuals depends on a specified co variate, replicating functionality from the nlme package.

3. Location-Scale Models: Applying a full SITAR-like model to the scale (sigma) parameter in addition to the location (mu).

Basic variance modelling: A simple example of modeling sigma directly:

  sigma_formula_manual = list(nlf(sigma ~ z) + lf(z ~ 1 + (1 | gr(id))))
  

Note that use can specify any external function in the linear predictor part of the formula 'lf()' in order to model nonlinear trajectories. For example, user can include splines2::nsk() in the model as follows:

  sigma_formula_manual = list(nlf(sigma ~ z) + 
  lf(z ~ 1 + splines2::nsk(age) + (1 | gr(id))))
  

Automatic Prior Assignment:
Priors for these linear predictors are assigned automatically for both mean and group-level random effects. The function uses priors specified via arguments like 'sigma_prior_beta', 'sigma_cov_prior_beta', 'sigma_prior_sd', etc. These are the same arguments otherwise used for setting priors on parameters defined by 'sigma_formula', 'sigma_formula_gr', and 'sigma_formula_gr_str'.

Advanced variance modelling: This approach models heteroscedasticity using an explicit variance function. The 'bsitar' package provides six different methods for variance modeling, five of which are implemented in the 'nlme' package. The variance modeling approach uses a helper function 'sigmavarfun' (short hand form, 'vf') that sets up the appropriate formulation for variance modeling within the 'bsitar' package.

For multivariate models, the 'sigmavarfun' is appropriately renamed to indicate response specific 'sigmavarfun' by adding '_response' as prefix. The documentation below provides a detailed overview of the available methods and their usage.

Available Methods (short hands in parentheses):

• 'varpower' ('vp'): Implements nlme::varPower().

• 'varconstpower' ('cp'): Implements nlme::varConstPower().

• 'varexp' ('ve'): Implements nlme::varExp() with 'form ~ x' as follows:
  'nlme::varExp(form ~ x)'

• 'fitted' ('fi'): Implements nlme::varExp() with 'form ~ fitted(.)' as follows:
  'nlme::varExp(form ~ fitted(.))'

• 'residual' ('re'): Implements nlme::varExp() with 'form ~ resid(.)' as follows:
  'nlme::varExp(form ~ resid(.))'

• 'mean' ('me'): A sixth method based on an example from the brms reference manual, which models the square root of the fitted values as follows:
  'nlme::varExp(form ~ sqrt(fitted(.)))'

Below are examples demonstrating how to use 'sigmavarfun' via the method sub-argument to specify each of the six variance models. To minimize the risk of typographical errors, we recommend using the shorthand form 'vf' instead of the full name 'sigmavarfun'.

All variance model specifications are provided through the sigma_formula_manual argument in bsitar(). In general, the syntax takes the form:

sigma_formula_manual = nlf(sigma ~ vf(...), method = ...) + lf(...)

Note that the predictor used within these variance functions should not be identical to x (see argument x, typically age). This is because the predictors for mu and sigma may require different transformations, which cannot be accommodated using a single shared variable.

Instead, define a separate predictor variable (even if it contains identical values). For example:

data$age_vf <- data$age

and then use age_vf within the variance function. This ensures that transformations applied to sigma do not interfere with those used for mu.

For method = 're', outcome in the nlf(sigma ~ vf()) is a placeholder that is internally replaced by the actual response variable such as height (See y for the response variables).

1. varpower:

  nlf(sigma ~ vf(param1, param2, predictor), method = 'vp') +
  lf(param1 + param2 ~ 1)

2. varConstPower:

  nlf(sigma ~ vf(param1, param2, param3, predictor), method = 'cp') + 
  lf(param1 + param2 + param3 ~ 1)

3. varExp:

  nlf(sigma ~ vf(param1, param2, predictor), method = 've') +
  lf(param1 + param2 ~ 1)

4. fitted:

  nlf(sigma ~ vf(param1, param2, identity()), method = 'fi') +
  lf(param1 + param2 ~ 1)

5. residual:

  nlf(sigma ~ vf(param1, param2, identity(), outcome), method = 're') + 
  lf(param1 + param2 ~ 1)

6. mean:

  nlf(sigma ~ vf(param1, param2, identity()), method = 'me') +
  lf(param1 + param2 ~ 1)

Internal Predictor Transformations: The function applies internal transformations based on the chosen method:

• For 'varpower' and 'varConstPower', the predictor is transformed to log(abs(predictor)).

• For 'varexp', the predictor is not transformed.

• For 'fitted' and 'residual', identity() is internally set to fitted(.).

• For 'mean', identity() is internally set to sqrt(fitted(.)).

Note that the default 'fitted' (see above, 4. fitted:) internally gets coded as method = 'fittedexp' (or short hand method = 'fe') which sets up the variance form similar to the varExp. To set up the variance form similar to the varpower, please use method = 'fittedpower' (or short hand method = 'fp'). Another form, which models the non zero values for the 'fitted', can be specified as method = 'fittedz' (or short hand method = 'fz').

Similar to the fitted (please above), the default method = 'residual' will set up the variance form similar to the varExp by internally coding the method argument as method = 'residualexp' (i.e., method = 're'). The power form can be set as method = 'residualpower' (or, method = 'rp').

Like fitted and residual (please above), the default method = 'mean' will set up the variance form similar to the varExp by internally coding the method argument as method = 'meanexp' (i.e., method = 'me'). The power form can be set as method = 'meanpower' (or, method = 'mp').

Covariates and Random Effects: The linear predictor, 'lf()', can be extended to include covariates and group-level random effects. For example, 'lf(param1 + param2 ~ 1)' can be expanded as follows:

  lf(param1 + param2 ~ 1 + covariate + (1 || gr(id, by = groupid)))

Automatic Prior Assignment: Priors for these linear predictors are assigned automatically for both mean and group-level random effects. The function uses priors specified via arguments like 'sigma_prior_beta', 'sigma_cov_prior_beta', 'sigma_prior_sd', etc. These are the same arguments otherwise used for setting priors on parameters defined by 'sigma_formula', 'sigma_formula_gr', and 'sigma_formula_gr_str'.

To disable this automatic prior assignment, user can add the argument prior = 'self' to the nlf() function. For example:

  nlf(..., method = 'vp', prior = 'self')

Note: If user disables the automatic prior assignment, it is advised to first get the required prior structure by calling the bsitar() with the 'get_prior = TRUE' argument. The relevant portions of this structure can then be edited and added back to the model via the 'add_self_priors' argument in bsitar().

For multivariate model, the sigma_formula_manual can be used set up different form and predictor variables for each outcome separately by using the list approach as shown below:
sigma_formula_manual = list( nlf(sigma ~ vf(param1, param2, identity()), method = 'fi') + lf(param1 + param2 ~ 1) , nlf(sigma ~ vf(param1, param2, identity()), method = 'fi') + lf(param1 + param2 ~ 1) ).
Note that formulation name 'nlf()' should be same across all outcomes. The appropriate response assignment is done internally. Also, parameter should not contain dots or underscores.

Location-Scale Models: Another important use case for sigma_formula_manual is modeling sigma in a location scale SITAR model, where the SITAR formula can be applied to the scale parameter (sigma) similar to the one used for modelling the location parameter mu. An example is:

nlf(sigma ~ ls(x, sigmaa, sigmab, sigmac, sigmas1, sigmas2, sigmas3, sigmas4), loop = FALSE) + lf(sigmaa ~ 1+(1 |110| gr(id, by = NULL))+(1 |330| gr(study, by = NULL))) + lf(sigmab ~ 1+(1 |110| gr(id, by = NULL))+(1 |330| gr(study, by = NULL))) + lf(sigmac ~ 1+(1 |110| gr(id, by = NULL))+(1 |330| gr(study, by = NULL))) + lf(sigmas1 + sigmas2 + sigmas3 + sigmas4 ~ 1).

Here, ls, which is an abbreviation for location-scale, is a placeholder and gets replaced by the the name of actual function that is used for modelling the scale (sigma) part of the model. For example, if the name of the mu function is SITARFun, then the name for the sigma function would be sigmaSITARFun and hence ls gets replaced as sigmaSITARFun. All functions and corresponding names are created internally.

The first argument x is again a placeholder for the actual predictor variable defined for the sigma modelling via the sigmax argument.In other words, the x gets replaced by the sigmax argument evaluated. For example, when sigmax = age, then the x gets replaced by age. Like ls, internal renaming of x is handled automatically.

The growth parameters sigmaa, sigmab, sigmac should match the input used for the sigmafixed and sigmarandom argument. In other words, either sigmafixed or sigmarandom must include the form for the corresponding growth parameter defined in the ls function. For example, when either or both sigmafixed and sigmarandom are defined as a+b+c, then all three growth parameters sigmaa, sigmab, sigmac should be part of the ls function. However, when only sunset of parameters are specified via the sigmafixed and sigmarandom form, say for example a+b, the ls function should exclude sigmac parameter.

Similarly, the number of spline parameters are based on the sigmadf argument. For example, when sigmadf = 4, then four spline parameters sigmas1, ..., sigmas4 are included, as shown above in the example.

The other relevant information that is passed to the ls function can be set via the sigmaknots, sigmaxoffset, sigmaxfun, and sigmabound arguments.

Note that for the location-scale model, priors must be set up manually using the add_self_priors argument. To see which priors are required, the user can run the code with get_priors = TRUE. Also note that the default initial values for location-scale model are random, or '0' depending on whether the init = 'random' or init = '0'.

Predictor for the distributional parameter sigma. See x for details. Ignored if sigma_formula_manual = NULL.

A factor variable uniquely identifying the groups (e.g., individuals) in the data frame. If NULL (default), then sigmaid is set same as id. For univariate_by and multivariate models, the sigmaid can be the same (typically) for all sub-models, or different for each sub-model (see the id argument for details on setting different arguments for sub-models). Ignored if sigma_formula_manual = NULL.

Degree of freedom for the spline function used for sigma. See df for details. Ignored if sigma_formula_manual = NULL.

Knots for the spline function used for sigma. See knots for details. Ignored if sigma_formula_manual = NULL.

Fixed effect formula for the sigma structure. See fixed for details. Ignored if sigma_formula_manual = NULL.

Random effect formula for the sigma structure. See random for details. Ignored if sigma_formula_manual = NULL. Currently not used even when sigma_formula_manual is specified.

Offset for the x in the sigma structure. See xoffset for details. Ignored if sigma_formula_manual = NULL.

Starting value for the b parameter in the sigma structure. See bstart for details. Ignored if sigma_formula_manual = NULL. Currently not used even when sigma_formula_manual is specified.

Starting value for the c parameter in the sigma structure. See cstart for details. Ignored if sigma_formula_manual = NULL. Currently not used even when sigma_formula_manual is specified.

Transformation of sigmax variable for sigma structure. See xfun for details. Ignored if sigma_formula_manual = NULL.

Bounds for the x in the sigma structure. See bound for details. Ignored if sigma_formula_manual = NULL.

Transformation applied to sigmaxoffset for sigmax variable (default TRUE). See sigmaxfun for details. The default sigmaxfunxoffset = TRUE sets its value to same as sigmaxfun. Users rarely need to specify sigmaxfunxoffset themselves. One potential application is when a user wants to turn it off by setting sigmaxfunxoffset = FALSE.

Note that sigmaxfunxoffset is called only when sigmaxoffset is a numeric values and not data based such as sigmaxoffset = mean. This is because even when sigmaxfunxoffset is not TRUE, sigmaxoffset is still automatically adjusted by sigmaxfun, as it is inferred from the transformed sigmax variable.

Formula for the distributional fixed effect parameter, sigma (default NULL). See sigma_formula for details.

Formula to set up the autocorrelation structure of residuals (default NULL). Allowed autocorrelation structures include:

• auto regressive moving average (arma) of order p and q, specified as autocor_formula = ~arma(p = 1, q = 1).

• auto regressive (ar) of order p, specified as autocor_formula = ~ar(p = 1).

• moving average (ma) of order q, specified as autocor_formula = ~ma(q = 1).

• unstructured (unstr) over time (and individuals), specified as autocor_formula = ~unstr(time, id).

See brms::brm() for further details on modeling the autocorrelation of residuals.

Note that for unstr autocorrelation structure, time is a placeholder which could be any name. Internally, a new variable time, or time_response for multivariate or univariate_by models, is created from the predictor variable x e.g., age.

This internal variable is required because the unstructured autocorrelation covariance matrix is indexed by time points treated as discrete, ordered indices. The underlying Stan implementation expects time to be represented as unique positive integers in increasing order (e.g., 1, 2, 3, ...) for each individual so that:

• Each observation can be mapped to a unique row/column in the N_time x N_time unstructured covariance matrix, where N_time is the number of distinct time points per individual.

• The autocorrelation structure is defined over a regular, ordered sequence of time indices, which is essential for constructing the covariance matrix in Stan.

• Individuals with different original time scales (e.g., age in years, months, or days) can share a common autocorrelation structure, since time is standardized to integer indices within each individual.

The predictor variable x is therefore transformed into this internal time variable such that it contains unique positive integers in increasing order, such as 1, 2, 3, for each individual identified by id. When the original x values do not satisfy this condition (e.g., they are non-integer ages like 0.5, 1.2, 3.7), the transformation automatically transforms them to meet this requirement during the creation of the internal time variable. No pre-processing is needed before using unstr autocorrelation.

Family distribution (default gaussian) and the link function (default identity). See brms::brm() for details on available distributions and link functions, and how to specify them. For univariate_by and multivariate models, the family can be the same for all sub-models (e.g., family = gaussian()) or different for each sub-model, such as family = list(gaussian(), student()), which sets gaussian distribution for the first sub-model and student_t distribution for the second. Note that the family argument is ignored if custom_family is specified (i.e., if custom_family is not NULL).

Specifies custom families (i.e., response distribution). Default is NULL. For details, see brms::custom_family(). Note that user-defined Stan functions must be exposed by setting expose_functions = TRUE.

Specifies custom formula. Default is NULL. For details, see brms::brmsformula(). Currently ignored.

Specifies custom prior Default is NULL. For details, see brms::prior(). Currently ignored. It is primarily designed to support setting custom prior for custom_formula. Note that custom_prior is different from the set_self_priors which is evaluated throughout the call.

Allows the passing of user-defined variables to be added to Stan's program blocks (default NULL). This is primarily useful when defining a custom_family. For more details on specifying stanvars, see brms::custom_family(). Note that custom_stanvars are passed directly without conducting any sanity checks.

Specify the grouping variable for group-level random effects. The group_arg should be a named list that may include groupvar, dist, cor, and by, as described below:

• groupvar specifies the grouping (subject) identifier. If NULL (default), it is automatically inferred from the id argument.

• dist specifies the distribution of the random effects (default gaussian). Currently, only the Gaussian distribution is supported (see [brms::brm()]).

• by allows estimation of separate variance-covariance structures (i.e., standard deviations and correlations) across levels of a factor variable (default NULL). The by variable must be a factor. For example, by = "sex" estimates separate variance-covariance structures for males and females.

• cor specifies the covariance (correlation) structure of the random effects. The default is unstructured (cor = un) for all model types (univariate, univariate_by, and multivariate).

  For univariate and univariate_by models, the alternative option is diagonal, which estimates only variance parameters while fixing correlations to zero.

  For multivariate models, available options are:

  • un: A fully unstructured covariance, where random effects across outcomes are jointly modeled with shared correlations.

  • diagonal: Variances are estimated separately for each outcome, with correlations fixed to zero.

  • un_s: Separate unstructured variance-covariance matrices are estimated for each outcome.

All components of group_arg are optional and default values are used when unspecified.

For univariate_by and multivariate models, only groupvar from group_arg is used, as additional model-specific options are handled via the corresponding arguments.

The group_arg argument is ignored when random effects are specified using vertical bar notation (see a_formula_gr) or when fitting hierarchical models with three or more levels (see a_formula_gr_str).

Specify arguments for modeling distributional-level random effects for sigma. The setup for sigma_group_arg follows the same approach as described for group-level random effects (see group_arg for details).

Set up the univariate-by-subgroup model fitting (default NULL) via a named list with the following elements:

• by (optional, character string): Specifies the factor variable used to define the sub-models (default NA).

• cor (optional, character string): Defines the correlation structure. Options include un (default) for a full unstructured variance-covariance structure and diagonal for a structure with only variance parameters (i.e., standard deviations) and no covariance (i.e., correlations set to zero).

• terms (optional, character string): Specifies the method for setting up the sub-models. Options are 'subset' (default) and 'weights'. See brms::`addition-terms` for more details.

Set up the multivariate model fitting (default NULL) using a named list with the following elements:

• mvar (logical, default FALSE): Indicates whether to fit a multivariate model.

• cor (optional, character string): Specifies the correlation structure for group-level random effects. Available options are:

  • "un" (default): Models a full unstructured correlation, where group-level random effects across response variables are drawn from a joint multivariate normal distribution with shared correlation parameters.

  • "diagonal": Estimates only the variance parameters for each sub-model, with the correlation parameters set to zero.

  • "un_s": Estimates unstructured variance-covariance parameters separately for each response variable.

• rescor (logical, default TRUE): Indicates whether to estimate the residual correlation between response variables.

• rcorr_by (character string, default NULL): The variable by which separate residual correlations between response variables are estimated. This must be a factor variable present in the data. Ignored when rescor = FALSE.

• rcorr_gr (character string, default NULL): The grouping variable (typically a level 2 variable like id) for which residual correlations between response variables are estimated. This must be a factor variable present in the data. Ignored when rescor = FALSE.

• rcorr_method (character string, default NULL): Specifies the method for estimating residual correlations: "lkj" (uses an LKJ distribution) or "cde" (uses the Cholesky decomposition of the covariance matrix with a uniform prior of -1, 1 for each off-diagonal element). If NULL, "lkj" is set as the default. Ignored when rescor = FALSE.

• rcorr_prior (numeric vector, default NULL): Used to specify the LKJ prior for each level of rcorr_by when rcorr_method = "lkj". The length of rcorr_prior should be either one (to use the same prior for all levels) or equal to the number of levels in rcorr_by. For rcorr_method = 'cde', uniform priors -1,1 are assigned to all correlation parameters. The argument rcorr_prior is ignored when rescor = FALSE.

Prior specification for the fixed effect parameter a (default: normal(ymean, ysd, autoscale = FALSE)). A robust alternative to the default setting is to use lm rather than ymean. This is because even when the predictor xvar is mean-centered, the intercept parameter a does not always correspond exactly to the overall sample mean. The @details section provides full details on prior specifications for the bsitar package.

Prior specification for the fixed effect parameter b (default: normal(0, 2.0, autoscale = FALSE)). The @details section provides full details on prior specifications for the bsitar package.

Prior specification for the fixed effect parameter c (default: normal(0, 1.0, autoscale = FALSE)). The @details section provides full details on prior specifications for the bsitar package.

Prior specification for the fixed effect parameter d (default: normal(0, 1.0, autoscale = FALSE)). For full details on prior specifications, please refer to a_prior_beta. Note: To set the scale for location-scale priors, users can specify the standard deviation (xsd) or the median absolute deviation (xmad) of the predictor variable x.

Prior specification for the fixed effect parameter s (spline coefficients) (default: normal(lm, lm, autoscale = FALSE)). For general guidelines on prior specifications, please refer to a_prior_beta. Key points specific to s include:

• Data-driven lm priors: When using lm for location-scale priors (e.g., normal(lm, lm)), the location parameter is derived from spline coefficients obtained via a simple linear model fit. The scale parameter is derived from the standard deviation of the spline design matrix.

• Hierarchical priors (sethp): The sethp argument (default: FALSE) enables hierarchical prior specifications. Setting sethp = TRUE modifies the prior from s ~ normal(0, lm) to s ~ normal(0, 'hp'), where hp ~ normal(0, lm). The scale for the hierarchical prior is inherited from s, but user can also explicitly pass a distribution. For example, normal(0, lm, sethp = cauchy) translates to s ~ normal(0, lm) and hp ~ cauchy(0, lm).

• Uniform distributions: The addrange argument can be used to symmetrically expand the prior range.

Scale Parameter Variants for lm When the scale parameter is set to lm, it defaults to lm1. The available variants calculate the scale differently (where y is the outcome, X is the design matrix, and lm_se is the standard error vector from the linear model fit):

• lm or lm1: sd(y) / sd(X)

• lm2: lm_se / sd(X)

• lm3: lm_se

Note: Location-scale priors (such as normal, student_t, and cauchy) generally perform well for spline coefficients. In most cases, the default normal(lm, lm) is highly effective and can be further refined by setting autoscale = TRUE depending on the specific data structure.

Prior specification for covariate effects included in the fixed effect parameter a (default: normal(0, 5.0, autoscale = FALSE)). The general approach to prior specification is similar to a_prior_beta, with the following differences:

• Restricted data-driven options: The options 'ymean', 'ymedian', 'ysd', and 'ymad' are not available for a_cov_prior_beta.

• lm for location: The lm option can be used for the location parameter, allowing covariate coefficients to be derived from a simple linear model fitted to the data. This option is available only for a_cov_prior_beta and not for covariates in other fixed or random effect parameters.

• Submodel-specific priors: Separate priors can be specified for individual sub models when fitting univariate_by or multivariate models; see a_prior_beta for details.

Prior specification for covariate effects included in the fixed effect parameter b (default: normal(0, 1.0, autoscale = FALSE)). See a_cov_prior_beta for details.

Prior specification for covariate effects included in the fixed effect parameter c (default: normal(0, 0.1, autoscale = FALSE)). See a_cov_prior_beta for details.

Prior specification for covariate effects included in the fixed effect parameter d (default: normal(0, 1.0, autoscale = FALSE)). See a_cov_prior_beta for details.

Prior specification for covariate effects included in the fixed effect parameter splines (default: normal(0, 10.0, autoscale = FALSE)). As described in s_formula, the SITAR model does not allow covariates in the spline design matrix. If covariates are specified (see s_formula), the approach to setting priors for the covariates in parameter s is the same as for a (see a_cov_prior_beta). For location-scale based priors, the option lm sets the location parameter based on spline coefficients obtained from fitting a simple linear model to the data.

Prior specification for the random effect parameter, a. (default: normal(0, 'ysd', autoscale = FALSE)). The prior is applied to the standard deviation (the square root of the variance), not the variance itself. The approach for setting the prior is similar to a_prior_beta, with the location parameter always set to zero. The lower bound is automatically set to 0 by brms::brm(). For univariate_by and multivariate models, priors can be the same or different for each sub model (see a_prior_beta).

Prior specification for the random effect parameter, b. (default: normal(0, 2.0, autoscale = FALSE)). See a_prior_sd for details.

Prior specification for the random effect parameter, c. (default: normal(0, 1.0, autoscale = FALSE)). See a_prior_sd for details.

Prior specification for the random effect parameter, d. (default: normal(0, 1.0, autoscale = FALSE)). See a_prior_sd for details.

Prior specification for covariate effects included in the random effect parameter a (default: normal(0, 5.0, autoscale = FALSE)). The approach is the same as described for a_cov_prior_beta, except that no pre-defined options (e.g., lm) are allowed.

Prior specification for covariate effects included in the random effect parameter b (default: normal(0, 1.0, autoscale = FALSE)). See a_cov_prior_sd for details.

Prior specification for covariate effects included in the random effect parameter c (default: normal(0, 0.1, autoscale = FALSE)). See a_cov_prior_sd for details.

Prior specification for covariate effects included in the random effect parameter d (default: normal(0, 1.0, autoscale = FALSE)). See a_cov_prior_sd for details.

Prior specification for the random effect parameter a when fitting a hierarchical model with three or more levels of hierarchy. (default NULL). See a_prior_sd for details on setting priors when needed.

Prior specification for the random effect parameter b when fitting a hierarchical model with three or more levels of hierarchy. (default NULL). See b_prior_sd for details on setting priors when needed.

Prior specification for the random effect parameter c when fitting a hierarchical model with three or more levels of hierarchy. (default NULL). See c_prior_sd for details on setting priors when needed.

Prior specification for the random effect parameter d when fitting a hierarchical model with three or more levels of hierarchy (default NULL). See a_prior_sd for details on setting priors when needed.

Prior specification for covariate effects in the random effect parameter a when fitting a hierarchical model with three or more levels of hierarchy (default NULL). The approach follows a_cov_prior_sd.

Prior specification for covariate effects in the random effect parameter b when fitting a hierarchical model with three or more levels of hierarchy (default NULL). The approach follows a_cov_prior_sd_str.

Prior specification for covariate effects in the random effect parameter c when fitting a hierarchical model with three or more levels of hierarchy (default NULL). The approach follows a_cov_prior_sd_str.

Prior specification for covariate effects in the random effect parameter d when fitting a hierarchical model with three or more levels of hierarchy (default NULL). The approach follows a_cov_prior_sd_str.

Prior specification for the fixed effect distributional parameter sigma (default normal(0, 1.0, autoscale = FALSE)). The general approach mirrors a_prior_beta.

Prior specification for covariate effects in the fixed effect distributional parameter sigma (default normal(0, 0.5, autoscale = FALSE)). The general approach mirrors a_cov_prior_beta.

Prior specification for the random effect distributional parameter sigma (default normal(0, 0.25, autoscale = FALSE)). The general approach mirrors a_prior_sd.

Prior specification for covariate effects in the random effect distributional parameter sigma (default normal(0, 0.15, autoscale = FALSE)). The general approach mirrors a_cov_prior_sd.

Prior specification for the random effect distributional parameter sigma when fitting a hierarchical model with three or more levels of hierarchy (default NULL). The general approach mirrors a_prior_sd_str.

Prior specification for covariate effects in the random effect distributional parameter sigma when fitting a hierarchical model with three or more levels of hierarchy (default NULL). The general approach mirrors a_cov_prior_sd_str.

Prior specification for the residual standard deviation parameter sigma (default normal(0, 'ysd', autoscale = FALSE)). Evaluated when both dpar_formula and sigma_formula are NULL. For location-scale distributions, the scale parameter can be set to the standard deviation (ysd) or median absolute deviation (ymad) of the outcome, or to the residual standard deviation from a linear mixed model ('lme_rsd') or linear model ('lm_rsd'). If nlme::lme() fails to converge, 'lm_rsd' is used automatically.

Prior specification for the fixed effect distributional parameter sigma when using dpar_formula (default normal(0, 'ysd', autoscale = FALSE)). Evaluated when sigma_formula is NULL. See rsd_prior_sigma for details.

Prior specification for covariate effects in the fixed effect distributional parameter sigma when using dpar_formula (default normal(0, 1.0, autoscale = FALSE)). Evaluated when sigma_formula is NULL.

Prior specification for autocorrelation parameters when fitting models with "arma", "ar", or "ma" autocorrelation structures (see autocor_formula). The only allowed distribution is uniform, bounded between -1 and 1 (default uniform(-1, 1, autoscale = FALSE)). For unstructured residual correlation structures, use autocor_prior_unstr_acor.

Prior specification for autocorrelation parameters when fitting models with the unstructured ("un") autocorrelation structure (see autocor_formula). The only allowed distribution is lkj (default lkj(1)). See gr_prior_cor for details on the lkj prior.

Prior specification for correlation parameters of group-level random effects (default lkj(1)). The only allowed distribution is lkj, specified via a single parameter eta.

Prior specification for correlation parameters of group-level random effects when fitting a hierarchical model with three or more levels of hierarchy (default lkj(1)). The general approach mirrors gr_prior_cor.

Prior specification for correlation parameters of distributional random effects sigma (default lkj(1)). The only allowed distribution is lkj (see gr_prior_cor). Note that brms::brm() does not currently allow different lkj priors for group-level and distributional random effects sharing the same group identifier.

Prior specification for correlation parameters of distributional random effects sigma when fitting a hierarchical model with three or more levels of hierarchy (default lkj(1)). The general approach mirrors sigma_prior_cor.

Prior specification for the residual correlation matrix in multivariate models (default lkj(1)). The only allowed distribution is lkj (see gr_prior_cor for details).

Initial values for the sampler. Options include:

• "random" (default): Stan randomly generates initial values for each parameter within a range defined by init_r (see below), or between -2 and 2 on the unconstrained scale if init_r = NULL.

• "0": All parameters are initialized to zero.

• "prior": Parameters are initialized from their specified priors.

• NULL: Initial values are taken from the parameter-specific init arguments defined below.

When init = NULL, fixed effects and variance–covariance parameters (with vcov_init_0 = FALSE) are initialized according to their individual settings. To force all parameters to use random initial values, set init = "random", which is internally translated to init = NULL with fully random initialization.

Positive real value specifying the range for random initial values (default 0.5). This argument is used only when init = "random". Note that the default in Stan is 2.0, which generates random initials between -2.0 and 2.0 on the unconstrained parameter space.

Initial specification for the fixed effect parameter a (default "random"). Available options include:

• "0": Initializes the parameter to zero.

• "random": Initializes with random values within a specified range.

• "prior": Uses values drawn from the prior distribution.

• "ymean": Initializes with the mean of the response variable.

• "ymedian": Initializes with the median of the response variable.

• "lm": Initializes with coefficients from a simple linear model fitted to the data.

Note that "ymean", "ymedian", and "lm" are available only for the fixed effect parameter a. For univariate_by and multivariate models, initial values can be the same across sub-models (e.g., a_init_beta = "0") or different for each sub-model (e.g., a_init_beta = list("0", "lm")).

Initial specification for the fixed effect parameter b (default "random"). See a_init_beta for details on available options.

Initial specification for the fixed effect parameter c (default "random"). See a_init_beta for details on available options.

Initial specification for the fixed effect parameter d (default "random"). See a_init_beta for details on available options.

Initial specification for the fixed effect parameter s (default "random"). Available options include:

• "0": Initializes the parameter to zero.

• "random": Initializes with random values within a specified range.

• "prior": Uses values drawn from the prior distribution.

• "lm": Initializes with coefficients from a simple linear model fitted to the data.

Initial specification for covariate effects in the fixed effect parameter a (default "random"). Available options include:

• "0": Initializes covariate effects to zero.

• "random": Initializes with random values within a specified range.

• "prior": Uses values drawn from the prior distribution.

• "lm": Initializes with coefficients from a simple linear model fitted to the data.

Note that "lm" is available only for a_cov_init_beta and not for covariates in other fixed effect parameters such as b, c, or d.

Initial specification for covariate effects in the fixed effect parameter b (default "random"). See a_cov_init_beta for details.

Initial specification for covariate effects in the fixed effect parameter c (default "random"). See a_cov_init_beta for details.

Initial specification for covariate effects in the fixed effect parameter d (default "random"). See a_cov_init_beta for details.

Initial specification for covariate effects in the fixed effect parameter s (default "lm"). See a_cov_init_beta for details. The option "lm" corresponds to spline coefficients obtained from a simple linear model fitted to the data. However, s_cov_init_beta acts only as a placeholder and is not evaluated, because covariates are not allowed for s; see s_formula for details.

Initial specification for the standard deviation of the group-level random effect parameter a (default "random"). Available options are:

• "random": Initializes with random values within a specified range.

• "prior": Uses values drawn from the prior distribution.

• "ysd": Uses the standard deviation (sd) of the response variable as the initial value.

• "ymad": Uses the median absolute deviation (mad) of the response variable as the initial value.

• "lme_sd_a": Uses the standard deviation of the random intercept from a linear mixed model (nlme::lme()) fitted to the data. If nlme::lme() fails to converge, "lm_sd_a" is used automatically.

• "lm_sd_a": Uses the square root of the residual variance from a simple linear model fitted to the data.

Note that "ysd", "ymad", "lme_sd_a", and "lm_sd_a" are available only for the random effect parameter a and not for other group-level random effects.

When fitting univariate_by and multivariate models, initial values can be shared across all sub-models or specified separately for each sub-model.

Initial specification for the standard deviation of the group-level random effect parameter b (default "random"). See a_init_sd for options and details.

Initial specification for the standard deviation of the group-level random effect parameter c (default "random"). See a_init_sd for options and details.

Initial specification for the standard deviation of the group-level random effect parameter d (default "random"). See a_init_sd for options and details.

Initial specification for covariate effects in the group-level random effect parameter a (default "random"). Available options include:

• "random": Random initialization.

• "prior": Uses values drawn from the prior distribution.

Initial specification for covariate effects in the group-level random effect parameter b (default "random"). See a_cov_init_sd for options and details.

Initial specification for covariate effects in the group-level random effect parameter c (default "random"). See a_cov_init_sd for options and details.

Initial specification for covariate effects in the group-level random effect parameter d (default "random"). See a_cov_init_sd for options and details.

Initial specification for the fixed effect distributional parameter sigma (default "random"). Available options include:

• "random": Random initialization.

• "prior": Uses values drawn from the prior distribution.

Initial specification for covariate effects in the fixed effect distributional parameter sigma (default "random"). See sigma_init_beta for options and details.

Initial specification for the standard deviation of the distributional random effect parameter sigma (default "random"). The approach is the same as for group-level random effects such as a (see a_init_sd for details).

Initial specification for covariate effects in the distributional random effect parameter sigma (default "random"). The approach mirrors a_cov_init_sd.

Initial specification for the correlation parameters of group-level random effects (default "random"). Allowed options are:

• "random": Random initialization.

• "prior": Uses values drawn from the prior distribution.

• Numeric vector: A vector of length equal to the number of distinct correlations (lower-triangular elements). For example, for three random effect parameters, one may use gr_init_cor = list(c(0.5, 0.5, 0.5)).

Note that when vcov_init_0 = TRUE, gr_init_cor is set to "0".

Initial specification for the correlation parameters of distributional random effects sigma (default "random"). Allowed options are:

• "random": Random initialization.

• "prior": Uses values drawn from the prior distribution.

Initial specification for the residual standard deviation parameter sigma (default "random"). Available options are:

• "0": Initializes the residual standard deviation to zero.

• "random": Random initialization of the residual standard deviation.

• "prior": Uses values drawn from the prior distribution.

• "lme_rsd": Uses the residual standard deviation from a linear mixed model (nlme::lme()) fitted to the data.

• "lm_rsd": Uses the square root of the residual variance from a simple linear model fitted to the data.

If nlme::lme() fails to converge, "lm_rsd" is used automatically. rsd_init_sigma is evaluated only when both dpar_formula and sigma_formula are NULL.

Initial specification for the distributional parameter sigma when using dpar_formula (default "random"). The approach and options mirror rsd_init_sigma. This argument is evaluated only when dpar_formula is not NULL.

Initial specification for covariate effects in the distributional parameter sigma (default "random"). Allowed options are "0", "random", and "prior".

Initial specification for autocorrelation parameters (see autocor_formula for details; default "random"). Allowed options are "0", "random", and "prior".

Initial specification for unstructured residual autocorrelation parameters (default "random"). Allowed options are "0", "random", and "prior". The approach mirrors gr_init_cor.

Initial specification for the residual correlation matrix in multivariate models (default "random"). Allowed options are "0", "random", and "prior".

Initial specification for standardized group-level random effect parameters (default "random"). These parameters arise from the non-centered parameterization used internally by brms::brm().

Logical indicator controlling initial values for variance (standard deviation) and covariance (correlation) parameters. When vcov_init_0 = TRUE, all variance–covariance parameters are initialized to zero, allowing custom initial specifications for fixed effects while keeping variance–covariance parameters at zero. When vcov_init_0 = FALSE (default), variance–covariance parameters are assigned random initial values unless a parameter-specific initial specification is provided (e.g., a_init_sd = "0"). This setting is ignored when global initial values are assigned via the init argument.

Named list or numeric value used to add a small amount of noise to initial values for population-level (fixed effect) parameters.

When specified as a numeric value, jitter_init_beta is interpreted as a percentage and must lie between 0 and 100. Internally, this percentage is converted to a proportion (percentage / 100) and passed as the amount argument to base::jitter(), with factor fixed at its default 1.

The default jitter_init_beta = NULL means no perturbation is applied, so identical initial values are used across chains. For mild perturbation, a value such as jitter_init_beta = 10 is often adequate.

Jitter is applied proportionally to the specified initial value, not as an absolute amount. For example, if the initial value is 100, setting jitter_init_beta = 10 yields perturbed values in the range 90–110, whereas an initial value of 10 yields values in the range 9–11.

When jitter_init_beta is provided as a named list, its elements are passed directly to base::jitter(). In addition to factor and amount, a percent argument may be supplied and is treated in the same way as a single numeric value. If both percent and factor are supplied, the effective perturbation is their product.

To use the defaults of base::jitter(), supply an empty list list(), which is then completed internally with factor = 1, amount = NULL. See base::jitter() for further details on factor and amount.

Named list or numeric value used to add a small amount of noise to initial values for the standard deviations of random-effect parameters. A pragmatic choice that has worked well in early testing is jitter_init_sd = 1, corresponding to a 1% perturbation. See jitter_init_beta for details on the available options and how perturbations are applied.

Named list or numeric value used to add a small amount of noise to initial values for the correlations of random-effect parameters. A pragmatic choice that has worked well in early testing is jitter_init_cor = 0.1, corresponding to a 0.1 (10%) perturbation. See jitter_init_beta for details on available options and how perturbations are applied.

Optional named list (default NULL) used to pass external objects into prior specifications for individual parameters (e.g., a_prior_beta). This is particularly useful when supplying long vectors or matrices as priors, which can be created in the R environment and then referenced via prior_data. For example, to pass vectors of location and scale parameters for 10 dummy-variable coefficients in the fixed effect parameter a, one can:

• Create named objects in R: prior_a_cov_location <- rnorm(10, mean = 0, sd = 1) prior_a_cov_scale <- rep(5, 10)

• Add them to prior_data: prior_data = list( prior_a_cov_location = prior_a_cov_location, prior_a_cov_scale = prior_a_cov_scale )

• Use these objects in the prior specification: a_cov_prior_beta = normal(prior_a_cov_location, prior_a_cov_scale)

Optional named list (default NULL) used to pass external objects into initialization arguments. The handling of init_data is identical to prior_data (see above).

Custom initialization object supplied as a named list. The list is passed directly to the underlying init argument without checking dimensions or name matching. If some parameters already have parameter-specific initial specifications (e.g., a_init_beta = "0"), init_custom is applied only to parameters without their own initial settings. To override this behavior and use init_custom for all parameters regardless of parameter-specific initials, set init = "custom".

Logical (default FALSE) indicating whether to expose Stan functions used in model fitting.

Logical (default FALSE) indicating whether to return the Stan code (see brms::stancode() for details).

Logical (default FALSE) indicating whether to return the Stan data (see brms::standata() for details).

Logical (default FALSE) indicating whether to return the model formula (see brms::brmsformula() for details).

Logical (default FALSE) indicating whether to return the Stan variables (see brms::stanvar() for details).

Logical or character (default FALSE) controlling prior retrieval (see brms::get_prior() for details). If get_priors = TRUE, priors corresponding to the final compiled model code are returned. To obtain basic default priors instead, use get_priors = "default".

Logical (default FALSE) indicating whether to return the priors as specified by the user.

Logical (default FALSE) indicating whether to return the initial values as specified by the user.

Logical (default FALSE) indicating whether to validate specified priors (see brms::validate_prior() for details).

Optional object (default NULL) used to manually specify priors. set_self_priors is passed directly to brms::brm() without additional checks.

Optional object (default NULL) used to append elements to an existing prior object. Intended for internal use only.

Optional object (default NULL) used to replace elements of an existing prior object. Intended for internal use only.

Optional object (default NULL) used to enforce the same prior structure across hierarchical levels. Intended for internal use only.

Optional specification (default NULL) for removing outliers. This should be a named list passed directly to sitar::velout() and sitar::zapvelout(). Intended for internal use only.

Optional formula defining variables that are unused in the model but should still be stored in the model frame, for example when they are required for post-processing.

Number of Markov chains (default 4).

Total number of iterations per chain, including warmup (default 2000).

Positive integer specifying the number of warmup (burn-in) iterations. Warmup iterations are also used for step-size adaptation and should not be used for inference. The number of warmup iterations must not exceed iter, and the default is iter / 2.

Positive integer specifying the thinning interval. Set thin > 1 to reduce memory usage and computation time when iter is large. Thinning is sometimes used when MCMC draws show high autocorrelation, evidenced by poor chain mixing (e.g., high rhat) despite good parameter recovery. Autocorrelation can be inspected using functions such as mcmc_acf() from the bayesplot package.

Number of cores used to execute chains in parallel. See brms::brm() for further details. Unlike brms::brm(), which defaults cores to getOption("mc.cores", 1), bsitar defaults cores to getOption("mc.cores", "optimize"), which attempts to optimize system resource usage. In this mode, the maximum number of cores used is the number of available cores minus one. If the available cores exceed the number of chains (see chains), cores is set equal to the number of chains.

Alternatively, cores may be set to getOption("mc.cores", "maximise"), which uses all available cores regardless of the number of chains, or specified as in brms::brm() via getOption("mc.cores", 1). This behavior can be set globally using options(mc.cores = x), where x is one of "optimize", "maximise", or 1. The cores argument can also be set directly as an integer (e.g., cores = 4).

Character string specifying the interface used to execute the Stan model. Available options are "rstan" (default) and "cmdstanr". The backend can also be set globally for the current R session via the "brms.backend" option. See brms::brm() for more details.

Number of threads used for within-chain parallelization. Unlike brms::brm(), which defaults threads to getOption("brms.threads", NULL) (no within-chain parallelization), bsitar defaults threads to getOption("brms.threads", "optimize") to better utilize modern multi-core systems. In this mode, the number of threads per chain is set to the number of available cores minus one.

An alternative is getOption("brms.threads", "maximise"), which sets the number of threads per chain equal to the total number of available cores. Users can also emulate the brms default via getOption("brms.threads", NULL). These behaviors can be set globally using options(brms.threads = x), where x is one of "optimize", "maximise", or NULL.

Alternatively, threads can be specified directly using threads = threading(x), where x is an integer. Additional arguments, such as grainsize and static, can also be passed via threading(). See brms::brm() for further details on within-chain parallelization.

Numeric vector of length two specifying the platform and device IDs for OpenCL-based GPU support during model fitting. If unsure, c(0, 0) usually corresponds to the default device. See brms::opencl() for guidance on identifying valid platform and device IDs. This parameter can also be set globally via the "brms.opencl" option.

Logical indicating whether normalization constants should be included in the Stan code (default TRUE). If FALSE, normalization constants are omitted, which may improve sampling efficiency, but requires Stan version 2.25 or higher. Note that normalize = FALSE disables some post-processing functionality (e.g., brms::bridge_sampler()). This behavior can be controlled globally via the "brms.normalize" option.

Character string specifying the estimation method. Available options are:

• "sampling" (default): MCMC sampling.

• "meanfield": Variational inference with independent normal approximations.

• "fullrank": Variational inference with a full-rank multivariate normal approximation.

• "fixed_param": Sampling with parameters fixed (no adaptation).

This argument can be set globally via the "brms.algorithm" option (see options).

Named list of sampler control parameters. Defaults are as in brms::brm(), except that max_treedepth is increased from 10 to 12 to better explore challenging posterior geometries common in nonlinear models. In contrast, adapt_delta retains the default value of 0.8 to avoid unnecessarily long sampling times. See brms::brm() for a complete description of available control parameters.

Logical. If TRUE, the Stan model is not created and compiled and the corresponding 'fit' slot of the brmsfit object will be empty. This is useful if you have estimated a brms-created Stan model outside of brms and want to feed it back into the package.

For internal use only.

Named list of arguments passed to the "pathfinder" algorithm, used to obtain pathfinder-based initial values for subsequent MCMC sampling. Currently, pathfinder_args is supported only when backend = "cmdstanr". If pathfinder_args is not NULL and backend = "rstan" is requested, the backend is automatically switched to "cmdstanr".

Logical (default FALSE) indicating whether to use initial values obtained from the "pathfinder" algorithm when fitting the final model via MCMC sampling. pathfinder_args is currently supported only when backend = "cmdstanr". If pathfinder_args is not NULL and backend = "rstan" is requested, the backend is automatically switched to "cmdstanr". The arguments passed to the "pathfinder" algorithm are supplied via pathfinder_args; if pathfinder_args is NULL, the default settings from cmdstanr are used.

A named list of objects containing data, which cannot be passed via argument data. Required for some objects used in autocorrelation structures to specify dependency structures as well as for within-group covariance matrices.

A data.frame object (default NULL). This is mainly for internal testing and not be used for routine model fitting.

A logical (default NULL) indicating whether to generate xyadj in the generated quantities block of the Stan.

A character string indicating whether to draw samples from the priors in addition to the posterior draws. Options are "no" (the default), "yes", and "only". These prior draws can be used for various purposes, such as calculating Bayes factors for point hypotheses via brms::hypothesis(). Note that improper priors (including the default improper priors used by brm) are not sampled. For proper priors, see brms::set_prior(). Also, prior draws for the overall intercept are not obtained by default for technical reasons. See brms::brmsformula() for instructions on obtaining prior draws for the intercept. If sample_prior is set to "only", draws will be taken solely from the priors, ignoring the likelihood, which allows user to generate draws from the prior predictive distribution. In this case, all parameters must have proper priors.

An object generated by brms::save_pars() that controls which parameters should be saved in the model. This argument does not affect the model fitting process itself but provides control over which parameters are retained in the final output.

A logical value indicating whether unused factor levels in the data should be dropped. The default is TRUE.

A list of additional arguments passed to rstan::stan_model when using the backend = "rstan" or backend = "cmdstanr". This allows customization of how models are compiled.

An integer specifying the frequency of printing every nth iteration. By default, NULL indicates that the refresh rate will be automatically set by brms::brm(). Setting refresh is especially useful when thin is greater than 1, in which case the refresh rate is recalculated as (refresh * thin) / thin.

A verbosity level between 0 and 2. When set to 1 (the default), most informational messages from the compiler and sampler are suppressed. Setting it to 2 suppresses even more messages. The sampling progress is still printed. To turn off all printing, set refresh = 0. Additionally, when using backend = "rstan", The user can prevent the opening of additional progress bars by setting open_progress = FALSE.

An integer or NA (default) specifying the seed for random number generation, ensuring reproducibility of results. If set to NA, Stan will randomly select the seed.

A character string or NULL (default). If provided, the Stan code for the model will be saved in a text file with the name corresponding to the string specified in save_model.

An instance of class brmsfit from a previous fit (default is NA). If a brmsfit object is provided, the compiled model associated with the fitted result is reused, and any arguments that modify the model code or data are ignored. It is generally recommended to use the update method for this purpose, rather than directly passing the fit argument.

Either NULL or a character string. If a character string is provided, the fitted model object is saved using saveRDS in a file named after the string supplied in file. The .rds extension is automatically added. If the specified file already exists, the existing model object is loaded and returned instead of refitting the model. To overwrite an existing file, user must manually remove the file or specify the file_refit argument. The file name is stored within the brmsfit object for later use.

Logical or a character string, specifying one of the compression algorithms supported by saveRDS. If the file argument is provided, this compression will be used when saving the fitted model object.

Modifies when the fit stored via the file argument is re-used. This can be set globally for the current R session via the "brms.file_refit" option (see options). The possible options are:

• "never" (default): The fit is always loaded if it exists, and fitting is skipped.

• "always": The model is always refitted, regardless of existing fits.

• "on_change": The model is refitted only if the model, data, algorithm, priors, sample_prior, stanvars, covariance structure, or similar parameters have changed.

If user believe a false positive occurred, the [brms::brmsfit_needs_refit()] can be used to investigate why a refit is deemed necessary. A refit will not be triggered for changes in additional parameters of the fit (e.g., initial values, number of iterations, control arguments). A known limitation is that a refit will be triggered if within-chain parallelization is switched on/off.

Logical; If TRUE, the future package is used for parallel execution of the chains. In this case, the cores argument will be ignored. The execution type is controlled via plan (see the examples section below). This argument can be set globally for the current R session via the "future" option.

Currently ignored. Placeholder for future development.

Currently ignored. Placeholder for future development.

A character string specifying the type of parameterization to use for drawing group-level random effects. Options are: 'ncp' for Non-Centered Parameterization (NCP), and 'cp' for Centered Parameterization (CP).

The NCP is generally recommended when the likelihood is weak (e.g., few observations per individual) and is the default approach (and only option) in brms::brm().

The CP parameterization is typically more efficient when a relatively large number of observations are available across individuals. We consider a 'relatively large number' as at least 10 repeated measurements per individual. If there are fewer than 10, NCP is used automatically. This behavior applies only when parameterization = NULL. To explicitly set CP parameterization, use parameterization = 'cp'.

Note that since brms::brm() does not support CP, the stancode generated by brms::brm() is edited internally before fitting the model using rstan::rstan() or "cmdstanr", depending on the chosen backend. Therefore, CP parameterization is considered experimental and may fail if the structure of the generated stancode changes in future versions of brms::brm().

An optional argument (logical, default FALSE) to indicate whether to print information collected during setting up the model formula priors, and initials. As an example, the user might be interested in knowing the response variables created for the sub model when fitting a univariate-by-subgroup model. This information can then be used in setting the desired order of options passed to each such model such as df, prior, initials etc.

Further arguments passed to brms::brm(). This may include additional arguments that are either passed directly to the underlying model fitting function, or used for internal purposes. Specifically, the ... can also be used to pass arguments used for testing and debugging, such as: match_sitar_a_form, sigmamatch_sitar_a_form, displayit, setcolh, setcolb, decomp and smat. Note that for all arguments (such as intercept) to correctly pass to the respective functions, use smat and not stype when setting up the spline arguments.

These internal arguments are typically not used in regular model fitting but can be relevant for certain testing scenarios or advanced customization. Users are generally not expected to interact with these unless working on debugging or testing specific features of the model fitting process.

An R object of class bsitar or a list of models.

More brmsfit objects.

The name of the criterion to be extracted from brmsfit objects.

If NULL (the default) will use model names derived from deparsing the call. Otherwise will use the passed values as model names.

A logical (default FALSE) indicating whether to check and add fit criterion if not already added to model objects

A named list to set up the arguments for the add_model_criterion() function. Ignored when check_criterion = FALSE or when each model objects has precomputed fit criteria.

A logical value or NULL (default FALSE). Controls whether Stan functions are exposed for post-processing.

• TRUE: Explicitly exposes Stan functions saved from the model fit. This is required when calculating fit criteria or bayes_R2 during model optimization.

• FALSE (default): Stan functions are not exposed.

• NULL: The setting is inherited from the original bsitar() model fit configuration.

Note: In the optimize_model() function, the default is NULL (inheriting behavior), whereas other post-processing functions default to FALSE.

A logical (default FALSE) to print some useful information.

An object of class bgmfit.

A character string containing the user-defined Stan function(s) in Stan code. If NULL (the default), the scode will be retrieved from the model.

A logical (default TRUE) to indicate whether to expose the functions and add them as an attribute to the model.

A character string (default NULL) to specify the model name. This parameter is for internal use only.

A logical (default TRUE) to specify whether to return the model object. If expose = TRUE, it is advisable to set returnobj = TRUE.

A logical (default FALSE) to indicate whether the exposed functions should be vectorized using base::Vectorize(). Note that currently, vectorize should be set to FALSE, as setting it to TRUE may not work as expected.

A logical argument (default FALSE) to specify whether to print information collected during the setup of the object(s).

A logical (default FALSE) to indicate whether to return the sigma functions. This parameter is for internal use only. Ignored

A character string (default NULL) to set up the the backend method ('rstan' or 'cmdstanr') for compiling and exposing stan functions. If NULL, then the backend is same as the backend method used for model fitting. Note that when backend = FALSE, then default backend will be set as 'rstan'. This is particularly useful because backend = 'cmdstanr' will fail on WSL system.

A character string to set up the path of installed CmdStan. If NULL (default)

The environment used for function evaluation. The default is NULL, which sets the environment to parent.frame(). Since most post-processing functions rely on brms, it is recommended to set envir = globalenv() or envir = .GlobalEnv, especially for derivatives like velocity curves.

Additional arguments passed to the rstan::expose_stan_functions() function. The "..." can be used to set the compiler, which can be either rstan::stanc() or rstan::stan_model(). You can also pass other compiler-specific arguments such as save_dso for rstan::stan_model(). Note that while both rstan::stanc() and rstan::stan_model() can be used as compilers before calling rstan::expose_stan_functions(), it is important to note that the execution time for rstan::stan_model() is approximately twice as long as rstan::stanc().

An object of class bgmfit.

An optional data frame for estimation. If NULL (default), newdata is retrieved from the model.

A character string (default NULL) to specify the response variable when processing posterior draws for univariate_by and multivariate models. See bsitar() for details on univariate_by and multivariate models.

Optional name of a predicted distributional parameter. If specified, expected predictions of this parameters are returned.

A positive integer indicating the number of posterior draws to use in estimation. If NULL (default), all draws are used.

An integer specifying the specific posterior draw(s) to use in estimation (default NULL).

Option to indicate whether or not to include individual/group-level effects in the estimation. When NA (default), individual-level effects are excluded, and population average growth parameters are computed. When NULL, individual-level effects are included in the computation, and the resulting growth parameters are individual-specific. In both cases (NA or NULL), continuous and factor covariates are appropriately included in the estimation. Continuous covariates are set to their means by default (see numeric_cov_at for details), while factor covariates remain unaltered, allowing for the estimation of covariate-specific population average and individual-specific growth parameters.

A flag indicating if new levels of group-level effects are allowed (defaults to FALSE). Only relevant if newdata is provided.

Indicates how to sample new levels for grouping factors specified in re_formula. This argument is only relevant if newdata is provided and allow_new_levels is set to TRUE. If "uncertainty" (default), each posterior sample for a new level is drawn from the posterior draws of a randomly chosen existing level. Each posterior sample for a new level may be drawn from a different existing level such that the resulting set of new posterior draws represents the variation across existing levels. If "gaussian", sample new levels from the (multivariate) normal distribution implied by the group-level standard deviations and correlations. This options may be useful for conducting Bayesian power analysis or predicting new levels in situations where relatively few levels where observed in the old_data. If "old_levels", directly sample new levels from the existing levels, where a new level is assigned all of the posterior draws of the same (randomly chosen) existing level.

A flag indicating if correlation structures originally specified via autocor should be included in the predictions. Defaults to TRUE.

An optional (named list) argument to specify the value of continuous covariate(s). The default NULL option sets the continuous covariate(s) to their mean. Alternatively, a named list can be supplied to manually set these values. For example, numeric_cov_at = list(xx = 2) will set the continuous covariate variable 'xx' to 2. The argument numeric_cov_at is ignored when no continuous covariates are included in the model.

An optional argument to specify the ids for the hierarchical model (default NULL). It is used only when the model is applied to data with three or more levels of hierarchy. For a two-level model, levels_id is automatically inferred from the model fit. For models with three or more levels, levels_id is inferred from the model fit under the assumption that hierarchy is specified from the lowest to the uppermost level, i.e., id followed by study, where id is nested within study. However, it is not guaranteed that levels_id is sorted correctly, so it is better to set it manually when fitting a model with three or more levels of hierarchy.

An optional argument (default NULL) to calculate (marginal/average) curves and growth parameters, such as APGV and PGV. If specified, it must be a named list indicating the over (typically a level 1 predictor, such as age), feby (fixed effects, typically a factor variable), and reby (typically NULL, indicating that parameters are integrated over the random effects). For example, avg_reffects = list(feby = 'study', reby = NULL, over = 'age').

An optional argument to specify the variable(s) that can be passed to the ipts argument (see below). This is useful when fitting location-scale models and measurement error models. If post-processing functions throw an error such as variable 'x' not found in either 'data' or 'data2', consider using aux_variables.

An optional argument to specify the variable(s) that can be passed to the marginaleffects::datagrid(). This is useful when fitting location-scale models and measurement error models. If post-processing functions throw an error such as variable 'x' not found in either 'data' or 'data2', consider using grid_add. Note that unlike aux_variables which are passed to the internal data functions such as 'get.newdata', the grid_add are passed to the marginaleffects::datagrid().

An integer specifying the number of points for interpolating the predictor variable (e.g., age) to generate smooth curves for predictions and plots. This value is used as the length.out argument for seq(), controlling the smoothness of distance and velocity curves without altering the predictor range.

NULL (the default)

Engages automatic behavior based on the dpar argument: it is internally set to 50 for the mean response (dpar = 'mu'), ensuring smooth curves, while for the distributional parameters (e.g., dpar = 'sigma'), no interpolation is performed.

An integer (e.g., 100)

Explicitly sets the number of interpolation points.

FALSE

Disables all interpolation, forcing predictions to be made only at the original data points of the predictor variable.

This argument affects the following post-processing functions:
fitted_draws(), predict_draws(), growthparameters(), plot_curves(), get_predictions(), get_comparisons(), and get_growthparameters().

An integer indicating whether to estimate the distance curve or its derivative (velocity curve). The default deriv = 0 is for the distance curve, while deriv = 1 is for the velocity curve.

A logical value specifying whether to estimate the velocity curve from the derivative function or by differentiating the distance curve. Set model_deriv = TRUE for functions that require the velocity curve, such as growthparameters() and plot_curves(). Set it to NULL for functions that use the distance curve (i.e., fitted values), such as loo_validation() and plot_ppc().

A logical value indicating whether only the estimate should be computed (TRUE), or whether the estimate along with SE and CI should be returned (FALSE, default). Setting summary to FALSE will increase computation time. Note that summary = FALSE is required to obtain correct estimates when re_formula = NULL.

A logical value to specify the summary options. If FALSE (default), the mean is used as the measure of central tendency and the standard deviation as the measure of variability. If TRUE, the median and median absolute deviation (MAD) are applied instead. Ignored if summary is FALSE.

A function applied to individual draws from the posterior distribution before computing summaries (default NULL). The argument transform_draws is derived from the marginaleffects::predictions() function and should not be confused with the transform argument from the deprecated brms::posterior_predict() function. It's important to note that for both marginaleffects::predictions() and marginaleffects::avg_predictions(), the transform_draws argument takes precedence over the transform argument. Note that when transform_draws = NULL, an attempt is made to automatically set transform_draws = 'exp' for dpar = 'sigma'. User can set transform_draws = FALSE to turn off this automatic assignment of 'exp' to the transform_draws. It is also important to set transform_draws = FALSE when computing the first derivative (velocity) for dpar = 'sigma'.

Either "response" or "linear". If "response", results are returned on the scale of the response variable. If "linear", results are returned on the scale of the linear predictor term, that is without applying the inverse link function or other transformations.

The percentiles to be computed by the quantile function. Only used if summary is TRUE.

An integer to set the predictor range (e.g., age) when executing the interpolation via ipts. By default, NULL sets the individual-specific predictor range. Setting xrange = 1 applies the same range for individuals within the same higher grouping variable (e.g., study). Setting xrange = 2 applies an identical range across the entire sample. Alternatively, a numeric vector (e.g., xrange = c(6, 20)) can be provided to set the range within the specified values.

A vector of length two or a character string 'range' to set the range of the predictor variable (x) within which growth parameters are searched. This is useful when there is more than one peak and the user wants to summarize the peak within a specified range of the x variable. The default value is xrange_search = NULL.

A logical value to specify whether or not to compute growth parameters on the fly. This is for internal use only and is mainly needed for compatibility across internal functions.

A character string specifying the method used when evaluating parms_eval. The default method is getPeak, which uses the sitar::getPeak() function from the sitar package. Alternatively, findpeaks uses the findpeaks function from the pracma package. This parameter is for internal use and ensures compatibility across internal functions.

A character string specifying the interpolation method (default NULL). The number of interpolation points is controlled by the ipts argument.

Available options:

• 'm1': Adapted from the iapvbs package (documented here). This method internally constructs the data frame based on the model configuration. Note: This method may fail if the model includes covariates (especially in univariate_by models).

• 'm2': Based on the JMbayes package (documented here). This method uses the exact data frame from the model fit (fit$data) and is generally more robust.

If idata_method = NULL, 'm2' is automatically selected. It is recommended to use 'm2' if 'm1' encounters errors with covariate-dependent models.

A logical argument (default FALSE) to specify whether to print information collected during the setup of the object(s).

A logical value indicating whether to return a fullframe object in which newdata is bound to the summary estimates. Note that fullframe cannot be used with summary = FALSE, and it is only applicable when idata_method = 'm2'. A typical use case is when fitting a univariate_by model. This option is mainly for internal use.

A named list (default NULL) to convert dummy variables into a factor variable. The list must include the following elements:

• factor.dummy: A character vector of dummy variables to be converted to factors.

• factor.name: The name for the newly created factor variable (default is 'factor.var' if NULL).

• factor.level: A vector specifying the factor levels. If NULL, levels are taken from factor.dummy. If factor.level is provided, its length must match factor.dummy.

A logical value or NULL (default FALSE). Controls whether Stan functions are exposed for post-processing.

• TRUE: Explicitly exposes Stan functions saved from the model fit. This is required when calculating fit criteria or bayes_R2 during model optimization.

• FALSE (default): Stan functions are not exposed.

• NULL: The setting is inherited from the original bsitar() model fit configuration.

Note: In the optimize_model() function, the default is NULL (inheriting behavior), whereas other post-processing functions default to FALSE.

A logical value (default NULL) indicating whether to use already exposed and saved Stan functions. This is typically set automatically based on the expose_functions argument from the bsitar() call. Manual specification of usesavedfuns is rarely needed and is intended for internal testing, as improper use can lead to unreliable estimates.

A logical value indicating whether to clear the exposed Stan functions from the environment (TRUE) or not (FALSE). If NULL, clearenvfuns is set based on the value of usesavedfuns: TRUE if usesavedfuns = TRUE, or FALSE if usesavedfuns = FALSE.

A list (default NULL) specifying function names. This is rarely needed, as required functions are typically retrieved automatically. A use case for funlist is when sigma_formula, sigma_formula_gr, or sigma_formula_gr_str use an external function (e.g., poly(age)). The funlist should include function names defined in the globalenv(). For functions needing both distance and velocity curves (e.g., plot_curves(..., opt = 'dv')), funlist must include two functions: one for the distance curve and one for the velocity curve.

A character string (default NULL) specifying the 'x' variable. Rarely used because xvar is inferred internally. A use case is when conflicting variables exist (e.g., sigma_formula) and user wants to set a specific variable as 'x'.

A character string (default NULL) specifying the 'x' variable that should be used for manual differentiation of the distance curve. Internally, the xvar is set as difx if specified. The argument difx is evaluated only when dpar = 'sigma', ignored otherwise. Note that argument xvar itself is rarely used because xvar is inferred internally. A use case is when conflicting variables exist (e.g., sigma_formula) and user wants to set a specific variable as 'x'.

A character string (default NULL) specifying the 'id' variable. Rarely used because idvar is inferred internally.

A character string (default NULL) indicating the variables names that are reverse transformed. Options are c("x", "y", "sigma"). The itransform is primarily used to get the xvar variable at original scale i.e., itransform = 'x'. To turn of all transformations, use itransform = "". when itransform = NULL, the appropriate transformation for xvar is selected automatically. Note that when no match for xvar is found in the data,frame, the itransform will be ignored within the calling function, 'prepare_transformations()'.

An indicator specifying how to handle the format and structure of the user-provided newdata before use. The available options are:

• NULL (default): Checks the data format and structure of newdata, and applies prepare_data2 and prepare_transformations as needed.

• TRUE: Returns newdata as-is, without any format checking or transformation.

• 0: Same as TRUE — returns newdata without checking the data format or applying prepare_data2 and prepare_transformations.

• 1: Checks the data format and, if needed, prepares it using prepare_data2, but does not apply prepare_transformations.

• 2: Applies prepare_transformations only, assuming the data format is already correct.

It is strongly recommended to leave both newdata = NULL and newdata_fixed = NULL (the defaults). In this case, the data used during model fitting is automatically retrieved, checked for the required format, and transformed as needed via prepare_data2 and prepare_transformations. The flags newdata_fixed = 0, 1, and 2 are intended primarily for internal use.

The environment used for function evaluation. The default is NULL, which sets the environment to parent.frame(). Since most post-processing functions rely on brms, it is recommended to set envir = globalenv() or envir = .GlobalEnv, especially for derivatives like velocity curves.

Additional arguments passed to the brms::fitted.brmsfit() function. For details on available options, please refer to brms::fitted.brmsfit().

An object of class bgmfit.

A character string (default NULL) to specify the response variable when processing posterior draws for univariate_by and multivariate models. See bsitar() for details on univariate_by and multivariate models.

Optional name of a predicted distributional parameter. If specified, expected predictions of this parameters are returned.

A positive integer indicating the number of posterior draws to use in estimation. If NULL (default), all draws are used.

An integer specifying the specific posterior draw(s) to use in estimation (default NULL).

An optional data frame for estimation. If NULL (default), newdata is retrieved from the model.

A data frame or named list for setting up a custom grid of predictor values to evaluate the quantities of interest. If NULL (default), no custom grid is used. The grid can be constructed using marginaleffects::datagrid(). If datagrid = list(), essential arguments such as model and newdata are inferred automatically from the respective arguments in the model fit.

Option to indicate whether or not to include individual/group-level effects in the estimation. When NA (default), individual-level effects are excluded, and population average growth parameters are computed. When NULL, individual-level effects are included in the computation, and the resulting growth parameters are individual-specific. In both cases (NA or NULL), continuous and factor covariates are appropriately included in the estimation. Continuous covariates are set to their means by default (see numeric_cov_at for details), while factor covariates remain unaltered, allowing for the estimation of covariate-specific population average and individual-specific growth parameters.

A named list of objects containing new data, which cannot be passed via argument newdata. Required for some objects used in autocorrelation structures, or stanvars.

A flag indicating if new levels of group-level effects are allowed (defaults to FALSE). Only relevant if newdata is provided.

Indicates how to sample new levels for grouping factors specified in re_formula. This argument is only relevant if newdata is provided and allow_new_levels is set to TRUE. If "uncertainty" (default), each posterior sample for a new level is drawn from the posterior draws of a randomly chosen existing level. Each posterior sample for a new level may be drawn from a different existing level such that the resulting set of new posterior draws represents the variation across existing levels. If "gaussian", sample new levels from the (multivariate) normal distribution implied by the group-level standard deviations and correlations. This options may be useful for conducting Bayesian power analysis or predicting new levels in situations where relatively few levels where observed in the old_data. If "old_levels", directly sample new levels from the existing levels, where a new level is assigned all of the posterior draws of the same (randomly chosen) existing level.

An integer to set the predictor range (e.g., age) when executing the interpolation via ipts. By default, NULL sets the individual-specific predictor range. Setting xrange = 1 applies the same range for individuals within the same higher grouping variable (e.g., study). Setting xrange = 2 applies an identical range across the entire sample. Alternatively, a numeric vector (e.g., xrange = c(6, 20)) can be provided to set the range within the specified values.

An integer (default 2) for rounding the estimates to the specified number of decimal places using base::round().

An optional (named list) argument to specify the value of continuous covariate(s). The default NULL option sets the continuous covariate(s) to their mean. Alternatively, a named list can be supplied to manually set these values. For example, numeric_cov_at = list(xx = 2) will set the continuous covariate variable 'xx' to 2. The argument numeric_cov_at is ignored when no continuous covariates are included in the model.

An optional argument to specify the variable(s) that can be passed to the ipts argument (see below). This is useful when fitting location-scale models and measurement error models. If post-processing functions throw an error such as variable 'x' not found in either 'data' or 'data2', consider using aux_variables.

An optional argument to specify the variable(s) that can be passed to the marginaleffects::datagrid(). This is useful when fitting location-scale models and measurement error models. If post-processing functions throw an error such as variable 'x' not found in either 'data' or 'data2', consider using grid_add. Note that unlike aux_variables which are passed to the internal data functions such as 'get.newdata', the grid_add are passed to the marginaleffects::datagrid().

An optional argument to specify the ids for the hierarchical model (default NULL). It is used only when the model is applied to data with three or more levels of hierarchy. For a two-level model, levels_id is automatically inferred from the model fit. For models with three or more levels, levels_id is inferred from the model fit under the assumption that hierarchy is specified from the lowest to the uppermost level, i.e., id followed by study, where id is nested within study. However, it is not guaranteed that levels_id is sorted correctly, so it is better to set it manually when fitting a model with three or more levels of hierarchy.

An optional argument (default NULL) to calculate (marginal/average) curves and growth parameters, such as APGV and PGV. If specified, it must be a named list indicating the over (typically a level 1 predictor, such as age), feby (fixed effects, typically a factor variable), and reby (typically NULL, indicating that parameters are integrated over the random effects). For example, avg_reffects = list(feby = 'study', reby = NULL, over = 'age').

A character string specifying the interpolation method (default NULL). The number of interpolation points is controlled by the ipts argument.

Available options:

• 'm1': Adapted from the iapvbs package (documented here). This method internally constructs the data frame based on the model configuration. Note: This method may fail if the model includes covariates (especially in univariate_by models).

• 'm2': Based on the JMbayes package (documented here). This method uses the exact data frame from the model fit (fit$data) and is generally more robust.

If idata_method = NULL, 'm2' is automatically selected. It is recommended to use 'm2' if 'm1' encounters errors with covariate-dependent models.

An integer specifying the number of points for interpolating the predictor variable (e.g., age) to generate smooth curves for predictions and plots. This value is used as the length.out argument for seq(), controlling the smoothness of distance and velocity curves without altering the predictor range.

NULL (the default)

Engages automatic behavior based on the dpar argument: it is internally set to 50 for the mean response (dpar = 'mu'), ensuring smooth curves, while for the distributional parameters (e.g., dpar = 'sigma'), no interpolation is performed.

An integer (e.g., 100)

Explicitly sets the number of interpolation points.

FALSE

Disables all interpolation, forcing predictions to be made only at the original data points of the predictor variable.

This argument affects the following post-processing functions:
fitted_draws(), predict_draws(), growthparameters(), plot_curves(), get_predictions(), get_comparisons(), and get_growthparameters().

An integer (default 123) that is passed to the estimation method to ensure reproducibility.

A logical value (default FALSE) indicating whether to perform parallel computations. If TRUE, posterior summaries are computed in parallel using future.apply::future_sapply().

A character string or a named list specifying the parallel plan when future = TRUE.

• Character string: Defaults to "multisession". Use "multicore" for forking (not supported on Windows).

• Named list: Useful for advanced plans like future.mirai::mirai_cluster(). The list must contain:

  • future_session: The planner function (e.g., mirai_cluster).

  • Additional named arguments passed to the planner (e.g., daemons for mirai::daemons()).

  Example: list(future_session = mirai_cluster, daemons = list(...)).

A list, numeric vector, or logical value (default TRUE). Controls how posterior draws are partitioned into smaller subsets for parallel computation. This helps manage memory and improve performance, particularly on Linux when using future::plan("multicore").

Input options:

• List: (e.g., list(1:6, 7:10)). Each element is a sequence of numbers passed to draw_ids to process in a separate chunk.

• Numeric vector of length 2: (e.g., c(100, 4)). The first element is the total number of draws, and the second is the number of splits. Indices are generated via parallel::splitIndices().

• Numeric vector of length 1: Specifies the total number of draws. Splits are calculated automatically.

• TRUE: Automatically creates splits based on ndraws and cores. If both ndraws and draw_ids are specified, draw_ids takes precedence. This is consistent with post-processing functions in the bsitar and brms packages.

• FALSE: Splitting is disabled.

Note: On Windows with future::plan("multisession"), background R processes may not free resources automatically. If needed, use installr::kill_all_Rscript_s() to terminate them.

A character string (default 'future') specifying the parallel computation backend. Options include:

• 'future': Uses future.apply::future_lapply() for execution.

• 'dofuture': Uses foreach::foreach() via the doFuture package.

A logical value (default NULL) indicating whether to re-expose internal Stan functions when future = TRUE. This is critical when future::plan() is set to "multisession", as compiled C++ functions cannot be exported across distinct R sessions.

• If NULL (default), it is automatically set to TRUE when the plan is "multisession".

• Explicitly setting this to TRUE is recommended for improved performance during parallel execution.

A logical (default FALSE) indicating whether to use the dtplyr package for summarizing the draws. This package uses data.table as a back-end. It is useful when the data has a large number of observations. For typical use cases, it does not make a significant performance difference. The usedtplyr argument is evaluated only when method = 'custom'.

A logical (default FALSE) to indicate whether to use the collapse package for summarizing the draws.

An integer specifying the number of cores for parallel execution.

• If NULL (default), the number of cores is set to future::availableCores() - 1.

• On non-Windows systems, this can be controlled globally via the mc.cores option.

A logical value indicating whether to return a fullframe object in which newdata is bound to the summary estimates. Note that fullframe cannot be used with summary = FALSE, and it is only applicable when idata_method = 'm2'. A typical use case is when fitting a univariate_by model. This option is mainly for internal use.

A logical indicating whether to call marginaleffects::comparisons() (if FALSE) or marginaleffects::avg_comparisons() (if TRUE). Default is FALSE.

A logical indicating whether to plot the comparisons using marginaleffects::plot_comparisons(). Default is FALSE.

A named list that can be used to pass the aesthetic mapping and facet arguments to the ggplot2 object (default NULL). The main use of this is to get overlay line plot instead of separate line plot for each id variable. This can also be used to remove a layer. For example, if user wants to remove the confidence intervals band from the lines i.e., geom_ribbon, then include rm_geom_ribbon = TRUE. Note the prefix 'rm_'. This is a general approach in which any layer name with prefix 'rm_' included in the mapping_facet will be removed. An example is shown below:
mapping_facet = list(group = 'id', colour = 'id', facet_wrap = c('study'), rm_geom_ribbon = TRUE).

A logical value to specify whether to show legends (TRUE) or not (FALSE). If NULL (default), the value of showlegends is internally set to TRUE if re_formula = NA, and FALSE if re_formula = NULL.

A named list specifying the level 1 predictor, such as age or time, used for estimating growth parameters in the current use case. The variables list is set via the esp argument (default value is 1e-6). If variables is NULL, the relevant information is retrieved internally from the model. Alternatively, users can define variables as a named list, e.g., variables = list('x' = 1e-6) where 'x' is the level 1 predictor. By default, variables = list('age' = 1e-6) in the marginaleffects package, as velocity is usually computed by differentiating the distance curve using the dydx approach. When using this default, the argument deriv is automatically set to 0 and model_deriv to FALSE. If parameters are to be estimated based on the model's first derivative, deriv must be set to 1 and variables will be defined as variables = list('age' = 0). Note that if the default behavior is used (deriv = 0 and variables = list('x' = 1e-6)), additional arguments cannot be passed to variables. In contrast, when using an alternative approach (deriv = 0 and variables = list('x' = 0)), additional options can be passed to the marginaleffects::comparisons() and marginaleffects::avg_comparisons() functions.

Conditional slopes

• Character vector (max length 4): Names of the predictors to display.

• Named list (max length 4): List names correspond to predictors. List elements can be:

  • Numeric vector

  • Function which returns a numeric vector or a set of unique categorical values

  • Shortcut strings for common reference values: "minmax", "quartile", "threenum"

• 1: x-axis. 2: color/shape. 3: facet (wrap if no fourth variable, otherwise cols of grid). 4: facet (rows of grid).

• Numeric variables in positions 2 and 3 are summarized by Tukey's five numbers ?stats::fivenum.

A numeric value indicating whether to estimate parameters based on the differentiation of the distance curve or the model's first derivative.

A logical value specifying whether to estimate the velocity curve from the derivative function or by differentiating the distance curve. Set model_deriv = TRUE for functions that require the velocity curve, such as growthparameters() and plot_curves(). Set it to NULL for functions that use the distance curve (i.e., fitted values), such as loo_validation() and plot_ppc().

A character string specifying whether to compute comparisons using the marginaleffects machinery (method = 'pkg') or via custom functions for efficiency (method = 'custom'). Default is 'pkg'. The 'custom' method is useful when testing hypotheses and should be used cautiously.

A character string specifying whether to compute comparisons using marginaleffects::predictions() (method_call = 'predictions') or marginaleffects::comparisons() (method_call = 'comparisons'). The 'method_call' is evaluated only when method = 'custom' ignored otherwise. Default method_call = NULL allows automatic selection of the 'method_call' depending on whether comparisons are for slopes or predictions.

A list, data.frame, or tibble of velocity returned by the marginaleffects functions (default NULL). This is only evaluated when method = 'custom'. The marginals can be the output from marginaleffects functions or posterior draws from marginaleffects::posterior_draws(). The marginals argument is primarily used for internal purposes only.

A character string (default FALSE) to indicate whether to return the original posterior draws for parameters. Options include:

• 'return': returns the original posterior draws,

• 'add': adds the original posterior draws to the outcome.

When pdrawso = TRUE, the default behavior is pdrawso = 'return'. Note that the posterior draws are returned before calling marginaleffects::posterior_draws().

A character string (default FALSE) to indicate whether to return the posterior draws for parameters. Options include:

• 'return': returns the posterior draws for parameters,

• 'add': adds the posterior draws to the outcome.

When pdrawsp = TRUE, the default behavior is pdrawsp = 'return'. The pdrawsp represent the parameter estimates for each of the posterior samples, and the summary of these are the estimates returned.

A character string (default FALSE) to indicate whether to return the posterior draws for parameter contrasts. Options include:

• 'return': returns the posterior draws for contrasts.

The summary of posterior draws for parameters is the default returned object. The pdrawsh represent the contrast estimates for each of the posterior samples, and the summary of these are the contrast returned.

A character string specifying the comparison type for growth parameter estimation. Options are 'difference' and 'differenceavg'. This argument sets up the internal function for estimating parameters using sitar::getPeak(), sitar::getTakeoff(), and sitar::getTrough() functions. These options are restructured according to the user-specified hypothesis argument.

string indicates the type (scale) of the predictions used to compute contrasts or slopes. This can differ based on the model type, but will typically be a string such as: "response", "link", "probs", or "zero". When an unsupported string is entered, the model-specific list of acceptable values is returned in an error message. When type is NULL, the first entry in the error message is used by default. See the Type section in the documentation below.

Aggregate unit-level estimates (aka, marginalize, average over). Valid inputs:

• FALSE: return the original unit-level estimates.

• TRUE: aggregate estimates for each term.

• Character vector of column names in newdata or in the data frame produced by calling the function without the by argument.

• Data frame with a by column of group labels, and merging columns shared by newdata or the data frame produced by calling the same function without the by argument.

• See examples below.

• For more complex aggregations, you can use the FUN argument of the hypotheses() function. See that function's documentation and the Hypothesis Test vignettes on the marginaleffects website.

numeric value between 0 and 1. Confidence level to use to build a confidence interval.

string or function. Transformation applied to unit-level estimates and confidence intervals just before the function returns results. Functions must accept a vector and return a vector of the same length. Support string shortcuts: "exp", "ln"

A function applied to individual draws from the posterior distribution before computing summaries (default NULL). The argument transform_draws is derived from the marginaleffects::predictions() function and should not be confused with the transform argument from the deprecated brms::posterior_predict() function. It's important to note that for both marginaleffects::predictions() and marginaleffects::avg_predictions(), the transform_draws argument takes precedence over the transform argument. Note that when transform_draws = NULL, an attempt is made to automatically set transform_draws = 'exp' for dpar = 'sigma'. User can set transform_draws = FALSE to turn off this automatic assignment of 'exp' to the transform_draws. It is also important to set transform_draws = FALSE when computing the first derivative (velocity) for dpar = 'sigma'.

• FALSE: Contrasts represent the change in adjusted predictions when one predictor changes and all other variables are held constant.

• TRUE: Contrasts represent the changes in adjusted predictions when all the predictors specified in the variables argument are manipulated simultaneously (a "cross-contrast").

logical, string or numeric: weights to use when computing average predictions, contrasts or slopes. These weights only affect the averaging in ⁠avg_*()⁠ or with the by argument, and not unit-level estimates. See ?weighted.mean

• string: column name of the weights variable in newdata. When supplying a column name to wts, it is recommended to supply the original data (including the weights variable) explicitly to newdata.

• numeric: vector of length equal to the number of rows in the original data or in newdata (if supplied).

• FALSE: Equal weights.

• TRUE: Extract weights from the fitted object with insight::find_weights() and use them when taking weighted averages of estimates. Warning: newdata=datagrid() returns a single average weight, which is equivalent to using wts=FALSE

specify a hypothesis test or custom contrast using a number , formula, string equation, vector, matrix, or function.

• Number: The null hypothesis used in the computation of Z and p (before applying transform).

• String: Equation to specify linear or non-linear hypothesis tests. Two-tailed tests must include an equal = sign. One-tailed tests must start with < or >. If the terms in coef(object) uniquely identify estimates, they can be used in the formula. Otherwise, use b1, b2, etc. to identify the position of each parameter. The ⁠b*⁠ wildcard can be used to test hypotheses on all estimates. When the hypothesis string represents a two-sided equation, the estimate column holds the value of the left side minus the right side of the equation. If a named vector is used, the names are used as labels in the output. Examples:

  • hp = drat

  • hp + drat = 12

  • b1 + b2 + b3 = 0

  • ⁠b* / b1 = 1⁠

  • ⁠<= 0⁠

  • ⁠>= -3.5⁠

  • b1 >= 10

• Formula: lhs ~ rhs | group

  • lhs

    • ratio (null = 1)

    • difference (null = 0)

    • Leave empty for default value

  • rhs

    • pairwise and revpairwise: pairwise differences between estimates in each row.

    • reference: differences between the estimates in each row and the estimate in the first row.

    • sequential: difference between an estimate and the estimate in the next row.

    • meandev: difference between an estimate and the mean of all estimates.

    • meanotherdev: difference between an estimate and the mean of all other estimates, excluding the current one.

    • poly: polynomial contrasts, as computed by the stats::contr.poly() function.

    • helmert: Helmert contrasts, as computed by the stats::contr.helmert() function. Contrast 2nd level to the first, 3rd to the average of the first two, and so on.

    • trt_vs_ctrl: difference between the mean of estimates (except the first) and the first estimate.

    • I(fun(x)): custom function to manipulate the vector of estimates x. The function fun() can return multiple (potentially named) estimates.

  • group (optional)

    • Column name of newdata. Conduct hypothesis tests withing subsets of the data.

  • Examples:

    • ~ poly

    • ~ sequential | groupid

    • ~ reference

    • ratio ~ pairwise

    • difference ~ pairwise | groupid

    • ~ I(x - mean(x)) | groupid

    • ⁠~ I(\(x) c(a = x[1], b = mean(x[2:3]))) | groupid⁠

• Matrix or Vector: Each column is a vector of weights. The the output is the dot product between these vectors of weights and the vector of estimates. The matrix can have column names to label the estimates.

• Function:

  • Accepts an argument x: object produced by a marginaleffects function or a data frame with column rowid and estimate

  • Returns a data frame with columns term and estimate (mandatory) and rowid (optional).

  • The function can also accept optional input arguments: newdata, by, draws.

  • This function approach will not work for Bayesian models or with bootstrapping. In those cases, it is easy to use get_draws() to extract and manipulate the draws directly.

• See the Examples section below and the vignette: https://marginaleffects.com/chapters/hypothesis.html

• Warning: When calling predictions() with type="invlink(link)" (the default in some models), hypothesis is tested and p values are computed on the link scale.

Numeric vector of length 2: bounds used for the two-one-sided test (TOST) of equivalence, and for the non-inferiority and non-superiority tests. For bayesian models, this report the proportion of posterior draws in the interval and the ROPE. See Details section below.

NULL or numeric value which determines the step size to use when calculating numerical derivatives: (f(x+eps)-f(x))/eps. When eps is NULL, the step size is 0.0001 multiplied by the difference between the maximum and minimum values of the variable with respect to which we are taking the derivative. Changing eps may be necessary to avoid numerical problems in certain models.

A character vector (default NULL) specifying the variables by which hypotheses should be tested (e.g., for post-draw comparisons). These variables should be a subset of the variables in the by argument of marginaleffects::comparisons().

A character vector (default NULL) specifying the values at which hypotheses should be tested. Useful for large estimates to limit testing to fewer rows.

A logical (default TRUE) to reformat the output from marginaleffects::comparisons() as a data.frame, with column names such as conf.low as Q2.5, conf.high as Q97.5, and dropping unnecessary columns like term, contrast, and predicted.

A character string (default NULL) specifying how to center estimates: either 'mean' or 'median'. This option sets the global options as follows: options("marginaleffects_posterior_center" = "mean") or options("marginaleffects_posterior_center" = "median"). These global options are restored upon function exit using base::on.exit().

A character string (default NULL) to specify the type of credible intervals: 'eti' for equal-tailed intervals or 'hdi' for highest density intervals. This option sets the global options as follows: options("marginaleffects_posterior_interval" = "eti") or options("marginaleffects_posterior_interval" = "hdi"), and is restored on exit using base::on.exit().

A named list (default NULL) to convert dummy variables into a factor variable. The list must include the following elements:

• factor.dummy: A character vector of dummy variables to be converted to factors.

• factor.name: The name for the newly created factor variable (default is 'factor.var' if NULL).

• factor.level: A vector specifying the factor levels. If NULL, levels are taken from factor.dummy. If factor.level is provided, its length must match factor.dummy.

A logical argument (default FALSE) to specify whether to print information collected during the setup of the object(s).

A logical value or NULL (default FALSE). Controls whether Stan functions are exposed for post-processing.

• TRUE: Explicitly exposes Stan functions saved from the model fit. This is required when calculating fit criteria or bayes_R2 during model optimization.

• FALSE (default): Stan functions are not exposed.

• NULL: The setting is inherited from the original bsitar() model fit configuration.

Note: In the optimize_model() function, the default is NULL (inheriting behavior), whereas other post-processing functions default to FALSE.

A logical value (default NULL) indicating whether to use already exposed and saved Stan functions. This is typically set automatically based on the expose_functions argument from the bsitar() call. Manual specification of usesavedfuns is rarely needed and is intended for internal testing, as improper use can lead to unreliable estimates.

A logical value indicating whether to clear the exposed Stan functions from the environment (TRUE) or not (FALSE). If NULL, clearenvfuns is set based on the value of usesavedfuns: TRUE if usesavedfuns = TRUE, or FALSE if usesavedfuns = FALSE.

A list (default NULL) specifying function names. This is rarely needed, as required functions are typically retrieved automatically. A use case for funlist is when sigma_formula, sigma_formula_gr, or sigma_formula_gr_str use an external function (e.g., poly(age)). The funlist should include function names defined in the globalenv(). For functions needing both distance and velocity curves (e.g., plot_curves(..., opt = 'dv')), funlist must include two functions: one for the distance curve and one for the velocity curve.

A character string (default NULL) specifying the 'x' variable. Rarely used because xvar is inferred internally. A use case is when conflicting variables exist (e.g., sigma_formula) and user wants to set a specific variable as 'x'.

A character string (default NULL) specifying the 'x' variable that should be used for manual differentiation of the distance curve. Internally, the xvar is set as difx if specified. The argument difx is evaluated only when dpar = 'sigma', ignored otherwise. Note that argument xvar itself is rarely used because xvar is inferred internally. A use case is when conflicting variables exist (e.g., sigma_formula) and user wants to set a specific variable as 'x'.

A character string (default NULL) specifying the 'id' variable. Rarely used because idvar is inferred internally.

A character string (default NULL) indicating the variables names that are reverse transformed. Options are c("x", "y", "sigma"). The itransform is primarily used to get the xvar variable at original scale i.e., itransform = 'x'. To turn of all transformations, use itransform = "". when itransform = NULL, the appropriate transformation for xvar is selected automatically. Note that when no match for xvar is found in the data,frame, the itransform will be ignored within the calling function, 'prepare_transformations()'.

An indicator specifying how to handle the format and structure of the user-provided newdata before use. The available options are:

• NULL (default): Checks the data format and structure of newdata, and applies prepare_data2 and prepare_transformations as needed.

• TRUE: Returns newdata as-is, without any format checking or transformation.

• 0: Same as TRUE — returns newdata without checking the data format or applying prepare_data2 and prepare_transformations.

• 1: Checks the data format and, if needed, prepares it using prepare_data2, but does not apply prepare_transformations.

• 2: Applies prepare_transformations only, assuming the data format is already correct.

It is strongly recommended to leave both newdata = NULL and newdata_fixed = NULL (the defaults). In this case, the data used during model fitting is automatically retrieved, checked for the required format, and transformed as needed via prepare_data2 and prepare_transformations. The flags newdata_fixed = 0, 1, and 2 are intended primarily for internal use.

The environment used for function evaluation. The default is NULL, which sets the environment to parent.frame(). Since most post-processing functions rely on brms, it is recommended to set envir = globalenv() or envir = .GlobalEnv, especially for derivatives like velocity curves.

Additional arguments passed to the brms::fitted.brmsfit() and brms::predict() functions.