Package 'beezdiscounting'

Description

The paper's F-test is exact when the random-effects variance estimate g-hat > 0 (Prop. 3.5's assumption); when g-hat = 0 the statistic is still reported, with a warning.

Usage

## S3 method for class 'beezdiscounting_linear'
anova(object, hypothesis = NULL, pairwise = FALSE, ...)

Arguments

Value

Tibble of class beezdiscounting_linear_anova.

See Also

beezdiscounting_linear-methods for the other S3 methods.

Description

Converts answers from 5.5 trial delay discounting from Qualtrics template

Usage

ans_dd(df)

Arguments

Value

A dataframe with the ResponseId, index, and response (ss or ll).

Examples

ans_dd(five.fivetrial_dd)

Description

Converts answers from 5.5 trial probability discounting from Qualtrics template

Usage

ans_pd(df)

Arguments

Value

A dataframe with the ResponseId, index, and response (sc or lu).

Examples

ans_pd(five.fivetrial_pd)

Description

Returns the data used for fitting (or newdata) as a tibble with three diagnostic columns appended:

• .fitted: subject-conditional fitted choice probability P(LL).

• .resid: raw (response) residual choice - .fitted.

• .std_resid: Pearson residual .resid / sqrt(p * (1 - p)).

Usage

## S3 method for class 'beezdiscounting_choice'
augment(x, newdata = NULL, ...)

Arguments

Value

A tibble with the same rows as the data plus .fitted, .resid, and .std_resid.

Description

Returns the data used for fitting (or newdata) as a tibble with three diagnostic columns appended:

• .fitted: subject-conditional fitted indifference proportion (clamped to ⁠(0, 1)⁠).

• .resid: raw residual y - .fitted on the response scale.

• .std_resid: Pearson (standardized) residual, i.e., .resid divided by the per-row response SD. For family = "gaussian" the SD is the constant sigma_e; for family = "sltb" it is the delta-method SLT SD sqrt(mu * (1 - mu) / (phi + 1)).

Usage

## S3 method for class 'beezdiscounting_tmb'
augment(x, newdata = NULL, ...)

Arguments

Value

A tibble with the same rows as the data plus .fitted, .resid, and .std_resid.

Examples

dd <- simulate_dd_ip(n_subjects = 20, seed = 1)
fit <- fit_dd_tmb(dd, verbose = 0)
head(augment(fit))

Description

S3 methods for objects of class beezdiscounting_linear returned by fit_dd_linear(). All quantities are closed form; nothing is re-optimized.

Usage

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

## S3 method for class 'beezdiscounting_linear'
summary(object, ...)

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

## S3 method for class 'beezdiscounting_linear'
tidy(x, effects = c("subject", "population"), ...)

## S3 method for class 'beezdiscounting_linear'
glance(x, ...)

## S3 method for class 'beezdiscounting_linear'
coef(object, ...)

## S3 method for class 'beezdiscounting_linear'
confint(object, parm = c("population", "subject"), level = NULL, ...)

## S3 method for class 'beezdiscounting_linear'
augment(x, ...)

## S3 method for class 'beezdiscounting_linear'
predict(object, type = "parameters", ...)

## S3 method for class 'beezdiscounting_linear'
nobs(object, ...)

## S3 method for class 'beezdiscounting_linear'
logLik(
  object,
  level = c("population", "subject"),
  scale = c("raw", "transformed"),
  ...
)

Arguments

Details

confint() uses parm to SELECT WHICH INTERVAL is returned, not to filter parameter names as in stats::confint(). logLik() reports the raw-scale (Jacobian-corrected) log-likelihood by default, which makes it comparable with this package's Gaussian NLS and TMB log-likelihoods.

Value

print() and print.summary() return their input invisibly; summary() returns an object of class summary.beezdiscounting_linear; tidy(), glance(), augment() and predict() return tibbles. augment() adds .fitted and .resid on the raw indifference-point scale and .std_resid, the transformed-scale residual ⁠y_lin - ln k_i⁠ divided by the model's error standard deviation (sqrt(sigma2) from the random-effects fit, or the pooled within-subject mean square when that component is unavailable); .std_resid is NA for points dropped under boundary = "drop". coef() returns the named vector of condition means mu; confint() returns a two-column matrix of lower/upper bounds; nobs() returns the number of usable transformed observations; logLik() returns a "logLik" object with df and nobs attributes. Methods needing the random-effects component error when the design was unbalanced and re is NULL.

See Also

fit_dd_linear(), anova.beezdiscounting_linear(), plot.beezdiscounting_linear()

Description

Calculates three area-under-the-curve (AUC) metrics for delay-discounting data: regular AUC (using raw delays), log10 AUC (using logarithmically scaled delays), and ordinal AUC (using ordinally scaled delays). The three differ only in how the delays are scaled before the area is computed. Metrics are computed separately for each id, so a data frame with several subjects returns one row per subject.

Usage

calc_aucs(dat)

Arguments

Value

A tibble with the following columns:

• id: The participant or group identifier.

• auc_regular: The regular AUC, calculated using the raw delay values.

• auc_log10: The log10 AUC, calculated using logarithmically transformed delay values.

• auc_ord: The ordinal AUC, calculated using ordinally scaled delay values.

Examples

# Example data
data <- data.frame(
  id = rep("P1", 6),
  x = c(1, 7, 30, 90, 180, 365),
  y = c(0.8, 0.5, 0.3, 0.2, 0.1, 0.05)
)

# Calculate AUC metrics for a single participant
calc_aucs(data)

Description

Computes the lower and upper bounds of a t-based confidence interval for a parameter estimate from its standard error, a specified significance level, and the degrees of freedom from the model.

Usage

calc_conf_int(estimate, std_error, model, alpha = 0.05)

Arguments

Value

A numeric vector of length two:

• First element: Lower bound of the confidence interval.

• Second element: Upper bound of the confidence interval.

Examples

# Example using a linear model
data <- data.frame(x = 1:10, y = c(2.3, 2.1, 3.7, 4.5, 5.1, 6.8, 7.3, 7.9, 9.2, 10.1))
lm_model <- lm(y ~ x, data = data)
calc_conf_int(estimate = 0.5, std_error = 0.1, model = lm_model, alpha = 0.05)

Description

Calculate scores, answers, and timing for 5.5 trial delay discounting from Qualtrics template

Usage

calc_dd(df)

Arguments

Value

A dataframe with k/ed50 values, answers, timing

Examples

calc_dd(five.fivetrial_dd)

Description

Calculate scores, answers, and timing for 5.5 trial probability discounting from Qualtrics template

Usage

calc_pd(df)

Arguments

Value

A dataframe with h/ep50 values, answers, timing

Examples

calc_pd(five.fivetrial_pd)

Description

Calculates the coefficient of determination ($R^2$) for a fitted model as one minus the ratio of the sum of squared errors (SSE) to the total sum of squares (SST).

Usage

calc_r2(model)

Arguments

Value

A numeric value representing the $R^2$ value of the model. Returns NA if the model is NULL.

Examples

# Example using a simple linear model
data <- data.frame(x = 1:10, y = c(1, 2, 3, 4, 5, 6, 7, 9, 10, 11))
lm_model <- lm(y ~ x, data = data)
calc_r2(lm_model)

Description

Checks a data set against the two Johnson and Bickel (2008) criteria for unsystematic delay-discounting data:

• Criterion 1: Any subsequent value of y exceeds the previous value by more than a specified proportion of the larger later reward (ll).

• Criterion 2: The last value of y is not at least a specified proportion less than the first value of y.

Usage

check_unsystematic(dat, ll = 1, c1 = 0.2, c2 = 0.1)

Arguments

Details

The criteria are applied separately to each id, so a data frame containing several subjects returns one row per subject. When an x (delay) column is present, each subject's points are ordered by delay before the criteria are applied.

Value

A tibble with the following columns:

• id: The unique identifier for the data set.

• c1_pass: Logical value indicating whether Criterion 1 was passed.

• c2_pass: Logical value indicating whether Criterion 2 was passed.

Examples

data <- tibble::tibble(
  id = c(rep("P1", 6)),
  x = c(1, 7, 30, 90, 180, 365), # delays
  y = c(0.9, 0.5, 0.3, 0.2, 0.1, 0.05) # indifference points
)
check_unsystematic(data, ll = 1, c1 = 0.2, c2 = 0.1)

Description

Returns the optimizer's flat named parameter vector: beta_k (one entry per fixed-effect design column, on the log-k scale), log_sigma_u, log_gamma, and beta0 (only when intercept = TRUE). This is the numeric escape hatch consumed by tooling such as car::deltaMethod.

Usage

## S3 method for class 'beezdiscounting_choice'
coef(object, ...)

Arguments

Value

Named numeric vector.

Description

Returns the optimizer's flat named parameter vector: beta_k (one entry per fixed-effect design column, on the log-k scale), log_sigma_u, and the auxiliary parameter (log_phi for family = "sltb", log_sigma_e for family = "gaussian"). This is the numeric escape hatch consumed by tooling such as car::deltaMethod.

Usage

## S3 method for class 'beezdiscounting_tmb'
coef(object, ...)

Arguments

Value

Named numeric vector.

Description

Equal-tailed quantile credible intervals on the report-space-transformed draws. Columns match confint.beezdiscounting_tmb(): term, estimate, conf.low, conf.high, level.

Usage

## S3 method for class 'beezdiscounting_brms'
confint(
  object,
  parm = NULL,
  level = 0.95,
  report_space = c("natural", "log10", "internal", "log"),
  ...
)

Arguments

Value

A tibble.

Description

Wald (Hessian-based) confidence intervals: ⁠estimate +/- z * se⁠ on the internal scale, then back-transformed to report_space.

Usage

## S3 method for class 'beezdiscounting_choice'
confint(
  object,
  parm = NULL,
  level = 0.95,
  report_space = c("internal", "natural"),
  ...
)

Arguments

Value

A tibble with columns term, estimate, conf.low, conf.high, level.

Description

Quantile credible intervals on the report-space-transformed draws; columns term, estimate, conf.low, conf.high, level. beta0 stays on the identity scale across report spaces, matching tidy.beezdiscounting_choice().

Usage

## S3 method for class 'beezdiscounting_choice_brms'
confint(
  object,
  parm = NULL,
  level = 0.95,
  report_space = c("natural", "log10", "internal", "log"),
  ...
)

Arguments

Value

A tibble.

Description

Wald (Hessian-based) confidence intervals: ⁠estimate +/- z * se⁠ on the internal (log-k) scale, then back-transformed to report_space.

Usage

## S3 method for class 'beezdiscounting_tmb'
confint(
  object,
  parm = NULL,
  level = 0.95,
  report_space = c("internal", "natural"),
  ...
)

Arguments

Value

A tibble with columns term, estimate, conf.low, conf.high, level.

Note

For a 2-RE fit the cor_re and log_sd_re rows are reported on their internal (atanh / log) scales and are NOT back-transformed by report_space = "natural" (only beta_k/log_s rows are); use VarCorr() to obtain the correlation and natural-log SDs.

Examples

dd <- simulate_dd_ip(n_subjects = 20, seed = 1)
fit <- fit_dd_tmb(dd, equation = "mazur", family = "sltb",
                  random_effects = k ~ 1, verbose = 0)
confint(fit)
confint(fit, parm = "k:(Intercept)", report_space = "natural")

Description

A dataset containing a set of fake delay discounting responses

Usage

dd_ip

Format

A data frame with 600 rows and 3 columns:

id

Subject identifier.

x

Delay until the larger, later reward (in days).

y

Indifference point, expressed as a proportion of the larger, later reward (these synthetic values are near, but not strictly within, the unit interval).

Description

loggamma ~ normal(1, 1) (choice sensitivity, typically 1-20) and the optional logit-scale bias b0 ~ normal(0, 1.5); logk defaults as in default_dd_priors() (autoscaled to the median delay when data is supplied; delay_var names the delay column).

Usage

default_dd_choice_priors(
  equation = c("mazur", "exponential"),
  intercept = FALSE,
  data = NULL,
  delay_var = "delay",
  factors = NULL,
  factor_interaction = FALSE,
  continuous_covariates = NULL,
  autoscale = !is.null(data)
)

Arguments

Value

A brmsprior data frame.

Description

Returns the default prior table used by fit_dd_brms(), so the defaults can be inspected, modified row-wise, and passed back via the fitter's prior argument. The logk location is the principled fix for k's delay-unit dependence: with autoscale = TRUE (the default whenever data is supplied) it centers k * median(delay) = 1 (the delay at which the Mazur curve crosses 0.5) via normal(-log(median(x)), 2.5); the static fallback is normal(-4.5, 2.5). The anchors used are attached as attr(, "autoscale_info"); numeric values are formatted with format(x, digits = 6, scientific = FALSE).

Usage

default_dd_priors(
  equation = c("mazur", "exponential", "green-myerson", "rachlin"),
  family = c("beta", "gaussian"),
  data = NULL,
  y_var = "y",
  x_var = "x",
  factors = NULL,
  factor_interaction = FALSE,
  continuous_covariates = NULL,
  autoscale = !is.null(data),
  random_effects = k ~ 1,
  covariance_structure = c("pdSymm", "pdDiag")
)

Arguments

Details

Other defaults: logs ~ normal(0, 0.5) (two-parameter equations; s is near 1 a priori), sd(logk) ~ student_t(3, 0, 1), phi ~ gamma(2, 0.1) (Beta precision; mean 20, far from brms's near-improper gamma(0.01, 0.01)), and sigma ~ student_t(3, 0, 0.25) for the Gaussian family (y is a proportion in the unit interval, so sd(y) <= 0.5).

Value

A brmsprior data frame, with attr(, "autoscale_info") when autoscaling was used.

Description

Bisection search over n_subjects for the smallest total N whose Monte Carlo power estimate from power_discounting() reaches target_power. The search accounts for Monte Carlo noise: at each evaluated N, replicates are added in batches (up to n_sim_max) until the Wilson interval for power lies wholly above or below the target; if it still straddles the target at the cap, the decision falls back to the point estimate and the result is flagged uncertain. The selected N and its lower neighbor are then re-evaluated with fresh replicates before minimality is claimed.

The returned n is an estimated minimum under Monte Carlo uncertainty rather than an exact bound. For grant-quality reporting, rerun power_discounting() at the returned n with a large n_sim (2000+) and report that estimate with its Monte Carlo confidence interval.

Monotonicity assumption. Bisection presumes that power is non-decreasing in n_subjects. Because every evaluated N is judged from its own independent replicates (and a convergence-conditioned denominator), a Monte Carlo fluctuation can make a lower N read "below" when its true power is above the target, so the search may step past a lower crossing that it never revisits. Evaluated N that contradict the assumption (a lower N reading "above" the selected N) demote the status to "uncertain"; N that were never evaluated cannot be checked. Widen n_sim/n_sim_max when the reported n matters. When the target is already met at n_range[1], that bound is likewise re-evaluated with fresh replicates before "at_lower_bound" is reported; if the second look does not clear the target the bound is treated as below and the bisection proceeds upward.

Usage

find_n_discounting(
  target_power = 0.8,
  effect = list(delta_k = NULL),
  design = list(),
  n_range = c(6, 200),
  n_sim = 200,
  n_sim_max = 4 * n_sim,
  alpha = 0.05,
  df = NULL,
  seed = NULL,
  equation = c("mazur", "exponential"),
  family = c("sltb", "gaussian"),
  random_effects = k ~ 1,
  multi_start = FALSE,
  verbose = TRUE,
  ...
)

Arguments

Value

An object of class beezdiscounting_power_n: a list with

n

Estimated smallest total n_subjects reaching target_power; NA when the confirmation pass contradicted the search (status = "unresolved").

target_power

As supplied.

status

"confirmed" (selected N re-confirmed above target and N - 1 below), "uncertain" (a decision relied on a point estimate, N - 1 also cleared the target on reconfirmation, or an evaluated lower N read above the target, so the returned N may not be minimal), "unresolved" (the selected N failed its own reconfirmation; n is NA), or "at_lower_bound" (the target was already met at n_range[1] on two independent looks; smaller N was not explored; widen n_range downward if that matters). These labels describe a heuristic Monte Carlo decision rule (repeated looks at ordinary Wilson intervals across several N) rather than a formal sequential error guarantee.

uncertain

Logical; TRUE when any search decision was made on a point estimate rather than a conclusive Wilson interval, or the status is not "confirmed"/"at_lower_bound".

evaluations

Tibble of every evaluation: n_subjects, n_sim_total, n_used, usable_fraction, power, ci_lower, ci_upper, decision. A warning fires when any evaluation had fewer than 95% usable fits.

alpha, df, effect, design, n_range, n_sim, n_sim_max, seed, settings, call

Echoed inputs and effective settings.

See Also

power_discounting() for the Monte Carlo engine.

Other power-analysis: power_discounting()

Examples

# Small search for demonstration (use larger n_sim for real planning)
res <- find_n_discounting(
  target_power = 0.8,
  effect = list(delta_k = log(8)),
  n_range = c(6, 20), n_sim = 4, n_sim_max = 4, seed = 1, verbose = FALSE
)
print(res)

Description

Fits a delay-discounting model to a data set with the specified equation and fitting method.

Usage

fit_dd(dat, equation, method)

Arguments

Value

A list object of class "fit_dd", containing:

• The fitted model(s).

• The original dataset (dat).

• The specified method (method).

Examples

data <- data.frame(
  id = rep(1:2, each = 6),
  x = rep(c(1, 7, 30, 90, 180, 365), 2),
  y = c(0.9, 0.5, 0.3, 0.2, 0.1, 0.05, 0.85, 0.55, 0.35, 0.15, 0.1, 0.05)
)
fit_dd(data, equation = "mazur", method = "two stage")

Description

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("squeeze", "zoib", "error"),
  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

Details

family = "beta" (default) uses Beta(link = "identity") with the mean squished into ⁠(1e-6, 1 - 1e-6)⁠, the closest brms analog of the TMB SLT-beta (family = "sltb" has no brms equivalent and errors with this pointer). Boundary observations (y exactly 0 or 1) are handled per boundary: "squeeze" (default) applies the Smithson-Verkuilen transform ⁠y* = (y (N - 1) + 0.5) / N⁠ to all responses (message reports the boundary count); "zoib" switches to zero_one_inflated_beta (statistically more honest, but it changes the estimand: k then describes interior responses only); "error" refuses to fit. 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).

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.

See Also

fit_dd_tmb(); default_dd_priors().

Description

Fits a trial-level binary choice model from two complementary perspectives. The structural model estimates the discount rate k directly from choices via the scale-invariant relative-value comparison ⁠logit P(LL) = beta0 + gamma * ((ll/ss) * D(k, delay) - 1)⁠, ⁠k = exp(X beta_k + sigma_u u)⁠. The descriptive model (Young 2018) forgoes a discount function and instead characterises choices via separate magnitude and delay sensitivities with optional correlated per-subject random slopes.

Usage

fit_dd_choice(
  data,
  mode = c("structural", "descriptive"),
  id_var = "id",
  ss_var = "ss_amount",
  ll_var = "ll_amount",
  delay_var = "delay",
  choice_var = "choice",
  equation = c("mazur", "exponential"),
  intercept = FALSE,
  predictors = NULL,
  random_slopes = TRUE,
  factors = NULL,
  factor_interaction = FALSE,
  continuous_covariates = NULL,
  start_values = NULL,
  tmb_control = list(iter_max = 1000, eval_max = 2000),
  multi_start = TRUE,
  verbose = 1,
  ...
)

Arguments

Details

The structural model (mode = "structural") parameterises choices through a classical discount function (Mazur hyperbolic or exponential) and shares the IP family's k/emmeans contract: k estimates and estimated marginal means are accessible via get_dd_param_emms() and get_dd_comparisons(). The descriptive model (mode = "descriptive") follows Young (2018): it regresses binary choice on log(ll_amount / ss_amount) (magnitude sensitivity) and log(delay + 1) (delay sensitivity) with no assumed discount function. When random_slopes = TRUE (default) each subject receives correlated random slopes on these two predictors; the primary inferential targets are the random-effect covariance (via nlme::VarCorr()) and per-subject slopes (via nlme::ranef()), not emmeans.

Value

An object of class beezdiscounting_choice.

See Also

simulate_dd_choice() for data generation; nlme::VarCorr() and nlme::ranef() for descriptive-mode random-effect output; get_dd_param_emms() and get_dd_comparisons() for structural-mode emmeans.

Examples

# Structural model: estimate a discount rate k from binary choices
sim_s <- simulate_dd_choice(n_subjects = 30, mode = "structural", seed = 1)
fit_s <- fit_dd_choice(sim_s, mode = "structural", equation = "mazur")
summary(fit_s)

# Descriptive (Young 2018) model: correlated per-subject magnitude/delay slopes
sim_d <- simulate_dd_choice(n_subjects = 30, mode = "descriptive", seed = 1)
fit_d <- fit_dd_choice(sim_d, mode = "descriptive")
summary(fit_d)
nlme::VarCorr(fit_d)

Description

The TMB structural likelihood of fit_dd_choice(): ⁠logit P(LL) = [b0] + gamma * ((ll/ss) * D(k, delay) - 1)⁠ with choice sensitivity gamma = exp(loggamma) and the discount rate k = exp(logk) carrying the subject random intercept. With bernoulli("logit") the brms nonlinear formula IS the logit, so the likelihood matches TMB exactly.

Usage

fit_dd_choice_brms(
  data,
  mode = c("structural", "descriptive"),
  id_var = "id",
  ss_var = "ss_amount",
  ll_var = "ll_amount",
  delay_var = "delay",
  choice_var = "choice",
  equation = c("mazur", "exponential"),
  intercept = FALSE,
  factors = NULL,
  factor_interaction = FALSE,
  continuous_covariates = NULL,
  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

Details

v1 implements mode = "structural" only: the descriptive (Young 2018) model is a plain logistic GLMM expressible directly with brms::brm().

Value

An object of class beezdiscounting_choice_brms. Coefficients are posterior medians on the estimation scale under the TMB names (one beta_k per design column, log_gamma, and beta0 when intercept = TRUE).

See Also

fit_dd_choice(); default_dd_choice_priors(); get_dd_param_emms() and get_dd_comparisons() for draws-based marginal means and contrasts of k.

Description

Linearizes ⁠D = 1/(1 + k t)⁠ to ln(1/D - 1) - ln(t) = ln(k) + e so that ln(k) has a closed-form per-subject estimator (the geometric mean of (1/D - 1)/t), a one-way random-effects model on the transformed scale has closed-form MLEs, and condition means can be compared with an exact F-test (anova.beezdiscounting_linear()); the test is exact when the random-effects variance estimate g-hat > 0, and when g-hat = 0 the statistic is still reported, with a warning. No numerical optimization is involved.

Usage

fit_dd_linear(
  data,
  y_var = "y",
  x_var = "x",
  id_var = "id",
  factors = NULL,
  response_scale = c("proportion", "percent", "amount"),
  ll = NULL,
  boundary = c("clamp", "drop", "error"),
  eps = NULL,
  conf_level = 0.95
)

Arguments

Details

Indifference points at exactly 0 or 1 are undefined under the transform; by default they are moved to eps and 1 - eps respectively (boundary = "clamp"). Only exact 0 and 1 are touched — interior points, however close to a bound, are used as observed. This is a package decision; the paper does not address boundary values. Use boundary = "drop" or "error" when you would rather not impute them.

Value

An object of class beezdiscounting_linear: a list with subjects (per-unit tibble of ln k estimates, intervals and log-likelihoods), re (closed-form random-effects MLEs mu, sigma2, g, or NULL if the design is unbalanced), design (factor name and levels), transform (boundary bookkeeping; the counts refer to the retained units), data (the long frame: id, x, y, the factor, and the derived columns condition, unit, y_lin, d_used, log_jac; factors may not name one of the derived columns, and other user columns are not carried), conf_level, and call.

References

Hinds, D., Tegge, A. N., Stein, J. S., LaConte, S. M., McClure, S. M., & Ferreira, M. A. R. (2026). To linearize or not to linearize: That is the Mazur delay discounting question. Journal of Mathematical Psychology, 130, 103006. doi:10.1016/j.jmp.2026.103006

See Also

beezdiscounting_linear-methods for the S3 surface, plot.beezdiscounting_linear(), anova.beezdiscounting_linear(), simulate_dd_linear(), and vignette("linearized-mazur").

Examples

fit <- fit_dd_linear(dd_ip)
fit

Description

Fits a 1-parameter discounting model (Mazur hyperbolic or exponential) with a random intercept on ⁠log k⁠, between-subject fixed effects, and either an SLT-beta or Gaussian observation family, using Template Model Builder for exact AD + Laplace approximation.

Usage

fit_dd_tmb(
  data,
  y_var = "y",
  x_var = "x",
  id_var = "id",
  equation = c("mazur", "exponential", "green-myerson", "rachlin"),
  family = c("sltb", "gaussian"),
  random_effects = k ~ 1,
  factors = NULL,
  factor_interaction = FALSE,
  continuous_covariates = NULL,
  ll = NULL,
  response_scale = c("proportion", "percent", "amount"),
  start_values = NULL,
  tmb_control = list(iter_max = 1000, eval_max = 2000),
  multi_start = TRUE,
  verbose = 1,
  covariance_structure = c("pdSymm", "pdDiag"),
  ...
)

Arguments

Value

An object of class beezdiscounting_tmb with components:

call

The matched call.

opt

Normalized optimizer result (par, objective, convergence, message).

model

List of coefficients, se, and variance_components.

sdr

TMB sdreport object (or NULL if SE computation failed).

hessian_pd

Logical positive-definiteness of the Hessian.

param_info

Model metadata (equation, family, dimensions, factor spec, parsed random effects).

formula_details

Fixed-effect design (X, rhs, contrasts).

subject_pars

Data frame of subject-level parameters. For a 1-RE fit (k ~ 1) the columns are ⁠id, u_i, k⁠; for a phi-target 2-RE fit (k + phi ~ 1) they are ⁠id, re_k, re_phi, k, phi⁠; for an s-target 2-RE fit (k + s ~ 1, GM/Rachlin) they are ⁠id, re_k, re_s, k, s⁠ where s is soft-clamped toward ⁠(0.05, 20)⁠.

loglik, AIC, BIC

Fit statistics.

converged, se_available

Convergence / SE-availability flags.

opt_warnings

Character vector of optimizer warnings.

data

The single filtered model frame (id/x/y + retained design columns), row-aligned with the design matrix.

data_all

The validated frame before complete-casing.

coercion_info

Scale-coercion/clamping audit list.

References

Young, M. E. (2017). Discounting: A practical guide to multilevel analysis of indifference data. Journal of the Experimental Analysis of Behavior, 108(1), 97-112. doi:10.1002/jeab.265

Kim, M., Koffarnus, M. N., & Franck, C. T. (2024). Thinking inside the bounds: Improved error distributions for indifference point data analysis and simulation via beta regression using common discounting functions. arXiv preprint arXiv:2404.18000.

Kim, M., Kaplan, B. A., Koffarnus, M. N., & Franck, C. T. (2025). Scale-location-truncated beta regression: Expanding beta regression to accommodate 0 and 1. arXiv preprint arXiv:2509.13167.

Examples

# Small two-subject long-format indifference-point data frame.
dd <- data.frame(
  id = rep(c("s1", "s2"), each = 5),
  x  = rep(c(7, 30, 180, 365, 730), times = 2),
  y  = c(0.95, 0.80, 0.45, 0.30, 0.15,
         0.90, 0.70, 0.40, 0.25, 0.10)
)
fit <- fit_dd_tmb(dd, equation = "mazur", family = "sltb", verbose = 0)
exp(fit$model$coefficients[["beta_k"]])  # population k

Description

Fitted values for a beezdiscounting_choice fit

Usage

## S3 method for class 'beezdiscounting_choice'
fitted(object, level = c("subject", "population"), ...)

Arguments

Value

Numeric vector of fitted choice probabilities, length nobs(object).

Description

Fitted values for a beezdiscounting_tmb fit

Usage

## S3 method for class 'beezdiscounting_tmb'
fitted(object, level = c("subject", "population"), ...)

Arguments

Value

Numeric vector of fitted indifference proportions, length nobs(object).

Examples

dd <- simulate_dd_ip(n_subjects = 20, seed = 1)
fit <- fit_dd_tmb(dd, verbose = 0)
head(fitted(fit))

Description

An example dataset containing four participants' data (two typical discounting patterns and two patterns suggesting potential misattention to the task).

Usage

five.fivetrial_dd

Format

Example Qualtrics output

Description

An example dataset containing four participants' data.

Usage

five.fivetrial_pd

Format

Example Qualtrics output

Description

Extract fixed effects from a structural choice model

Usage

## S3 method for class 'beezdiscounting_choice'
fixef(object, ...)

Arguments

Value

Named numeric vector (identical to coef()).

Description

Extract fixed effects from a TMB discounting model

Usage

## S3 method for class 'beezdiscounting_tmb'
fixef(object, ...)

Arguments

Value

Named numeric vector (identical to coef()).

Description

Generate fake MCQ data

Usage

generate_data_mcq(n_ids = 100, n_items = 27, seed = 1234, prop_na = 0)

Arguments

Value

Dataframe of subjectid, questionid, and response

Examples

generate_data_mcq(n_ids = 2, n_items = 27, prop_na = .01)

Description

Computes factor-level contrasts of the discount rate k from a fitted beezdiscounting_tmb model (or a structural beezdiscounting_choice fit). k is linear in the fixed-effect coefficients on the natural-log scale (log k = X beta_k), so each contrast is a linear combination of beta_k with a Wald standard error from the beta_k block of TMB::sdreport()'s fixed-effect covariance. Contrasts are reported on the log10 scale and, optionally, as multiplicative ratios (ratio = exp(est_log)). The returned container mirrors the beezdemand beezdemand_comparison shape, so tidy.beezdiscounting_comparison() gives a flat, cross-backend frame. Bayesian fits (beezdiscounting_brms, beezdiscounting_choice_brms) compute per-draw contrasts over the same reference grid: posterior medians with quantile credible intervals and post.prob in place of adjusted p-values (no multiplicity adjustment; the joint posterior already encodes contrast dependence).

Usage

get_dd_comparisons(
  fit,
  compare_specs = NULL,
  contrast_type = c("pairwise", "trt.vs.ctrl"),
  contrast_by = NULL,
  adjust = "holm",
  at = NULL,
  ci_level = 0.95,
  report_ratios = TRUE,
  ...
)

Arguments

Value

A beezdiscounting_comparison object: a list with a single element k, itself a list with emmeans (the EMM tibble from get_dd_param_emms()), contrasts_log10 (columns contrast, estimate, std.error, statistic, df, conf.low, conf.high, p.value), and (if report_ratios) contrasts_ratio (columns contrast, ratio, conf.low, conf.high, p.value). When contrast_by is active, the contrast tables gain leading by-column(s) before contrast. Attributes backend, adjustment_method, compare_specs_used, contrast_type_used, and contrast_by_used describe the call.

See Also

tidy.beezdiscounting_comparison() for the cross-backend frame.

Examples

sim <- simulate_dd_ip(
  n_subjects = 30, n_conditions = 2, delta_k = c(0, log(3)), seed = 1
)
fit <- fit_dd_tmb(sim, factors = "condition", verbose = 0)
get_dd_comparisons(fit, contrast_type = "pairwise")

Description

Computes estimated marginal means (EMMs) of the discount rate k from a fitted beezdiscounting_tmb model (or a structural beezdiscounting_choice fit, which shares the log k design). EMMs are computed on the log k scale (linear in the fixed-effect coefficients) using the averaging-matrix reference grid, then back-transformed with exp() so that k = exp(k_log). Standard errors use the beta_k block of TMB::sdreport()'s fixed-effect covariance; intervals are Wald on the log scale and exponentiated. Bayesian fits (beezdiscounting_brms, beezdiscounting_choice_brms) use the same reference grid with draws-based summaries: posterior medians and equal-tailed quantile intervals, no delta method.

Usage

get_dd_param_emms(fit, factors_in_emm = NULL, at = NULL, ci_level = 0.95, ...)

Arguments

Value

A tibble with columns level, k, k_log, std.error, conf.low, conf.high. k_log is the marginal mean on the log k scale; k = exp(k_log); std.error is the SE of k_log; the intervals are on the k (natural) scale.

Examples

sim <- simulate_dd_ip(
  n_subjects = 30, n_conditions = 2, delta_k = c(0, log(3)), seed = 1
)
fit <- fit_dd_tmb(sim, factors = "condition", verbose = 0)
get_dd_param_emms(fit)

Description

Get an internal instrument design table

Usage

get_lookup_table(items = NULL, instrument = NULL)

Arguments

Value

Data frame with the complete item design. MCQ tables: questionid, magnitude, kindiff, k_rank, ss_amount, ll_amount, delay (days). PDQ table: questionid, block, h_rank, overall_rank (the item's position in the pooled overall ladder used by the overall_h extension in score_pdq()), sc_amount, lu_amount, prob, theta (odds against), hindiff.

Examples

get_lookup_table()
get_lookup_table(items = 21)
get_lookup_table(instrument = "pdq")

Description

The dd glance columns with backend = "brms" and logLik/AIC/BIC = NA_real_ by design, plus the Bayesian comparison and diagnostic currency (elpd_loo, p_loo, looic, rhat_max, ess_bulk_min, num_divergent).

Usage

## S3 method for class 'beezdiscounting_brms'
glance(x, ...)

Arguments

Value

A one-row tibble.

Description

Returns a one-row tibble of model-level summary statistics suitable for use with broom::glance() workflows and multi-model comparison tables.

Usage

## S3 method for class 'beezdiscounting_choice'
glance(x, ...)

Arguments

Value

A one-row tibble with columns model_class, backend, mode, equation, nobs, n_subjects, n_random_effects, converged, logLik, AIC, BIC.

Description

Glance at a beezdiscounting_choice_brms model

Usage

## S3 method for class 'beezdiscounting_choice_brms'
glance(x, ...)

Arguments

Value

A one-row tibble.

Description

Returns a one-row tibble of model-level summary statistics suitable for use with broom::glance() workflows and multi-model comparison tables.

Usage

## S3 method for class 'beezdiscounting_tmb'
glance(x, ...)

Arguments

Value

A one-row tibble with columns model_class, backend, equation, family, nobs, n_subjects, n_random_effects, converged, logLik, AIC, BIC.

Description

Reshape MCQ data long to wide

Usage

long_to_wide_mcq(dat, q_col = "questionid", ans_col = "response")

Arguments

Value

Wide format data frame

Examples

long_to_wide_mcq(mcq27)

Description

Reshape MCQ data from long to wide (as used in the 21- and 27-Item Monetary Choice Questionnaire Automated Scorer)

Usage

long_to_wide_mcq_excel(dat, subj_col = "subjectid", ans_col = "response")

Arguments

Value

Wide format MCQ data that can be used in the Excel Automated Scorers

Examples

long_to_wide_mcq_excel(generate_data_mcq(2))

Description

Reshapes long-form Monetary Choice Questionnaire (MCQ) responses into the per-trial smaller-sooner versus larger-later choice frame consumed by fit_dd_choice(), joining each questionid to the canonical item design (immediate amount, delayed amount, delay) bundled in the lookup table (see get_lookup_table(); Kirby, Petry, & Bickel (1999) for the 27-item MCQ, or Kirby & Maraković (1996, Table 1) for the 21-item MCQ).

Usage

mcq_to_choice(
  responses,
  items = 27,
  id_var = "subjectid",
  question_var = "questionid",
  response_var = "response"
)

Arguments

Details

Unknown or non-coercible question ids raise an error rather than silently producing unmatched rows. Question ids 1-21 are valid in both the 21- and 27-item designs, so passing the wrong items does not error; it silently returns the wrong amounts/delays for those ids. Make sure items matches the instrument actually administered. Ragged input is allowed (subjects need not have all items) and NA responses are preserved (they are complete-cased by fit_dd_choice()). For the strict scorer see score_mcq().

Value

A tibble with columns id (character), ss_amount, ll_amount, delay (days), and choice (0/1, 1 = chose LL), in the input row order. Ready to pass to fit_dd_choice().

See Also

fit_dd_choice(), score_mcq(), get_lookup_table(), mcq27_to_choice()

Examples

ch <- mcq_to_choice(mcq27, items = 27)
head(ch)
# feeds directly into the structural choice model (requires TMB):
# fit_dd_choice(ch, mode = "structural", equation = "mazur")

Description

A dataset containing two illustrative participants' 21-item MCQ (Kirby & Marakovic, 1996) responses: subject 1 responds like the median subject in Kirby & Marakovic (1996) (switching between the eighth and ninth ladder items, overall k of about 0.007) and subject 2 always chooses the larger delayed reward. These responses are constructed for illustration only, not drawn from an empirical sample.

Usage

mcq21

Format

Long-form data.frame with columns: subjectid, questionid, response.

Source

Constructed in data-raw/mcq21_data.R.

Description

A dataset containing two participants' data (same data as in the paper by Kaplan et al., 2016)

Usage

mcq27

Format

Long-form data.frame with columns: subjectid, questionid, response.

Description

Reshapes long-form Monetary Choice Questionnaire (MCQ) responses into the per-trial smaller-sooner versus larger-later choice frame consumed by fit_dd_choice(), joining each questionid to the canonical item design (immediate amount, delayed amount, delay) bundled in the lookup table (see get_lookup_table(); Kirby, Petry, & Bickel (1999) for the 27-item MCQ, or Kirby & Maraković (1996, Table 1) for the 21-item MCQ).

Usage

mcq27_to_choice(
  responses,
  id_var = "subjectid",
  question_var = "questionid",
  response_var = "response"
)

Arguments

Details

Unknown or non-coercible question ids raise an error rather than silently producing unmatched rows. Question ids 1-21 are valid in both the 21- and 27-item designs, so passing the wrong items does not error; it silently returns the wrong amounts/delays for those ids. Make sure items matches the instrument actually administered. Ragged input is allowed (subjects need not have all items) and NA responses are preserved (they are complete-cased by fit_dd_choice()). For the strict scorer see score_mcq().

Value

A tibble with columns id (character), ss_amount, ll_amount, delay (days), and choice (0/1, 1 = chose LL), in the input row order. Ready to pass to fit_dd_choice().

See Also

fit_dd_choice(), score_mcq27(), get_lookup_table(), mcq_to_choice()

Examples

ch <- mcq27_to_choice(mcq27)
head(ch)
# feeds directly into the structural choice model (requires TMB):
# fit_dd_choice(ch, mode = "structural", equation = "mazur")

Description

A dataset containing two illustrative participants' Probability Discounting Questionnaire (Madden, Petry, & Johnson, 2009) responses: subject 1 chooses the guaranteed reward through the fifth h rank of every block and the risky reward after (block h values of about 1.22), and subject 2 always chooses the risky reward. These responses are constructed for illustration only, not drawn from an empirical sample.

Usage

pdq

Format

Long-form data.frame with columns: subjectid, questionid, response (0 = guaranteed, 1 = risky).

Source

Constructed in data-raw/pdq_data.R.

Description

Reshapes long-form Probability Discounting Questionnaire (PDQ) responses into a per-trial certain-versus-risky choice frame, joining each questionid to the canonical item design (certain amount, probabilistic amount, probability, and odds against theta = (1 - p)/p) from Madden, Petry, & Johnson (2009; see get_lookup_table()).

Usage

pdq_to_choice(
  responses,
  id_var = "subjectid",
  question_var = "questionid",
  response_var = "response"
)

Arguments

Details

Unknown or non-coercible question ids raise an error rather than silently producing unmatched rows. Ragged input is allowed (subjects need not have all items) and NA responses are preserved. For the strict scorer see score_pdq().

Value

A tibble with columns id (character), sc_amount, lu_amount, prob, theta (odds against winning), and choice (0/1, 1 = chose the risky reward), in the input row order.

See Also

score_pdq(), get_lookup_table(), mcq_to_choice()

Examples

ch <- pdq_to_choice(pdq)
head(ch)

Description

Plots the delay-discounting data with the fitted model overlaid.

Usage

plot_dd(
  fit_dd_object,
  xlabel = "Delay",
  ylabel = "Indifference Point",
  title = "",
  logx = TRUE
)

Arguments

Value

A ggplot object representing the fitted model and data.

Examples

data <- data.frame(
  id = rep(1:2, each = 6),
  x = rep(c(1, 7, 30, 90, 180, 365), 2),
  y = c(0.9, 0.5, 0.3, 0.2, 0.1, 0.05, 0.85, 0.55, 0.35, 0.15, 0.1, 0.05)
)
fit <- fit_dd(data, equation = "mazur", method = "mean")
plot_dd(fit)

Description

QQ plot of the estimated subject random-effect deviates against a normal reference (the standard check on the Gaussian random-effects assumption). Methods exist for the TMB indifference-point (fit_dd_tmb()) and choice (fit_dd_choice()) models. Bayesian fits are intentionally excluded; use brms::pp_check() and MCMC diagnostics there instead. The plotted deviates are shrunken (empirical-Bayes) estimates, so the normal reference is approximate.

Usage

## S3 method for class 'beezdiscounting_tmb'
plot_qq(object, which = NULL, ...)

## S3 method for class 'beezdiscounting_choice'
plot_qq(object, which = NULL, ...)

Arguments

Value

A ggplot2::ggplot object.

Description

Visualize a fitted fit_dd_brms() model. type = "population" (default) draws the posterior-median discount curve with a credible band for the population conditional mean (random effects set to zero; a credible band rather than a posterior-predictive band) over the observed (boundary-squeezed) indifference points; type = "individual" adds per-subject median curves; type = "parameters" shows the posterior subject-k caterpillar; type = "resid" plots standardized residuals against fitted values.

Usage

## S3 method for class 'beezdiscounting_brms'
plot(
  x,
  type = c("population", "individual", "parameters", "resid"),
  ids = NULL,
  at = NULL,
  n_points = 200,
  x_trans = c("log10", "linear"),
  show_observed = TRUE,
  ...
)

Arguments

Value

A ggplot2::ggplot object.

Description

Visualize a fitted fit_dd_choice() model. For a structural fit, type = "population" (default) draws the implied discount curve from the estimated discount rate, type = "individual" adds the per-subject curves, and type = "parameters" shows the subject-k distribution. For both structural and descriptive (Young 2018) fits, type = "calibration" plots the observed choice proportion against the fitted P(choose LL), binned by fitted probability (a likelihood-geometry-agnostic goodness-of-fit view). Descriptive fits have no structural discount rate, so the curve/parameter types are unavailable.

Usage

## S3 method for class 'beezdiscounting_choice'
plot(
  x,
  type = NULL,
  ids = NULL,
  at = NULL,
  n_points = 200,
  x_trans = c("log10", "linear"),
  show_observed = TRUE,
  ...
)

Arguments

Value

A ggplot2::ggplot object.

Examples

ch <- simulate_dd_choice(n_subjects = 12, seed = 1)
fit <- fit_dd_choice(ch, equation = "mazur")
plot(fit)
plot(fit, type = "calibration")

Description

Visualize a fitted fit_dd_choice_brms() model. type = "population" (default) draws the posterior implied discount curve with a credible band; type = "individual" adds per-subject (posterior-median) curves; type = "parameters" shows the posterior subject-k caterpillar; and type = "calibration" plots observed choice proportion against fitted P(choose LL), binned by fitted probability.

Usage

## S3 method for class 'beezdiscounting_choice_brms'
plot(
  x,
  type = c("population", "individual", "calibration", "parameters"),
  ids = NULL,
  at = NULL,
  n_points = 200,
  x_trans = c("log10", "linear"),
  show_observed = TRUE,
  ...
)

Arguments

Value

A ggplot2::ggplot object.

Description

Forest plot of the pairwise (or treatment-vs-control) contrasts from get_dd_comparisons(). type = "ratio" (default) shows each k ratio on a log axis with a reference line at 1; type = "difference" shows the difference in log10 k with a reference line at 0. A contrast is flagged when its interval excludes the null (a backend-agnostic rule that applies to both the frequentist (TMB) and Bayesian (brms) backends).

Usage

## S3 method for class 'beezdiscounting_comparison'
plot(x, type = c("ratio", "difference"), ...)

Arguments

Value

A ggplot2::ggplot object.

Examples

sim <- simulate_dd_ip(n_subjects = 30, n_conditions = 2, delta_k = c(0, 0.5), seed = 1)
fit <- fit_dd_tmb(sim, equation = "mazur", factors = "condition")
cmp <- get_dd_comparisons(fit)
plot(cmp)

Description

Visualize a fitted fit_dd_linear() model. type = "population" (default) draws the hyperbola implied by each condition's geometric-mean discount rate exp(mu) over the observed indifference points (one curve per condition when the fit has a factor); type = "individual" adds each subject's own closed-form hyperbola ⁠1 / (1 + k_i x)⁠ as thin lines; type = "transformed" shows the linearization itself, ln(1/D - 1) against ln(delay) per subject with a slope-1 line at that subject's ⁠ln k⁠; type = "parameters" shows the subject ⁠ln k⁠ estimates with their t-intervals (by condition, with the condition mean mu and its interval from confint() overlaid, when the fit has a factor; otherwise ordered as a caterpillar), or the same quantities as k on a log10 axis with k_scale = "log10"; type = "resid" plots the standardized residual on the transformed scale against ln(delay).

Usage

## S3 method for class 'beezdiscounting_linear'
plot(
  x,
  type = c("population", "individual", "transformed", "parameters", "resid"),
  ids = NULL,
  n_points = 200,
  x_trans = c("log10", "linear"),
  show_observed = TRUE,
  k_scale = c("ln", "log10"),
  ...
)

Arguments

Details

The method shares its core arguments (type, ids, n_points, x_trans, show_observed) with plot.beezdiscounting_tmb(); it adds k_scale and omits at (a linearized fit has no covariates or reference grid, so there is nothing to condition on). x_trans, n_points and show_observed apply to the "population" and "individual" curves only; ids applies to "individual" (default: every subject) and "transformed" (default: the first 12 subjects, with a message when the fit has more).

Unlike the mixed-model tiers there is no shrinkage: every subject's curve is its own closed-form estimate, and the per-subject pictures do not need the random-effects component. "population" requires it (the condition means mu); when the design was unbalanced and re is NULL, "population" errors, "individual" draws the subject curves only, and "parameters" omits the condition means with a message.

Two pictures deliberately differ from the mixed-model tiers because this model is fitted on the linearized scale. "parameters" defaults to ⁠ln k⁠ (the scale of coef(), confint() and anova()'s effect size); k_scale = "log10" gives the sibling tiers' picture, k on a log10 axis. "resid" plots the transformed-scale residual against ln(delay), the model's regressor with its slope fixed at 1, so a trend across delay may indicate delay-dependent lack of fit; the sibling tiers plot residuals against fitted values instead. k_scale affects only "parameters".

The "transformed" abscissa is ln(delay) by construction (a log10 axis would break the slope-1 reference), so x_trans is ignored there. Indifference points at exactly 0 or 1 are shown at their observed value on the raw scale, and at their boundary-handled value d_used on the transformed scale; points dropped under boundary = "drop" do not appear on the transformed scale. "resid" uses augment()'s .std_resid, the transformed-scale residual divided by the model's error standard deviation.

Value

A ggplot2::ggplot object.

See Also

fit_dd_linear(), beezdiscounting_linear-methods, vignette("linearized-mazur").

Examples

sim <- simulate_dd_linear(
  n_subjects = 10, delays = c(7, 30, 90, 365),
  mu = c(A = -6, B = -4.5), sigma2 = 2, g = 8, seed = 1
)
fit <- fit_dd_linear(sim, factors = "condition")
plot(fit)
plot(fit, type = "individual")
plot(fit, type = "transformed", ids = c("A_1", "B_1"))
plot(fit, type = "parameters")

Description

Visualize a fitted fit_dd_tmb() model. type = "population" draws the population (random-effects-at-zero) discount curve over the observed indifference points; type = "individual" adds the per-subject curves (the shrinkage picture). For a fit with factors the population curve is drawn once per factor-level combination; use at to condition on specific levels or covariate values. type = "parameters" shows the distribution of the subject-specific discount rate k (log scale); type = "resid" plots standardized (Pearson) residuals against fitted values.

Usage

## S3 method for class 'beezdiscounting_tmb'
plot(
  x,
  type = c("population", "individual", "parameters", "resid"),
  ids = NULL,
  at = NULL,
  n_points = 200,
  x_trans = c("log10", "linear"),
  show_observed = TRUE,
  ...
)

Arguments

Value

A ggplot2::ggplot object.

Examples

sim <- simulate_dd_ip(n_subjects = 12, seed = 1)
fit <- fit_dd_tmb(sim, equation = "mazur")
plot(fit)
plot(fit, type = "individual")

Description

Plot Proportion of Guaranteed Choices by h Rank

Usage

## S3 method for class 'prop_sc_output'
plot(
  x,
  ...,
  pt_shape = 21,
  pt_fill = "white",
  pt_size = 3,
  title = "Proportion of guaranteed choices by h rank",
  xlab = "h value rank",
  ylab = "Proportion of guaranteed choices"
)

Arguments

Value

A ggplot object.

Examples

plot(prop_sc(pdq))

Description

Plots the proportion of SIR/SS choices by k value using the output of the prop_ss function.

Usage

## S3 method for class 'prop_ss_output'
plot(
  x,
  ...,
  pt_shape = 21,
  pt_fill = "white",
  pt_size = 3,
  title = "Proportion of SIR/SS choices by k value",
  xlab = "k value rank",
  ylab = "Proportion of SS choices"
)

Arguments

Value

A ggplot object.

Examples

plot(prop_ss(mcq27))

Description

Boxplots the MCQ scores for different metrics (small_k, medium_k, large_k, geomean_k, overall_k). Log-transformed k values (trans = "log"/"ln") are detected and the y-axis label adjusted accordingly.

Usage

## S3 method for class 'score_mcq_output'
plot(x, ..., xlab = "Metric", alpha = 0.3)

Arguments

Value

A ggplot object showing the boxplot of MCQ scores.

Examples

dat21 <- data.frame(subjectid = 1:5, questionid = rep(1:21, 5), response = 1)
plot(score_mcq(dat21, items = 21))

Description

Boxplots the MCQ-27 scores for different metrics (small_k, medium_k, large_k, geomean_k, overall_k). Log-transformed k values (trans = "log"/"ln") are detected and the y-axis label adjusted accordingly.

Usage

## S3 method for class 'score_mcq27_output'
plot(x, ..., xlab = "Metric", alpha = 0.3)

Arguments

Value

A ggplot object showing the boxplot of MCQ-27 scores.

Examples

plot(score_mcq27(mcq27))

Description

Boxplot of PDQ h metrics (block1_h, block2_h, block3_h, mean_h, geomean_h, and the pooled overall_h extension), handling the same log transformations as plot.score_mcq_output().

Usage

## S3 method for class 'score_pdq_output'
plot(x, ..., xlab = "Metric", alpha = 0.3)

Arguments

Value

A ggplot object showing the boxplot of PDQ scores.

Examples

plot(score_pdq(pdq))

Description

Estimates statistical power to detect a between-subject difference in the discount rate (a shift delta_k on natural-log k) by simulation: each replicate (1) simulates a two-condition between-subject indifference-point dataset with simulate_dd_ip() under assumed population parameters plus the effect, (2) refits it with fit_dd_tmb(), and (3) tests the log-k condition contrast with a Wald test at level alpha, referred to a t distribution with df degrees of freedom (see the df argument). Power is the proportion of usable fits (converged, positive-definite Hessian, finite standard error) that reject.

Because the power estimate is a proportion from finitely many replicates, it is reported with a Wilson score confidence interval (power_mc_ci). Both a p-value verdict (p < alpha) and a confidence-interval verdict (Wald CI excludes 0) are recorded per replicate; they use the same standard error and reference distribution, so they coincide by construction, and both rates are returned.

This mirrors beezdemand::power_demand(). The design asymmetry is intentional: this package's simulator models a between-subject condition (subjects are split across conditions round-robin), while beezdemand's models a within-subject condition.

Usage

power_discounting(
  n_subjects,
  effect = list(delta_k = NULL),
  design = list(),
  n_sim = 500,
  alpha = 0.05,
  df = NULL,
  seed = NULL,
  equation = c("mazur", "exponential"),
  family = c("sltb", "gaussian"),
  random_effects = k ~ 1,
  multi_start = FALSE,
  verbose = TRUE,
  ...
)

Arguments

Details

A replicate whose fit fails (non-convergence, non-positive-definite Hessian, unusable standard error, or an error) is excluded from the power denominator and reported through the ⁠n_*⁠ counts and ⁠$replicates$status⁠. It is never counted as "no effect detected", which would bias power in an unpredictable direction. A warning is issued when fewer than 95% of replicates are usable.

The v1 scope is a single fixed-effect delta_k under the package's existing between-subject simulator. Effects on s or phi, power for derived measures (ED50, AUC), and arbitrary designs are out of scope; see vignette("power-analysis").

Value

An object of class beezdiscounting_power: a list with

power

Estimated power: proportion of usable replicates whose Wald CI excludes 0 (equal to hit_rate_ci). NA if no replicate was usable.

power_mc_ci

Wilson 95% confidence interval on power, reflecting Monte Carlo uncertainty from n_used replicates.

hit_rate_p

Proportion of usable replicates with p < alpha.

hit_rate_ci

Proportion of usable replicates whose Wald CI excludes 0 (the same decision rule as hit_rate_p, since both use the same SE and t reference; both reported).

n_sim

Total replicates attempted.

n_converged

Replicates whose fit converged.

n_hessian_pd

Replicates with a positive-definite Hessian.

n_used

Replicates entering the power denominator (converged, positive-definite Hessian, finite SE).

alpha

Nominal test level.

df

Degrees of freedom of the t reference distribution actually used (n_subjects - 2 unless overridden).

effect

The validated effect specification (name and delta).

target_term

The tested coefficient ("k:conditionC2").

design

The merged design list actually used.

n_subjects

As supplied.

replicates

Tibble with one row per replicate: sim, status ("ok", "nonconverged", "hessian_not_pd", "se_unusable", "error"), converged, hessian_pd, estimate, se, statistic, p_value, ci_lower, ci_upper, hit_p, hit_ci, and message (error text, if any). Estimates are on the natural-log-k scale of the simulated delta_k.

seed

As supplied.

settings

List of equation, family, multi_start, and the deparsed random-effects formula.

call

The matched call.

See Also

find_n_discounting() to search for the smallest adequate sample size; fit_dd_tmb() for the model being refit; simulate_dd_ip() for the data-generating process.

Other power-analysis: find_n_discounting()

Examples

# Quick exploratory run (use n_sim >= 500 for real planning)
res <- power_discounting(
  n_subjects = 8,
  effect = list(delta_k = log(2)),
  n_sim = 4, seed = 1, verbose = FALSE
)
print(res)

Description

type = "response" summarizes brms::posterior_epred() draws (the indifference-proportion mean) to the posterior median with equal-tailed intervals; level = "subject" conditions on the subject random effects, "population" sets them to zero. type = "parameters" returns the per-subject k summaries.

Usage

## S3 method for class 'beezdiscounting_brms'
predict(
  object,
  newdata = NULL,
  type = c("response", "parameters"),
  level = "subject",
  probs = c(0.025, 0.975),
  draws = FALSE,
  ...
)

Arguments

Value

A tibble (or draws matrix).

Description

Predict from a structural choice discounting model

Usage

## S3 method for class 'beezdiscounting_choice'
predict(
  object,
  newdata = NULL,
  type = c("prob", "parameters"),
  level = "subject",
  ...
)

Arguments

Value

• type = "parameters": the subject-parameter tibble (one row per subject).

• type = "prob": newdata as a tibble plus a .prob column of fitted choice probabilities.

Description

Predict P(LL) from a beezdiscounting_choice_brms model

Usage

## S3 method for class 'beezdiscounting_choice_brms'
predict(
  object,
  newdata = NULL,
  type = c("response", "parameters"),
  level = "subject",
  probs = c(0.025, 0.975),
  draws = FALSE,
  ...
)

Arguments

Value

A tibble (or draws matrix).

Description

Predict from a TMB mixed-effects discounting model

Usage

## S3 method for class 'beezdiscounting_tmb'
predict(
  object,
  newdata = NULL,
  type = c("response", "parameters"),
  level = "subject",
  ...
)

Arguments

Value

• type = "parameters": the subject-parameter tibble (one row per subject).

• type = "response", level = "subject" only: newdata as a tibble plus a .fitted column.

• type = "response", "population" requested: newdata plus predict.fixed (and predict.id when "subject" is also requested), matching the nlme::predict.lme(level = 0:1) column-name convention so nlme-based plotting code runs unchanged.

Examples

dd <- simulate_dd_ip(n_subjects = 20, seed = 1)
fit <- fit_dd_tmb(dd, equation = "mazur", family = "sltb",
                  random_effects = k ~ 1, verbose = 0)

# Subject-conditional fitted values
head(predict(fit, type = "response"))

# Population-mean curve at specific delays (no id column needed)
nd <- data.frame(x = c(7, 30, 180, 365))
predict(fit, newdata = nd, level = "population")

# Per-subject parameters
head(predict(fit, type = "parameters"))

Description

Print a structural choice discounting fit

Usage

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

Arguments

Value

Invisibly returns x.

Description

Print a TMB mixed-effects discounting fit

Usage

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

Arguments

Value

Invisibly returns x.

Examples

dd <- simulate_dd_ip(n_subjects = 20, seed = 1)
fit <- fit_dd_tmb(dd, equation = "mazur", family = "sltb",
                  random_effects = k ~ 1, verbose = 0)
print(fit)

Description

Print a structural choice discounting model summary

Usage

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

Arguments

Value

Invisibly returns x.

Description

Print a TMB discounting model summary

Usage

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

Arguments

Value

Invisibly returns x.

Description

The PDQ analog of prop_ss(): pools all subjects' responses and reports, for each of the 10 h ranks (one item per block per rank), the proportion choosing the smaller guaranteed reward. 1 - prop_sc is the risky-choice proportion at that rank.

Usage

prop_sc(dat)

Arguments

Value

Dataframe with proportion of guaranteed choices at each h rank

Examples

prop_sc(pdq)

Description

Calculate proportion of SIR/SS responses at each k value

Usage

prop_ss(dat, items = 27)

Arguments

Details

items must match the instrument actually administered. Question ids 1-21 are valid in both the 21- and 27-item designs, so passing the wrong items does not error; it silently pools responses into the wrong k-rank rows. If the observed question ids do not exactly match the requested design, prop_ss() warns.

Value

Dataframe with proportion of SIR/SS responses at each k rank

Examples

prop_ss(mcq27)
dat21 <- data.frame(subjectid = 1, questionid = 1:21, response = 1)
prop_ss(dat21, items = 21)

Description

Structural: id, the standardized random-intercept deviate u_i (such that ⁠log k_i = X beta + sigma_u * u_i⁠), and the resolved per-subject discount rate k. Descriptive (Young 2018): id plus the per-subject magnitude and delay random slopes (b_mag, b_delay); pooled (random_slopes = FALSE) returns an id-only frame.

Usage

## S3 method for class 'beezdiscounting_choice'
ranef(object, ...)

Arguments

Value

Data frame: structural id/u_i/k; descriptive id/b_mag/ b_delay (or id-only when pooled).

Description

Extract subject-level random effects from a TMB discounting model

Usage

## S3 method for class 'beezdiscounting_tmb'
ranef(object, ...)

Arguments

Value

Data frame keyed by id. For a 1-RE fit (k ~ 1): the standardized random-intercept deviate u_i (such that ⁠log k_i = X beta + sigma_u * u_i⁠) and the resolved per-subject discount rate k (no phi or s column; both are population-level). For a phi-target 2-RE fit (k + phi ~ 1): ⁠(re_k, re_phi)⁠ offsets plus per-subject k and phi. For an s-target 2-RE fit (k + s ~ 1, GM/Rachlin only): ⁠(re_k, re_s)⁠ offsets plus per-subject k and s (soft-clamped toward ⁠(0.05, 20)⁠).

Description

Residuals for a beezdiscounting_choice fit

Usage

## S3 method for class 'beezdiscounting_choice'
residuals(
  object,
  type = c("response", "pearson"),
  level = c("subject", "population"),
  ...
)

Arguments

Value

Numeric vector of residuals, length nobs(object).

Description

Residuals for a beezdiscounting_tmb fit

Usage

## S3 method for class 'beezdiscounting_tmb'
residuals(
  object,
  type = c("response", "pearson"),
  level = c("subject", "population"),
  ...
)

Arguments

Value

Numeric vector of residuals, length nobs(object).

Examples

dd <- simulate_dd_ip(n_subjects = 20, seed = 1)
fit <- fit_dd_tmb(dd, verbose = 0)
head(residuals(fit))
head(residuals(fit, type = "pearson"))

Description

Extracts parameter estimates, fit statistics, and confidence intervals from a fitted delay-discounting model.

Usage

results_dd(fit_dd_object)

Arguments

Value

A tibble containing the following columns:

• id: The participant or group ID (if applicable).

• term: The model parameter (e.g., k).

• estimate: The estimated value of the parameter.

• std.error: The standard error of the parameter estimate.

• statistic: The t-statistic for the parameter estimate.

• p.value: The p-value for the parameter estimate.

• conf_low: The lower bound of the 95% confidence interval.

• conf_high: The upper bound of the 95% confidence interval.

• R2: The coefficient of determination ($R^2$).

Examples

data <- data.frame(
  id = rep(1:2, each = 6),
  x = rep(c(1, 7, 30, 90, 180, 365), 2),
  y = c(0.9, 0.5, 0.3, 0.2, 0.1, 0.05, 0.85, 0.55, 0.35, 0.15, 0.1, 0.05)
)
fit <- fit_dd(data, equation = "mazur", method = "two stage")
results_dd(fit)

Description

Score 5.5 trial delay discounting from Qualtrics template

Usage

score_dd(df)

Arguments

Details

Currently assumes the attending questions are present and labeled "Attend-LL" and "Attend-SS"

Value

A dataframe with id, indexes, response, k value, and effective delay 50.

Examples

score_dd(five.fivetrial_dd)

Description

Score an MCQ (27- or 21-item)

Usage

score_mcq(
  dat = dat,
  items = 27,
  impute_method = "none",
  round = 6,
  random = FALSE,
  trans = "none",
  return_data = FALSE,
  verbose = FALSE
)

Arguments

Details

Each subject's data must satisfy a strict contract: exactly one row per canonical question id (items of them; no duplicates, no unknown ids, none missing) and responses coded 0, 1, or NA (numeric, logical, or character/factor values that coerce to 0/1). Malformed input errors rather than silently mis-scoring. Contrast with mcq_to_choice()'s lenient, ragged contract, which accepts partial per-subject coverage.

Value

If return_data = FALSE (default), a summary data frame with one row per subject. If return_data = TRUE, a list with results (that summary data frame) and data (the input data plus a newresponse column reflecting any imputation).

Examples

score_mcq(mcq27, items = 27)
dat21 <- data.frame(subjectid = 1, questionid = 1:21, response = 1)
score_mcq(dat21, items = 21)

Description

Score 27-item MCQ

Usage

score_mcq27(
  dat = dat,
  impute_method = "none",
  round = 6,
  random = FALSE,
  trans = "none",
  return_data = FALSE,
  verbose = FALSE
)