Skip to contents

Fit a Bayesian Mixed-Effects Discounting Model via brms

Source: R/dd-brms.R
fit_dd_brms.Rd

Fits the four discounting equations of fit_dd_tmb() ("mazur", "exponential", "green-myerson", "rachlin") as Bayesian nonlinear mixed-effects models, with the discount rate estimated on the natural-log scale (logk, subject random intercept; k ~ 1 by default, or k + phi ~ 1 to add a per-subject precision random effect) and the shape exponent (logs) population-level for the two-parameter equations.

Usage

fit_dd_brms(
  data,
  y_var = "y",
  x_var = "x",
  id_var = "id",
  equation = c("mazur", "exponential", "green-myerson", "rachlin"),
  family = c("beta", "gaussian"),
  boundary = c("error", "squeeze", "zoib"),
  random_effects = k ~ 1,
  covariance_structure = c("pdSymm", "pdDiag"),
  factors = NULL,
  factor_interaction = FALSE,
  continuous_covariates = NULL,
  ll = NULL,
  response_scale = c("proportion", "percent", "amount"),
  prior = NULL,
  autoscale_priors = TRUE,
  chains = 4,
  iter = 2000,
  warmup = floor(iter/2),
  thin = 1,
  cores = getOption("mc.cores", 1L),
  seed = NA,
  backend = getOption("brms.backend", "rstan"),
  control = list(adapt_delta = 0.95),
  init = c("prior_center", "tmb", "random", "0"),
  sample_prior = "no",
  loo = TRUE,
  file = NULL,
  file_refit = getOption("brms.file_refit", "on_change"),
  verbose = 1,
  ...
)

Arguments

data

Long-format data frame (one row per subject-delay).

y_var, x_var, id_var

Column names (canonical defaults), as in fit_dd_tmb().

equation

Discounting equation.

family

"beta" or "gaussian"; "sltb" errors with guidance.

boundary

Boundary handling for the beta family: "error" (default), "squeeze", or "zoib" (see Details). Ignored for family = "gaussian".

random_effects

k ~ 1 (single random intercept on log k, the default) or k + phi ~ 1 (adds a per-subject precision random effect; beta family only). The latter mirrors fit_dd_tmb(random_effects = k + phi ~ 1): phi becomes a predicted distributional parameter on the log link carrying a subject random intercept.

covariance_structure

For k + phi ~ 1, the (log k, log phi) covariance: "pdSymm" (default; correlated, the TMB parity choice) or "pdDiag" (independent). Ignored for k ~ 1.

factors, factor_interaction, continuous_covariates

Fixed-effect design on logk, as in fit_dd_tmb().

ll, response_scale

Response coercion, as in fit_dd_tmb().

prior

Optional brmsprior; user rows override the defaults.

autoscale_priors

Anchor the logk prior to the median delay (see default_dd_priors()). For equation = "rachlin" this also fits on the normalised delay x / median(x) (recorded as param_info$delay_scale) so the prior is delay-unit invariant; reported draws are back-transformed to data-unit k, but the raw brmsfit draws of the logk intercept stay normalised.

chains, iter, warmup, thin, cores, seed, backend, control, sample_prior

MCMC settings passed to brms::brm().

init

"prior_center" (default), "tmb" (a quiet fit_dd_tmb() pre-fit supplies the centers, with prior_center fallback on failure; the beta family maps to the TMB sltb pre-fit), "random", or "0"; or a list/function passed through to brms::brm().

loo

Compute and store brms::loo() at fit time.

file, file_refit

Passed to brms::brm() for fit caching.

verbose

0 (silent), 1 (messages), 2 (full Stan output).

...

Passed through to brms::brm().

Value

An object of class beezdiscounting_brms: model (posterior medians/SDs on the estimation scale under TMB names beta_k/log_s, plus variance_components), brmsfit, subject_pars (id, k, k_lower, k_upper, plus phi/phi_lower/phi_upper for k + phi ~ 1), converged (Rhat < 1.01, no divergences, bulk ESS >= 400), mcmc_info, loo, data, param_info, priors, autoscale_info.

Details

family = "beta" (default) uses Beta(link = "identity") with the mean squished into (1e-6, 1 - 1e-6). This is a different likelihood from the TMB scale-location-truncated beta (family = "sltb", which has no brms equivalent and errors with this pointer): the ordinary beta density is undefined at exactly 0 and 1, so its estimates are not expected to match fit_dd_tmb(family = "sltb") whenever responses sit at or near the boundaries. Boundary observations (y exactly 0 or 1) are handled per boundary: "error" (default) refuses to fit and reports how many there are; "squeeze" applies the Smithson-Verkuilen transform y* = (y (N - 1) + 0.5) / N to every response (an explicit opt-in: it moves each value toward 0.5 and floors the response at 0.5 / N, so when the fitted curve approaches that floor – long delays, steep discounting – the squeezed values pull k downward); "zoib" switches to zero_one_inflated_beta (changes the estimand – k then describes interior responses only). The fraction of exact-boundary and near-boundary (y <= 0.01 or y >= 0.99) responses is stored in param_info$boundary_info and printed by summary(). family = "gaussian" matches fit_dd_tmb(family = "gaussian") wherever the TMB template's mu clamp into [1e-6, 1 - 1e-6] does not bind (everywhere except extreme decay underflow).

See also

fit_dd_tmb(); default_dd_priors().