Added plot_diagnostics(), a new convenience function for generating diagnostic plots from fitted bsitar models. The function provides a unified interface for residual diagnostics, MCMC diagnostics, and posterior predictive checks, including residual-versus-fitted plots, residual Q-Q plots, trace plots, autocorrelation plots, density overlays, R-hat, effective sample size diagnostics, and posterior predictive summaries. It also supports different fitted-value scales through set_draws = "epred", "linpred", or "prediction", along with optional arguments for multivariate responses and distributional parameters such as resp, re_formula, category, and dpar. When multiple compatible plots are requested, the function can combine them automatically into a single patchwork display while still returning the individual plot objects for further customization.
Added compare_models(), a new helper function for comparing one or more fitted models using common Bayesian model fit criteria such as LOO and WAIC. The function provides a convenient interface for evaluating relative model performance across candidate models, making it easier to assess expected out-of-sample predictive fit and support model selection workflows within the package. By allowing users to compare multiple models in a single step, compare_models() streamlines routine model comparison and improves consistency in reporting fit-based evaluation results.
The functions marginal_draws(), marginal_comparison(), and growthparameters_comparison() have been renamed to get_predictions(), get_comparisons(), and get_growthparameters(), respectively, to better reflect their roles and to harmonise the naming scheme across the package. In particular, the earlier names with the marginal_ prefix suggested that these functions are used only for marginal inference, whereas they in fact support both marginal and conditional inferences. The functions with the marginal_* prefix are now soft-deprecated and will be removed in the next release. Therefore, it is recommended that users update their code to use the correspondingget_* prefix functions instead.
Changed the default behavior of random initial values so that Stan’s randomly generated initial values are now controlled by init_r. When init_r = NULL, Stan draws initial values uniformly from −2,2 on the unconstrained scale. When init_r is specified (for example, init_r = 0.5, the default for init = "random"), initial values are drawn uniformly from −0.5,0.5 on the unconstrained scale. Users should note that this change in initialization may affect convergence, model fit, and posterior estimates, especially in sensitive or weakly identified models, or when adaptation is difficult, when using init = "random". To recover the previous random-initialization behavior, set init_r = NULL. The default setting remains init = NULL, as before. In this case, user specified initial values are assigned to the fixed-effect parameters (a, b, and c, as well as the spline coefficients), while random initial values are assigned to the variance-covariance parameters, drawn uniformly from −2,2 on the unconstrained scale.
Changed default priors. All model parameters — including population-average effects, standard deviations of individual-level random effects, and residual standard deviation — are now assigned normal()-based priors. The vignette Bayesian_SITAR_model_An_Introduction has been updated to reflect the default priors for each parameter. Users should be aware that this change may affect convergence, model fit, and posterior estimates, particularly in sensitive or weakly identified models. It is strongly recommended that users carefully review and set priors appropriately for their specific model.
A major refactoring of the internal code to streamline the integration of modelling the distributional parameter sigma, aligning it more closely with the location parameter mu. This update marks the first step in a series of improvements aimed at enabling robust modelling of location-scale models. Additionally, the core functions are significantly rewritten to enhance the speed and efficiency of model fitting. Note: This release should be considered experimental. Future versions may introduce further significant changes. Additionally, documentation has not yet been fully updated to reflect these modifications. Due to the above changes, previously saved model objects will need to be refitted for post-processing to work correctly. No other changes are required from the user’s perspective—only that the model should be re-run.
Support for explicit modelling of scale parameter, sigma (distributional model)
Like location parameter mu, the scale parameter sigma, which models the residual variance structure, can be estimated as a function of covariates to accommodate heteroskedasticity. This flexibility allows researchers to specify complex variance relationships that better reflect the underlying data-generating process. The updates include support for prior specifications setting custom initial values. Note that the post-processing support remains limited and experimental The documentation provides comprehensive details on modeling the distributional parameter sigma. Briefly, the available modeling options are:
Constant Variance: By default, sigma is treated as a constant parameter across all observations, assuming homoskedastic residuals.
Covariate-Dependent Variance: The distributional parameter can be modeled as a function of predictors using the formula syntax sigma ~ covariates. This approach enables heteroskedastic modeling where residual variance changes systematically with covariates.
Variance Function Specifications: Common variance structures include exponential relationships (sigma ~ exp(age)), power functions (sigma ~ pow(abs(fitted_values), power)), and linear combinations of multiple predictors. The bsitar package provide six different types of variance function that are adapted from the nlme package.
Random Effects Integration
At higher hierarchical levels, sigma can incorporate random effects to model between-group variance heterogeneity. This is particularly valuable in multilevel growth models where different clusters exhibit varying degrees of residual scatter.
Prior Specification
The bsitar allows full control on specifying appropriate prior distributions for sigma parameters typically include half-normal, half-Cauchy, or exponential priors for group level standard deviation parameters that respect the positive constraint while providing regularization for model stability.
Added support for estimating model-based individual growth parameters, such as age at peak growth velocity (APGV), peak growth velocity (PGV) and size at peak growth velocity (SPGV). For details, refer to the function modelbased_growthparameters(). Note multivariate models are not supported yet.
Support has been added to allow the use of external functions, such as splines::ns(), for modeling the distributional parameter sigma. Users can now specify functions for both the fixed sigma_formula and random sigma_formula_gr effects. Further, different functions can be used in the fixed and random effects formulas. For example. For example. sigma_formula = ~ 1 + splines::ns(age, df = 3), and sigma_formula_gr = ~ 1 + stats::poly(age, degree = 2)
Added support for tag feature implementing parameter specific prior sensitivity analysis in priorsense package via brms.
An new feature that allows for optimizing the number and /or the placements of knots via the argumentknots_selection argument. This feature explores an extended candidate set (usually df + 4) to select an optimal subset of knots minimizing a specified information criterion such as AIC, BIC or the residual error from Cross-validation procedure. It supports flexible strategies via select (knots, degrees of freedom, or both), with options for returning and plotting all scores during df selection. Diagnostic outputs and plotting types can be selected, with options to return or print these results. This feature enhances model flexibility and diagnostic capability for improved spline regression performance.
Added hypothesis_test() that provides comprehensive hypothesis testing framework for the Bayesian SITAR models. The hypothesis_test() not only estimates the contrasts but also support testing for the ROPE (Region of Practical Equivalence) and pd (probability of direction) results in a unified interface.
The efficiency of the post-processing function has been improved (average improvement in speed ~ 2x).
The default stype is rcs
The function block for the multivariate model did not render properly.
The random initial values for the residual correlation parameter of the multivariate model resulted in an empty list warning.
The bsitar now supports three different types of splines: 'rcs', 'nsp', and 'nsk'.
While 'rcs' constructs the spline design matrix using the truncated power basis, both 'nsp'
and 'nsk' implement a B-spline-based natural cubic spline basis. The truncated power basis method,
often referred to as Harrell's method, is implemented in the rcspline.eval() function of the
Hmisc package. The B-spline-based implementations of 'nsp' and 'nsk' are the same as
those described in the splines2 package. Previously, only 'rcs' was available. Now, the
default method is 'nsp'.
The bsitar package allows for fitting the SITAR model with a four-parameter formulation (a + b + c + d),
where a fourth parameter, d, is added. The parameter d, the age slope, allows the adult part of the growth
curve to vary in slope. The age slope d represents the regression coefficient of y on x.
Like the sitar package, the bsitar offers two parameterizations: one in which x is adjusted for the random
effects b and c i.e., (x-b)*exp(c), and the other in which the unadjusted x is used. This is controlled
by the argument d_adjusted provided in the bsitar::bsitar() function. When d_adjusted = TRUE, the first
version is applied (x adjusted for random effects b and c), which means that individual developmental age,
rather than chronological age, is used in the slope regression. This makes parameter d more sensitive to the timing
of puberty in individuals. When d_adjusted = FALSE (the default), the unadjusted x is used. The default choice,
d_adjusted = FALSE, is primarily to match the behavior of the sitar package, which sets d.adjusted = FALSE.
Added support for computing and comparing growth curves using the marginaleffects package as
the back-end (see marginal_draws(), marginal_comparison(), and growthparameters_comparison()).
This enables the use of the computational flexibility offered by the marginaleffects package to
estimate various quantities of interest, such as adjusted growth curves (distance and velocity), and
growth parameters like age at peak growth velocity. All three functions support parallel computation via
the future and doFuture packages.
The optimize_model() function now allows users to specify custom functions in optimize_x and optimize_y
when optimizing the Bayesian SITAR model. For example, it is now possible to use
optimize_x = list(function(x) log(x + 3/4)). Thanks to Tim Cole for suggesting this feature.
This update greatly enhances the flexibility of optimize_model() and enables users to search for a
range of optimal x and y transformations.
An experimental feature has been added to use $pathfinder() based initial values for the MCMC sampling
$sample() (via the argument pathfinder_init = TRUE, default is FALSE). The arguments for
$pathfinder() can be specified as a named list using pathfinder_args. Note that this feature is
only available when backend = 'cmdstanr'.
Added a new vignette comparing growth curves and growth parameters obtained from the frequentist
(sitar package) and Bayesian (bsitar package) versions of the SITAR model applied to the
height data.
The prior distribution for each parameter has been changed from student_t() to normal().
The prior distribution for all parameters, including regression coefficients as well as the standard
deviation (sd) for the group-level random effects and the distributional parameter (sigma), has been
changed to normal(). Previously, the distribution for regression coefficients and the sd for the
group-level random effects was student_t(), while the distribution for the sd of the distributional
parameter (sigma) was exponential(). Note that the same location and scale parameters used earlier
for student_t() are now used for the normal() distribution. For example, prior specified
earlier as student_t(3, 0, 15) (3 degree of freedom, location 0 and scale) has been revised as
normal(0, 15).
The default setting for initial values is now random except for the population average parameters:
size (a_init_beta), timing (b_init_beta), intensity (c_init_beta), and spline coefficients
(s_init_beta). For size and spline coefficient parameters, the initial values are derived from the
linear regression fit and specified as a_init_beta = lm and s_init_beta = lm. The initial
values for both timing and intensity parameters are set to '0', i.e., b_init_beta = 0 and
c_init_beta = 0.
Improved documentation
The sigma_cov_init_beta = random argument was setting incorrect initial values for the covariates
included in the sigma formula. The initial values for the Intercept (sigma_init_beta) were
incorrectly applied to the covariates as well.
Users no longer need to set the environment to globalenv(), i.e., envir = globalenv(), for
post-processing functions to work properly. The environment is now automatically set to match the
environment of the exposed functions. It is important to note that manually setting the environment
(via the envir argument) may result in errors. The envir argument is now primarily for internal
use, which is needed during tests.
Minor corrections and changes have been made to improve the efficiency of the R code.
add_model_criterion() function to compute fit criteria such as "loo", "waic", "kfold",
"loo_subsample", "bayes_R2" (Bayesian R-squared), "loo_R2" (LOO-adjusted R-squared), and "marglik"
(log marginal likelihood). The computed fit criteria are added to the model object for later use including
comparison of models. The add_model_criterion() is a wrapper around the add_criterion()
function available from the brms packagebsitar::bsitar()```received options ```file```, file_refit, and ``file_compress to save and retreive fitted objects. See brms::brm help file for details.