bsitar {bsitar}R Documentation

Fit Bayesian SITAR Growth Curve Model

Description

The bsitar() function fits the Bayesian version of the Super Imposition by Translation and Rotation (SITAR) model. The SITAR model is a nonlinear mixed-effects model that has been widely used to summarize growth processes (such as height and weight) from early childhood through adulthood.

The frequentist version of the SITAR model can be fit using the already available R package, sitar (Cole 2022). However, the bsitar package offers an enhanced Bayesian implementation that improves modeling capabilities. In addition to univariate analysis (i.e., modeling a single outcome), bsitar also supports:

The Bayesian implementation in bsitar provides a more flexible and robust framework for growth curve modeling compared to the frequentist approach.

Usage

bsitar(
  x,
  y,
  id,
  data,
  df = 4,
  knots = NA,
  fixed = a + b + c,
  random = a + b + c,
  xoffset = mean,
  bstart = xoffset,
  cstart = 0,
  xfun = NULL,
  yfun = NULL,
  bound = 0.04,
  stype = nsp,
  terms_rhs = NULL,
  a_formula = ~1,
  b_formula = ~1,
  c_formula = ~1,
  d_formula = ~1,
  s_formula = ~1,
  a_formula_gr = ~1,
  b_formula_gr = ~1,
  c_formula_gr = ~1,
  d_formula_gr = ~1,
  a_formula_gr_str = NULL,
  b_formula_gr_str = NULL,
  c_formula_gr_str = NULL,
  d_formula_gr_str = NULL,
  d_adjusted = FALSE,
  sigma_formula = NULL,
  sigma_formula_gr = NULL,
  sigma_formula_gr_str = NULL,
  sigma_formula_manual = NULL,
  sigmax = NULL,
  sigmadf = 4,
  sigmaknots = NA,
  sigmafixed = a + b + c,
  sigmarandom = "",
  sigmaxoffset = mean,
  sigmabstart = sigmaxoffset,
  sigmacstart = 0,
  sigmaxfun = NULL,
  sigmabound = 0.04,
  dpar_formula = NULL,
  autocor_formula = NULL,
  family = gaussian(),
  custom_family = NULL,
  custom_stanvars = NULL,
  group_arg = list(groupvar = NULL, by = NULL, cor = un, cov = NULL, dist = gaussian),
  sigma_group_arg = list(groupvar = NULL, by = NULL, cor = un, cov = NULL, dist =
    gaussian),
  univariate_by = list(by = NA, cor = un, terms = subset),
  multivariate = list(mvar = FALSE, cor = un, rescor = TRUE),
  a_prior_beta = normal(lm, ysd, autoscale = TRUE),
  b_prior_beta = normal(0, 1.5, autoscale = FALSE),
  c_prior_beta = normal(0, 0.5, autoscale = FALSE),
  d_prior_beta = normal(0, 1, autoscale = TRUE),
  s_prior_beta = normal(lm, lm, autoscale = TRUE),
  a_cov_prior_beta = normal(0, 5, autoscale = FALSE),
  b_cov_prior_beta = normal(0, 1, autoscale = FALSE),
  c_cov_prior_beta = normal(0, 0.1, autoscale = FALSE),
  d_cov_prior_beta = normal(0, 1, autoscale = FALSE),
  s_cov_prior_beta = normal(lm, lm, autoscale = TRUE),
  a_prior_sd = normal(0, ysd, autoscale = FALSE),
  b_prior_sd = normal(0, 1, autoscale = FALSE),
  c_prior_sd = normal(0, 0.25, autoscale = FALSE),
  d_prior_sd = normal(0, 1, autoscale = TRUE),
  a_cov_prior_sd = normal(0, 5, autoscale = FALSE),
  b_cov_prior_sd = normal(0, 1, autoscale = FALSE),
  c_cov_prior_sd = normal(0, 0.1, autoscale = FALSE),
  d_cov_prior_sd = normal(0, 1, autoscale = FALSE),
  a_prior_sd_str = NULL,
  b_prior_sd_str = NULL,
  c_prior_sd_str = NULL,
  d_prior_sd_str = NULL,
  a_cov_prior_sd_str = NULL,
  b_cov_prior_sd_str = NULL,
  c_cov_prior_sd_str = NULL,
  d_cov_prior_sd_str = NULL,
  sigma_prior_beta = normal(0, 1, autoscale = FALSE),
  sigma_cov_prior_beta = normal(0, 0.5, autoscale = FALSE),
  sigma_prior_sd = normal(0, 0.25, autoscale = FALSE),
  sigma_cov_prior_sd = normal(0, 0.15, autoscale = FALSE),
  sigma_prior_sd_str = NULL,
  sigma_cov_prior_sd_str = NULL,
  rsd_prior_sigma = normal(0, ysd, autoscale = FALSE),
  dpar_prior_sigma = normal(0, ysd, autoscale = FALSE),
  dpar_cov_prior_sigma = normal(0, 1, autoscale = FALSE),
  autocor_prior_acor = uniform(-1, 1, autoscale = FALSE),
  autocor_prior_unstr_acor = lkj(1),
  gr_prior_cor = lkj(1),
  gr_prior_cor_str = lkj(1),
  sigma_prior_cor = lkj(1),
  sigma_prior_cor_str = lkj(1),
  mvr_prior_rescor = lkj(1),
  init = NULL,
  init_r = NULL,
  a_init_beta = lm,
  b_init_beta = 0,
  c_init_beta = 0,
  d_init_beta = random,
  s_init_beta = lm,
  a_cov_init_beta = random,
  b_cov_init_beta = random,
  c_cov_init_beta = random,
  d_cov_init_beta = random,
  s_cov_init_beta = random,
  a_init_sd = random,
  b_init_sd = random,
  c_init_sd = random,
  d_init_sd = random,
  a_cov_init_sd = random,
  b_cov_init_sd = random,
  c_cov_init_sd = random,
  d_cov_init_sd = random,
  sigma_init_beta = random,
  sigma_cov_init_beta = random,
  sigma_init_sd = random,
  sigma_cov_init_sd = random,
  gr_init_cor = random,
  sigma_init_cor = random,
  rsd_init_sigma = random,
  dpar_init_sigma = random,
  dpar_cov_init_sigma = random,
  autocor_init_acor = random,
  autocor_init_unstr_acor = random,
  mvr_init_rescor = random,
  r_init_z = random,
  vcov_init_0 = FALSE,
  jitter_init_beta = NULL,
  jitter_init_sd = NULL,
  jitter_init_cor = NULL,
  prior_data = NULL,
  init_data = NULL,
  init_custom = NULL,
  verbose = FALSE,
  expose_function = FALSE,
  get_stancode = FALSE,
  get_standata = FALSE,
  get_formula = FALSE,
  get_stanvars = FALSE,
  get_priors = FALSE,
  get_priors_eval = FALSE,
  get_init_eval = FALSE,
  validate_priors = FALSE,
  set_self_priors = NULL,
  add_self_priors = NULL,
  set_replace_priors = NULL,
  set_same_priors_hierarchy = FALSE,
  outliers = NULL,
  unused = NULL,
  chains = 4,
  iter = 2000,
  warmup = floor(iter/2),
  thin = 1,
  cores = getOption("mc.cores", "optimize"),
  backend = getOption("brms.backend", "rstan"),
  threads = getOption("brms.threads", "optimize"),
  opencl = getOption("brms.opencl", NULL),
  normalize = getOption("brms.normalize", TRUE),
  algorithm = getOption("brms.algorithm", "sampling"),
  control = list(adapt_delta = 0.8, max_treedepth = 15),
  empty = FALSE,
  rename = TRUE,
  pathfinder_args = NULL,
  pathfinder_init = FALSE,
  sample_prior = "no",
  save_pars = NULL,
  drop_unused_levels = TRUE,
  stan_model_args = list(),
  refresh = NULL,
  silent = 1,
  seed = 123,
  save_model = NULL,
  fit = NA,
  file = NULL,
  file_compress = TRUE,
  file_refit = getOption("brms.file_refit", "never"),
  future = getOption("future", FALSE),
  parameterization = "ncp",
  ...
)

Arguments

x

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, you can specify x = list(x1) or simply x = x1.

y

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.

id

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

data

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

df

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.

knots

A numeric vector specifying the knots for the natural cubic spline design matrix (default NULL). Note that you 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).

fixed

A character string specifying the fixed effects structure (default 'a+b+c'). For univariate_by and multivariate models, you 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'.

random

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

xoffset

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

bstart

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

cstart

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

xfun

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

yfun

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

bound

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

stype

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

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

  • 'nks': 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 default is 'nsp'.

The 'rcs' method uses a truncated power basis, whereas 'nks' and 'nsp' are B-spline-based methods. Unlike splines2::nsp() and splines2::nsk(), which normalize the spline basis by default, 'nks' 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().

terms_rhs

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 used in fitting univariate_by models.

a_formula

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.

b_formula

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.

c_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.

d_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.

s_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.

a_formula_gr

Formula for the random effect parameter, a (default ~ 1). Similar to a_formula, users can specify different formulas when fitting univariate_by and multivariate models. The formula can include continuous and/or factor variables, including their interactions as covariates (see a_formula for details). In addition to setting up the design matrix for the random effect parameter a, users can define the group identifier and the correlation structure for random effects using the vertical bar || notation. For example, to include only an intercept for the random effects a, b, and c, you can specify:

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

To specify the group identifier (e.g., id) and an unstructured correlation structure, use the vertical bar notation:

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

Here, i within the vertical bars is a placeholder, and a common identifier (e.g., i) shared across the random effect formulas will model them as unstructured correlated random effects. For more details on this vertical bar approach, please see [brms::brm()].

An alternative approach to specify the group identifier and correlation structure is through the group_by argument. To achieve the same setup as described above with the vertical bar approach, users can define the formula part as:

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

and use group_by as group_by = list(groupvar = id, cor = un), where id specifies the group identifier and un sets the unstructured correlation structure. See the group_by argument for more details.

b_formula_gr

Formula for the random effect parameter, b (default ~ 1). Similar to a_formula_gr, user can specify different formulas when fitting univariate_by and multivariate models. The formula can include continuous and/or factor variable(s), including their interactions as covariates (see a_formula_gr for details). In addition to setting up the design matrix for the random effect parameter b, the user can set up the group identifier and the correlation structure for random effects via the vertical bar || approach. For example, consider only an intercept for the random effects a, b, and c specified as a_formula_gr = ~1, b_formula_gr = ~1 and c_formula_gr = ~1. To specify the group identifier (e.g., id) and an unstructured correlation structure, the formula argument can be specified as:
a_formula_gr = ~ (1|i|id)
b_formula_gr = ~ (1|i|id)
c_formula_gr = ~ (1|i|id)
where i within the vertical bars || is just a placeholder. A common identifier (i.e., i) shared across random effect formulas are modeled as unstructured correlated. For more details on the vertical bar approach, please see brms::brm().

c_formula_gr

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

d_formula_gr

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

a_formula_gr_str

Formula for the random effect parameter, a (default NULL), used when fitting a hierarchical model with three or more levels of hierarchy. For example, a model applied to data that includes repeated measurements (level 1) on individuals (level 2), which are further nested within growth studies (level 3).

For a_formula_gr_str argument, only the vertical bar approach (see a_formula_gr) can be used to define the group identifiers and correlation structure. An example of setting up a formula for a three-level model with random effect parameters a, b, and c is as follows:
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 example, |i| and |i2| set up unstructured correlation structures for the random effects at the individual and study levels, respectively. Note that |i| and |i2| must be distinct, as random effect parameters cannot be correlated across different levels of hierarchy.

Additionally, users can specify models with any number of hierarchical levels and include covariates in the random effect formula.

b_formula_gr_str

Formula for the random effect parameter, b (default NULL), used when fitting a hierarchical model with three or more levels of hierarchy. For details, see a_formula_gr_str.

c_formula_gr_str

Formula for the random effect parameter, c (default NULL), used when fitting a hierarchical model with three or more levels of hierarchy. For details, see a_formula_gr_str.

d_formula_gr_str

Formula for the random effect parameter, d (default NULL), used when fitting a hierarchical model with three or more levels of hierarchy. For details, see a_formula_gr_str.

d_adjusted

A logical indicator to adjust the scale of the predictor variable x when fitting the model with the random effect parameter d. The coefficient of parameter d is estimated as a linear function of x, i.e., d * x. If FALSE (default), the original x is used. When d_adjusted = TRUE, x is adjusted for the timing (b) and intensity (c) parameters as x - b) * exp(c) i.e., d * ((x-b)*exp(c)). The adjusted scale of x reflects individual developmental age rather than chronological age. This makes d more sensitive to the timing of puberty in individuals. See sitar::sitar() function for details.

sigma_formula

Formula for the fixed effect distributional parameter, sigma. The sigma_formula sets up the fixed effect design matrix, which may include continuous and/or factor variables (and their interactions) as covariates for the distributional parameter. In other words, setting up the covariates for sigma_formula follows the same approach as for other fixed parameters, such as a (see a_formula for details). Note that sigma_formula estimates the sigma parameter on the log scale. By default, sigma_formula is NULL, as the brms::brm() function itself models sigma as a residual standard deviation (RSD) parameter on the link scale. The sigma_formula, along with the arguments sigma_formula_gr and sigma_formula_gr_str, allows sigma to be estimated as a random effect. The setup for fixed and random effects for sigma is similar to the approach used for other parameters such as a, b, and c.

It is important to note that an alternative way to set up the fixed effect design matrix for the distributional parameter sigma is to use the dpar_formula argument. The advantage of dpar_formula over sigma_formula is that it allows users to specify both linear and nonlinear formulations using the brms::lf() and brms::nlf() syntax. These functions offer more flexibility, such as centering the predictors and enabling or disabling cell mean centering by excluding the intercept via 0 + formulation. However, a disadvantage of the dpar_formula approach is that random effects cannot be included for sigma.

sigma_formula and dpar_formula cannot be specified together. When either sigma_formula or dpar_formula is used, the default estimation of RSD by brms::brm() is automatically turned off.

Users can specify an external function, such as poly, but only with a single argument (the predictor), i.e., poly(age). Additional arguments must be specified externally. For example, to set the degree of the polynomial to 3, a copy of the poly function can be created and modified as follows:
mypoly = poly; formals(mypoly)[['degree']] <- 3; mypoly(age).

sigma_formula_gr

Formula for the random effect parameter, sigma (default NULL). See a_formula_gr for details. Similar to sigma_formula, external functions such as poly can be used. For further details, please refer to the description of sigma_formula.

sigma_formula_gr_str

Formula for the random effect parameter, sigma, when fitting a hierarchical model with three or more levels of hierarchy. See a_formula_gr_str for details. As with sigma_formula, external functions such as poly can be used. For further details, please refer to the description of sigma_formula.

sigma_formula_manual

Formula for the random effect parameter, sigma, provided as a character string that explicitly uses the brms::nlf() and brms::lf() functions (default NULL). An example is:
nlf(sigma ~ z) + lf(z ~ 1 + age + (1 + age |55| gr(id, by = NULL))).

Another use case for sigma_formula_manual is modeling a location-scale model in the SITAR framework, where the same SITAR formula can be used to model the scale (sigma). An example is:

nlf(sigma ~ sigmaSITARFun(logage, 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, sigmaSITARFun (and all other required sub-functions) are defined through the sigmax, sigmadf, sigmaknots, sigmafixed, sigmarandom, sigmaxoffset, sigmaxfun, and sigmabound arguments. Ensure the sigma_formula_manual code matches the sigmaSITARFun function created by these arguments.

Note that for sigma_formula_manual, 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. Additionally, no initial values are defined, so initial values for these parameters should be set to either 0 or random.

sigmax

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

sigmadf

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

sigmaknots

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

sigmafixed

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

sigmarandom

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.

sigmaxoffset

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

sigmabstart

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.

sigmacstart

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.

sigmaxfun

Transformation function for x in the sigma structure. See xfun for details. Ignored if sigma_formula_manual = NULL.

sigmabound

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

dpar_formula

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

autocor_formula

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

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

  • autoregressive (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 structure of residuals.

family

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

custom_family

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.

custom_stanvars

Allows the preparation and 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.

group_arg

Specify arguments 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 subject identifier. If groupvar = NULL (default), groupvar is automatically assigned based on the id argument.

  • dist specifies the distribution from which the random effects are drawn (default gaussian). Currently, gaussian is the only available distribution (as per the brms::brm() documentation).

  • by can be used to estimate a separate variance-covariance structure (i.e., standard deviation and correlation parameters) for random effect parameters (default NULL). If specified, the variable used for by must be a factor variable. For example, by = 'sex' estimates separate variance-covariance structures for males and females.

  • cor specifies the covariance (i.e., correlation) structure for random effect parameters. The default covariance is unstructured (cor = un) for all model types (i.e., univariate, univariate_by, and multivariate). The alternative correlation structure available for univariate and univariate_by models is diagonal, which estimates only the variance parameters (standard deviations), while setting the covariance (correlation) parameters to zero. For multivariate models, options include un, diagonal, and un_s. The un structure models a full unstructured correlation, meaning that the group-level random effects across response variables are drawn from a joint multivariate normal distribution with shared correlation parameters. The cor = diagonal option estimates only variance parameters for each sub-model, while setting the correlation parameters to zero. The cor = un_s option allows for separate estimation of unstructured variance-covariance parameters for each response variable.

Note that it is not necessary to define all or any of these options (groupvar, dist, cor, or by), as they will automatically be set to their default values if unspecified. Additionally, only groupvar from the group_arg argument is passed to the univariate_by and multivariate models, as these models have their own additional options specified via the univariate_by and multivariate arguments. Lastly, the group_arg is ignored when random effects are specified using the vertical bar || approach (see a_formula_gr for details) or when fitting a hierarchical model with three or more levels of hierarchy (see a_formula_gr_str for details).

sigma_group_arg

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

univariate_by

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.

multivariate

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. 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.

a_prior_beta

Specify priors for the fixed effect parameter, a. (default normal(lm, ysd, autoscale = TRUE)). The following key points are applicable for all prior specifications. For full details, see brms::prior():

  • Allowed distributions: normal, student_t, cauchy, lognormal, uniform, exponential, gamma, and inv_gamma (inverse gamma).

  • For each distribution, upper and lower bounds can be set via lb and ub (default NA).

  • Location-scale based distributions (such as normal, student_t, cauchy, and lognormal) have an autoscale option (default FALSE). This option multiplies the scale parameter by a numeric value. While brms typically uses a scaling factor of 1.0 or 2.5, the bsitar package allows any real number to be used (e.g., autoscale = 5.0).

  • For location-scale distributions, fxl (function location) and fxs (function scale) are available to apply transformations to the location and scale parameters. For example, setting normal(2, 5, fxl = 'log', fxs = 'sqrt') translates to normal(log(2), sqrt(5)).

  • fxls (function location scale) transforms both location and scale parameters. The transformation applies when both parameters are involved, as in the log-transformation for normal priors: log_location = log(location / sqrt(scale^2 / location^2 + 1)), log_scale = sqrt(log(scale^2 / location^2 + 1)). This can be specified as a character string or a list of functions.

  • For strictly positive distributions like exponential, gamma, and inv_gamma, the lower bound (lb) is automatically set to zero.

  • For uniform distributions, the option addrange widens the prior range symmetrically. For example, uniform(a, b, addrange = 5) adjusts the range to uniform(a-5, b+5).

  • For exponential distributions, the rate parameter is evaluated as the inverse of the specified value. For instance, exponential(10.0) is internally treated as exponential(1.0 / 10.0) = exponential(0.1).

  • Users do not need to specify each option explicitly, as missing options will automatically default to their respective values. For example, a_prior_beta = normal(location = 5, scale = 1) is equivalent to a_prior_beta = normal(5, 1).

  • For univariate_by and multivariate models, priors can either be the same for all submodels (e.g., a_prior_beta = normal(5, 1)) or different for each submodel (e.g., a_prior_beta = list(normal(5, 1), normal(10, 5))).

  • For location-scale distributions, the location parameter can be specified as the mean (ymean) or median (ymedian) of the response variable, and the scale parameter can be specified as the standard deviation (ysd) or median absolute deviation (ymad). Alternatively, coefficients from a simple linear model can be used (e.g., lm(y ~ age)).

    Example prior specifications include: a_prior_beta = normal(ymean, ysd), a_prior_beta = normal(ymedian, ymad), a_prior_beta = normal(lm, ysd).

    Note that options such as ymean, ymedian, ysd, ymad, and lm are available only for the fixed effect parameter a, not for other parameters like b, c, or d.

b_prior_beta

Specify priors for the fixed effect parameter, b. The default prior is normal(0, 1.5, autoscale = FALSE). For full details on prior specification, please refer to a_prior_beta.

  • Allowed distributions include normal, student_t, cauchy, lognormal, uniform, exponential, gamma, and inv_gamma.

  • You can set upper and lower bounds (lb, ub) as needed (default is NA).

  • The autoscale option controls scaling of the prior’s scale parameter. By default, this is set to FALSE.

  • Further customization and transformations can be applied, similar to the a_prior_beta specification.

c_prior_beta

Specify priors for the fixed effect parameter, c. The default prior is normal(0, 0.5, autoscale = FALSE). For full details on prior specification, please refer to a_prior_beta.

  • Allowed distributions include normal, student_t, cauchy, lognormal, uniform, exponential, gamma, and inv_gamma.

  • Upper and lower bounds (lb, ub) can be set as necessary (default is NA).

  • The autoscale option is also available for scaling the prior's scale parameter (default FALSE).

  • Similar to a_prior_beta, further transformations or customization can be applied.

d_prior_beta

Specify priors for the fixed effect parameter, d. The default prior is normal(0, 1.0, autoscale = FALSE). For full details on prior specification, please refer to a_prior_beta.

  • Allowed distributions include normal, student_t, cauchy, lognormal, uniform, exponential, gamma, and inv_gamma.

  • The option to set upper and lower bounds (lb, ub) is available (default is NA).

  • autoscale allows scaling of the prior’s scale parameter and is FALSE by default.

  • For more advanced transformations or customization, similar to a_prior_beta, these options are available.

s_prior_beta

Specify priors for the fixed effect parameter, s (i.e., spline coefficients). The default prior is normal(0, 'lm', autoscale = TRUE). The general approach is similar to the one described for other fixed effect parameters (see a_prior_beta for details). Key points to note:

  • When using location-scale based priors with 'lm' (e.g., s_prior_beta = normal(lm, 'lm')), the location parameter is set from the spline coefficients obtained from the simple linear model fit, and the scale parameter is based on the standard deviation of the spline design matrix. The location parameter is typically set to 0 (default), and autoscale is set to TRUE.

  • For location-scale based priors, the option sethp (logical, default FALSE) is available to define hierarchical priors. Setting sethp = TRUE alters the prior setup to use hierarchical priors: s ~ normal(0, 'lm') becomes s ~ normal(0, 'hp'), where 'hp' is defined as hp ~ normal(0, 'lm'). The scale for the hierarchical prior is automatically taken from the s parameter, and it can also be defined using the same sethp option. For example, s_prior_beta = normal(0, 'lm', sethp = cauchy) will result in s ~ normal(0, 'lm'), hp ~ cauchy(0, 'lm') .

  • For uniform priors, you can use the option addrange to symmetrically expand the prior range.

It is observed that location-scale based prior distributions (such as normal, student_t, and cauchy) typically work well for spline coefficients.

a_cov_prior_beta

Specify priors for the covariate(s) included in the fixed effect parameter, a (default normal(0, 5.0, autoscale = FALSE)). The approach for specifying priors is similar to a_prior_beta, with a few differences:

  • The options 'ymean', 'ymedian', 'ysd', and 'ymad' are not allowed for a_cov_prior_beta.

  • The 'lm' option for the location parameter allows the covariate coefficient(s) to be obtained from a simple linear model fit to the data. Note that the 'lm' option is only allowed for a_cov_prior_beta and not for covariates in other fixed or random effect parameters.

  • Separate priors can be specified for submodels when fitting univariate_by and a_prior_beta models (see a_prior_beta for details).

b_cov_prior_beta

Specify priors for the covariate(s) included in the fixed effect parameter, b (default normal(0, 1.0, autoscale = FALSE)). See a_cov_prior_beta for details.

c_cov_prior_beta

Specify priors for the covariate(s) included in the fixed effect parameter, c (default normal(0, 0.1, autoscale = FALSE)). See a_cov_prior_beta for details.

d_cov_prior_beta

Specify priors for the covariate(s) included in the fixed effect parameter, d (default normal(0, 1.0, autoscale = FALSE)). See a_cov_prior_beta for details.

s_cov_prior_beta

Specify priors for the covariate(s) included in the fixed effect parameter, s (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.

a_prior_sd

Specify priors 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 submodel (see a_prior_beta).

b_prior_sd

Specify priors for the random effect parameter, b. (default normal(0, 1.0, autoscale = FALSE)). See a_prior_sd for details.

c_prior_sd

Specify priors for the random effect parameter, c. (default normal(0, 0.25, autoscale = FALSE)). See a_prior_sd for details.

d_prior_sd

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

a_cov_prior_sd

Specify priors for the covariate(s) 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.

b_cov_prior_sd

Specify priors for the covariate(s) included in the random effect parameter, b. (default normal(0, 1.0, autoscale = FALSE)). See a_cov_prior_sd for details.

c_cov_prior_sd

Specify priors for the covariate(s) included in the random effect parameter, c. (default normal(0, 0.1, autoscale = FALSE)). See a_cov_prior_sd for details.

d_cov_prior_sd

Specify priors for the covariate(s) included in the random effect parameter, d. (default normal(0, 1.0, autoscale = FALSE)). See a_cov_prior_sd for details.

a_prior_sd_str

Specify priors for the random effect parameter, a, when fitting a hierarchical model with three or more levels of hierarchy. (default NULL). The approach is the same as described for a_prior_sd.

b_prior_sd_str

Specify priors for the random effect parameter, b, when fitting a hierarchical model with three or more levels of hierarchy. (default NULL). The approach is the same as described for a_prior_sd_str.

c_prior_sd_str

Specify priors for the random effect parameter, c, when fitting a hierarchical model with three or more levels of hierarchy. (default NULL). The approach is the same as described for a_prior_sd_str.

d_prior_sd_str

Specify priors for the random effect parameter, d, when fitting a hierarchical model with three or more levels of hierarchy. (default NULL). The approach is the same as described for a_prior_sd_str.

a_cov_prior_sd_str

Specify priors for the covariate(s) included in the random effect parameter, a, when fitting a hierarchical model with three or more levels of hierarchy. (default NULL). The approach is the same as described for a_cov_prior_sd.

b_cov_prior_sd_str

Specify priors for the covariate(s) included in the random effect parameter, b, when fitting a hierarchical model with three or more levels of hierarchy. (default NULL). The approach is the same as described for a_cov_prior_sd_str.

c_cov_prior_sd_str

Specify priors for the covariate(s) included in the random effect parameter, c, when fitting a hierarchical model with three or more levels of hierarchy. (default NULL). The approach is the same as described for a_cov_prior_sd_str.

d_cov_prior_sd_str

Specify priors for the covariate(s) included in the random effect parameter, d, when fitting a hierarchical model with three or more levels of hierarchy. (default NULL). The approach is the same as described for a_cov_prior_sd_str.

sigma_prior_beta

Specify priors for the fixed effect distributional parameter, sigma. (default normal(0, 1.0, autoscale = FALSE)). The approach is similar to that for a_prior_beta.

sigma_cov_prior_beta

Specify priors for the covariate(s) included in the fixed effect distributional parameter, sigma. (default normal(0, 0.5, autoscale = FALSE)). Follows the same approach as a_cov_prior_beta.

sigma_prior_sd

Specify priors for the random effect distributional parameter, sigma. (default normal(0, 0.25, autoscale = FALSE)). Same approach as a_prior_sd.

sigma_cov_prior_sd

Specify priors for the covariate(s) included in the random effect distributional parameter, sigma. (default normal(0, 0.15, autoscale = FALSE)). Follows the same approach as a_cov_prior_sd.

sigma_prior_sd_str

Specify priors for the random effect distributional parameter, sigma, when fitting a hierarchical model with three or more levels of hierarchy. (default NULL). Same approach as a_prior_sd_str.

sigma_cov_prior_sd_str

Specify priors for the covariate(s) included in the random effect distributional parameter, sigma, when fitting a hierarchical model with three or more levels of hierarchy. (default NULL). Follows the same approach as a_cov_prior_sd_str.

rsd_prior_sigma

Specify priors 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 based distributions, user can specify standard deviation (ysd) or the median absolute deviation (ymad) of outcome as the scale parameter. Also, residual standard deviation from the linear mixed model (nlme::lme()) or the linear model (base::lm()) fitted to the data. These are specified as 'lme_rsd' and 'lm_rsd', respectively. Note that if nlme::lme() fails to converge, the option 'lm_rsd' is set automatically. The argument rsd_prior_sigma is evaluated when both dpar_formula and sigma_formula are set to NULL.

dpar_prior_sigma

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

dpar_cov_prior_sigma

Specify priors for the covariate(s) included in the fixed effect distributional parameter sigma. (default normal(0, 1.0, autoscale = FALSE)). Evaluated when sigma_formula is NULL.

autocor_prior_acor

Specify priors for the autocorrelation parameters when fitting a model 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 the unstructured residual correlation structure, use autocor_prior_unstr_acor.

autocor_prior_unstr_acor

Specify priors for the autocorrelation parameters when fitting a model 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 setting up the lkj prior.

gr_prior_cor

Specify priors for the correlation parameter(s) of group-level random effects (default lkj(1)). The only allowed distribution is lkj, specified via a single parameter eta (see brms::prior() for details).

gr_prior_cor_str

Specify priors for the correlation parameter(s) of group-level random effects when fitting a hierarchical model with three or more levels of hierarchy (default lkj(1)). Same as gr_prior_cor.

sigma_prior_cor

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

sigma_prior_cor_str

Specify priors for the correlation parameter(s) of distributional random effects sigma when fitting a hierarchical model with three or more levels of hierarchy (default lkj(1)). Same as sigma_prior_cor.

mvr_prior_rescor

Specify priors for the residual correlation parameter when fitting a multivariate model (default lkj(1)). The only allowed distribution is lkj (see gr_prior_cor for details).

init

Initial values for the sampler. Options include:

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

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

  • 'prior': Initializes parameters based on the specified prior.

  • NULL (default): Initial values are provided by the corresponding init arguments defined below.

init_r

A positive real value that defines the range for the random generation of initial values (default NULL). This argument is used only when init = 'random'.

a_init_beta

Initial values 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 the coefficients from a simple linear model fitted to the data.

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

b_init_beta

Initial values for the fixed effect parameter, b (default 'random'). See a_init_beta for details on available options.

c_init_beta

Initial values for the fixed effect parameter, c (default 'random'). See a_init_beta for details on available options.

d_init_beta

Initial values for the fixed effect parameter, d (default 'random'). See a_init_beta for details on available options.

s_init_beta

Initial values 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 the coefficients from a simple linear model fitted to the data.

a_cov_init_beta

Initial values for the covariate(s) included in the fixed effect parameter, a (default 'random'). Available options include:

  • '0': Initializes the covariates to zero.

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

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

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

Note that the 'lm' option is only available for a_cov_init_beta and not for covariates in other parameters such as b, c, or d.

b_cov_init_beta

Initial values for the covariate(s) included in the fixed effect parameter, b (default 'random'). See a_cov_init_beta for details.

c_cov_init_beta

Initial values for the covariate(s) included in the fixed effect parameter, c (default 'random'). See a_cov_init_beta for details.

d_cov_init_beta

Initial values for the covariate(s) included in the fixed effect parameter, d (default 'random'). See a_cov_init_beta for details.

s_cov_init_beta

Initial values for the covariate(s) included in the fixed effect parameter, s (default 'lm'). See a_cov_init_beta for details. The option 'lm' sets the spline coefficients obtained from a simple linear model fitted to the data. However, note that s_cov_init_beta serves as a placeholder and is not evaluated, as covariates are not allowed for the s parameter. For more details on covariates for s, refer to s_formula.

a_init_sd

Initial value 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': Sets the standard deviation (sd) of the response variable as the initial value.

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

  • 'lme_sd_a': Sets the initial value based on the standard deviation of the random intercept obtained from a linear mixed model (nlme::lme()) fitted to the data. If nlme::lme() fails to converge, the option 'lm_sd_a' will be used automatically.

  • 'lm_sd_a': Sets the square root of the residual variance obtained from a simple linear model applied to the data as the initial value.

Note that the options '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.

Additionally, when fitting univariate_by and multivariate models, the user can set the same initial values for all sub-models, or different initial values for each sub-model.

b_init_sd

Initial value for the standard deviation of the group-level random effect parameter, b (default 'random'). Refer to a_init_sd for available options and details.

c_init_sd

Initial value for the standard deviation of the group-level random effect parameter, c (default 'random'). Refer to a_init_sd for available options and details.

d_init_sd

Initial value for the standard deviation of the group-level random effect parameter, d (default 'random'). Refer to a_init_sd for available options and details.

a_cov_init_sd

Initial values for the covariate(s) included in the random effect parameter a (default 'random'). Available options include:

  • 'random': Random initialization.

  • 'prior': Uses prior distribution values.

b_cov_init_sd

Initial values for the covariate(s) included in the random effect parameter b (default 'random'). Refer to a_cov_init_sd for available options and details.

c_cov_init_sd

Initial values for the covariate(s) included in the random effect parameter c (default 'random'). Refer to a_cov_init_sd for available options and details.

d_cov_init_sd

Initial values for the covariate(s) included in the random effect parameter d (default 'random'). Refer to a_cov_init_sd for available options and details.

sigma_init_beta

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

  • 'random': Random initialization.

  • 'prior': Uses prior distribution values.

sigma_cov_init_beta

Initial values for the covariate(s) included in the fixed effect distributional parameter sigma (default 'random'). Refer to sigma_init_beta for available options and details.

sigma_init_sd

Initial value for the standard deviation of the distributional random effect parameter sigma (default 'random'). The approach is the same as described earlier for the group-level random effect parameters such as a (See a_init_sd for details).

sigma_cov_init_sd

Initial values for the covariate(s) included in the distributional random effect parameter sigma (default 'random'). The approach is the same as described for a_cov_init_sd (See a_cov_init_sd for details).

gr_init_cor

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

  • 'random': Random initialization.

  • 'prior': Uses prior distribution values.

sigma_init_cor

Initial values for the correlation parameters of distributional random effects parameter sigma (default 'random'). Allowed options are:

  • 'random': Random initialization.

  • 'prior': Uses prior distribution values.

rsd_init_sigma

Initial values for the residual standard deviation parameter, sigma (default 'random'). Options available are:

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

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

  • 'prior': Initializes the residual standard deviation based on prior distribution values.

  • 'lme_rsd': Sets the initial value based on the standard deviation of residuals obtained from the linear mixed model (nlme::lme()) fitted to the data.

  • 'lm_rsd': Sets the initial value as the square root of the residual variance from the simple linear model fitted to the data.

Note that if nlme::lme() fails to converge, the option 'lm_rsd' is set automatically. The argument rsd_init_sigma is evaluated when both dpar_formula and sigma_formula are set to NULL.

dpar_init_sigma

Initial values for the distributional parameter sigma (default 'random'). The approach and available options are the same as described for rsd_init_sigma. This argument is evaluated only when dpar_formula is not NULL.

dpar_cov_init_sigma

Initial values for the covariate(s) included in the distributional parameter sigma (default 'random'). Allowed options are '0', 'random', and 'prior'.

autocor_init_acor

Initial values for the autocorrelation parameter (see autocor_formula for details). Allowed options are '0', 'random', and 'prior' (default 'random').

autocor_init_unstr_acor

Initial values for unstructured residual autocorrelation parameters (default 'random'). Allowed options are '0', 'random', and 'prior'. The approach for setting initials for autocor_init_unstr_acor is the same as for gr_init_cor.

mvr_init_rescor

Initial values for the residual correlation parameter when fitting a multivariate model (default 'random'). Allowed options are '0', 'random', and 'prior'.

r_init_z

Initial values for the standardized group-level random effect parameters (default 'random'). These parameters are part of the Non-Centered Parameterization (NCP) approach used in the brms::brm().

vcov_init_0

A logical (default FALSE) to set initial values for variance (standard deviation) and covariance (correlation) parameters to zero. This allows for setting custom initial values for the fixed effects parameters while keeping the variance-covariance parameters at zero.

jitter_init_beta

A proportion (between 0 and 1) to perturb the initial values for fixed effect parameters. The default is NULL, which means that the same initial values are used across all chains. A sensible option might be jitter_init_beta = 0.1, which mildly perturbs the initial values. Note that the jitter is applied as a proportion of the specified initial value, not an absolute amount. For example, if the initial value is 100, setting jitter_init_beta = 0.1 means the perturbed initial value will be within the range 90 to 110. Conversely, if the initial value is 10, the perturbed value will fall within the range 9 to 11.

jitter_init_sd

A proportion (between 0 and 1) to perturb the initial values for the standard deviation of random effect parameters. The default is NULL, which means the same initial values are used across all chains. A reasonable option might be jitter_init_sd = 0.01, which was found to work well during early testing.

jitter_init_cor

A proportion (between 0 and 1) to perturb the initial values for the correlation parameters of random effects. The default is NULL, which means the same initial values are used across all chains. An option of setting jitter_init_cor = 0.001 was found to be effective during early testing.

prior_data

An optional argument (a named list, default NULL) that can be used to pass information to the prior arguments for each parameter (e.g., a_prior_beta). The prior_data is particularly helpful when passing a long vector or matrix as priors. These vectors and matrices can be created in the R framework and then passed using the prior_data. For example, to pass a vector of location and scale parameters when setting priors for covariate coefficients (with 10 dummy variables) included in the fixed effects parameter a, the following steps can be used:

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

  • Specify these objects in the prior_data list: prior_data = list(prior_a_cov_location = prior_a_cov_location, prior_a_cov_scale = prior_a_cov_scale).

  • Use the prior_data objects to set up the priors: a_cov_prior_beta = normal(prior_a_cov_location, prior_a_cov_scale).

init_data

An optional argument (a named list, default NULL) that can be used to pass information to the initial arguments. The approach is identical to how prior_data is handled (as described above).

init_custom

Specify a custom initialization object (a named list). The named list is directly passed to the init argument without verifying the dimensions or name matching. If initial values are set for some parameters via parameter-specific arguments (e.g., a_init_beta = 0), init_custom will only be passed to those parameters that do not have initialized values. To override this behavior and use all of init_custom values regardless of parameter-specific initials, set init = 'custom'.

verbose

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.

expose_function

An optional argument (logical, default FALSE) to indicate whether to expose the Stan function used in model fitting.

get_stancode

An optional argument (logical, default FALSE) to retrieve the Stan code (see [brms::stancode()] for details).

get_standata

An optional argument (logical, default FALSE) to retrieve the Stan data (see [brms::standata()] for details).

get_formula

An optional argument (logical, default FALSE) to retrieve the model formula (see [brms::brmsformula()] for details).

get_stanvars

An optional argument (logical, default FALSE) to retrieve the Stan variables (see [brms::stanvar()] for details).

get_priors

An optional argument (logical, default FALSE) to retrieve the priors (see [brms::get_prior()] for details).

get_priors_eval

An optional argument (logical, default FALSE) to retrieve the priors specified by the user.

get_init_eval

An optional argument (logical, default FALSE) to retrieve the initial values specified by the user.

validate_priors

An optional argument (logical, default FALSE) to validate the specified priors (see [brms::validate_prior()] for details).

set_self_priors

An optional argument (default NULL) to manually specify the priors. set_self_priors is passed directly to [brms::brm()] without performing any checks.

add_self_priors

An optional argument (default NULL) to append part of the prior object. This is for internal use only.

set_replace_priors

An optional argument (default NULL) to replace part of the prior object. This is for internal use only.

set_same_priors_hierarchy

An optional argument (default NULL) to replace part of the prior object. This is for internal use only.

outliers

An optional argument (default NULL) to remove outliers. This should be a named list passed directly to [sitar::velout()] and [sitar::zapvelout()] functions. This is for internal use only.

unused

An optional formula defining variables that are unused in the model but should still be stored in the model's data frame. Useful when variables are needed during post-processing.

chains

The number of Markov chains (default 4).

iter

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

warmup

A positive integer specifying the number of warmup (aka burn-in) iterations. This also specifies the number of iterations used for stepsize adaptation, so warmup draws should not be used for inference. The number of warmup iterations should not exceed iter, and the default is iter/2.

thin

A positive integer specifying the thinning interval. Set thin > 1 to save memory and computation time if iter is large. Thinning is often used in cases with high autocorrelation of MCMC draws. An indication of high autocorrelation is poor mixing of chains (i.e., high rhat values) despite the model recovering parameters well. A useful diagnostic to check for autocorrelation of MCMC draws is the mcmc_acf function from the bayesplot package.

cores

Number of cores to be used when executing the chains in parallel. See brms::brm() for details. Unlike brms::brm(), which defaults the cores argument to cores=getOption("mc.cores", 1), the default cores in the bsitar package is cores=getOption("mc.cores", 'optimize'), which optimizes the utilization of system resources. The maximum number of cores that can be deployed is calculated as the maximum number of available cores minus 1. When the number of available cores exceeds the number of chains (see chains), then the number of cores is set equal to the number of chains.

Another option is to set cores as getOption("mc.cores", 'maximise'), which sets the number of cores to the maximum number of cores available on the system regardless of the number of chains specified. Alternatively, the user can specify cores in the same way as brms::brm() with getOption("mc.cores", 1).

These options can be set globally using options(mc.cores = x), where x can be 'optimize', 'maximise', or 1. The cores argument can also be directly specified as an integer (e.g., cores = 4).

backend

A character string specifying the package to be used when executing the Stan model. The available options are "rstan" (the default) or "cmdstanr". The backend can also be set globally for the current R session using the "brms.backend" option. See brms::brm() for more details.

threads

Number of threads to be used in within-chain parallelization. Note that unlike the brms::brm() which sets the threads argument as getOption("brms.threads", NULL) implying that no within-chain parallelization is used by default, the bsitar package, by default, sets threads as getOption("brms.threads", 'optimize') to utilize the available resources from the modern computing systems. The number of threads per chain is set as the maximum number of cores available minus 1. Another option is to set threads as getOption("brms.threads", 'maximise') which set the number threads per chains same as the maximum number of cores available. User can also set the threads similar to the brms i.e., getOption("brms.threads", NULL). All these three options can be set globally as options(brms.threads = x) where x can be 'optimize', 'maximise' or NULL. Alternatively, the number of threads can be set directly as threads = threading(x) where X is an integer. Other arguments that can be passed to the threads are grainsize and the static. See brms::brm() for further details on within-chain parallelization.

opencl

The platform and device IDs of the OpenCL device to use for GPU support during model fitting. If you are unsure about the IDs of your OpenCL device, c(0,0) is typically the default that should work. For more details on how to find the correct platform and device IDs, refer to brms::opencl(). This parameter can also be set globally for the current R session using the "brms.opencl" option.

normalize

Logical flag indicating whether normalization constants should be included in the Stan code (default is TRUE). If set to FALSE, normalization constants are omitted, which may increase sampling efficiency. However, this requires Stan version >= 2.25. Note that setting normalize = FALSE will disable some post-processing functions, such as brms::bridge_sampler(). This option can be controlled globally via the brms.normalize option.

algorithm

A character string specifying the estimation method to use. Available options are:

  • "sampling" (default): Markov Chain Monte Carlo (MCMC) method.

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

  • "fullrank": Variational inference with a multivariate normal distribution.

  • "fixed_param": Sampling from fixed parameter values.

This parameter can be set globally via the "brms.algorithm" option (see options for more details).

control

A named list to control the sampler's behavior. The default settings are the same as those in brms::brm(), with one exception: the max_treedepth has been increased from 10 to 12 to better explore the typically challenging posterior geometry in nonlinear models. However, the adapt_delta, which is often increased for nonlinear models, retains its default value of 0.8 to avoid unnecessarily increasing sampling time. For full details on control parameters and their default values, refer to brms::brm().

empty

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.

rename

For internal use only.

pathfinder_args

A named list of arguments passed to the 'pathfinder' algorithm. This is used to set 'pathfinder'-based initial values for the 'MCMC' sampling. Note that 'pathfinder_args' currently only works when backend = "cmdstanr". If pathfinder_args is not NULL and the user specifies backend = "rstan", the backend will automatically be changed to cmdstanr.

pathfinder_init

A logical value (default FALSE) indicating whether to use initial values from the 'pathfinder' algorithm when fitting the final model (i.e., 'MCMC' sampling). Note that 'pathfinder_args' currently works only when backend = "cmdstanr". If pathfinder_args is not NULL and the user specifies backend = "rstan", the backend will automatically switch to cmdstanr. The arguments passed to the 'pathfinder' algorithm are specified via 'pathfinder_args'; if 'pathfinder_args' is NULL, the default arguments from 'cmdstanr' will be used.

sample_prior

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 you to generate draws from the prior predictive distribution. In this case, all parameters must have proper priors.

save_pars

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.

drop_unused_levels

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

stan_model_args

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.

refresh

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.

silent

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", you can prevent the opening of additional progress bars by setting open_progress = FALSE.

seed

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.

save_model

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.

fit

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.

file

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, you must manually remove the file or specify the file_refit argument. The file name is stored within the brmsfit object for later use.

file_compress

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.

file_refit

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 you believe a false positive occurred, you can use [brms::brmsfit_needs_refit()] 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.

future

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.

parameterization

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().

...

Further arguments passed to brms::brm(). This can 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, match_sitar_d_form, sigmamatch_sitar_a_form, displayit, setcolh, setcolb.

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.

Details

The SITAR is a shape-invariant nonlinear mixed-effects growth curve model that fits a population average (i.e., mean) curve to the data and aligns each individual's growth trajectory to the underlying population curve via a set of (typically) three random effects: size, timing, and intensity. Additionally, a slope parameter can be included as a random effect to estimate the variability in adult growth rate (see sitar::sitar() for details).

The concept of a shape-invariant model (SIM) was first introduced by Lindstrom (1995), and later used by Beath (2007) to model infant growth data (birth to 2 years). The current version of the SITAR model was developed by Cole et al. (2010) and has been extensively used for modeling growth data (see Nembidzane et al. 2020 and Sandhu 2020).

The frequentist version of the SITAR model can be fit using the already available R package, sitar (Cole 2022). The framework of the Bayesian implementation of the SITAR model in the bsitar package is similar to the sitar package, with the main difference being that sitar uses splines::ns() to construct the B-splines based natural cubic spline design matrix, whereas bsitar implements a different strategy to create natural cubic splines. The bsitar offers three different types of splines: nsp, nsk, and rcs. Both nsp and nsk use the B-splines basis to generate the natural cubic spline design matrix as implemented in splines2::nsp() and splines2::nsk(), whereas rcs is based on the truncated power basis approach (see Harrell and others (2001) and Harrell Jr. (2022) for details) to construct the spline design matrix. While all approaches produce the same growth curves, the model-estimated spline coefficients differ from each other.

Like the sitar package (Cole et al. 2010), the bsitar package fits the SITAR model with (usually) three random effects: size (parameter a), timing (parameter b), and intensity (parameter c). Additionally, there is a slope parameter (parameter d) that models the variability in the adult slope of the growth curve (see sitar::sitar() for details).

Note that the author of the sitar package (Cole et al. 2010) enforces the inclusion of the d parameter as a random effect only, excluding it from the fixed structure of the model. However, the bsitar package allows inclusion of the d parameter in both the fixed and/or random effects structures of the SITAR model.

For the three-parameter version of the SITAR model (default), the fixed effects structure (i.e., population average trajectory) is specified as fixed = 'a+b+c', and the random effects structure, capturing the deviation of individual trajectories from the population average curve, is specified as random = 'a+b+c'.

The bsitar package offers flexibility in model specification. For example:

The sitar package internally depends on the brms package (see Bürkner 2022; Bürkner 2021), which fits a wide range of hierarchical linear and nonlinear regression models, including multivariate models. The brms package itself depends on Stan for full Bayesian inference (see Stan Development Team 2023; Gelman et al. 2015). Like brms, the bsitar package allows flexible prior specifications based on user's knowledge of growth processes (e.g., timing and intensity of growth spurts).

The brms package uses a combination of normal and student_t distributions for regression coefficients, group-level random effects, and the distributional parameter (sigma), while rstanarm uses normal distributions for regression coefficients and group-level random effects, but sets exponential for the distributional parameter (sigma). By default, bsitar uses normal distributions for all parameters, including regression coefficients, standard deviations of group-level random effects, and the distributional parameter. Additionally, bsitar provides flexibility in choosing scale parameters for location-scale distributions (such as normal and student_t).

The bsitar package also allows three types of model specifications: 'univariate', 'univariate_by', and 'multivariate':

The bsitar package offers full flexibility in specifying predictors, degrees of freedom for design matrices, priors, and initial values. The package also allows users to specify options in a user-friendly manner (e.g., univariate_by = sex is equivalent to univariate_by = 'sex').

Value

An object of class brmsfit, bsitar, which contains the posterior draws, model coefficients, and other useful information related to the model fitting. This object includes details such as the fitted model, the data used, prior distributions, and any other relevant outputs from the Stan model fitting process. The resulting object can be used for further analysis, diagnostics, and post-processing, including model summary statistics, predictions, and visualizations.

Note

The package is under continuous development, and new models, post-processing features, and improvements are being actively worked on. Keep an eye on future releases for additional functionality and updates to enhance model fitting, diagnostics, and analysis capabilities.

Author(s)

Satpal Sandhu satpal.sandhu@bristol.ac.uk

References

Beath KJ (2007). “Infant growth modelling using a shape invariant model with random effects.” Statistics in Medicine, 26(12), 2547–2564. doi:10.1002/sim.2718, Type: Journal article.

Bürkner P (2021). “Bayesian Item Response Modeling in R with brms and Stan.” Journal of Statistical Software, 100(5), 1–54. doi:10.18637/jss.v100.i05.

Bürkner P (2022). brms: Bayesian Regression Models using Stan. R package version 2.18.0, https://CRAN.R-project.org/package=brms.

Cole T (2022). sitar: Super Imposition by Translation and Rotation Growth Curve Analysis. R package version 1.3.0, https://CRAN.R-project.org/package=sitar.

Cole TJ, Donaldson MDC, Ben-Shlomo Y (2010). “SITAR—a useful instrument for growth curve analysis.” International Journal of Epidemiology, 39(6), 1558–1566. ISSN 0300-5771, doi:10.1093/ije/dyq115, tex.eprint: https://academic.oup.com/ije/article-pdf/39/6/1558/18480886/dyq115.pdf.

Gelman A, Lee D, Guo J (2015). “Stan: A Probabilistic Programming Language for Bayesian Inference and Optimization.” Journal of Educational and Behavioral Statistics, 40(5), 530-543. doi:10.3102/1076998615606113.

Harrell FE, others (2001). Regression modeling strategies: with applications to linear models, logistic regression, and survival analysis, volume 608. Springer.

Harrell Jr. FE (2022). Hmisc: Harrell Miscellaneous. R package version 4.7-2, https://hbiostat.org/R/Hmisc/.

Lindstrom MJ (1995). “Self-modelling with random shift and scale parameters and a free-knot spline shape function.” Statistics in Medicine, 14(18), 2009-2021. doi:10.1002/sim.4780141807, https://pubmed.ncbi.nlm.nih.gov/8677401/.

Nembidzane C, Lesaoana M, Monyeki KD, Boateng A, Makgae PJ (2020). “Using the SITAR Method to Estimate Age at Peak Height Velocity of Children in Rural South Africa: Ellisras Longitudinal Study.” Children, 7(3), 17. ISSN 2227-9067, doi:10.3390/children7030017, https://www.mdpi.com/2227-9067/7/3/17.

Sandhu SS (2020). Analysis of longitudinal jaw growth data to study sex differences in timing and intensity of the adolescent growth spurt for normal growth and skeletal discrepancies. Thesis, University of Bristol.

Stan Development Team (2023). Stan Reference Manual version 2.31. https://mc-stan.org/docs/reference-manual/.

See Also

brms::brm() brms::brmsformula() brms::prior()

Examples




# Below, we fit a SITAR model to a subset of the Berkley height data, 
# specifically the data for 70 girls between the ages of 8 and 18.  
# This subset is used as an example in the vignette for the 'sitar' package.
#
# The original Berkley height data contains repeated growth measurements for
# 66 boys and 70 girls (ages 0-21). For this example, we use a subset of the 
# data for 70 girls aged 8 to 18 years.
#
# For details on the full Berkley height dataset, refer to 'sitar' package
# documentation (help file: ?sitar::berkeley). Further details on the subset
# of the data used here can be found in the vignette ('Fitting_models_with_SITAR', 
# package = 'sitar').

# Load the 'berkeley_exdata' that has been pre-saved
berkeley_exdata <- getNsObject(berkeley_exdata)

# Fit frequentist SITAR model with df = 3 using the sitar package 

model_ml <- sitar::sitar(x = age, y = height, id = id, 
                          df = 3, 
                          data = berkeley_exdata, 
                          xoffset = 'mean',
                          fixed = 'a+b+c', 
                          random = 'a+b+c',
                          a.formula = ~1, 
                          b.formula = ~1, 
                          c.formula = ~1
                          )


# Fit Bayesian SITAR model 

# To avoid time-consuming model estimation, the Bayesian SITAR model fit has 
# been saved as an example fit ('berkeley_exfit'). This model was fit using 
# 2 chains (2000 iterations per chain) with thinning set to 5 for memory  
# efficiency. Users are encouraged to refit the model using default settings 
# (4 chains, 2000 iterations per chain, thin = 1) as suggested by the Stan team.
# Note that with thinning set to 5 (thin = 5), only one fifth of total draws 
# will be saved and hence the effective sample size is expected to be small.

# Check if the pre-saved model 'berkeley_exfit' exists
# berkeley_exfit <- bsitar:::berkeley_exfit

berkeley_exfit <- getNsObject(berkeley_exfit)
 
if(exists('berkeley_exfit')) {
  model <- berkeley_exfit
} else {
  # Fit model with default priors
  # Refer to the documentation for prior on each parameter
  model <- bsitar(x = age, y = height, id = id, 
                  df = 3, 
                  data = berkeley_exdata,
                  xoffset = 'mean', 
                  fixed = 'a+b+c', 
                  random = 'a+b+c',
                  a_formula = ~1, 
                  b_formula = ~1, 
                  c_formula = ~1, 
                  threads = brms::threading(NULL),
                  chains = 2, cores = 2, iter = 2000, thin = 5)
                  
}

# Generate model summary
summary(model)

# Compare model summary with the frequentist SITAR model
print(model_ml)

# Check model fit via posterior predictive checks using plot_ppc.
# This function is based on pp_check from the 'brms' package.
plot_ppc(model, ndraws = 100)

# Plot distance and velocity curves using plot_conditional_effects.
# This function works like conditional_effects from the 'brms' package,
# with the added option to plot velocity curves.

# Distance curve
plot_conditional_effects(model, deriv = 0)

# Velocity curve
plot_conditional_effects(model, deriv = 1)

# Plot distance and velocity curves along with parameter estimates using 
# plot_curves (similar to plot.sitar from the sitar package).
plot_curves(model, apv = TRUE)

# Compare plots with the frequentist SITAR model
plot(model_ml)



[Package bsitar version 0.3.2 Index]