Package {BayesRTMB}

BayesRTMB: Bayesian Modeling with RTMB

Description

The 'BayesRTMB' package provides an intuitive, R-based Bayesian modeling environment using the 'RTMB' package as its computational backend. It allows users to perform Maximum A Posteriori (MAP) estimation, Markov Chain Monte Carlo (MCMC) sampling, and Automatic Differentiation Variational Inference (ADVI) for hierarchical and complex non-linear models through a unified and user-friendly interface.

Details

The standard workflow of the package consists of three main steps:

1. **Model Definition**: Use [rtmb_code()] to define the model structure (likelihood, prior distributions, transformations, etc.) using standard R syntax. 2. **Model Instantiation**: Pass the data and the generated model code to [rtmb_model()] to create an 'RTMB_Model' R6 object. 3. **Estimation**: Call the appropriate method on the model object: * '$optimize()': For MAP estimation and Laplace approximation. * '$sample()': For MCMC sampling using the No-U-Turn Sampler (NUTS). * '$variational()': For approximate Bayesian inference using ADVI.

For more detailed examples and tutorials, please refer to the package vignettes.

Author(s)

Maintainer: Hiroshi Shimizu simizu706@gmail.com

See Also

Useful links: * [https://github.com/norimune/BayesRTMB](https://github.com/norimune/BayesRTMB) * Report bugs at [https://github.com/norimune/BayesRTMB/issues](https://github.com/norimune/BayesRTMB/issues)

Automatic Differentiation Variational Inference (ADVI)

Description

Automatic Differentiation Variational Inference (ADVI)

Usage

ADVI_method(
  model,
  par_list,
  pl_full,
  iter = 3000,
  tol_rel_obj = 0.001,
  window_size = 100,
  num_samples = 1000,
  alpha = 0.01,
  laplace = FALSE,
  print_freq = 500,
  method = c("meanfield", "fullrank", "hybrid"),
  update_progress = NULL,
  update_interval = 100
)

Arguments

Value

A list containing 'fit', 'random_fit', 'elbo_history', 'elbo_final', 'rel_obj_final', and 'converged'.

Big Five Personality Traits Data

Description

A dataset containing responses to 20 items measuring the Big Five personality traits. This dataset was originally collected by the package author. The items are arranged in a cyclical order representing the five factors: Extraversion, Neuroticism, Openness, Agreeableness, and Conscientiousness.

Usage

BigFive

Format

A data frame with 10 rows and 20 columns:

BF1

Extraversion item 1

BF2

Neuroticism item 1

BF3

Openness item 1

BF4

Agreeableness item 1

BF5

Conscientiousness item 1

BF6

Extraversion item 2

BF7

Neuroticism item 2

BF8

Openness item 2

BF9

Agreeableness item 2

BF10

Conscientiousness item 2

BF11

Extraversion item 3

BF12

Neuroticism item 3

BF13

Openness item 3

BF14

Agreeableness item 3

BF15

Conscientiousness item 3

BF16

Extraversion item 4

BF17

Neuroticism item 4

BF18

Openness item 4

BF19

Agreeableness item 4

BF20

Conscientiousness item 4

Source

Originally collected by the package author.

Classic fit object

Description

An R6 class representing the results of a classical (frequentist) estimation.

Super class

BayesRTMB::RTMB_Fit_Base -> Classic_Fit

Public fields

Methods

Public methods

• Classic_Fit$get_point_estimate()

• Classic_Fit$new()

• Classic_Fit$robust_se()

• Classic_Fit$compute_robust()

• Classic_Fit$bootstrap()

• Classic_Fit$compute_bootstrap()

• Classic_Fit$AIC()

• Classic_Fit$BIC()

• Classic_Fit$WAIC()

• Classic_Fit$diagnose()

• Classic_Fit$print()

• Classic_Fit$logLik()

• Classic_Fit$EAP()

• Classic_Fit$MAP()

• Classic_Fit$summary()

• Classic_Fit$anova()

• Classic_Fit$lsmeans()

• Classic_Fit$.calc_contrast()

• Classic_Fit$.get_lsmeans_df()

• Classic_Fit$.construct_par_list()

• Classic_Fit$clone()

Method get_point_estimate()

Get point estimate for a target parameter.

Usage

Classic_Fit$get_point_estimate(target, ...)

Arguments

Returns

Matrix, array, vector, or scalar point estimate.

Method new()

Create a new 'Classic_Fit' object.

Usage

Classic_Fit$new(
  model,
  par_vec = NULL,
  par = NULL,
  objective = NULL,
  log_lik = NULL,
  restricted_log_lik = NULL,
  convergence = NULL,
  sd_rep = NULL,
  df_fixed = NULL,
  random_effects = NULL,
  df_transform = NULL,
  df_generate = NULL,
  opt_history = NULL,
  transform = NULL,
  generate = NULL,
  se_samples = NULL,
  par_unc = NULL,
  vcov_unc = NULL,
  ci_method = "wald",
  laplace = TRUE,
  map = NULL,
  test_results = list(),
  view = NULL,
  fit = NULL,
  vcov = NULL,
  df_method = "auto",
  idx_fix_active = NULL,
  show_df = TRUE,
  rss = NULL,
  df_residual = NULL,
  ...
)

Arguments

Method robust_se()

Compute robust standard errors (sandwich estimator).

Usage

Classic_Fit$robust_se(
  cluster = NULL,
  type = c("HC3", "HC0", "HC1", "CR1", "CR0"),
  inplace = FALSE,
  ...
)

Arguments

Returns

Self.

Method compute_robust()

(Deprecated) Use robust_se() instead.

Usage

Classic_Fit$compute_robust(...)

Arguments

Returns

Self.

Method bootstrap()

Compute nonparametric bootstrap standard errors and confidence intervals. Currently implemented only for mediation models. Bootstrap refits mediation equations with base R lm.fit()/glm.fit() when possible, and falls back to RTMB refits for unsupported families.

Usage

Classic_Fit$bootstrap(n_boot = 1000, seed = NULL, inplace = FALSE, ...)

Arguments

Returns

A Classic_Fit object.

Method compute_bootstrap()

Compute nonparametric bootstrap standard errors and confidence intervals. Currently implemented only for mediation models. Bootstrap refits mediation equations with base R lm.fit()/glm.fit() when possible, and falls back to RTMB refits for unsupported families.

Usage

Classic_Fit$compute_bootstrap(n_boot = 1000, seed = NULL, inplace = FALSE, ...)

Arguments

Returns

A Classic_Fit object.

Method AIC()

Get the AIC of the fitted model.

Usage

Classic_Fit$AIC()

Method BIC()

Get the BIC of the fitted model.

Usage

Classic_Fit$BIC()

Method WAIC()

WAIC is not defined for classical fits.

Usage

Classic_Fit$WAIC()

Method diagnose()

Run basic diagnostics for the classical fit.

Usage

Classic_Fit$diagnose(...)

Arguments

Returns

A 'diagnose_BayesRTMB' object.

Method print()

Print the fit results.

Usage

Classic_Fit$print(...)

Arguments

Method logLik()

Get the Log-Likelihood of the fitted model.

Usage

Classic_Fit$logLik()

Method EAP()

EAP estimates are not available for Classic_Fit.

Usage

Classic_Fit$EAP(...)

Arguments

Method MAP()

MAP estimates are not available for Classic_Fit.

Usage

Classic_Fit$MAP(...)

Arguments

Method summary()

Display a summary of the estimation results.

Usage

Classic_Fit$summary(view = NULL, digits = 5, max_rows = 10, ranef = FALSE)

Arguments

Method anova()

Perform ANOVA (Wald F-tests / Chisq-tests) on the fitted model.

Usage

Classic_Fit$anova(method = c("reml", "ls"), type = 3)

Arguments

Returns

A data frame containing the ANOVA table.

Method lsmeans()

Calculate Least Squares Means (Marginal Means) and contrasts.

Usage

Classic_Fit$lsmeans(
  specs = NULL,
  pairwise = FALSE,
  simple = NULL,
  adjust = "holm",
  protect = FALSE
)

Arguments

Returns

A data frame containing the marginal means or contrasts.

Method .calc_contrast()

(Internal) Calculate metrics for a contrast.

Usage

Classic_Fit$.calc_contrast(L, specs)

Arguments

Returns

A data frame with estimate, SE, df, t-value, and p-value.

Method .get_lsmeans_df()

(Internal) Get representative DF for lsmeans.

Usage

Classic_Fit$.get_lsmeans_df(specs)

Arguments

Returns

Degrees of freedom.

Method .construct_par_list()

(Internal) Construct a list of parameters from the fit.

Usage

Classic_Fit$.construct_par_list(fit)

Arguments

Returns

A named list of parameters.

Method clone()

The objects of this class are cloneable with this method.

Usage

Classic_Fit$clone(deep = FALSE)

Arguments

Define parameter dimensions and types

Description

Define parameter dimensions and types

Usage

Dim(dim = 1, type = NULL, lower = NULL, upper = NULL, random = FALSE)

Arguments

MAP fit object

Description

MAP fit object

MAP fit object

Details

An R6 class storing optimization results from maximum a posteriori (MAP) estimation.

Super class

BayesRTMB::RTMB_Fit_Base -> map_fit

Public fields

Methods

Public methods

• MAP_Fit$get_point_estimate()

• MAP_Fit$new()

• MAP_Fit$ranef()

• MAP_Fit$draws()

• MAP_Fit$summary()

• MAP_Fit$print()

• MAP_Fit$generated_quantities()

• MAP_Fit$WAIC()

• MAP_Fit$diagnose()

• MAP_Fit$profile()

• MAP_Fit$clone()

Method get_point_estimate()

Get point estimate for a target parameter.

Usage

MAP_Fit$get_point_estimate(target, ...)

Arguments

Returns

Matrix, array, vector, or scalar point estimate.

Method new()

Create a new 'MAP_Fit' object.

Usage

MAP_Fit$new(
  model,
  par_vec = NULL,
  par = NULL,
  objective = NULL,
  log_ml = NULL,
  convergence = NULL,
  sd_rep = NULL,
  df_fixed = NULL,
  random_effects = NULL,
  df_transform = NULL,
  df_generate = NULL,
  opt_history = NULL,
  transform = NULL,
  generate = NULL,
  se_samples = NULL,
  par_unc = NULL,
  ci_method = "wald",
  laplace = TRUE,
  map = NULL,
  vcov_unc = NULL,
  marginal_vars = NULL,
  laplace_random_vars = NULL,
  idx_fix_active = NULL,
  show_df = TRUE,
  view = NULL,
  fallback_needed = NULL
)

Arguments

Method ranef()

Return random effect estimates as a named list.

Usage

MAP_Fit$ranef()

Returns

A named list of random effect estimates.

Method draws()

Extract samples from the asymptotic posterior distribution.

Usage

MAP_Fit$draws(
  pars = NULL,
  inc_random = FALSE,
  inc_transform = TRUE,
  inc_generate = TRUE,
  ...
)

Arguments

Returns

An array of samples [iterations, 1, parameters].

Method summary()

Summarize MAP estimates.

Usage

MAP_Fit$summary(
  pars = NULL,
  max_rows = 10,
  digits = 5,
  ranef = FALSE,
  view = NULL
)

Arguments

Returns

A summary object, typically a data frame.

Method print()

Print a brief summary of the fitted object.

Usage

MAP_Fit$print(pars = NULL, max_rows = 10, digits = 5, ...)

Arguments

Returns

The object itself, invisibly.

Method generated_quantities()

Compute generated quantities from the MAP estimate.

Usage

MAP_Fit$generated_quantities(code)

Arguments

Returns

The 'MAP_Fit' object itself (invisibly). Results are added or updated in the 'generate' list and 'df_generate'.

Method WAIC()

Compute an approximate WAIC from sampling-based uncertainty propagation.

Usage

MAP_Fit$WAIC()

Returns

A 'waic_BayesRTMB' object.

Method diagnose()

Run basic diagnostics for the MAP fit.

Usage

MAP_Fit$diagnose(...)

Arguments

Returns

A 'diagnose_BayesRTMB' object.

Method profile()

Calculate Profile Likelihood confidence intervals for specific parameters.

Usage

MAP_Fit$profile(
  pars = NULL,
  level = 0.95,
  trace = FALSE,
  digits = 5,
  show_plot = FALSE,
  quiet = FALSE,
  jacobian = "none",
  ...
)

Arguments

Returns

A data frame containing the profile-based confidence intervals, with the raw profile objects stored in the "profiles" attribute.

Method clone()

The objects of this class are cloneable with this method.

Usage

MAP_Fit$clone(deep = FALSE)

Arguments

MCMC fit object

Description

MCMC fit object

MCMC fit object

Details

An R6 class storing posterior samples and related information from MCMC estimation.

Bayes factors are computed as ratios of marginal likelihoods estimated by bridge sampling. The comparison model can be constructed by fixing parameters with 'fixed = list(...)', or supplied as an already fitted 'MCMC_Fit' object via 'comparison_fit'. Bayes factors are computed only by marginal-likelihood model comparison.

Super class

BayesRTMB::RTMB_Fit_Base -> mcmc_fit

Public fields

Methods

Public methods

• MCMC_Fit$get_point_estimate()

• MCMC_Fit$new()

• MCMC_Fit$print()

• MCMC_Fit$draws()

• MCMC_Fit$summary()

• MCMC_Fit$unconstrain_draws()

• MCMC_Fit$log_prob()

• MCMC_Fit$bridgesampling()

• MCMC_Fit$WAIC()

• MCMC_Fit$diagnose()

• MCMC_Fit$bayes_factor()

• MCMC_Fit$transformed_draws()

• MCMC_Fit$generated_quantities()

• MCMC_Fit$resolve_switching()

• MCMC_Fit$clone()

Method get_point_estimate()

Get point estimate for a target parameter.

Usage

MCMC_Fit$get_point_estimate(target, chains = NULL, best_chains = NULL)

Arguments

Returns

Matrix, array, vector, or scalar point estimate.

Method new()

Create a new 'MCMC_Fit' object.

Usage

MCMC_Fit$new(
  model,
  fit,
  random_fit,
  eps,
  accept,
  treedepth,
  laplace,
  posterior_mean,
  max_treedepth = NULL,
  pd_error_count = NULL
)

Arguments

Method print()

Print a brief summary of the fitted object.

Usage

MCMC_Fit$print(...)

Arguments

Returns

The object itself, invisibly.

Method draws()

Extract posterior draws for selected parameters.

Usage

MCMC_Fit$draws(
  pars = NULL,
  chains = NULL,
  best_chains = NULL,
  inc_random = FALSE,
  inc_transform = TRUE,
  inc_generate = TRUE
)

Arguments

Returns

Posterior draws.

Method summary()

Summarize posterior draws.

Usage

MCMC_Fit$summary(
  pars = NULL,
  chains = NULL,
  best_chains = NULL,
  max_rows = 10,
  digits = 2,
  inc_random = FALSE,
  inc_transform = TRUE,
  inc_generate = TRUE
)

Arguments

Returns

A summary object.

Method unconstrain_draws()

Transform posterior draws to the unconstrained scale.

Usage

MCMC_Fit$unconstrain_draws()

Returns

Posterior draws on the unconstrained scale.

Method log_prob()

Evaluate log-probability values.

Usage

MCMC_Fit$log_prob(safe = FALSE)

Arguments

Returns

Numeric vector of log-probability values.

Method bridgesampling()

Estimate the marginal likelihood by bridge sampling.

Usage

MCMC_Fit$bridgesampling(
  method = "normal",
  use_neff = TRUE,
  seed = NULL,
  max_iter = 100
)

Arguments

Returns

Bridge sampling result.

Method WAIC()

Compute WAIC from pointwise generated log likelihood.

Usage

MCMC_Fit$WAIC(...)

Arguments

Returns

A 'waic_BayesRTMB' object.

Method diagnose()

Run basic diagnostics for the MCMC fit.

Usage

MCMC_Fit$diagnose(...)

Arguments

Returns

A 'diagnose_BayesRTMB' object.

Method bayes_factor()

Calculate the Bayes factor by marginal-likelihood model comparison.

Usage

MCMC_Fit$bayes_factor(
  fixed = NULL,
  comparison_fit = NULL,
  bs_method = "normal",
  error_threshold = 0.2,
  ...
)

Arguments

Returns

A list of class bayes_factor_rtmb containing Bayes factor results.

Method transformed_draws()

Compute transformed parameters from posterior draws.

Usage

MCMC_Fit$transformed_draws(tran_fn = NULL)

Arguments

Returns

Transformed parameter draws.

Method generated_quantities()

Compute generated quantities from posterior draws.

Usage

MCMC_Fit$generated_quantities(code)

Arguments

Returns

The 'MCMC_Fit' object itself (invisibly). Results are appended to the 'generate_fit' field.

Method resolve_switching()

Resolve label switching in posterior draws.

Usage

MCMC_Fit$resolve_switching(
  target,
  linked = NULL,
  overwrite = TRUE,
  scalar_fns = list()
)

Arguments

Returns

Relabeled draws or updated object.

Method clone()

The objects of this class are cloneable with this method.

Usage

MCMC_Fit$clone(deep = FALSE)

Arguments

Base class for RTMB Fit objects

Description

Base class for RTMB Fit objects

Base class for RTMB Fit objects

Details

An R6 base class providing common methods for Bayesian and MAP inference results.

Public fields

Methods

Public methods

• RTMB_Fit_Base$get_point_estimate()

• RTMB_Fit_Base$estimate()

• RTMB_Fit_Base$EAP()

• RTMB_Fit_Base$MAP()

• RTMB_Fit_Base$rotate()

• RTMB_Fit_Base$fa_rotate()

• RTMB_Fit_Base$clone()

Method get_point_estimate()

Abstract method to get a point estimate for a target parameter.

Usage

RTMB_Fit_Base$get_point_estimate(target, ...)

Arguments

Returns

A numeric array or matrix of the point estimate.

Method estimate()

Get point estimates for parameters, transformed parameters, and generated quantities.

Usage

RTMB_Fit_Base$estimate(
  pars = NULL,
  type = c("mean", "EAP", "marginal_map", "joint_map", "MAP"),
  component = c("all", "parameters", "transform", "generate"),
  chains = NULL,
  best_chains = NULL,
  drop = TRUE,
  ...
)

Arguments

Returns

A named list of point estimates, or a single value if 'drop = TRUE'.

Method EAP()

Calculate Expected A Posteriori (EAP) estimates from posterior samples.

Usage

RTMB_Fit_Base$EAP(
  pars = "parameters",
  chains = NULL,
  best_chains = NULL,
  drop = FALSE,
  ...
)

Arguments

Returns

A named list of EAP estimates.

Method MAP()

Calculate Maximum A Posteriori (MAP) estimates.

Usage

RTMB_Fit_Base$MAP(
  pars = "parameters",
  chains = NULL,
  best_chains = NULL,
  type = c("marginal", "joint"),
  drop = FALSE,
  ...
)

Arguments

Returns

A named list of MAP estimates.

Method rotate()

Rotate sampled parameters.

Usage

RTMB_Fit_Base$rotate(target, reference = NULL, linked = NULL)

Arguments

Returns

The updated object invisibly.

Method fa_rotate()

Rotate factor loadings and optional factor scores.

Usage

RTMB_Fit_Base$fa_rotate(
  target = "L",
  linked = NULL,
  scores = NULL,
  rotate = "promax",
  ...
)

Arguments

Returns

The updated object invisibly.

Method clone()

The objects of this class are cloneable with this method.

Usage

RTMB_Fit_Base$clone(deep = FALSE)

Arguments

RTMB model object

Description

An R6 class representing a Bayesian model built with RTMB. This class stores model components and provides methods for building the automatic differentiation object, optimizing the posterior, and drawing posterior samples.

Public fields

Methods

Public methods

• RTMB_Model$new()

• RTMB_Model$prepare_init()

• RTMB_Model$get_par_list()

• RTMB_Model$get_ad_obj()

• RTMB_Model$to_unconstrained()

• RTMB_Model$to_constrained()

• RTMB_Model$unconstrained_vector_to_list()

• RTMB_Model$constrained_vector_to_list()

• RTMB_Model$calculate_satterthwaite_df()

• RTMB_Model$calculate_reml_satterthwaite_df()

• RTMB_Model$build_ad_obj()

• RTMB_Model$optimize()

• RTMB_Model$classic()

• RTMB_Model$sample()

• RTMB_Model$variational()

• RTMB_Model$print_code()

• RTMB_Model$print_log_prob()

• RTMB_Model$print_transform()

• RTMB_Model$print_generate()

• RTMB_Model$get_n_obs()

• RTMB_Model$fixed_model()

• RTMB_Model$clone()

Method new()

Create a new 'RTMB_Model' object.

Usage

RTMB_Model$new(
  data,
  par_list,
  log_prob,
  transform = NULL,
  generate = NULL,
  par_names = NULL,
  init = NULL,
  view = NULL,
  code = NULL
)

Arguments

Method prepare_init()

Prepare and format initial values for the model parameters.

Usage

RTMB_Model$prepare_init(init_arg)

Arguments

Returns

A flat numeric vector of constrained initial values.

Method get_par_list()

Get the current parameters as a named list.

Usage

RTMB_Model$get_par_list(init = NULL)

Arguments

Returns

A named list of constrained parameters.

Method get_ad_obj()

Get the RTMB automatic differentiation object.

Usage

RTMB_Model$get_ad_obj(...)

Arguments

Returns

A TMB objective object (ad_obj).

Method to_unconstrained()

Convert a list of constrained parameters to unconstrained space.

Usage

RTMB_Model$to_unconstrained(par_list_constrained)

Arguments

Returns

A named list of parameters in unconstrained space.

Method to_constrained()

Convert a list of unconstrained parameters to constrained (natural) space.

Usage

RTMB_Model$to_constrained(par_list_unconstrained)

Arguments

Returns

A named list of parameters in constrained space.

Method unconstrained_vector_to_list()

Convert a flat unconstrained vector to a named list.

Usage

RTMB_Model$unconstrained_vector_to_list(vec)

Arguments

Returns

A named list of unconstrained parameters.

Method constrained_vector_to_list()

Convert a flat constrained vector to a named list.

Usage

RTMB_Model$constrained_vector_to_list(vec)

Arguments

Returns

A named list of constrained parameters.

Method calculate_satterthwaite_df()

Calculate Satterthwaite degrees of freedom for parameters.

Usage

RTMB_Model$calculate_satterthwaite_df(
  ad_obj,
  idx_fix_active = NULL,
  L_u_total = NULL,
  opt_par = NULL,
  max_df = NULL,
  silent = FALSE,
  return_sensitivities = FALSE
)

Arguments

Returns

A numeric vector of estimated degrees of freedom (length = L_u_total). Inf for random effects.

Method calculate_reml_satterthwaite_df()

Calculate Satterthwaite degrees of freedom for integrated fixed effects (REML).

Usage

RTMB_Model$calculate_reml_satterthwaite_df(
  ad_obj,
  opt_par,
  beta_idx,
  max_df = NULL,
  silent = FALSE,
  return_sensitivities = FALSE
)

Arguments

Returns

A numeric vector of estimated degrees of freedom for the fixed effects.

Method build_ad_obj()

Build the RTMB automatic differentiation object.

Usage

RTMB_Model$build_ad_obj(
  init = NULL,
  laplace = FALSE,
  jacobian_target = "none",
  map = NULL,
  .marginal_vars = NULL
)

Arguments

Returns

An RTMB objective object (ad_obj).

Method optimize()

Perform Maximum Likelihood or Maximum A Posteriori (MAP) estimation.

Usage

RTMB_Model$optimize(
  laplace = TRUE,
  init = NULL,
  num_estimate = 1,
  control = list(),
  optimizer = "nlminb",
  method = "BFGS",
  map = NULL,
  fixed = NULL,
  se_method = c("wald", "sampling", "none"),
  num_samples = 1000,
  seed = 123,
  marginal = "none",
  df_method = "auto",
  view = NULL,
  ...
)

Arguments

Returns

A fitted 'MAP_Fit' object.

Method classic()

Perform frequentist inference (REML/ML) with model-appropriate degrees of freedom.

Usage

RTMB_Model$classic(
  df_method = "auto",
  se_method = c("wald", "sampling"),
  num_samples = 1000,
  seed = 123,
  view = NULL,
  map = NULL,
  fixed = NULL,
  ...
)

Arguments

Returns

A 'Classic_Fit' object.

Method sample()

Draw posterior samples from the model.

Usage

RTMB_Model$sample(
  sampling = 1000,
  warmup = 1000,
  chains = 4,
  thin = 1,
  seed = sample.int(1e+06, 1),
  delta = 0.8,
  max_treedepth = 10,
  parallel = FALSE,
  laplace = FALSE,
  init = NULL,
  init_jitter = 0.1,
  save_csv = NULL,
  map = NULL,
  fixed = NULL
)

Arguments

Returns

A fitted 'MCMC_Fit' object.

Method variational()

Run Automatic Differentiation Variational Inference (ADVI).

Usage

RTMB_Model$variational(
  iter = 3000,
  tol_rel_obj = 0.005,
  window_size = 100,
  num_samples = 1000,
  num_estimate = 4,
  alpha = 0.01,
  laplace = FALSE,
  print_freq = 1000,
  method = c("meanfield", "fullrank", "hybrid"),
  parallel = FALSE,
  seed = sample.int(1e+06, 1),
  init = NULL,
  save_csv = NULL,
  map = NULL,
  fixed = NULL
)

Arguments

Returns

A fitted 'VB_Fit' object containing posterior samples and diagnostic information.

Method print_code()

Print model code or model structure.

Usage

RTMB_Model$print_code()

Returns

The object itself, invisibly.

Method print_log_prob()

Print the internal log-probability function.

Usage

RTMB_Model$print_log_prob()

Returns

The object itself, invisibly.

Method print_transform()

Print the internal transformation function.

Usage

RTMB_Model$print_transform()

Returns

The object itself, invisibly.

Method print_generate()

Print the internal generation function.

Usage

RTMB_Model$print_generate()

Returns

The object itself, invisibly.

Method get_n_obs()

Automatically detect and count the number of independent data points (observations).

Usage

RTMB_Model$get_n_obs()

Returns

Integer; total number of observations, or NULL if model code is missing.

Method fixed_model()

Create a model with fixed parameters.

Usage

RTMB_Model$fixed_model(fixed = list(), silent = FALSE)

Arguments

Returns

A new RTMB_Model object.

Method clone()

The objects of this class are cloneable with this method.

Usage

RTMB_Model$clone(deep = FALSE)

Arguments

VB fit object

Description

VB fit object

VB fit object

Details

An R6 class storing posterior samples and related information from Automatic Differentiation Variational Inference (ADVI).

Super class

BayesRTMB::RTMB_Fit_Base -> advi_fit

Public fields

Methods

Public methods

• VB_Fit$get_point_estimate()

• VB_Fit$new()

• VB_Fit$print()

• VB_Fit$draws()

• VB_Fit$summary()

• VB_Fit$plot_elbo()

• VB_Fit$plot_trajectory()

• VB_Fit$transformed_draws()

• VB_Fit$generated_quantities()

• VB_Fit$WAIC()

• VB_Fit$diagnose()

• VB_Fit$clone()

Method get_point_estimate()

Get point estimate for a target parameter.

Usage

VB_Fit$get_point_estimate(target, chains = NULL, best_chains = NULL)

Arguments

Returns

Matrix, array, vector, or scalar point estimate.

Method new()

Create a new 'VB_Fit' object.

Usage

VB_Fit$new(
  model,
  fit,
  random_fit,
  elbo_history,
  laplace,
  posterior_mean,
  ELBO,
  rel_obj_vals,
  best_chain,
  mu_history
)

Arguments

Method print()

Print a brief summary of the fitted object.

Usage

VB_Fit$print(...)

Arguments

Returns

The object itself, invisibly.

Method draws()

Extract posterior draws for selected parameters.

Usage

VB_Fit$draws(
  pars = NULL,
  chains = NULL,
  best_chains = NULL,
  inc_random = FALSE,
  inc_transform = TRUE,
  inc_generate = TRUE,
  best_only = FALSE
)

Arguments

Returns

A 3D array of posterior draws '[iterations, chains, parameters]'.

Method summary()

Summarize posterior draws.

Usage

VB_Fit$summary(
  pars = NULL,
  max_rows = 10,
  digits = 2,
  inc_random = FALSE,
  inc_transform = TRUE,
  inc_generate = TRUE
)

Arguments

Returns

A data frame containing the summarized posterior statistics.

Method plot_elbo()

Plot the ELBO history to diagnose convergence.

Usage

VB_Fit$plot_elbo(tail_n = 1000, ests = NULL, type = "l", ...)

Arguments

Returns

The object itself, invisibly.

Method plot_trajectory()

Plot the parameter trajectory from the final optimization window.

Usage

VB_Fit$plot_trajectory(pars = NULL, type = "l", ...)

Arguments

Returns

The object itself, invisibly.

Method transformed_draws()

Compute transformed parameters from posterior draws.

Usage

VB_Fit$transformed_draws(tran_fn = NULL)

Arguments

Returns

The 'VB_Fit' object itself, invisibly.

Method generated_quantities()

Compute generated quantities from posterior draws.

Usage

VB_Fit$generated_quantities(code)

Arguments

Returns

The 'VB_Fit' object itself (invisibly). Results are appended to the 'generate_fit' field.

Method WAIC()

Compute WAIC from pointwise generated log likelihood.

Usage

VB_Fit$WAIC(...)

Arguments

Returns

A 'waic_BayesRTMB' object.

Method diagnose()

Run basic diagnostics for the variational fit.

Usage

VB_Fit$diagnose(...)

Arguments

Returns

A 'diagnose_BayesRTMB' object.

Method clone()

The objects of this class are cloneable with this method.

Usage

VB_Fit$clone(deep = FALSE)

Arguments

Calculate Bayes Factor

Description

Compare two models by calculating the Bayes Factor based on their marginal likelihoods.

Usage

bayes_factor(logml1, logml2, error_threshold = 0.2)

Arguments

Value

An object of class 'bayes_factor' containing the Bayes factor, log Bayes factor, approximate estimation error, and interpretation.

Examples

  # Compare two models using Bayes Factor
  data(debate, package = "BayesRTMB")
  fit1 <- rtmb_lm(sat ~ talk, data = debate)
  fit2 <- rtmb_lm(sat ~ talk + perf, data = debate)
  map1 <- fit1$optimize()
  map2 <- fit2$optimize()
  bf <- bayes_factor(map1, map2)
  print(bf)

Bernoulli log-probability mass function with logit parameterization

Description

Bernoulli log-probability mass function with logit parameterization

Usage

bernoulli_logit_lpmf(x, eta, sum = TRUE)

Arguments

Value

The sum of the log-probability.

Bernoulli log-probability mass function

Description

Bernoulli log-probability mass function

Usage

bernoulli_lpmf(x, prob, sum = TRUE)

Arguments

Value

The sum of the log-probability.

Beta log-probability density function

Description

Beta log-probability density function

Usage

beta_lpdf(x, a, b, sum = TRUE)

Arguments

Value

The sum of the log-density.

Beverage Preference Data

Description

A dataset containing preference ratings for various beverages. This data was used to demonstrate a random effect model in metric multidimensional unfolding.

Usage

beverage

Format

A data frame (or matrix) with [INSERT NUMBER OF ROWS] rows and [INSERT NUMBER OF COLUMNS] columns.

AJ

Apple Juice preference rating

BT

Black Tea preference rating

CF

Coffee preference rating

CL

Cola preference rating

GT

Green Tea preference rating

OJ

Orange Juice preference rating

OT

Oolong Tea preference rating

Source

Adachi, K. (2000). A Random Effect Model in Metric Multidimensional Unfolding. The Japanese Journal of Behaviormetrics, 27(1), 12-23. doi:10.2333/jbhmk.27.12

Binomial log-probability mass function with logit parameterization

Description

Binomial log-probability mass function with logit parameterization

Usage

binomial_logit_lpmf(x, size, eta, sum = TRUE)

Arguments

Value

The sum of the log-probability.

Binomial log-probability mass function

Description

Binomial log-probability mass function

Usage

binomial_lpmf(x, size, prob, sum = TRUE)

Arguments

Value

The sum of the log-probability.

Best-Worst Categorical Logit Log-PMF (RTMB optimized)

Description

Best-Worst Categorical Logit Log-PMF (RTMB optimized)

Usage

bw_categorical_logit_lpmf(x, U, lambda = 1)

Arguments

Value

Log-probability (scalar advector).

Categorical log-probability mass function with logit parameterization

Description

Categorical log-probability mass function with logit parameterization

Usage

categorical_logit_lpmf(x, eta, sum = TRUE)

Arguments

Value

The sum of the log-probability.

Categorical log-probability mass function

Description

Categorical log-probability mass function

Usage

categorical_lpmf(x, prob, sum = TRUE)

Arguments

Value

The sum of the log-probability.

Cauchy log-probability density function

Description

Cauchy log-probability density function

Usage

cauchy_lpdf(x, location, scale, sum = TRUE)

Arguments

Value

The sum of the log-density.

Centered / Centered matrix multivariate normal log-probability density function

Description

Centered / Centered matrix multivariate normal log-probability density function

Usage

centered_multi_normal_lpdf(x, sigma, K = NULL)

Arguments

Value

The log-density.

Centered triangular multivariate normal log-probability density function

Description

Centered triangular multivariate normal log-probability density function

Usage

centered_tri_multi_normal_lpdf(x, sigma)

Arguments

Value

The log-density.

Calculate Conditional Effects

Description

Calculate and visualize the predicted values (marginal effects) of a variable, potentially conditional on the levels of another variable (interaction).

Usage

conditional_effects(fit, effect, prob = 0.95, sd_multiplier = 1, ...)

Arguments

Value

An object of class 'ce_rtmb' containing the predicted values and their credible intervals.

Examples

  data(debate, package = "BayesRTMB")
  fit <- rtmb_lm(sat ~ talk * perf, data = debate)
  map_fit <- fit$optimize()
  ce <- conditional_effects(map_fit, effect = "talk:perf")
  plot(ce)
  summary(ce)

Calculate conditional effects for MCMC fit objects

Description

Calculate conditional effects for MCMC fit objects

Usage

## S3 method for class 'mcmc_fit'
conditional_effects(
  fit,
  effect,
  prob = 0.95,
  sd_multiplier = 1,
  resolution = 100,
  ...
)

Arguments

Value

A 'ce_rtmb' object.

Debate Simulation Data

Description

A simulated dataset containing individual satisfaction scores and related variables from a 3-person group debate experiment. This dummy data was created by the package author to demonstrate hierarchical (multilevel) modeling.

Usage

debate

Format

A data frame with 300 rows and 6 columns:

group

Group identifier. Each group consists of 3 members.

sat

Self-reported satisfaction level of the individual after the debate.

talk

Amount of talk or utterances by the individual during the debate.

perf

Evaluated performance of the group debate (group-level variable).

skill

Evaluated conversation skill of the individual.

cond

Experimental condition indicator (e.g., 0 = control, 1 = treatment).

Source

Simulated dummy data created by the package author.

Dirichlet log-probability density function

Description

Dirichlet log-probability density function

Usage

dirichlet_lpdf(x, alpha)

Arguments

Value

The sum of the log-density.

Euclidean distance

Description

Euclidean distance

Usage

distance(x, y, eps = 1e-08)

Arguments

Value

The Euclidean distance between x and y.

Probability Distributions for RTMB Models

Description

This page summarizes the probability distributions available within the 'rtmb_code' block. These functions are designed to be used with the Stan-like sampling syntax ('~'), which internally adds the log-density to the model's total log-posterior.

Details

Syntax Styles: In 'rtmb_code', you can specify distributions in two ways:

• Sampling Syntax (Recommended): y ~ normal(mu, sigma)

• Explicit Function Call: lp <- lp + normal_lpdf(y, mu, sigma)

Continuous Distributions (LPDF):

• normal(mean, sd): Normal distribution.

• lognormal(meanlog, sdlog): Lognormal distribution.

• exponential(rate): Exponential distribution.

• cauchy(location, scale): Cauchy distribution.

• student_t(df, mu, sigma): Student's t-distribution.

• gamma(shape, rate): Gamma distribution.

• inverse_gamma(shape, scale): Inverse-gamma distribution.

• beta(a, b): Beta distribution.

Discrete Distributions (LPMF):

• bernoulli(prob) / bernoulli_logit(eta): Binary outcomes.

• binomial(size, prob) / binomial_logit(size, eta): Binomial outcomes.

• poisson(mean): Poisson count data.

• neg_binomial_2(mu, size): Negative binomial (mean/dispersion parameterization).

• ordered_logistic(eta, cutpoints): Ordered categorical outcomes.

• sequential_logistic(eta, cutpoints): Sequential ordered categorical outcomes.

Multivariate and Matrix Distributions:

• multi_normal(mean, Sigma): Standard multivariate normal distribution.

• multi_student_t(df, mean, Sigma): Multivariate Student's t-distribution.

• multi_cauchy(mean, Sigma): Multivariate Cauchy distribution (Student-t with df=1).

• lkj_corr(eta): LKJ prior for correlation matrices.

• dirichlet(alpha): Dirichlet distribution for simplexes.

• lower_tri_normal(mean, sd): Normal distribution for elements of a lower-triangular matrix.

• centered_tri_multi_normal(sigma): Multivariate normal for centered triangular matrices (used in identification constraints).

• sufficient_multi_normal_fa(S_mat, N, y_bar, mu, psi, Lambda): Factor analysis likelihood using sufficient statistics (highly efficient for large sample sizes).

Vectorization: Most univariate distributions are vectorized. If y and mu are vectors, y ~ normal(mu, sigma) will calculate the sum of log-densities for all elements efficiently.

Internal function to convert an environment to a list while maintaining original order

Description

Internal function to convert an environment to a list while maintaining original order

Usage

env_to_ordered_list(env, orig_list, code_ast = NULL)

Basic Effective Sample Size for a single chain or pooled chains

Description

Basic Effective Sample Size for a single chain or pooled chains

Usage

ess_basic(sims)

Arguments

Value

A numeric value.

Calculate Bulk Effective Sample Size

Description

Calculate Bulk Effective Sample Size

Usage

ess_bulk(sims)

Arguments

Value

A numeric value.

Calculate Tail Effective Sample Size (at 2.5% and 97.5% quantiles)

Description

Calculate Tail Effective Sample Size (at 2.5% and 97.5% quantiles)

Usage

ess_tail95(sims)

Arguments

Value

A numeric value.

Exponential log-probability density function

Description

Exponential log-probability density function

Usage

exponential_lpdf(x, rate, sum = TRUE)

Arguments

Value

The sum of the log-density.

Factor analysis multivariate normal log-probability density function

Description

Woodbury matrix identity is used for efficient computation.

Usage

fa_multi_normal_lpdf(x, mu, Lambda, psi, sum = TRUE)

Arguments

Value

The sum of the log-density.

Smooth absolute value function

Description

Smooth absolute value function

Usage

fabs(x, eps = 1e-08)

Arguments

Value

The smooth absolute value of x.

Multivariate normal log-probability density function parameterized by Cholesky factor of correlation matrix (FIML)

Description

Multivariate normal log-probability density function parameterized by Cholesky factor of correlation matrix (FIML)

Usage

fiml_multi_normal_CF_lpdf(x, mean, sd, CF_Omega)

Arguments

Value

The sum of the log-density using Full Information Maximum Likelihood.

Factor analysis multivariate normal log-probability density function (FIML)

Description

Factor analysis multivariate normal log-probability density function (FIML)

Usage

fiml_multi_normal_fa_lpdf(x, mu, Lambda, psi)

Arguments

Value

The sum of the log-density.

Gamma log-probability density function

Description

Gamma log-probability density function

Usage

gamma_lpdf(x, shape, rate, sum = TRUE)

Arguments

Value

The sum of the log-density.

Gaussian Process Log-Density (Squared Exponential Kernel)

Description

Calculates the log-density of a Gaussian Process with a Squared Exponential (RBF) kernel.

Usage

gaussian_process_lpdf(
  y,
  x,
  mean = 0,
  magnitude = 1,
  smoothing = 1,
  noise = 0.01,
  sum = TRUE
)

Arguments

Value

Log-density value.

Generate Random Initial Values

Description

Generates random initial values for model parameters by drawing from a uniform distribution on the unconstrained scale and then transforming them to the constrained scale.

Usage

generate_random_init(pl_full, par_list, range = 2)

Arguments

Value

A numeric vector of initial values where 'NA' elements are replaced with randomly generated constrained values.

Half-Normal log-probability density function

Description

Half-Normal log-probability density function

Usage

half_normal_lpdf(x, sd, sum = TRUE)

Arguments

Value

The sum of the log-density.

Internal function to search AST and inject namespace to package functions

Description

Internal function to search AST and inject namespace to package functions

Usage

inject_namespace(expr, pkg = "BayesRTMB")

Arguments

Inverse logit function

Description

Inverse logit function

Usage

inv_logit(x)

Arguments

Value

The inverse logit of x.

Inverse-gamma log-probability density function

Description

Inverse-gamma log-probability density function

Usage

inverse_gamma_lpdf(x, shape, scale, sum = TRUE)

Arguments

Value

The sum of the log-density.

Calculate Item Response Curve / Category Response Curve

Description

Calculate Item Response Curve / Category Response Curve

Usage

item_curve(x, ...)

Arguments

Item Response Curve for RTMB_Fit_Base

Description

Item Response Curve for RTMB_Fit_Base

Usage

## S3 method for class 'RTMB_Fit_Base'
item_curve(x, theta_seq = seq(-4, 4, length.out = 100), items = NULL, ...)

Arguments

Examples

  fit <- rtmb_irt(data = BigFive[, 1:5], model = "2PL", type = "ordered")
  map_fit <- fit$optimize()
  ic <- item_curve(map_fit)
  plot(ic)

Calculate Item Information Function

Description

Calculate Item Information Function

Usage

item_info(x, ...)

Arguments

Examples

  fit <- rtmb_irt(data = BigFive[, 1:5], model = "2PL", type = "ordered")
  map_fit <- fit$optimize()
  ii <- item_info(map_fit)
  plot(ii)

Item Information Function for RTMB_Fit_Base

Description

Item Information Function for RTMB_Fit_Base

Usage

## S3 method for class 'RTMB_Fit_Base'
item_info(x, theta_seq = seq(-4, 4, length.out = 100), items = NULL, ...)

Arguments

Laplace log-probability density function

Description

Laplace log-probability density function

Usage

laplace_lpdf(x, location = 0, scale = 1, sum = TRUE)

Arguments

Value

The sum of the log-density.

LKJ correlation log-probability density function for Cholesky factors

Description

LKJ correlation log-probability density function for Cholesky factors

Usage

lkj_CF_corr_lpdf(CF_Omega, eta = 1)

Arguments

Value

The log-density.

LKJ correlation log-probability density function

Description

LKJ correlation log-probability density function

Usage

lkj_corr_lpdf(Omega, eta = 1)

Arguments

Value

The log-density.

Log of one minus x

Description

Log of one minus x

Usage

log1m(x)

Arguments

Value

The natural logarithm of 1 - x.

Log of one minus exponential of x

Description

Log of one minus exponential of x

Usage

log1m_exp(x)

Arguments

Value

The natural logarithm of 1 - exp(x).

Log of one plus exponential of x

Description

Log of one plus exponential of x

Usage

log1p_exp(x)

Arguments

Value

The natural logarithm of 1 + exp(x).

Log determinant of a Cholesky factor

Description

Log determinant of a Cholesky factor

Usage

log_det_chol(L)

Arguments

Value

The log determinant of the corresponding covariance matrix.

Log mixture of two probabilities

Description

Log mixture of two probabilities

Usage

log_mix(theta, lp1, lp2)

Arguments

Value

The log-probability of the mixture.

Log-softmax function

Description

Log-softmax function

Usage

log_softmax(x)

Arguments

Value

A numeric vector of log-softmax probabilities.

Log-sum-exp function

Description

Log-sum-exp function

Usage

log_sum_exp(x, y = NULL)

Arguments

Value

The log of the sum of the exponentials.

Log-sum-exp function for matrices (row-wise)

Description

Log-sum-exp function for matrices (row-wise)

Usage

log_sum_exp_matrix(M)

Arguments

Value

A numeric vector of row-wise log-sum-exp.

Logit function

Description

Logit function

Usage

logit(x)

Arguments

Value

The logit of x.

Logit-Normal log-probability density function

Description

Logit-Normal log-probability density function

Usage

logit_normal_lpdf(x, mu, sd, sum = TRUE)

Arguments

Value

The sum of the log-density.

Lognormal log-probability density function

Description

Lognormal log-probability density function

Usage

lognormal_lpdf(x, meanlog, sdlog, sum = TRUE)

Arguments

Value

The sum of the log-density.

Lower-triangular normal log-probability density function

Description

Lower-triangular normal log-probability density function

Usage

lower_tri_normal_lpdf(x, mean = 0, sd = 1)

Arguments

Value

The log-density calculated only for the non-zero lower triangular elements.

Least Squares Means (Marginal Means)

Description

Calculate least squares means (also known as marginal means or predicted means) for categorical factors in a fitted model. It also supports pairwise comparisons and simple main effects.

Usage

lsmeans(object, specs, ...)

Arguments

Make Best and Worst Responses from Best-Worst Pair Indices

Description

Converts a matrix of Best-Worst pair indices ('Y_dif') into separate 'Best' and 'Worst' data frames. The returned values are positions within each row of 'sets', not the item labels stored in 'sets'.

Usage

make_bw_from_ydif(Y_dif, sets)

restore_bw_from_ydif(Y_dif, sets)

Arguments

Value

A list with two data frames: 'Best' and 'Worst'.

Reconstruct an Observation-Level Random-Effect Design Matrix

Description

Converts the transposed block random-effect design matrix used internally by 'rtmb_glmer()' into an observation-level matrix. This is kept as a small utility for users who want to work from the block matrix returned by 'make_glmer_re_terms()'.

Usage

make_glmer_Z_matrix(Zt, group_idx, N = length(group_idx))

Arguments

Value

A matrix with 'N' rows and one column per random-effect coefficient.

Prepare GLMM Formula Components

Description

Builds the response, fixed-effect model matrix, and lme4-style random-effect design components used by 'rtmb_glmer()'. This is the user-facing helper that reproduces what the wrapper places in the generated 'setup' block.

Usage

make_glmer_re_terms(
  formula,
  data,
  family = "gaussian",
  resid_group = NULL,
  resid_time = NULL,
  within = NULL,
  factors = NULL,
  missing = "listwise"
)

Arguments

Value

A list containing 'Y', 'X', 'trials', 'offset', 'N', fixed-effect metadata, and random-effect terms.

Create Initial Values for Multidimensional Unfolding

Description

Create Initial Values for Multidimensional Unfolding

Usage

make_init_mdu(
  Y,
  D,
  distance = c("squared", "euclidean"),
  alpha = c("random", "fix"),
  distance_eps = 1e-04,
  min_sigma = 0.1
)

Arguments

Value

A named list of initial values.

Make Best-Worst Pair Indices from Best and Worst Responses

Description

Converts separate 'Best' and 'Worst' response matrices into 'Y_dif' pair indices for Best-Worst MDU models. The returned index is the position of the ordered '(best, worst)' pair among all 'C * (C - 1)' position pairs generated from each row of 'sets'.

Usage

make_ydif_from_bw(Best, Worst, sets)

Arguments

Value

An integer matrix of pair indices ('Y_dif').

Maximum A Posteriori (MAP) Estimate

Description

Maximum A Posteriori (MAP) Estimate

Usage

map_est(z)

Arguments

Value

A numeric value.

Mathematical and Matrix Utility Functions for RTMB Models

Description

This page summarizes the utility functions provided to assist in writing efficient and numerically stable model code within 'rtmb_code'. These functions are specifically designed to be compatible with RTMB's Automatic Differentiation (AD).

Details

1. Link and Inverse Functions:

• logit(x): Computes the logit transformation log(x/(1-x)).

• inv_logit(x): Computes the inverse logit (logistic) transformation 1/(1+exp(-x)).

2. Numerical Stability Utilities: These functions use specialized algorithms to prevent overflow/underflow in AD calculations.

• log_sum_exp(x): Safely computes log(sum(exp(x))) using the log-sum-exp trick.

• log1p_exp(x): Computes log(1 + exp(x)) stably.

• log1m_exp(x): Computes log(1 - exp(x)) stably for x < 0.

• log_softmax(x): Computes the log of the softmax function.

• fabs(x): A smooth version of the absolute value function abs(x) to ensure differentiability at zero.

3. Matrix and Vector Transformations: Used to convert unconstrained vectors into structured matrices (e.g., for factor analysis or identification constraints).

• sum_to_zero(x): Transforms a vector of length K-1 into a vector of length K that sums to zero.

• to_lower_tri(x, M, D): Fills a matrix of size M x D with elements from vector x in a lower-triangular fashion.

• to_centered_matrix(x, R, C): Creates an R x C matrix where each column sums to zero.

• to_centered_tri(x, R, C): Creates an R x C matrix with column-wise sum-to-zero constraints on the lower elements (useful for identification in factor analysis).

4. Linear Algebra for AD:

• log_det_chol(L): Calculates the log-determinant of a covariance matrix from its Cholesky factor L.

• quad_form_chol(x, L): Computes the quadratic form x^T Sigma^(-1) x using the Cholesky factor L.

• distance(x, y): Computes the Euclidean distance between two vectors with a small epsilon for stability.

Generic Mixture log-probability density function

Description

Generic Mixture log-probability density function

Usage

mixture_lpdf(x, pi_w, lpdf_list, sum = TRUE)

Arguments

Value

The sum of the log-density.

Model Code Wrapper for RTMB

Description

Model Code Wrapper for RTMB

Usage

model_code(expr, env = parent.frame())

Arguments

Value

A standard R function object taking (dat, par).

Multivariate Cauchy log-probability density function

Description

Multivariate Cauchy log-probability density function

Usage

multi_cauchy_lpdf(x, mean, Sigma, sum = TRUE)

Arguments

Value

The sum of the log-density.

Multivariate normal log-probability density function parameterized by Cholesky factor of correlation matrix

Description

Multivariate normal log-probability density function parameterized by Cholesky factor of correlation matrix

Usage

multi_normal_CF_lpdf(x, mean, sd, CF_Omega, sum = TRUE)

Arguments

Value

The sum of the log-density.

Multivariate normal log-probability density function

Description

Multivariate normal log-probability density function

Usage

multi_normal_lpdf(x, mean, Sigma, sum = TRUE)

Arguments

Value

The sum of the log-density.

Multivariate Student-t log-probability density function

Description

Multivariate Student-t log-probability density function

Usage

multi_student_t_lpdf(x, df, mean, Sigma, sum = TRUE)

Arguments

Value

The sum of the log-density.

Multinomial log-probability mass function

Description

Multinomial log-probability mass function

Usage

multinomial_lpmf(x, size, prob)

Arguments

Value

The sum of the log-probability.

Negative binomial log-probability mass function (alternative parameterization)

Description

Negative binomial log-probability mass function (alternative parameterization)

Usage

neg_binomial_2_lpmf(x, mu, size, sum = TRUE)

Arguments

Value

The sum of the log-probability.

Negative binomial log-probability mass function

Description

Negative binomial log-probability mass function

Usage

neg_binomial_lpmf(x, size, prob, sum = TRUE)

Arguments

Value

The sum of the log-probability.

Normal log-probability density function

Description

Normal log-probability density function

Usage

normal_lpdf(x, mean, sd, sum = TRUE)

Arguments

Value

The sum of the log-density (scalar) or a vector of log-densities.

Normal Mixture log-probability density function

Description

Normal Mixture log-probability density function

Usage

normal_mixture_lpdf(x, pi_w, mean, sd, sum = TRUE)

Arguments

Value

The sum of the log-density.

Ordered logistic log-probability mass function

Description

Ordered logistic log-probability mass function

Usage

ordered_logistic_lpmf(x, eta, cutpoints, sum = TRUE)

Arguments

Value

The sum of the log-probability.

Parameter Types and Constraints in RTMB Models

Description

When declaring parameters in the parameters block using Dim(), you can specify a type to automatically apply structural constraints. These constraints ensure that parameters remain within their valid mathematical space during estimation.

Details

Standard Types:

• "vector", "matrix", "array": Standard unconstrained containers. Use lower and upper for element-wise bounds.

Ordering Constraints:

• "ordered": A vector where $x_1 < x_2 < ... < x_K$.

• "positive_ordered": A vector where $0 < x_1 < x_2 < ... < x_K$.

Constrained Vectors:

• "simplex": A vector where all elements are >= 0 and sum(x) = 1.

• "centered": A vector of length K where sum(x) = 0. (Estimated using K-1 degrees of freedom).

Correlation and Covariance Matrices:

• "corr_matrix": A symmetric, positive-definite correlation matrix (diagonal elements are 1).

• "cov_matrix": A symmetric, positive-definite covariance matrix.

• "CF_corr": The Cholesky factor of a correlation matrix.

• "CF_cov": The Cholesky factor of a covariance matrix (diagonal elements are positive).

Special Structural Matrices:

• "lower_tri": A lower-triangular matrix.

• "positive_lower_tri": A lower-triangular matrix with positive diagonal elements.

• "centered_matrix": A matrix where each column sums to zero.

• "centered_tri": A triangular matrix with column-wise sum-to-zero constraints (often used for identification in factor analysis).

See Also

rtmb_code, math_functions, distributions

Code block for parameter definitions

Description

Defines the list of parameters to be passed to 'rtmb_model'. By writing in the format 'variable_name = Dim(...)' within a block enclosed in '', evaluation is delayed, enabling strict error checking during model construction.

Usage

parameters_code(expr)

Arguments

Value

A lazily evaluated list of parameters.

Plot method for ce_rtmb class (Base R)

Description

Plot method for ce_rtmb class (Base R)

Usage

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

Arguments

Plot marginal means with error bars

Description

Plot method for rtmb_lsmeans objects.

Usage

## S3 method for class 'rtmb_lsmeans'
plot(
  x,
  y,
  error_bar = "se",
  col = NULL,
  main = NULL,
  xlab = NULL,
  ylab = NULL,
  ...
)

Arguments

Plot autocorrelation for one variable across chains

Description

Draw autocorrelation plots for a selected parameter across all chains.

Usage

plot_acf(x, var_idx = 1)

Arguments

Value

No return value. This function is called for its side effect of plotting.

Plot conditional effects

Description

Calculate and plot conditional effects for a fitted model.

Usage

plot_conditional_effects(fit, effect, ...)

Arguments

Value

Invisibly returns the plotted object of class ce_rtmb.

Plot posterior densities for MCMC samples

Description

Draw kernel density plots for each parameter across chains.

Usage

plot_dens(x, mono = FALSE)

Arguments

Value

No return value. This function is called for its side effect of plotting.

Plot parameter estimates and credible intervals (Forest Plot)

Description

Plot parameter estimates and credible intervals (Forest Plot)

Usage

plot_forest(
  x,
  prob = 0.95,
  point_estimate = c("median", "EAP", "mean", "MAP", "marginal_map", "joint_map")
)

Arguments

Value

No return value.

Plot item/category response curves

Description

Calculate and plot item/category response curves for a fitted IRT model.

Usage

plot_item_curve(fit, ...)

Arguments

Value

Invisibly returns the plotted object of class rtmb_item_curve.

Plot item information functions

Description

Calculate and plot item information functions for a fitted IRT model.

Usage

plot_item_info(fit, ...)

Arguments

Value

Invisibly returns the plotted object of class rtmb_item_info.

Plot least-squares marginal means

Description

Calculate and plot marginal means for a fitted model.

Usage

plot_lsmeans(fit, specs, ...)

Arguments

Value

Invisibly returns the plotted object.

Plot Multidimensional Unfolding Configuration

Description

Draws an item-person map for multidimensional unfolding models. Item locations are shown as labels, optional blue circles show item alpha/radius values, and the respondent locations are summarized by a two- dimensional kernel density contour.

Usage

plot_mdu(
  delta,
  theta = NULL,
  item_alpha = NULL,
  phi = NULL,
  dims = c(1, 2),
  radius = NULL,
  signs = c(1, 1),
  item_labels = NULL,
  show_phi = TRUE,
  show_density = TRUE,
  circle_scale = 1,
  alpha = 0.2,
  contour_n = 60,
  distance = c("auto", "squared", "euclidean"),
  point_estimate = c("EAP", "MAP", "mean", "marginal_map", "joint_map"),
  main = "MDU Configuration",
  xlab = NULL,
  ylab = NULL,
  ...
)

Arguments

Value

Invisibly returns a list with plotted 'delta', 'theta', and 'item_alpha'.

Plot pairs for posterior samples

Description

Draw a scatterplot matrix to examine posterior correlations between parameters.

Usage

plot_pairs(x, pars = NULL)

Arguments

Value

No return value.

Plot test information function

Description

Calculate and plot test information function for a fitted IRT model.

Usage

plot_test_info(fit, ...)

Arguments

Value

Invisibly returns the plotted object of class rtmb_test_info.

Plot MCMC trace plots

Description

Draw trace plots for each parameter across chains.

Usage

plot_trace(x, mono = FALSE)

Arguments

Value

No return value. This function is called for its side effect of plotting.

Poisson log-probability mass function

Description

Poisson log-probability mass function

Usage

poisson_lpmf(x, mean, sum = TRUE)

Arguments

Value

The sum of the log-probability.

Positive centered triangular multivariate normal log-probability density function

Description

Positive centered triangular multivariate normal log-probability density function

Usage

positive_centered_tri_multi_normal_lpdf(x, sigma = 1)

Arguments

Value

The log-density.

Positive lower-triangular normal log-probability density function

Description

Positive lower-triangular normal log-probability density function

Usage

positive_lower_tri_normal_lpdf(x, mean = 0, sd = 1)

Arguments

Value

The log-density calculated for the non-zero lower triangular elements.

Print method for bayes_factor objects

Description

Print method for bayes_factor objects

Usage

## S3 method for class 'bayes_factor'
print(x, digits = 4, ...)

Arguments

Print method for bayes_factor_rtmb objects

Description

Print method for bayes_factor_rtmb objects

Usage

## S3 method for class 'bayes_factor_rtmb'
print(x, digits = 4, ...)

Arguments

Print method for ce_rtmb class (automatically calls plot)

Description

Print method for ce_rtmb class (automatically calls plot)

Usage

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

Arguments

Print simple effects

Description

Print simple effects

Usage

## S3 method for class 'ce_simple'
print(x, digits = 3, ...)

Arguments

print for summary_BayesRTMB class

Description

print for summary_BayesRTMB class

Usage

## S3 method for class 'summary_BayesRTMB'
print(x, digits = NULL, ...)

Arguments

Specify a flat prior

Description

'prior_flat()' specifies that no additional prior density is added by the wrapper. This is the only prior type allowed for 'classic()'.

Usage

prior_flat()

Value

A list with class '"rtmb_prior"' and 'type = "flat"'.

Specify a JZS (Jeffrey-Zellner-Siow) prior for t-tests

Description

Specify a JZS (Jeffrey-Zellner-Siow) prior for t-tests

Usage

prior_jzs(r = 0.707)

Arguments

Value

A list with class "rtmb_prior"

Specify normal/exponential priors for MAP and Bayesian inference

Description

'prior_normal()' specifies normal priors for location parameters and exponential priors for scale parameters. It is intended for MAP and Bayesian inference, not for 'classic()'.

Usage

prior_normal(
  Intercept_sd = 10,
  mu_sd = 10,
  b_sd = 10,
  sigma_rate = 5,
  tau_rate = 5,
  ...
)

Arguments

Value

A list with class '"rtmb_prior"' and 'type = "normal"'.

Specify a Regularized Horseshoe prior for continuous shrinkage

Description

Specify a Regularized Horseshoe prior for continuous shrinkage

Usage

prior_rhs(expected_vars = 3, slab_scale = 2, slab_df = 4, ...)

Arguments

Value

A list with class "rtmb_prior"

Specify a Spike-and-Slab prior for variable selection

Description

Specify a Spike-and-Slab prior for variable selection

Usage

prior_ssp(ssp_ratio = 0.25, max_beta = 1, ...)

Arguments

Value

A list with class "rtmb_prior"

Specify a flat prior

Description

Deprecated alias of 'prior_flat()'.

Usage

prior_uniform()

Value

A list with class '"rtmb_prior"' and 'type = "flat"'.

Specify a weakly informative prior

Description

Specify a weakly informative prior

Usage

prior_weak(sd_ratio = 0.5, max_beta = 1, ...)

Arguments

Value

A list with class "rtmb_prior"

Quadratic form using a Cholesky factor

Description

Quadratic form using a Cholesky factor

Usage

quad_form_chol(x, L)

Arguments

Value

The value of the quadratic form.

Quadratic form with a diagonal matrix

Description

Quadratic form with a diagonal matrix

Usage

quad_form_diag(m, v)

Arguments

Value

The value of the quadratic form.

Calculate 95% Quantiles

Description

Calculate 95% Quantiles

Usage

quantile95(x)

Arguments

Value

A numeric vector of length 2.

Calculate Rank-normalized Split-R-hat

Description

Calculate Rank-normalized Split-R-hat

Usage

r_hat(sims)

Arguments

Value

A numeric value.

Restore MCMC Fit from CSV

Description

Restore MCMC Fit from CSV

Usage

read_mcmc_csv(model, name, dir = "BayesRTMB_mcmc", chains = 4, laplace = FALSE)

Arguments

Value

An MCMC_Fit object.

Define an RTMB Model with Stan-like Syntax

Description

rtmb_code is the primary interface for defining Bayesian or Frequentist models within the RTMB framework. It allows you to organize model logic into functional blocks while utilizing a Stan-like ~ operator for specifying priors and likelihoods.

Usage

rtmb_code(...)

Arguments

Details

Key Blocks and Logic:

• setup: Pre-compilation data processing.

• parameters: Declaration of estimated parameters. See parameter_types for available constraints.

• transform: Definition of intermediate variables. Use math_functions for stable AD calculations.

• model: Likelihood and Priors. Refer to distributions for available sampling statements.

• generate: Post-hoc calculation of predictions or diagnostics.

Automatic Differentiation (AD) Note: All code within parameters, transform, and model is automatically differentiated by RTMB. To ensure numerical stability, use provided utility functions like log_sum_exp or inv_logit instead of raw algebraic implementations.

Value

A list of unevaluated code blocks.

See Also

rtmb_model for building the model object. parameter_types for constraining parameters. distributions for likelihood functions. math_functions for numerical utilities.

Examples

# Simulate data for a linear regression
set.seed(123)
N <- 100
x <- rnorm(N)
y <- 2.0 + 1.5 * x + rnorm(N, mean = 0, sd = 1)
data_list <- list(N = N, x = x, y = y)

# Define the model using rtmb_code
code <- rtmb_code(
  setup = {
    # Center the predictor variable (executed once)
    x_centered <- x - mean(x)
  },
  parameters = {
    # Define parameters and their constraints
    alpha = Dim(1)
    beta  = Dim(1)
    sigma = Dim(1, lower = 0)
  },
  transform = {
    # Calculate the linear predictor
    mu <- alpha + beta * x_centered
  },
  model = {
    # Priors
    alpha ~ normal(0, 10)
    beta  ~ normal(0, 10)
    sigma ~ exponential(1)

    # Likelihood (Vectorized)
    y ~ normal(mu, sigma)
  },
  generate = {
    # Calculate generated quantities
    y_pred <- mu

    # Must return a named list
    list(y_pred = y_pred)
  }
)

# Create the model object
mod <- rtmb_model(data = data_list, code = code)

# Fit the model using MAP estimation
map_res <- mod$optimize()
map_res$summary(pars = c("alpha", "beta", "sigma"))

# The generated quantity 'y_pred' can also be summarized
map_res$summary("y_pred", max_rows = 5)

Fit a Correlation Model using RTMB

Description

'rtmb_corr' fits a correlation model to estimate means, standard deviations, and correlation structures. It supports simple correlation, multilevel correlation, and classical frequentist estimation.

Usage

rtmb_corr(
  x = NULL,
  data = NULL,
  ID = NULL,
  covariates = NULL,
  method = c("pearson", "spearman", "reml"),
  prior = prior_flat(),
  y_range = NULL,
  init = NULL,
  fixed = NULL,
  missing = c("listwise", "fiml", "pairwise"),
  WAIC = FALSE
)

Arguments

Examples

  # Estimate the correlation between two variables in the debate dataset
  data(debate, package = "BayesRTMB")

  fit_corr <- rtmb_corr(cbind(sat, perf), data = debate)

  
  mcmc_corr <- fit_corr$sample(sampling = 500, warmup = 500, chains = 2)
  mcmc_corr$summary()

  bf_corr <- mcmc_corr$bayes_factor(fixed = list(corr = 0))
  print(bf_corr)
  

RTMB-based Factor Analysis Wrapper

Description

Performs exploratory factor analysis (EFA) using RTMB. It supports standard rotation methods (e.g., varimax, promax) as well as regularized factor analysis using a Spike-and-Slab Prior (SSP) for estimating sparse loading matrices.

Usage

rtmb_fa(
  data,
  nfactors = 1,
  rotate = NULL,
  score = FALSE,
  prior = prior_flat(),
  y_range = NULL,
  init = NULL,
  fixed = NULL,
  missing = c("listwise", "fiml"),
  WAIC = FALSE
)

Arguments

Examples

  # Prepare a subset of the BigFive dataset for factor analysis
  data(BigFive, package = "BayesRTMB")
  fa_data <- BigFive[, 1:10]

  # --- 1. Standard Exploratory Factor Analysis (1 Factor) ---
  fit_fa1 <- rtmb_fa(data = fa_data, nfactors = 1)

  # Maximum A Posteriori (MAP) estimation
  map_fa1 <- fit_fa1$optimize()
  map_fa1$summary()



  # --- 2. Factor Analysis with Rotation and Factor Scores (2 Factors) ---
  # Extract 2 factors, apply Promax rotation during model fitting, and calculate factor scores
  fit_fa2 <- rtmb_fa(data = fa_data, nfactors = 2, rotate = "promax", score = TRUE)

  # Setting se_method = "sampling" enables standard errors and 95% CIs
  # for transformed and generated quantities, such as factor scores and post-hoc rotations.
  # (It uses multivariate normal sampling from the unconstrained parameter space)
  map_fa2 <- fit_fa2$optimize(se_method = "sampling")
  map_fa2$summary()

  # Post-hoc rotation using the fa_rotate() method
  # You can also apply other rotation methods (e.g., "varimax" from the stats package)
  # to the unrotated loading matrix ("L") after estimation.
  map_fa2$fa_rotate(target = "L", rotate = "varimax")

  # The post-hoc rotated loadings are automatically stored with the method's
  # suffix (e.g., "L_varimax")
  map_fa2$summary("L_varimax")



  # --- 3. Regularized Factor Analysis using Spike-and-Slab Prior (SSP) ---
  # Specifying 'rotate = "ssp"' enables sparse loading matrix estimation
  fit_ssp <- rtmb_fa(data = fa_data, nfactors = 2, rotate = "ssp")

  map_ssp <- fit_ssp$optimize()
  map_ssp$summary()

RTMB-based GLM wrapper function (no random effects)

Description

RTMB-based GLM wrapper function (no random effects)

Usage

rtmb_glm(
  formula,
  data,
  family = "gaussian",
  prior = prior_flat(),
  y_range = NULL,
  init = NULL,
  fixed = NULL,
  gmc = NULL,
  centering = NULL,
  view = NULL,
  factors = NULL,
  contrasts = "treatment",
  missing = c("listwise", "fiml"),
  WAIC = FALSE,
  ...
)

Arguments

Examples

  # --- 1. Linear Regression (rtmb_lm) ---
  # Fit a linear regression model using the debate dataset
  data(debate, package = "BayesRTMB")
  fit_lm <- rtmb_lm(sat ~ talk + perf, data = debate)
  map_lm <- fit_lm$optimize()
  map_lm$summary()

  # --- 2. Generalized Linear Model (rtmb_glm) ---
  # Fit a logistic regression model using the debate dataset
  data(debate, package = "BayesRTMB")
  fit_glm <- rtmb_glm(cond ~ talk + sat, data = debate, family = "bernoulli")
  map_glm <- fit_glm$optimize()
  map_glm$summary()

  # --- 3. Generalized Linear Mixed Model (rtmb_glmer) ---
  # Fit a linear mixed-effects model using the debate dataset
  data(debate, package = "BayesRTMB")
  fit_glmer <- rtmb_glmer(talk ~ cond + (1 | group), data = debate, family = "gaussian")

  # MAP estimation using Laplace approximation for random effects
  map_glmer <- fit_glmer$optimize(laplace = TRUE)
  map_glmer$summary()

  # MCMC sampling (chains and iterations reduced for faster execution)
  
  mcmc_glmer <- fit_glmer$sample(sampling = 500, warmup = 500, chains = 2)
  mcmc_glmer$summary()
  

  # --- 4. Linear Mixed Model (rtmb_lmer) ---
  # A convenient wrapper for Gaussian mixed models (identical to rtmb_glmer with family="gaussian")
  fit_lmer <- rtmb_lmer(sat ~ talk + (1 | group), data = debate)
  map_lmer <- fit_lmer$optimize()
  map_lmer$summary()

  # --- 5. Regularized Regression (Variable Selection) ---
  # You can apply regularization to the fixed effects to shrink noise variables towards zero.
  # Use prior = prior_rhs() for the Regularized Horseshoe prior,
  # or prior_ssp() for the Spike-and-Slab prior.
  # Note: When using regularization, you must specify 'y_range' (the theoretical minimum and maximum
  # values of the response variable) to automatically set up the required weakly informative priors.

  # Fit a linear regression using debate predictors with the Horseshoe prior
  fit_rhs <- rtmb_lm(
    sat ~ talk + perf + skill,
    data = debate,
    prior = prior_rhs(),
    y_range = c(1, 5)
  )
  map_rhs <- fit_rhs$optimize()
  # Summarize only the fixed effects (slopes)
  map_rhs$summary("b")

  # Fit a linear regression with the Spike-and-Slab prior
  fit_ssp <- rtmb_lm(
    sat ~ talk + perf + skill,
    data = debate,
    prior = prior_ssp(),
    y_range = c(1, 5)
  )
  map_ssp <- fit_ssp$optimize()
  map_ssp$summary("b")

RTMB-based GLMM wrapper function

Description

RTMB-based GLMM wrapper function

Usage

rtmb_glmer(
  formula,
  data,
  family = "gaussian",
  laplace = FALSE,
  prior = prior_flat(),
  y_range = NULL,
  init = NULL,
  fixed = NULL,
  gmc = NULL,
  centering = NULL,
  cwc = NULL,
  view = NULL,
  within = NULL,
  factors = NULL,
  contrasts = "treatment",
  sigma_by = NULL,
  resid_corr = NULL,
  resid_time = NULL,
  resid_group = NULL,
  generate = NULL,
  missing = c("listwise", "fiml"),
  WAIC = FALSE,
  .force_sum = FALSE
)

Arguments

Examples

  # --- 1. Linear Regression (rtmb_lm) ---
  # Fit a linear regression model using the debate dataset
  data(debate, package = "BayesRTMB")
  fit_lm <- rtmb_lm(sat ~ talk + perf, data = debate)
  map_lm <- fit_lm$optimize()
  map_lm$summary()

  # --- 2. Generalized Linear Model (rtmb_glm) ---
  # Fit a logistic regression model using the debate dataset
  data(debate, package = "BayesRTMB")
  fit_glm <- rtmb_glm(cond ~ talk + sat, data = debate, family = "bernoulli")
  map_glm <- fit_glm$optimize()
  map_glm$summary()

  # --- 3. Generalized Linear Mixed Model (rtmb_glmer) ---
  # Fit a linear mixed-effects model using the debate dataset
  data(debate, package = "BayesRTMB")
  fit_glmer <- rtmb_glmer(talk ~ cond + (1 | group), data = debate, family = "gaussian")

  # MAP estimation using Laplace approximation for random effects
  map_glmer <- fit_glmer$optimize(laplace = TRUE)
  map_glmer$summary()

  # MCMC sampling (chains and iterations reduced for faster execution)
  
  mcmc_glmer <- fit_glmer$sample(sampling = 500, warmup = 500, chains = 2)
  mcmc_glmer$summary()
  

  # --- 4. Linear Mixed Model (rtmb_lmer) ---
  # A convenient wrapper for Gaussian mixed models (identical to rtmb_glmer with family="gaussian")
  fit_lmer <- rtmb_lmer(sat ~ talk + (1 | group), data = debate)
  map_lmer <- fit_lmer$optimize()
  map_lmer$summary()

  # --- 5. Regularized Regression (Variable Selection) ---
  # You can apply regularization to the fixed effects to shrink noise variables towards zero.
  # Use prior = prior_rhs() for the Regularized Horseshoe prior,
  # or prior_ssp() for the Spike-and-Slab prior.
  # Note: When using regularization, you must specify 'y_range' (the theoretical minimum and maximum
  # values of the response variable) to automatically set up the required weakly informative priors.

  # Fit a linear regression using debate predictors with the Horseshoe prior
  fit_rhs <- rtmb_lm(
    sat ~ talk + perf + skill,
    data = debate,
    prior = prior_rhs(),
    y_range = c(1, 5)
  )
  map_rhs <- fit_rhs$optimize()
  # Summarize only the fixed effects (slopes)
  map_rhs$summary("b")

  # Fit a linear regression with the Spike-and-Slab prior
  fit_ssp <- rtmb_lm(
    sat ~ talk + perf + skill,
    data = debate,
    prior = prior_ssp(),
    y_range = c(1, 5)
  )
  map_ssp <- fit_ssp$optimize()
  map_ssp$summary("b")

RTMB-based IRT (Item Response Theory) Wrapper

Description

Performs Item Response Theory modeling (1PL, 2PL, 3PL, and Graded Response Model) using RTMB. Missing values (NA) in the data are automatically removed internally for efficient computation.

Usage

rtmb_irt(
  data,
  model = c("2PL", "1PL", "3PL"),
  type = c("binary", "ordered"),
  prior = prior_flat(),
  init = NULL,
  fixed = NULL,
  view = NULL,
  missing = c("listwise", "fiml"),
  WAIC = FALSE
)

Arguments

Examples

  # --- 1. Binary Data (e.g., correct/incorrect answers) ---
  # Simulate binary response data for 100 persons and 5 items
  set.seed(123)
  bin_data <- matrix(rbinom(500, size = 1, prob = 0.6), nrow = 100, ncol = 5)
  colnames(bin_data) <- paste0("Item", 1:5)

  # Introduce some missing values (NA) to demonstrate automatic handling
  #bin_data[sample(1:500, 10)] <- NA

  # Fit a 2-Parameter Logistic (2PL) model
  fit_2pl <- rtmb_irt(data = bin_data, model = "2PL", type = "binary")

  # Maximum A Posteriori (MAP) estimation
  map_2pl <- fit_2pl$optimize()
  map_2pl$summary()

  # --- 2. Ordered Data (e.g., Likert scale) ---
  # Simulate ordered response data (categories 1 to 5) with a true latent structure
  set.seed(123)
  N <- 100
  J <- 5
  K <- 4

  # True parameters
  theta <- rnorm(N, 0, 1)          # Person abilities (latent traits)
  a <- runif(J, 0.8, 2.0)          # Item discriminations

  # Item-specific thresholds (must be strictly increasing)
  b <- matrix(NA, nrow = J, ncol = K - 1)
  for (j in 1:J) {
    b[j, ] <- sort(c(-1.5, 0, 1.5) + rnorm(3, 0, 0.2))
  }

  ord_data <- matrix(NA, nrow = N, ncol = J)
  colnames(ord_data) <- paste0("Item", 1:J)

  # Generate responses based on the Graded Response Model
  for (i in 1:N) {
    for (j in 1:J) {
      # Latent continuous response (eta + standard logistic noise)
      eta <- a[j] * theta[i]
      y_star <- eta + rlogis(1)

      # Apply thresholds to determine the observed category (1 to 4)
      category <- 1
      for (k in 1:(K - 1)) {
        if (y_star > b[j, k]) category <- k + 1
      }
      ord_data[i, j] <- category
    }
  }

  # Fit a Graded Response Model (2PL for ordered data)
  fit_ord <- rtmb_irt(data = ord_data, model = "2PL", type = "ordered")
  fit_ord$print_code()

  map_ord <- fit_ord$optimize()
  map_ord$summary()

RTMB-based Linear Regression wrapper function

Description

RTMB-based Linear Regression wrapper function

Usage

rtmb_lm(
  formula,
  data,
  prior = prior_flat(),
  y_range = NULL,
  init = NULL,
  fixed = NULL,
  gmc = NULL,
  centering = NULL,
  view = NULL,
  factors = NULL,
  contrasts = "treatment",
  missing = c("listwise", "fiml"),
  WAIC = FALSE,
  ...
)

Arguments

Examples

  # --- 1. Linear Regression (rtmb_lm) ---
  # Fit a linear regression model using the debate dataset
  data(debate, package = "BayesRTMB")
  fit_lm <- rtmb_lm(sat ~ talk + perf, data = debate)
  map_lm <- fit_lm$optimize()
  map_lm$summary()

  # --- 2. Generalized Linear Model (rtmb_glm) ---
  # Fit a logistic regression model using the debate dataset
  data(debate, package = "BayesRTMB")
  fit_glm <- rtmb_glm(cond ~ talk + sat, data = debate, family = "bernoulli")
  map_glm <- fit_glm$optimize()
  map_glm$summary()

  # --- 3. Generalized Linear Mixed Model (rtmb_glmer) ---
  # Fit a linear mixed-effects model using the debate dataset
  data(debate, package = "BayesRTMB")
  fit_glmer <- rtmb_glmer(talk ~ cond + (1 | group), data = debate, family = "gaussian")

  # MAP estimation using Laplace approximation for random effects
  map_glmer <- fit_glmer$optimize(laplace = TRUE)
  map_glmer$summary()

  # MCMC sampling (chains and iterations reduced for faster execution)
  
  mcmc_glmer <- fit_glmer$sample(sampling = 500, warmup = 500, chains = 2)
  mcmc_glmer$summary()
  

  # --- 4. Linear Mixed Model (rtmb_lmer) ---
  # A convenient wrapper for Gaussian mixed models (identical to rtmb_glmer with family="gaussian")
  fit_lmer <- rtmb_lmer(sat ~ talk + (1 | group), data = debate)
  map_lmer <- fit_lmer$optimize()
  map_lmer$summary()

  # --- 5. Regularized Regression (Variable Selection) ---
  # You can apply regularization to the fixed effects to shrink noise variables towards zero.
  # Use prior = prior_rhs() for the Regularized Horseshoe prior,
  # or prior_ssp() for the Spike-and-Slab prior.
  # Note: When using regularization, you must specify 'y_range' (the theoretical minimum and maximum
  # values of the response variable) to automatically set up the required weakly informative priors.

  # Fit a linear regression using debate predictors with the Horseshoe prior
  fit_rhs <- rtmb_lm(
    sat ~ talk + perf + skill,
    data = debate,
    prior = prior_rhs(),
    y_range = c(1, 5)
  )
  map_rhs <- fit_rhs$optimize()
  # Summarize only the fixed effects (slopes)
  map_rhs$summary("b")

  # Fit a linear regression with the Spike-and-Slab prior
  fit_ssp <- rtmb_lm(
    sat ~ talk + perf + skill,
    data = debate,
    prior = prior_ssp(),
    y_range = c(1, 5)
  )
  map_ssp <- fit_ssp$optimize()
  map_ssp$summary("b")

RTMB-based Linear Mixed Model (LMM) wrapper function

Description

RTMB-based Linear Mixed Model (LMM) wrapper function

Usage

rtmb_lmer(
  formula,
  data,
  laplace = TRUE,
  prior = prior_flat(),
  y_range = NULL,
  init = NULL,
  fixed = NULL,
  gmc = NULL,
  centering = NULL,
  cwc = NULL,
  view = NULL,
  sigma_by = NULL,
  factors = NULL,
  contrasts = "treatment",
  resid_corr = NULL,
  resid_time = NULL,
  resid_group = NULL,
  within = NULL,
  missing = c("listwise", "fiml"),
  WAIC = FALSE,
  ...
)

Arguments

Value

RTMB_Model object

Examples

  # --- 1. Linear Regression (rtmb_lm) ---
  # Fit a linear regression model using the debate dataset
  data(debate, package = "BayesRTMB")
  fit_lm <- rtmb_lm(sat ~ talk + perf, data = debate)
  map_lm <- fit_lm$optimize()
  map_lm$summary()

  # --- 2. Generalized Linear Model (rtmb_glm) ---
  # Fit a logistic regression model using the debate dataset
  data(debate, package = "BayesRTMB")
  fit_glm <- rtmb_glm(cond ~ talk + sat, data = debate, family = "bernoulli")
  map_glm <- fit_glm$optimize()
  map_glm$summary()

  # --- 3. Generalized Linear Mixed Model (rtmb_glmer) ---
  # Fit a linear mixed-effects model using the debate dataset
  data(debate, package = "BayesRTMB")
  fit_glmer <- rtmb_glmer(talk ~ cond + (1 | group), data = debate, family = "gaussian")

  # MAP estimation using Laplace approximation for random effects
  map_glmer <- fit_glmer$optimize(laplace = TRUE)
  map_glmer$summary()

  # MCMC sampling (chains and iterations reduced for faster execution)
  
  mcmc_glmer <- fit_glmer$sample(sampling = 500, warmup = 500, chains = 2)
  mcmc_glmer$summary()
  

  # --- 4. Linear Mixed Model (rtmb_lmer) ---
  # A convenient wrapper for Gaussian mixed models (identical to rtmb_glmer with family="gaussian")
  fit_lmer <- rtmb_lmer(sat ~ talk + (1 | group), data = debate)
  map_lmer <- fit_lmer$optimize()
  map_lmer$summary()

  # --- 5. Regularized Regression (Variable Selection) ---
  # You can apply regularization to the fixed effects to shrink noise variables towards zero.
  # Use prior = prior_rhs() for the Regularized Horseshoe prior,
  # or prior_ssp() for the Spike-and-Slab prior.
  # Note: When using regularization, you must specify 'y_range' (the theoretical minimum and maximum
  # values of the response variable) to automatically set up the required weakly informative priors.

  # Fit a linear regression using debate predictors with the Horseshoe prior
  fit_rhs <- rtmb_lm(
    sat ~ talk + perf + skill,
    data = debate,
    prior = prior_rhs(),
    y_range = c(1, 5)
  )
  map_rhs <- fit_rhs$optimize()
  # Summarize only the fixed effects (slopes)
  map_rhs$summary("b")

  # Fit a linear regression with the Spike-and-Slab prior
  fit_ssp <- rtmb_lm(
    sat ~ talk + perf + skill,
    data = debate,
    prior = prior_ssp(),
    y_range = c(1, 5)
  )
  map_ssp <- fit_ssp$optimize()
  map_ssp$summary("b")

RTMB-based Log-linear analysis (Poisson regression)

Description

Performs Bayesian or Frequentist log-linear analysis (Poisson regression) on a contingency table or raw data.

Usage

rtmb_loglinear(
  formula,
  data,
  prior = prior_flat(),
  fixed = NULL,
  WAIC = FALSE,
  ...
)

Arguments

Value

An 'RTMB_Model' object.

Examples

  # Create a contingency table
  tab <- matrix(c(10, 20, 30, 40), nrow = 2)
  dimnames(tab) <- list(A = c("A1", "A2"), B = c("B1", "B2"))
  
  # Fit a log-linear model (independence model: ~ A + B)
  fit_log <- rtmb_loglinear(~ A + B, data = tab)
  
  # MAP estimation
  map_log <- fit_log$optimize()
  map_log$summary()

Fit a Latent Rank Theory (LRT) Model

Description

Fits a Latent Rank Theory model, which is a mixture model with ordered ranks and Gaussian Process smoothing on the mean profiles.

Usage

rtmb_lrt(
  formula,
  k = 3,
  data = NULL,
  rank_coords = NULL,
  covariance = c("diagonal", "diagonal_equal", "full", "full_equal", "full_equal_corr"),
  magnitude = NULL,
  smoothing = NULL,
  noise = 0.01,
  prob_smoothing = FALSE,
  link = c("ordered", "sequential"),
  prior = prior_flat(),
  y_range = NULL,
  fixed = NULL,
  two_stage = FALSE,
  WAIC = FALSE,
  ...
)

Arguments

Value

A RTMB_Model object.

RTMB-based Multidimensional Unfolding Wrapper

Description

Fits a multidimensional unfolding model for preference/rating data. Rows are persons or observations and columns are stimuli/items. The model represents both row scores ('theta') and item locations ('delta') in a shared D-dimensional space.

Usage

rtmb_mdu(
  data,
  ndim = 2,
  distance = c("squared", "euclidean"),
  alpha = c("random", "fix"),
  method = c("rating", "Best", "Best-Worst", "MDS"),
  sets = NULL,
  prior = prior_flat(),
  y_range = NULL,
  init = NULL,
  fixed = NULL,
  view = NULL,
  distance_eps = 1e-04,
  missing = c("listwise", "fiml"),
  WAIC = FALSE
)

Arguments

Value

An 'RTMB_Model' object.

Examples

  # Simulate rating data for Multidimensional Unfolding (MDU)
  set.seed(123)
  N <- 50  # Number of persons
  M <- 10  # Number of items
  D <- 2   # Number of dimensions
  
  # True person and item coordinates in a 2D space
  theta <- matrix(rnorm(N * D), N, D)
  delta <- matrix(rnorm(M * D), M, D)
  
  # Generate distance-like ratings (smaller means more preferred)
  Y <- matrix(NA, N, M)
  for(i in 1:N) {
    for(j in 1:M) {
      Y[i, j] <- sum((theta[i,] - delta[j,])^2) + rnorm(1, 0, 0.5)
    }
  }

  # Fit a 2-dimensional MDU model
  fit_mdu <- rtmb_mdu(Y, ndim = 2, distance = "squared", method = "rating")
  
  # MAP estimation
  map_mdu <- fit_mdu$optimize()
  map_mdu$summary()
  
  # Note: MDU models have many parameters, so MCMC sampling might take time.
  
  # mcmc_mdu <- fit_mdu$sample(sampling = 500, warmup = 500, chains = 2)
  

RTMB-based Mediation Analysis Wrapper

Description

'rtmb_mediation' performs mediation analysis by simultaneously estimating multiple GLM regression equations. It automatically identifies mediation paths and calculates indirect, direct, and total effects.

Usage

rtmb_mediation(
  formula,
  data,
  family = "gaussian",
  prior = prior_flat(),
  y_range = NULL,
  fixed = NULL,
  view = NULL,
  WAIC = FALSE,
  ...
)

Arguments

Details

The function identifies mediation paths by looking for variables that are responses in one equation and predictors in another. Indirect effects are calculated as the product of coefficients along these paths (a * b).

Uncertainty Estimation: When using '$optimize(ci_method = "sampling")', the function provides asymmetric confidence intervals for indirect effects based on the distribution of products, which is more accurate than the standard Sobel test (Delta Method).

Value

An 'RTMB_Model' object.

Examples

  # Simulate mediation data
  set.seed(123)
  N <- 100
  X <- rnorm(N)
  # M is influenced by X
  M <- 0.5 * X + rnorm(N, 0, 0.5)
  # Y is influenced by both X and M
  Y <- 0.3 * X + 0.8 * M + rnorm(N, 0, 0.5)
  dat <- data.frame(X, M, Y)

  # Fit a mediation model
  # The formula list specifies the two regression equations
  fit_med <- rtmb_mediation(list(M ~ X, Y ~ X + M), data = dat)
  
  # Maximum A Posteriori (MAP) estimation
  map_med <- fit_med$optimize()
  # The summary automatically calculates indirect, direct, and total effects
  map_med$summary()

Mixture Model Wrapper for RTMB

Description

Provides a user-friendly interface for fitting Gaussian mixture models with optional covariates on class membership probabilities and various covariance structures.

Usage

rtmb_mixture(
  formula,
  k = 2,
  data = NULL,
  covariance = c("diagonal", "diagonal_equal", "full", "full_equal", "full_equal_corr"),
  prior = prior_flat(),
  y_range = NULL,
  fixed = NULL,
  WAIC = FALSE,
  ...
)

Arguments

Value

A RTMB_Model object.