Feature release covering everything since 0.2.0. Headline changes:
fit_demand_tmb()) with automatic
differentiation, a Laplace approximation, multi-start optimization,
optional estimation of k, factor-expanded / multi-block (pdBlocked)
and continuous within-subject random-slope structures, and the full
post-hoc surface (EMMs, contrasts, subject-level parameters, parametric
bootstrap CIs, diagnostics).power_demand(), find_n_demand()) for
within- and between-subject designs.fit_demand_fixed() and the zben
Pmax/Omax are computed numerically; see "Bug fixes that can change
estimates" for exactly which outputs can differ from 0.2.0 and how to pin
the old numbers.check_demand_model() no longer
reports failed checks as passing.simulate_hurdle_data(part2 = "snd")), a
breaking predict.beezdemand_hurdle() default (type = "demand"), and
the harmonized get_demand_comparisons() / EMM API across backends.The subsections below give the per-change detail, oldest at the bottom.
power_demand() estimates statistical power for detecting a condition
difference in Q0 or alpha, by simulating from the package's
mixed-effects demand model and refitting each replicate with
fit_demand_tmb(). Supports both a two-condition within-subject design
(default) and, via design_type = "between", a two-arm between-subject
design testing the group difference; the between design composes the two
arms from the same simulator (one condition per arm), uses df = n - 2,
and refits per-subject intercept random effects that match its
data-generating process exactly. Reports the power estimate with a Wilson
Monte Carlo confidence interval, p-value and CI-exclusion hit rates, and
convergence diagnostics; non-usable fits are excluded from the denominator
and surfaced, never counted as misses.find_n_demand() searches for the smallest n_subjects reaching a target
power via bisection, adding replicates adaptively where the Monte Carlo
verdict is ambiguous and re-confirming the selected N before reporting.
Also takes design_type.pwr::pwr.t.test() (one- and two-sample), monotonicity,
and seed reproducibility are verified in the test suite; see
vignette("power-analysis") for scope and validity notes.The fixes in this subsection correct wrong numbers rather than add features, so
outputs can differ under the stated conditions. Single-subject fits and paths
that were already correct are unchanged. To reproduce the old numbers exactly,
pin the previous release: remotes::install_version("beezdemand", "0.2.0").
The comparison is against 0.2.0 except where an item says otherwise: the TMB
tier, the hurdle tier and the power functions are new in 0.3.0, so for those the
change is against the development version rather than against a release.
calc_group_metrics() on an NLME fit now averages over the factor cells
that were observed, with equal weight per cell. It previously let
emmeans average over the full factorial grid, so a design with an empty
cell (say, no Female subjects at site B) extrapolated the additive model
into that cell before taking the geometric mean of Q0 and alpha. The TMB
method already used observed cells only; the two backends now agree on any
design fit in log space, and the policy is written out under
"Marginalisation policy" in both methods' help. Unconditioned numbers
change only when a factor cell has no subjects. In the same method, a
factor level passed in
at on a collapse_levels fit was silently ignored (the restriction was
keyed by the original factor name, which the collapsed model does not
have) and the result was the unconditioned one; at is now applied to
each parameter's collapsed cells, and a requested cell with no data is an
error rather than an extrapolation.
Hurdle marginal P(zero) now integrates over the fitted normal random
effect by default. predict() and plot() for beezdemand_hurdle fits
used marginal_method = "kde", a kernel density of the shrunken subject
intercepts, which understates the random-effect spread and pulls the
population curve toward the conditional one; the default is now
"normal", the distribution the likelihood itself integrates over, and
the integral runs over the whole real line instead of a truncated
interval whose weights were not renormalised. "kde" and "empirical"
remain available as descriptive summaries. The hurdle tier is new in
0.3.0, so this differs only from the development version.
fit_demand_tmb() now fixes k at 2 by default (estimate_k = FALSE).
The TMB tier is new in 0.3.0 and has never been released, so this changes
results only for users who installed the development version; nothing on CRAN
is affected. A free k is identified only by the curvature that appears as
consumption approaches its floor: the response depends on k through
k * (exp(-alpha * Q0 * price) - 1), and while alpha * Q0 * price stays
small that term is linear in price with slope k * alpha, so only the product
is identified. On the package's own apt data with
equation = "exponentiated" the old default walked that ridge to k = 1.8e13
with alpha = 2e-16, ending non-converged with a non-positive-definite
Hessian; at the new default the same call converges. k = 2 is the convention
of Hursh & Silberberg (2008) and the default of fit_demand_fixed(). Being a
convention rather than an estimate, it deserves a sensitivity fit at a second
k. Pass estimate_k = TRUE for the old behaviour, or k for another
constant. Two consequences for development-version users: update() replays a
stored call, so updating a fit whose original call omitted estimate_k now
fixes k (a call that named estimate_k = TRUE is unaffected), and a
tmb_control$warm_start vector saved from a free-k fit no longer matches the
parameter count of a default fit.
check_demand_model() screens free-k TMB fits for that failure.
The boundary slot of the returned diagnostics was previously a hard-coded
empty list for beezdemand_tmb fits. It now carries k_identification,
which flags an implausible k, a decay exponent alpha * Q0 * price that
never leaves its linear regime over the observed prices, a degenerate alpha,
or log_k resting on a user-supplied optimizer bound; a non-positive-definite
Hessian with a free k is reported as the weaker "may not be identified".
The finding is repeated in summary()'s notes and by print() on the
diagnostics object. The check is a heuristic screen rather than a formal
identification test, so no flag means only that nothing was detected.
equation = "zben" numerical Pmax/Omax could return a non-global
maximum that depended on the observed price grid. The zben expenditure
curve on the back-transformed scale can have two local maxima; the engine's
single optimize() call converged to whichever one its bracketing happened
to find, so two fits of the same curve observed through different price
ranges could report different Pmax (e.g. Q0 = 19, alpha = 0.012:
17.5 vs the true 3.64). zben now searches an analytic domain that provably
contains every stationary point of its expenditure curve (all maxima lie
below 4 * max(log10(Q0), 1e-3) / (alpha * Q0)), scanning a dense
log-price grid and refining the grid-local maxima with optimize(), so
the observed price range no longer enters the model Pmax/Omax
(method "numerical_optimize_analytic_domain"; the domain-expansion
search remains only as a fallback when that bound cannot be formed).
Condition under which output differs: zben fits whose expenditure curve is
bimodal (about 5 % of a broad random sweep of (Q0, alpha)) or whose
observed price range excluded the higher peak. The same grid-then-refine
search is used by the hurdle numerical fallbacks and the unconditional
hurdle Pmax/Omax; smooth unimodal curves give the same answer to
optimizer tolerance. Equations with closed forms are unchanged.
residuals(<beezdemand_hurdle>, type = "pearson") is now the Part-II
standardized log-scale residual (log(y) - mu_i) / sigma_e (NA at
zeros). It previously divided the raw-scale response residual by the
log-scale sigma_e, giving values that scaled with the consumption unit.
augment(<cp_model_nls>) returns model-scale residuals. For
equation = "exponential" (fit to log10(y)) .resid was y minus the
log10-scale fitted value; it is now residuals(model) on the model scale,
matching .fitted.
simulate_hurdle_data(n_random_effects = 3, part2 = "koff") now draws
alpha_i = exp(log(alpha) + c_i), the model src/HurdleDemand3RE.h
fits, instead of the additive alpha + c_i (which could generate
increasing curves). Two-random-effect simulations are unchanged
(c_i = 0).
run_hurdle_monte_carlo() also compares the fitted rho_bc_raw to the
partial-correlation raw value that generates the requested rho_bc, not
to atanh(rho_bc).
Non-converged fit_demand_tmb() fits now warn (class
beezdemand_tmb_convergence_warning) at fit time regardless of verbose,
and tidy(), confint(), vcov() and predict() re-warn before
returning numbers from such a fit; previously only a non-PD Hessian was
flagged, so a fit that hit the iteration limit with a PD Hessian gave
silent p-values and intervals.
Multi-start is now the default fitting protocol in fit_demand_fixed()
(TICKET-047). Previously each subject was fit from a single
production-heuristic starting value; a subject whose start led to a
failed or at-a-bound fit was simply reported as non-converged. Now, any
subject whose production-heuristic fit is not strict-converged
(converged_strict: optimizer convergence AND finite coefficients/
objective AND not sitting on a user-supplied bound) is automatically
re-fit from several additional sampled starting values (8 for
2-parameter equations hs/koff/simplified with a fixed k, 32 when
k = "fit"), and the best strict-converged result is kept. Condition
under which output differs from 0.2.0: only subjects whose
production-heuristic fit previously failed to converge or landed on a
bound (some previously non-converged/NA rows may now report a
converged fit). Subjects whose production fit was already
strict-converged are never refit and are byte-identical (this is
guaranteed by construction as well as tested). A sampled starting value is
only ever accepted as a rescue if it is BOTH strict-converged AND
domain-valid (natural-scale Q0 > 0 and Alpha > 0); a sampled start
that only "succeeds" by landing in a domain-invalid region (e.g.
negative alpha) is never preferred over leaving the subject
non-converged. fit_demand_fixed(..., multistart = FALSE) or S = 1
restores the exact legacy single-start behavior. FitCurves() itself is
completely unchanged.
Essential value (EV) in boot_demand() and get_demand_param_emms().
The TMB and NLME tiers computed EV = 1 / (100 * alpha) for every equation
form, silently dropping the k^1.5 term for k-bearing forms and applying a
spurious factor of 100 to the k-free SND form. Both now mirror analyze.R
exactly: EV = 1 / (100 * alpha * k^1.5) for k-bearing forms (exponential /
exponentiated, with k fixed or fitted) and EV = 1 / alpha for the SND
("simplified" / zben) form. Condition under which output differs: every
EV point estimate, CI, and bootstrap draw from these two functions; Pmax,
Omax, Qmax, and elasticity_at_pmax are unaffected. FitCurves() /
fit_demand_fixed() EV values were already correct and do not change. For
NLME fits built with a custom_model_formula the equation form is unknown,
so get_demand_param_emms() now reports EV/LCL_EV/UCL_EV as NA with
a warning instead of applying a guessed formula.
residuals(fit, scale = "natural") was scale-mixed for equation = "zben" in the TMB tier. zben's y_var is the caller-supplied
LL4-transformed response (ll4(y, lambda = 4)); the natural-scale residual
subtracted the natural-scale fitted value from that still-LL4-transformed
y_var instead of from ll4_inv(y_var), so the "natural" residual was
really LL4(y) - Q_hat_natural. Condition under which output differs:
equation = "zben" fits only, and only residuals(fit, scale = "natural"); scale = "model" (the default), fitted(), predict(), and
augment() were already correct. All other equation forms are unchanged
(#18).
equation = "zben" per-subject and group Pmax/Omax used the SND
closed form in the TMB tier. get_subject_pars() (both the fit-time and
the expanded = TRUE computation), calc_group_metrics(), and
boot_demand() passed zben fits through model_type = "snd", i.e.
Pmax = 1 / (alpha * Q0). That closed form assumes SND's (Q0, alpha)
coupling, which does not hold for zben's LL4-scale exponential decay;
stored Pmax was up to ~3.4x too small at the median and rank-inverted
relative to the true expenditure-maximizing price. All four call sites now
numerically optimize expenditure on the back-transformed (ll4_inv()-ed)
natural-scale curve via beezdemand_calc_pmax_omax(model_type = "zben").
That numerical search also adaptively expands beyond the observed price
domain (up to 6 decades) when the true expenditure-maximizing price lies
outside it, rather than silently returning the domain edge as Pmax (the
fix as first landed still had this domain-truncation defect; two curves
observed through different price ranges could report different Pmax for
an otherwise-identical fitted curve). get_subject_pars() and
calc_group_metrics() gain a pmax_at_bound field/column, TRUE only
when the search's expansion cap was hit without finding the true
(interior) maximum, in which case Pmax/Omax are a lower-bound
estimate rather than a converged value; boot_demand() does not gain a
pmax_at_bound column (its output schema is unchanged) but emits a
cli::cli_warn() naming the count of affected bootstrap draws when any
occur among a requested Pmax/Omax/Qmax/elasticity_at_pmax
statistic. Condition under which output differs: equation = "zben" fits
only; Pmax, Omax, and (where reported) Qmax and
elasticity_at_pmax from get_subject_pars(), calc_group_metrics(),
and boot_demand(). All other equation forms are unchanged (#19).
FitCurves() / fit_demand_fixed() batch fitting (legacy fixed-effect
engine). Two interacting defects in the per-subject loop: (1) the default
Q0/alpha start values were computed once from the first subject's data and
then reused (sticky) for every subsequent subject instead of being
recomputed per subject, silently making batch estimates order- and
scale-dependent; (2) when the entire nls fallback chain (wrapnlsr →
nlxb → nls2 brute-force) failed for a subject, an unguarded
fit$m$Rmat() dereference on the resulting try-error crashed the whole
batch call instead of recording that subject as non-converged. Condition
under which output differs from 0.2.0: any batch call (2+ subjects) where
a subject after the first has a different consumption scale, or where any
subject's nls fallback chain fully fails. A batch that previously crashed
now returns one row per subject with the failing subject(s) flagged
non-converged; a batch that previously "succeeded" but silently used a
wrong start value for subjects 2..N may now report different (correct,
order-invariant) estimates for those subjects. Single-subject calls are
unchanged (they never had a sticky prior subject to inherit from).
FitCurves() / fit_demand_fixed() with k = "fit" and
param_space = "log10". The k start value's log10() transform was
applied inside the per-subject loop instead of once before it, so the
transform compounded on every iteration: subject 1 started from
log10(K) (correct), subject 2 from log10(log10(K)) (silently wrong for
typical K < 10), and subject 3+ from log10(<negative>) = NaN,
producing non-converged rows for later subjects. Condition under which
output differs from 0.2.0: batch calls with k = "fit",
param_space = "log10", and 2 or more subjects: subject 1's estimates
are unchanged, subjects 2+ now converge (previously mis-started or
NaN-started). Natural-space and fixed/individual/shared-k paths are
unaffected.
FitCurves() / fit_demand_fixed() reported unverified fallback
endpoints as estimates, with no way to tell a genuine fit from a stalled
one. When wrapnlsr failed and the chain fell back to nlxb, the old
code re-fit that endpoint with nls2::nls2(..., algorithm = "brute-force")
and a single-point start, which is a snapshot rather than a fit; it always
"succeeds" and reports whatever point nlxb stalled at, including points
with a singular Jacobian, as Notes = "wrapnls failed to converge, reverted to nlxb" (indistinguishable from a genuine rescue). Separately,
a numerically converged fit could still land on a physiologically
impossible point (Q0 <= 0 or Alpha <= 0, e.g. for flat or otherwise
degenerate data) and be reported as Notes = "converged" with no flag.
The fallback endpoint is now verified with a genuine iterative refit
(stats::nls(algorithm = "port")) that must itself report
convInfo$isConv; an endpoint that fails this verification is recorded as
a non-converged row (Notes = "endpoint unverified: fallback refit did not converge") instead of being reported as an estimate; a verified
rescue is recorded as Notes = "wrapnls failed to converge; nlxb endpoint verified by port refit". Results now carry converged (the optimizer's
own isConv verdict) and converged_strict (isConv AND finite
coefficients/objective AND not at a user-supplied bound) columns. A
"converged" fit with non-positive Q0 and/or Alpha now also raises a
warning() naming the subject and which parameter is non-positive.
Domain validity is signaled only by that warning: Notes is never
modified and converged_strict is never demoted for a domain-invalid
estimate (it is only reachable in param_space = "natural"; the log10
parameterization's 10^x back-transform is always positive), so a single
subject that converges on the first wrapnlsr attempt keeps byte-identical
Notes/converged/converged_strict regardless of domain validity.
fit_demand_fixed()$results$converged now derives from
converged_strict instead of grepping Notes for failure keywords,
including for domain-invalid-but-numerically-converged fits, which are
therefore reported as converged = TRUE (flagged only by the warning,
not excluded from downstream success counts). Default bounds are
unchanged (still c(-Inf, -Inf)/c(Inf, Inf) unless lobound/hibound
are supplied; this release does not add default non-negativity bounds).
Condition under which output differs from 0.2.0: any subject whose
wrapnlsr fit fails and falls back to nlxb (now either genuinely
verified or reported as non-converged, not a raw snapshot), and any
subject reported "converged" with a non-positive Q0 or Alpha (now also
raises a warning; Notes, converged, and converged_strict are
otherwise unaffected). Subjects that converge cleanly on the first
wrapnlsr attempt are unchanged.
With the brute-force refit gone, nls2 is no longer used anywhere in the
package and has been dropped from Imports. Removing that block also
removed a crash: with constrainq0 set only alpha is free, and the
old code indexed fit$m$Rmat()[2, 2] on a 1x1 matrix outside every
try(), so a single non-fittable subject (for example one with a lone
non-zero consumption point) aborted the whole FitCurves() /
fit_demand_fixed() call with "subscript out of bounds". Such a
subject is now reported as a non-converged row and the other subjects
are fit normally; a regression test covers the individual koff case.
GetValsForSim() (used by SimulateDemand()'s Koffarnus et al., 2015
simulation workflow) misaligned or dropped per-price residuals. Residual
columns are keyed to the global price set (unique(dat$x) order);
resid(fit) is returned in the fitted subject's own row order. The old
code assigned residuals to price columns by POSITION
(dfres[i, 4:NCOL(dfres)] <- resid(fit)), which either errored ("replacement
has N items, need M") when a subject was missing a price row, or (when a
subject had a full but differently-ordered price grid) silently placed
residuals under the wrong price column with no error. Since sdindex
(per-price residual SD, which directly controls simulated variance) is
computed from these columns, a subject whose row order didn't match
unique(dat$x) order silently corrupted sdindex without any warning.
Residuals are now matched to price columns by price VALUE
(adf$x, tolerating a missing price as NA, which sdindex's existing
na.rm = TRUE already handles) instead of by position. A subject whose
fit fails now also raises a warning() naming the subject and cause
instead of silently contributing an all-NA row. Condition under which
output differs from 0.2.0: any call where a subject's row order (within
that subject) doesn't match unique(dat$x) order, or a subject is missing
one or more price rows (previously miscomputed sdindex or crashed
outright; now computed correctly, with NA for missing price cells).
Subjects with a complete price grid in canonical order are unchanged.
Separately, SimulateDemand() now works in a fresh R session that has
never touched the RNG (RunOneSim() previously read .Random.seed
unconditionally, which does not exist until the RNG has been used once).
run_hurdle_monte_carlo() summary statistics can change. See
"Silent-failure fixes" below (TICKET-062): converged replicates with a
non-positive-definite Hessian are no longer counted as valid Monte Carlo
evidence in $summary's bias/coverage calculations.
fit_demand_tmb() rescues an nlminb false-convergence exit. R's
nlminb() reports convergence codes 0 and 1 only and puts the PORT status
in the message, so "false convergence (8)" arrives as code 1. When that
happens the fitter now restarts nlminb from the stalled point, tries
L-BFGS-B, and, for a fixed-k fit, warm-starts from a free-k refit; a
candidate replaces the stalled result only if it reports code 0, its
gradient is small (tmb_control$rescue_grad_tol, default 0.01) and its
Hessian is positive definite, and the lowest-NLL accepted candidate wins.
fit$opt$rescued_from / fit$opt$rescue_method record a rescue and
summary() notes it; tmb_control = list(rescue = FALSE) restores the
loud failure. Fits that converge first time are untouched. The worked
example (apt_full, exponential, gender, fixed k = 2, where the stalled
point has a gradient of order 1e11 and an indefinite Hessian) is in
vignette("convergence-guide").
NLME summary() and tidy() gain df_method. nlme's containment rule
assigns the observation-level residual df to between-subject terms, whose
effective sample size is the number of subjects, so those p-values were
anticonservative and the help could only warn about it. Both methods now
carry a df column (nlme's containment df by default, bit-identical to
before) and accept df_method = "between", which gives each
between-subject term n_subjects - rank(X_between) from that parameter's
subject-level design and recomputes its p-value from the unchanged t
statistic; intercepts and within-subject terms keep containment df. The
printed summary shows the df in use. For column parity,
tidy.beezdemand_tmb() reports df = Inf on its asymptotic z rows.
glance.cp_model_lmer() reports converged. The mixed-effects
cross-price wrapper stored lme4's convergence status but glance() did not
expose it; a converged column (logical, NA on a failed fit or an object
saved before the field existed) is now appended so a batch of fits can be
screened without printing each one, matching broom::glance.nls()'s
isConv.
A hurdle fit records when its random-effects covariance fell back to a
diagonal approximation. .hurdle_chol_or_fallback() warned once at fit
time and left no trace; fit$re_cov_fallback now stores the flag,
summary() and print() carry a note, and marginal predict() output
gets a re_cov_fallback attribute when its own draws used the fallback.
check_demand_model() on a nested-grouping NLME fit mis-read
nlme::VarCorr(). The variance table for random = ~ 1 | outer/inner
interleaves group-header rows with per-level parameter rows whose names
repeat, and indexing by row name turned the headers into NA and reported
the outer level's variance for every level, so the inner level's variance
never appeared and the near-zero flag could be NA. Rows are now walked by
position, header rows dropped, and entries named <level>:<term> when more
than one grouping level is present; single-level fits keep the bare term
names.
NLME residual diagnostics are now guarded like the hurdle and fixed
ones. A residuals() failure or an empty residual vector is reported as a
failed computation (classed warning plus a "could not be computed" issue)
instead of passing as "no outliers". The hurdle and fixed checks also treat
a non-numeric .resid column as a failed computation.
The fixes in this subsection change status/diagnostic output (warnings,
issue lists) rather than point estimates; fits that were correct before
still return the same numbers. Two exceptions, both scoped to NLME
param_space = "natural" fits: the get_demand_param_emms() bullet below
fixes a wrong-by-orders-of-magnitude back-transformation, and the
get_demand_comparisons() bullet (TICKET-075) changes what a natural-space
contrast reports (a difference, not a 10^-exponentiated ratio); param_space = "log10" fits (the default) are unaffected by both.
boot_demand() and confint(method = "simulate") silently repaired an
indefinite covariance. .tmb_parametric_draws() refused a non-finite
fixed-effect covariance but clamped genuinely negative eigenvalues of a
finite one to zero, so draws came from a different, rank-deficient
distribution than the requested asymptotic posterior, with no condition
raised. A materially indefinite covariance is now refused with a classed
error (beezdemand_indefinite_vcov_error); eigenvalues negative only to
within numerical tolerance are still clamped, so a positive-semidefinite
singular covariance keeps working, as does a well-formed covariance whose
parameters differ wildly in scale. "Materially" is judged both against the
covariance's largest absolute eigenvalue and against that of its
correlation-scaled form, so a negative variance in a small-scale parameter
is not masked by a large-scale one. Affected fits are those whose Hessian is
not positive definite; use check_demand_model(fit) and fit$hessian_pd,
and note that confint(method = "wald") rests on the same curvature.
get_demand_param_trends() did not gate on NLME convergence. It was the
only NLME inference surface that never called the convergence guard, so
trends from a fit whose apVar could not be inverted were returned with no
warning. It now warns once, like every other NLME surface.
check_demand_model()$random_effects$variances reported standard
deviations, not variances, for beezdemand_tmb fits (and print()
labelled them "variance"). The log10-scale SDs now live in a new sd_log10
element -- that is the one matching summary(fit)$variance_components -- and
variances holds their squares, so the field name and the printed label are
true. sd_internal_log (raw natural-log-scale SDs, used by the near-zero
degeneracy check) is unchanged, as are all near-zero predicates and the
issues they raise. The NLME path now takes nlme::VarCorr()'s Variance
column explicitly. No estimate changes; code reading $variances from a TMB
fit and expecting an SD should read $sd_log10.
anova(fit, test = "Wald") grouped columns under the wrong model term
when the design had a factor plus a covariate (or two factors). The
term-assignment helper indexed term labels with the raw assign vector;
the intercept's 0 dropped an element and shifted every label by one, so
the default group_by = "auto" table reported e.g. Q0 ~ grp as the joint
test of one factor column with the covariate. Fixed; group_by = "term"
and test = "LRT" were unaffected.
TMB and hurdle inference surfaces now honor hessian_pd. When
TMB::sdreport() reports a non-positive-definite Hessian
(fit$hessian_pd == FALSE), sdr$cov.fixed is a pseudo-inverse of an
indefinite matrix, so standard errors, confidence intervals, p-values, and
parametric draws computed from it are unreliable even though the point
estimates are unaffected. vcov(), confint(), anova() (single-fit
Wald test), get_demand_param_emms(), get_demand_comparisons(), and
boot_demand() now each emit one classed warning
(beezdemand_hessian_not_pd_warning / beezdemand_warning) the first time
they consume such a covariance, for both the beezdemand_tmb and
beezdemand_hurdle classes. Previously only summary()'s print method and
check_demand_model() surfaced this; a user calling confint() or
boot_demand() directly received unreliable intervals with no indication
anything was wrong. Values are unchanged; healthy (PD-Hessian) fits emit no
new conditions.
NLME inference surfaces now honor .check_nlme_convergence().
summary(), glance(), and check_demand_model() already gated on
whether an NLME fit's final apVar inverted cleanly; get_demand_param_emms(),
get_demand_comparisons(), calc_group_metrics(), confint(),
get_subject_pars(), tidy(), get_individual_coefficients(), and
anova() (per compared model) now each emit one classed warning
(beezdemand_nlme_convergence_warning / beezdemand_warning) on a
non-converged fit instead of computing inference silently.
calc_group_metrics() also replaces a blanket
suppressWarnings(suppressMessages(...)) around its internal Q0/alpha
EMM calls with targeted muffling of two specifically-matched benign
conditions, so a real warning raised inside those calls (e.g. the
estimate-column-guess fallback, which flags a possibly-wrong Pmax/Omax)
now reaches the caller instead of being silently dropped.
get_demand_param_trends() now warns once, naming every dropped
(parameter, covariate) combination and its cause, instead of silently
shrinking the returned table. Values are unchanged; converged fits emit
no new conditions. See also the get_demand_param_emms() /
param_space = "natural" bullet below, which fixes a related but
distinct back-transformation bug on the same NLME surface.
get_demand_param_emms() back-transformed with 10^ unconditionally,
giving wrong *_natural columns (and EV) for NLME fits made with
param_space = "natural". fit_demand_mixed(..., param_space = "natural") fits Q0/alpha directly on the natural scale (supported for
equation_form = "simplified"/"exponentiated"), so the raw emmeans
summary is already natural-scale; exponentiating it again with 10^
inflated Q0_natural/alpha_natural by orders of magnitude and
propagated into EV, with no warning. get_demand_param_emms() now
resolves param_space the same way every other NLME surface does and,
for a natural-space fit, uses the emmeans summary directly for the
*_natural columns and fills *_param_log10 with log10() of those
values (keeping the same column set across both spaces). Because a
natural-space fit is an unconstrained parameterization, a Wald CI bound
(or, rarely, the point estimate) can be non-positive; *_param_log10 is
NA (not NaN, and without a raw "NaNs produced" warning) wherever the
corresponding *_natural value is <= 0; *_natural itself is never
affected. Condition under which output differs from 0.2.0: beezdemand_nlme
fits made with param_space = "natural", only in get_demand_param_emms()
(and anything built on its EV branch). param_space = "log10" fits (the
default) and the TMB tier (already space-aware) are unaffected.
get_demand_comparisons() (NLME) had the same param_space = "natural"
gap as the bullet above, in $contrasts_ratio. (TICKET-075.)
emmeans::contrast()'s estimate/CI in $contrasts_log10 are on the
fit's internal scale (log10 for param_space = "log10", natural for
param_space = "natural"); $contrasts_ratio always computed
ratio_estimate = 10^estimate to turn a log10-scale difference into a
multiplicative fold-change, which is meaningless for an already-natural-scale
difference. For param_space = "natural" fits, $contrasts_ratio now
reports the difference again (same column names/shape:
ratio_estimate/LCL_ratio/UCL_ratio) instead of exponentiating it a
second time; the returned object's new contrasts_ratio_scale attribute
("ratio" or "difference") says which content a given call got.
$contrasts_log10 itself was already correct (unaffected); only
$contrasts_ratio's content for natural-space fits changes, from a
previously-meaningless number to a documented, correct one.
param_space = "log10" fits (the default) are unaffected.
check_demand_model() no longer reports a failed internal check as a
passing one. The fixed and hurdle residual sub-checks, and the NLME
random-effects sub-check, converted an internal augment() / VarCorr()
error (or a missing/all-NA .resid column) into the same clean-looking
"no outliers found" / "nothing near zero" result the check returns when it
actually ran and found nothing, so a report that never examined
residuals or random-effect variances printed "No issues detected"
indistinguishably from one that genuinely checked and passed. Each
sub-check now sets an explicit computation_failed flag, raises one
classed warning naming the cause, and check_demand_model() adds a
"...could not be computed" issue instead of silently passing. Mirrors the
pattern the TMB tier's residual check already used. Healthy fits are
byte-identical; only fits where one of these internal checks errors are
affected.
fit_demand_tmb() failure paths now name their causes. Three gaps on
the failure/error-reporting side of the TMB fitter:
message()d
only at verbose >= 2; optimizer-sentinel causes in opt$message were
never re-read once a start was rejected). It now appends a Causes:
section naming at least one underlying cause per failed start,
regardless of verbosity.TMB::sdreport() failure now warns regardless of verbose
(previously gated on verbose >= 1, so a verbose = 0 fit failed
completely silently); the warning includes the fallback attempt's
message when it differs from the first; the same classed warning now
also fires on ADREPORT/variance-component extraction failure (a
previously fully silent path); and hessian_pd = NA (meaning
"unknowable because sdreport failed") is now explained via tidy()'s
hessian_warning attribute in addition to summary()'s print note.equation = "exponential"'s
zero-consumption filter (e.g. all-zero y) now aborts immediately with
an informative message naming the equation and the dropped-row count,
instead of proceeding to "fit" a 0-observation model, reporting a
spurious "Converged (NLL = 0.00)", and then crashing during SE
extraction with a cryptic no 'dimnames' attribute for array (the only
output a verbose = 0 caller saw). Healthy fits are unaffected; no
numeric output changes.print()/summary() for beezdemand_hurdle fits now surface a
false-converged 3-random-effect fit prominently. The 3RE spec
(random_effects = c("zeros", "q0", "alpha")) can converge according to
nlminb()'s reported code while the Hessian is not positive definite on
real purchase task data (weak identification of the alpha random effect,
not a broken spec); fit$opt$message/fit$opt$convergence and
fit$hessian_pd already existed but were easy to miss. Both print methods
now show a warning block (quoting the optimizer message) whenever
converged is FALSE or hessian_pd is FALSE, naming the recommended
stability check: refit with random_effects = c("zeros", "q0") and
compare empirical-Bayes subject parameters.
summary.beezdemand_hurdle()'s notes field (previously computed but
never printed by print.summary.beezdemand_hurdle()) now reaches the
console. Converged, PD-Hessian fits print byte-identically to before. No
change to converged/hessian_pd semantics, EB parameter extraction, or
the TMB templates/likelihoods.
confint(fit, method = "simulate") and boot_demand() on a TMB fit whose
sdreport covariance contains non-finite values (non-PD Hessian /
non-converged fit, platform-dependent) now stop with a classed,
informative error (beezdemand_nonfinite_vcov_error) pointing at
check_demand_model() and fit$hessian_pd, instead of an opaque
eigen() "infinite or missing values in 'x'" failure.
Cross-price method contracts. tidy(), glance() and augment() on a
cp_model_nls / cp_model_lm / cp_model_lmer object whose fit failed
(model = NULL) now return typed empties with the same columns as the
successful call (zero rows for tidy()/augment(), one all-NA row for
glance()) instead of six different shapes. tidy.cp_model_lmer() accepts
effects = "ran_vals" and keeps "random" as an alias (it used to pass
argument matching and then error inside broom.mixed). fit_cp_linear(type = "mixed") stores lme4's convergence verdict (converged,
convergence_messages), warns once with class
beezdemand_cp_lmer_nonconverged_warning when it is bad, and print() /
summary() show it; print.cp_model_nls() shows a non-converged winning
fit. The summary print now describes qalone as the asymptote reached as
the alternative's price grows (the old text called it the zero-price value).
calc_group_metrics() on an NLME fit recorded a covariate value it had
not conditioned on. A multi-value continuous at entry was forwarded
whole to emmeans (the grid expanded over every value) while
conditioned_on reported only the first. The NLME method now warns and
uses the first value, as the TMB method already did. Both methods'
documentation now states the marginalisation grid each backend uses and
when the two can differ.
predict.beezdemand_tmb() rebuilt its design under the contrasts in force
at call time rather than at fit time. Changing options("contrasts")
after fitting a model with factors kept the column count but changed the
basis, so population- and subject-level predictions silently changed. The
rebuilt fixed-effect design is now pinned to the fitted contrasts
attribute and verified column-by-column (aborting loudly on a mismatch),
the same rule the EMM grid already used. The factor-expanded random-effect
design stores the fit-time contrasts on the parsed block and reuses them,
and the anova() term-map fallback does the same. Default-contrast sessions
are unaffected. In the same area, predict() on a collapse_levels fit
rejected newdata in the original shape because the fitted factor columns
are internal (<factor>_Q0 / <factor>_alpha); those columns are now
recreated from the original factor using the training data's level map.
Hurdle random-effects covariance chol() failure silently substituted an
uncorrelated diagonal Sigma at five sites (the RE-transform helpers, the
live fit_demand_hurdle() inline path for 2- and 3-RE models, and
.compute_marginal_demand()'s Monte Carlo draws). When the assembled
covariance was not positive definite (near-boundary rhos, overflow, or
tanh(raw) rounding to exactly +/-1), subject-level effects and marginal
demand curves were silently computed from zeroed correlations while the
reported correlations/summary rho estimates continued to show the fitted
values. All five sites now share .hurdle_chol_or_fallback(), which emits
one classed warning (beezdemand_hurdle_chol_fallback_warning) when the
fallback fires. Condition under which output differs: only calls where the
assembled Sigma is not positive definite (rare given the partial-correlation
parameterization); healthy (PD) fits are bit-identical and silent
(TICKET-061).
run_hurdle_monte_carlo() discarded per-replicate diagnostics and
counted non-PD-Hessian fits as valid Monte Carlo evidence. Failed and
non-converged replicates collapsed to NULL with no record of why; the
return value carried no per-replicate status at all. Separately, a
replicate that converged with hessian_pd = FALSE passed the same filter
as a clean fit and contributed its (unreliable) estimate and SE to the
bias/coverage summary. The return value now includes $diagnostics (one
row per replicate: sim_id, status ("error", "nonconverged",
"converged_non_pd", "converged_hessian_unavailable", or "clean"),
converged, hessian_pd, opt_convergence, opt_message) and
$n_hessian_not_pd/$n_hessian_unavailable; $estimates gains a
hessian_pd column. hessian_pd = NA (i.e. sdreport() itself failed) is
kept distinct from an explicit hessian_pd = FALSE, since they are
different conditions, though both are excluded from $summary the same
way. A classed warning (beezdemand_hurdle_mc_hessian_excluded_warning)
fires naming both excluded counts (e.g. "1 non-PD, 1 Hessian unavailable")
when either happens. This changes $summary output for any prior run
that had converged-but-non-PD-Hessian or Hessian-unavailable replicates
(TICKET-062).
fit_cp_nls() discarded all backend convergence diagnostics, and
summary.cp_model_nls()/confint.cp_model_nls() reported SEs, p-values,
and CIs with no convergence gate. A maxiter-capped or otherwise
non-converged nlsLM/nls.multstart fit produced a clean-looking
coefficient table with zero warnings. fit_cp_nls() now returns a
convergence field (isConv, finIter, stopCode, stopMessage for the
winning backend, read from model$convInfo; populated for nls/nlsLM
fits and for nlsr::wrapnlsr() fits too when it returns a plain
nls-class object, its usual successful case; isConv = NA only when the
winning fit carries no convInfo at all). summary() and confint() now
emit a classed warning
(beezdemand_cp_nls_nonconverged_warning) when isConv is explicitly
FALSE. Separately, nlsLM_fit is now NULL (never the caught error
condition) in the branch where wrapnlsr won after nlsLM failed
(TICKET-065).
Requested plot annotations, extractor rows, and a summary CI section
disappeared with no condition when their computation errored.
plot_expenditure() (hurdle and TMB): with show_pmax = TRUE /
show_omax = TRUE explicitly requested, a calc_group_metrics() error
silently omitted the annotation; now warns
(beezdemand_plot_annotation_warning) naming the cause, matching the
existing pattern used by plot_compare(). coef.beezdemand_fixed():
subject rows for a stored try-error/NULL fit, or one whose coef()
call failed, were dropped with no way to tell "subject absent" from
"subject failed"; now warns (beezdemand_fixed_coef_omitted_warning)
naming the dropped ids. augment.cp_model_nls() /
augment.cp_model_lm() / augment.cp_model_lmer(): the documented
.fitted/.resid/.fixed columns vanished on a fitted()/
residuals()/predict() error or a length mismatch with no indication;
now warn (beezdemand_cp_augment_omitted_warning) naming the omitted
column(s). summary.cp_model_nls(): a genuine nlstools::confint2()
error (as opposed to nlstools simply not being installed, which stays
silent) silently dropped the confidence-interval section; now warns
(beezdemand_cp_summary_ci_omitted_warning). Healthy paths are
unaffected and remain silent throughout (TICKET-068).
FitCurves(equation = "linear") (and fit_demand_fixed(equation = "linear")) now degrades a per-subject fit failure to an NA-parameter row
with an informative Notes message instead of crashing. The linear
extractor (ExtractCoefs.linear()) was missing the try-error guard its
nonlinear sibling (ExtractCoefs()) already has; on a failed wrapnlsr
fit it dereferenced the resulting try-error immediately
(coef(fit)[c("l", "b", "a")]), which raised $ operator is invalid for atomic vectors and aborted the entire batch with no per-subject failure
record; reproducible even for a single unfittable subject called alone.
The extraction is now wrapped end-to-end, so a failure later in the
extraction (summary(), nlstools::confint2(), deviance()) produces
an NA-parameter row the same way a failed initial coef() call does.ExtraF() now reports which group's per-group fit failed
("ExtraF: unable to fit group '<name>': ...") instead of an opaque
no applicable method for 'predict' applied to an object of class "try-error". GetSharedK()'s shared-k search (start-value grid
construction and the final nlxb() fit, neither of which was previously
guarded) is now wrapped so any internal failure reaches its designed
sentinel return ("Unable to find a shared k.") instead of escaping raw;
FitCurves(k = "share") can therefore actually reach its documented
fallback to GetK() with a warning, which previously could not fire
because the sentinel path was unreachable (the final nlxb() call was
never wrapped in try(), so the inherits(fit, "try-error") check after
it could never be true). A shared-k group set that drops below 2 usable
groups (after the existing <3-row-per-group drop) now returns an
informative sentinel instead of proceeding into a nonsensical
single-group contrast.simulate_hurdle_data(part2 = "snd") (TICKET-044). Adds an SND
positive-part generator to the hurdle simulator, matching
src/HurdleDemand3RE_SND.h / src/HurdleDemand2RE_SND.h exactly: a
log-linear (no k) mean with lognormal errors on the positive part, and
the same zero-inflation logistic as the existing (now part2 = "koff")
generator. Random-effect correlations for part2 = "snd" are specified
via rho_ab_raw/rho_ac_raw/rho_bc_raw, mirroring the TMB model's own
raw-parameter coefficients exactly (rho_ab = tanh(rho_ab_raw), rho_ac = tanh(rho_ac_raw), rho_bc via the LKJ-Cholesky partial-correlation
transform), so a fitted fit_demand_hurdle(part2 = "snd") model's own
coefficients can be plugged directly into the simulator for a parametric
bootstrap or recovery study. part2 defaults to "koff", and its output
is byte-identical to previous releases.
simulate_hurdle_data(seed = ) and run_hurdle_monte_carlo(seed = ) no
longer overwrite the caller's RNG stream: the global .Random.seed is
restored on exit (or removed if none existed), matching what
power_demand(), boot_demand() and the TMB parametric-draw helpers
already guaranteed. Simulated outputs for a given seed are unchanged.
fit_demand_tmb() now treats a continuous within-subject covariate as a
random slope in its own right (dose-response demand). Specify it with a numeric
random-effects term, e.g.
random_effects = nlme::pdSymm(Q0 + alpha ~ dose_c), and pair it with
continuous_covariates = "dose_c" for the population (fixed) dose slope.
Each subject's intensity and elasticity then change with the covariate at
their own rate, with far fewer parameters than the equivalent multi-level
factor random slope.get_subject_pars() and predict(type = "parameters") expose the
per-subject slope deviations (q0_<term> / alpha_<term>, matching
ranef()) and accept an at = argument to evaluate per-subject Q0/alpha
at a chosen covariate value (e.g. at = c(dose_c = 0)); the two surfaces now
agree and no longer return NA.summary(), VarCorr(), and tidy(effects = "ran_pars") label the
continuous-term variance component and its intercept correlation by the
covariate name rather than a positional index.check_demand_model() flags a near-singular intercept/slope
covariance.predict.beezdemand_hurdle() now defaults to type = "demand" (the
marginal expectation (1 - p0) * E[Y | Y > 0]) instead of
type = "response" (the conditional positive mean E[Y | Y > 0]).
Rationale: observed consumption includes zeros, so scoring predictions
against raw data with the conditional mean systematically overstates
hurdle error wherever the probability of zero consumption is large (a
documented cross-validation erratum). This is a deliberate breaking
change from 0.2.0: the marginal expectation is the statistically
correct default for prediction scoring, and a clean flip (with
communication) was judged better than keeping the misleading default for another
release. Omitting type emits a once-per-session message naming the
change; every type's computation is unchanged, so
type = "response" restores the old behavior exactly. See the new
"Scoring predictions" section in ?predict.beezdemand_hurdle.TMB::sdreport()/summary.sdreport() SE
extraction (one per offending block). They are replaced by at most ONE
classed warning per fit (beezdemand_sdreport_warning, with parent
class beezdemand_warning) pointing at $hessian_pd and summary()
diagnostics. Matching is on the warning's sqrt() call rather than the
locale-dependent message text, and unrelated warnings pass through
untouched. The existing "Standard error computation failed" warning now
carries the same classes. No fit results, SEs, or hessian_pd logic
changed.Bug fix. lambertW() now accepts vector input. The Halley-iteration
convergence check combined possibly-length-n operands with &&, which
errors under R >= 4.3 ("length > 1 in coercion to logical(1)"); it now uses a
vectorized all(... &). Scalar results are unchanged, so existing scalar
callers are unaffected. This hardens every analytic Pmax/Omax path that
vectorizes over price.
Bug fix. get_empirical_measures() now computes the breakpoint BP0 and
the empirical Pmaxe order-invariantly by sorting each series by price before
walking consumption. Previously a subject whose rows were not in
ascending-price order could receive an incorrect breakpoint; results are
unchanged for already price-ordered input (the typical case). The deprecated
legacy path GetEmpirical() received the same per-subject sort, so the two
functions stay in agreement on unsorted input.
Consistency fix. The analytic Pmax/Omax existence threshold is now a
strict interior maximum across all equations. Both hurdle variants previously
returned a finite Pmax at the exact boundary k = e; they now return NA
there, matching the Hursh & Silberberg path (already strict at
k = e / ln(10)). At the boundary the two stationary points merge into a
tangent inflection, so no strict interior maximum exists. The practical effect
is measure-zero. The deprecated numerical helper calc_omax_pmax() now
classifies k = e the same way: it takes its bounded-range fallback (with
the usual warning and note) instead of the analytic root path.
Reporting change (broom convention). summary() and tidy() now report the
Wald statistic and p.value on the estimation scale (log/log10; logit for
hurdle participation rows) for every report_space, including the default
"natural"; only estimate and
std.error are back-transformed. Previously the natural scale recomputed the
test from the back-transformed estimate/SE, which is degenerate for log-scale
parameters: the statistic reduces to 1/(c·SE), independent of (and dropping
the sign of) the estimate, and for factor effects tests an impossible null
(ratio = 0 rather than ratio = 1). Keeping the test on the estimation scale
matches broom/emmeans/glmmTMB: for NLME the natural-scale test is nlme's native
containment-t (DF-aware) test; for the TMB and hurdle tiers it is the pnorm()
z-test. This changes the default reported statistic/p.value for the core
demand parameters (Q0/alpha/k); estimate and std.error are unchanged.
Bug fix. summary(fit_nlme)$converged now reflects the same operational
convergence gate as glance(fit_nlme)$converged (positive-definite apVar
and no terminal error; TICKET-020). It previously hard-coded TRUE, so a
summary could report convergence on a fit that glance() correctly flagged as
unusable for inference. Any diagnostic message now appears in the summary's
notes.
compare_specs_used attribute on
get_demand_comparisons(fit_nlme, ...) output now records the user's
requested comparison spec (or "all fitted factors" when none is given),
matching the TMB backend. It previously stored the last parameter's
internally-derived per-parameter formula, which under asymmetric
collapse_levels was both misleading (it could show a collapsed
intercept-only ~1, or an internal collapsed column name) and dependent on
the order of param.VarCorr(fit_tmb) now places random-effect correlations on the
correct rows for multi-block pdBlocked fits. The Corr entries were
previously positioned using indices local to each covariance block, which was
correct only when the correlated (pdSymm) block came first; a correlated
block in any later position was placed on the wrong row. Positions are now
derived from the same block map used to compute summary(fit)$correlations,
adding each block's global row offset. Single-block fits are unaffected.New method. get_subject_pars() now has a beezdemand_nlme method,
closing the last cross-backend asymmetry (every other demand-model extractor
already had one). It returns the same columns as
get_subject_pars.beezdemand_tmb(): wide form id, b_i, c_i (if alpha
has random effects), Q0, alpha, Pmax, Omax; expanded form adds the
within-subject factor column(s) with one row per (subject, factor-level)
cell. Q0, alpha, Pmax, Omax are on the natural scale.
Honors the same expanded contract as the TMB method: NULL (default)
auto-detects within-id variation and expands when present; TRUE always
attempts expansion; FALSE returns the wide shape and warns (with affected
subjects' Q0/alpha/Pmax/Omax set to NA) when within-id variation
exists. Per-subject parameters are reconstructed from the fixed effects plus
each subject's random-effect deviations and back-transformed from the model's
param_space (log10 by default) to the natural scale.
The b_i / c_i random-effect aliases are reported on the natural-log
linear-predictor scale to match the TMB method (the log10 deviations are
multiplied by log(10) for param_space = "log10"). The full
per-coefficient random effects remain available via ranef(). The
backend-agnostic expansion scaffold is currently duplicated from the TMB
method; unifying the two into a shared internal is a planned follow-up.
Bug fix. glance(fit_nlme)$converged no longer flips to FALSE merely
because nlme emitted iteration-level convergence warnings (false convergence,
singular precision matrix, step-halving, iteration limit, ...) during its
PNLS-LME alternation. nlme prints these routinely while iterating even when the
fit settles to a usable optimum, so the old heuristic produced false negatives
on fits that were perfectly usable for inference (positive-definite apVar,
valid emmeans/contrasts/predictions, no terminal error). converged now
reports the operational gate: apVar is positive-definite and there is no
terminal error.
glance(fit_nlme) gains two additive NLME-only columns: final_fit_ok
(the canonical usable-for-inference gate; identical to converged) and
fit_warned (a diagnostic flag, TRUE when iteration-level convergence
warnings were emitted; informational only, does not gate converged).
Breaking change (semantic). Some fits previously flagged
converged = FALSE (despite a usable apVar) now correctly read TRUE. Code
that scripted around the old false-FALSE behavior should gate on
final_fit_ok (the canonical check) and inspect fit_warned for iteration
diagnostics. check_demand_model(fit_nlme) likewise no longer reports
"Model did not converge" for warned-but-usable fits. The genuine-failure path
is unchanged: a non-positive-definite apVar (or a terminal error) still reads
converged = FALSE.
Breaking (NLME, narrow): when contrast_by conditions on a factor that
was asymmetrically collapsed via collapse_levels, the nested
$contrasts_log10 / $contrasts_ratio by-column(s) now use the
user-requested original factor name (e.g. age_group) instead of the
collapse-mapped name (e.g. age_group_alpha). This completes the TICKET-016
harmonization: the NLME nested object now matches the TMB backend and the
flat tidy() output, all of which already used the original name. The
by-column values (the collapsed level labels) are unchanged. Only code
that reads the suffixed nested by-column directly
(comps$alpha$contrasts_log10$age_group_alpha) is affected; replace with the
original name (...$age_group). tidy() output and the uncollapsed-fit case
are unaffected.
Two robustness guards were added to the NLME contrast_by path: supplying two
by-variables that resolve to the same effective column (e.g.
c("age_group", "age_group_alpha") under collapse) now aborts, mirroring the
TMB backend; and a factor name that would collide with a reserved contrast
column during the nested rename (e.g. a factor literally named estimate) now
aborts with a clear message rather than failing cryptically (this collision
guard is NLME-specific because it arises from the rename introduced here).
get_demand_comparisons() now supports contrast_by on the TMB backend,
completing the NLME/TMB harmonization begun in TICKET-016. Within each
observed combination of the by-level(s), pairwise (or trt.vs.ctrl)
contrasts are computed over the remaining factors, with p-value adjustment
applied per by-cell. Results match the NLME backend in shape, direction,
by-cell labels, and message UX (numerics differ by design: TMB uses
asymptotic z, NLME t). A by-cell of a single contrast reproduces the
corresponding at = -filtered call exactly.
Nested contrast tables ($contrasts_log10, $contrasts_ratio) and the flat
tidy() frame gain leading by-column(s) using the user-requested original
factor name (e.g. age_cut, not the collapse-mapped age_cut_alpha). The
tidy() schema is unchanged (the canonical 9 columns) when contrast_by is
inactive. print() shows the by-column(s) before contrast.
Both backends now populate a contrast_by_map attribute (a per-parameter
named map from the original by-name to the effective, possibly
collapse-mapped, column). contrast_by resolution is soft per parameter
(a by-variable absent from a parameter's design is skipped) whereas
compare_specs resolution remains strict (an unresolvable name aborts).
Behavior change (NLME): supplying a contrast_by factor that is not in
compare_specs now aborts loudly with a parameter-scoped message, on both
backends. Previously the NLME backend silently returned an empty contrasts
table (with a $contrasts_log10_error note). No released code relied on the
silent-empty path.
get_demand_comparisons() now returns a classed beezdemand_comparison
object on both the NLME and TMB backends, and a new
tidy.beezdemand_comparison() method gives a backend-agnostic flat
contrasts frame (param, contrast, estimate, std.error, statistic,
df, conf.low, conf.high, p.value) with identical columns regardless
of backend. Estimates and CIs are reported on the log10 scale on both
backends; tidy(res, exponentiate = TRUE) returns base-invariant ratios.
(Per-backend inference is unchanged: NLME reports a t statistic with finite
df, TMB an asymptotic z with df = Inf.)
get_demand_comparisons() now compares both Q0 and alpha by default
on both backends (the TMB backend previously returned Q0 only).
The TMB backend gains compare_specs, a formal at argument,
report_ratios, factor-level contrast ordering (previously
data-appearance order, which could flip signs when input rows were
reordered), and equal-weight marginalization over omitted factors
(averaging across the full crossing of their levels, emmeans' default
weights = "equal", matching the NLME backend). get_demand_param_emms()
for TMB fits likewise marginalizes when factors_in_emm names a subset of
the fitted factors, rather than erroring.
Behavior change (NLME): the default p-value adjustment is now "holm"
(was "tukey"); pass adjust = "tukey" to retain the previous default.
Rationale: cross-backend reproducibility and the base-R pairwise default.
Deprecation (NLME): get_demand_comparisons(params_to_compare = ) is
deprecated in favor of param. Supplying both is an error.
TMB API (new in this release, so no released code is affected):
get_demand_comparisons.beezdemand_tmb() renames p_adjust to adjust
(no alias) and validates it against stats::p.adjust.methods; emmeans-only
methods (e.g. "tukey", "sidak") are rejected.
contrast_by (by-grouped contrasts) is now supported on the TMB backend
(TICKET-032, see below), completing the backend harmonization.
Correctness fixes (post-review). (1) The NLME backend now validates
compare_specs against the union of the model's fitted factors (originals
and per-parameter collapsed columns) and errors on names not in that union,
matching the TMB boundary check; cross-parameter aliases that pass the
boundary but cannot resolve for a given parameter (e.g. ~ age_group_alpha
with param = "Q0") abort with a parameter-scoped message rather than
silently producing an intercept-only EMM. Factors whose collapsed column has
fewer than 2 levels for a parameter still return empty contrasts for that
parameter without error (the existing intentional behavior). (2) Under
asymmetric collapse_levels, naming the original factor in
compare_specs/factors_in_emm (e.g. ~ age_group) now resolves to that
parameter's collapsed column on the TMB backend (age_group_Q0 /
age_group_alpha), as it already did on NLME, instead of silently returning
zero contrasts; a name that cannot be resolved for the parameter is rejected
with an error. (3) TMB EMM and contrast reference grids are now built using
the fitted design's contrasts, so changing the global
options("contrasts") between fitting and calling no longer silently changes
the estimates.
fit_demand_tmb() fits are now substantially smaller on disk by default
(often >80% smaller via saveRDS() on large datasets). The full
covariance matrix of all ADREPORT'd quantities ($sdr$cov), which is read
by no method, is no longer materialized; $sdr$cov is a scalar NA unless
the new store_report_cov = TRUE argument is supplied. Standard errors,
cov.fixed, variance components, and all inference (coef(), vcov(),
tidy(), confint(), get_demand_param_emms(), get_demand_comparisons(),
boot_demand()) are unchanged. Pass store_report_cov = TRUE to restore
the prior behavior.calc_group_metrics() gains a beezdemand_nlme method, giving
fit_demand_mixed() fits the same population-level metric summary
(Pmax, Omax, Qmax, elasticity_at_pmax) already available for
beezdemand_tmb fits. Fixed-effect Q0 and alpha estimated marginal
means are marginalized parameter-first (a geometric mean across the
reference grid; continuous covariates at their training mean and factor
levels equally weighted, unless conditioned via at), then the scalar
metrics are derived through the shared Pmax/Omax engine. Returns the same
flat list (Pmax, Omax, Qmax, elasticity_at_pmax, method,
conditioned_on) as the TMB method.boot_demand() computes parametric-bootstrap confidence intervals on
derived demand metrics (Pmax, Omax, Qmax, EV,
elasticity_at_pmax) for beezdemand_tmb fits. Draws of the
fixed-effect parameter vector are taken from the joint asymptotic
posterior, mapped to per-condition (Q0, alpha, k) through the model
design, and summarized by empirical quantiles. Returns one row per
(statistic, condition) (statistic, condition, estimate,
conf.low, conf.high, level); the per-cell point estimate matches
calc_group_metrics(fit, at = cell). Uncertainty in k is propagated
when k is estimated. NLME fits and nonparametric resampling are
planned for a follow-up.get_subject_pars() on beezdemand_tmb fits now auto-detects fits
with within-id-varying design columns (factor-expanded random
effects, within-id continuous covariates, or multi-block
pdBlocked specs) and runs the expansion machinery by default.
The returned shape depends on the kind of within-id variation:
for fits with within-id factors, rows are expanded across factor
levels (one row per (subject, factor-level) cell with per-cell
Q0, alpha, Pmax, Omax); for fits whose only within-id
variation is in numeric covariates, numerics are conditioned at
the subject's mean and the return is one row per subject with
finite (non-NA) Q0/alpha. Previously the default returned the
wide one-row-per-subject shape with NA in Q0, alpha, Pmax,
and Omax for affected subjects, which left users with no next step. The new default
signature is expanded = NULL (auto-detect); pass expanded = TRUE
or expanded = FALSE for explicit override. For fits without
within-id variation the behavior is unchanged (the auto-detect
path resolves to the wide shape).get_subject_pars(fit, expanded = FALSE) on a fit with within-id
variation now emits a one-line warning to flag that the returned
Q0 / alpha / Pmax / Omax columns are NA. Pre-change this
case was silent.expanded = FALSE explicitly (and suppressWarnings() the new
one-line warning) or check any(is.na(fit$subject_pars$Q0))
directly. The fit-time 4-line warning at fit time
(R/tmb-demand.R) is unchanged.coef() on beezdemand_tmb fits gains
type = c("internal", "subject", "combined", "fixed"). The default
("internal") is unchanged: it still returns the raw optimizer
coefficient vector, so fixef() and tooling that dispatches via
coef() (e.g. car::deltaMethod, multcomp::glht) are unaffected.
type = "subject" (alias "combined") returns the per-subject
parameter tibble (get_subject_pars(fit), auto-detecting within-id
factor expansion); type = "fixed" returns a one-row tibble of the
fixed-effect coefficients. Supplying report_space through ... is an
error (no scale conversion in coef()).confint() on beezdemand_tmb fits gains
method = c("wald", "simulate"). The default ("wald") is
unchanged. "simulate" draws R parametric Monte Carlo samples
(default R = 1000, with an optional seed) from the joint asymptotic
Gaussian posterior N(coef(fit), vcov(fit)) and reports per-coefficient
empirical quantiles. It is diagnostic: the simulated intervals are
asymptotically Wald-equivalent on per-coefficient CIs (useful as a
side-by-side check on the Gaussian approximation), do not improve on
Wald at boundary cases, and carry no positivity guarantee on the
internal scale. No new package dependency is added.fit_demand_tmb() now validates that every column of the fixed-effect
design matrix is constant within each id. When a factor or continuous
covariate varies within subject, subject_pars$Q0, $alpha, $Pmax,
and $Omax are set to NA_real_ for affected subjects and a
cli::cli_warn() names the offending columns. Previously the function
silently returned row-order-dependent values. New
validate_subject_pars = TRUE argument provides an escape hatch for
users who have reasoned about the behavior. Factor-expanded random
slopes landed in TICKET-011 Phase 2 (single-block) and Phase 3
(multi-block); the silent first-observed-row fallback was retired
in Phase 5A in favor of NA-on-within-id-variation plus an opt-in
expanded argument on get_subject_pars().get_demand_param_emms() on a beezdemand_tmb fit now honors
continuous_covariates even when no factors are present. Previously,
the early-return for factor-less models ignored the at argument and
always returned the intercept.get_demand_param_emms() on a beezdemand_tmb fit now raises a clear
error when factors_in_emm drops any fitted factor. Previously, the
shorter reference row was silently recycled against the full
coefficient vector, either producing wrong numbers or crashing
downstream with a generic "non-conformable arguments" error. Proper
marginalization over omitted factors lands in TICKET-011 Phase 5.fit_demand_tmb() now accepts formula-based random_effects
arguments. The default signature is random_effects = Q0 + alpha ~ 1,
equivalent to the previous c("q0", "alpha") default. Single-parameter
Q0 ~ 1, pdMat objects, lists of pdMat, and nlme::pdBlocked are
all parsed and attached to the fit object's
$param_info$random_effects_parsed as a canonical block representation.fit_demand_tmb() gains covariance_structure = c("pdSymm", "pdDiag")
matching the fit_demand_mixed() argument of the same name.random_effects (e.g. c("q0", "alpha"))
are soft-deprecated via lifecycle::deprecate_soft() and internally
translated to the equivalent formula. A hard deprecation follows in
0.4.0.Q0 + alpha ~ condition, pdBlocked(list(...))) are now fully
supported. Template generalization to a Z-matrix-driven covariance
landed in TICKET-011 Phase 2 (single-block factor expansion) and
Phase 3 (multi-block pdBlocked).R/random-effects-utils.R
(.classify_re_input, .normalize_re_input, .validate_re_input,
.re_is_phase1_fittable, .deprecate_character_re,
.re_shape_summary, .re_parsed_to_character) factor out the
formula/pdMat parsing so both fit_demand_mixed() and
fit_demand_tmb() consume the same canonical representation.Adversarial review surfaced two more silent wrong-answer paths in the TMB post-fit layer; both are sister bugs to Phase 0.2 / 0.3 fixes and land here as the foundation for the Phase 2 factor-RE work.
get_demand_comparisons() on a beezdemand_tmb fit now consumes the
same conditioned reference grid as get_demand_param_emms(). New
internal helper .tmb_build_emm_ref_grid() in R/tmb-methods.R
ensures both functions honor at (factor-level filters AND
continuous-covariate value overrides) and factors_in_emm. Before
this fix the wrapper forwarded ... to emms but rebuilt its own
contrast grid from the unfiltered training data, producing off-grid
contrasts and "NA" labels when at filtered factor levels.calc_group_metrics() on a beezdemand_tmb fit now warns and
returns a conditioned_on field when continuous covariates are
present. The numeric output is unchanged (intercept-only Q0 /
alpha, i.e. covariates held at 0) but the warning matches the
convention from predict(type = "demand") so users cannot
silently misread reference-intercept metrics as population means.
summary.beezdemand_tmb() propagates the warning through the
standard summary path. Phase 5 will replace warn-and-label with
explicit conditioning via the .tmb_build_emm_ref_grid() helper.fit_demand_tmb() now fits formula-based random effects with
factor-expanded slopes (e.g. pdDiag(Q0 + alpha ~ condition) or
pdSymm(Q0 + alpha ~ condition)). Single-block pdDiag and pdSymm
with arbitrary RHS terms are accepted. Multi-block pdBlocked /
list() of pdMats lands in Phase 3 (also in this release).src/MixedDemand.h rewritten with a block-aware DATA interface
(Z_q0, Z_alpha, block-structure metadata) and a generalized
per-block Cholesky loop. pdSymm blocks of size > 2 use the
Lewandowski-Kurowicka-Joe Cholesky construction; the d == 2 case
reduces exactly to the previous tanh(rho_bc_raw) parameterization,
so existing intercept-only fits produce bit-identical loglik /
coefficients (verified to 1e-10 on apt)..tmb_build_z_matrices() and
.tmb_build_block_map() in R/tmb-demand.R consume the canonical
block representation and emit the design matrices and metadata the
template needs..tmb_compute_subject_pars() generalized into a per-block
reconstruction. For factor-expanded fits, the Phase 2 implementation
used first-observed-row X and Z for subject-level Q0 / alpha;
Phase 5A in this release supersedes that with NA-on-within-id-
variation plus an opt-in expanded = TRUE argument on
get_subject_pars() for per-(subject, condition) rows..simulate_within_subject_demand() and parity tests
confirm TMB Laplace approximation agrees with NLME's iterative
algorithm to within ~1% on the loglik across all four target specs.fit_demand_tmb() now accepts multi-block nlme::pdBlocked(list(...))
and bare list(pdMat, pdMat, ...) random_effects specifications.
The motivating use case is the load-bearing M1 spec from in-house
manuscript work:
pdBlocked(list(pdSymm(Q0+alpha~1), pdDiag(Q0+alpha~condition-1))),
which combines a correlated subject-baseline block with an
uncorrelated subject-by-condition slopes block. The intercepts-only
alternative inverts the cigarette Q0 direction on that data; the
multi-block spec recovers it.n_blocks > 1; Phase 3 lifts the gate that rejected those shapes.
The renamed gate helper .re_is_phase3_fittable() accepts any number
of single-grouping-level blocks of class pdDiag or pdSymm. Other
pdMat classes (pdCompSymm, pdIdent, pdLogChol, ...) remain
unsupported pending a triggering use case; use fit_demand_mixed()
for those.pdBlocked(list(pdSymm(Q0+alpha~1)))
(single-block-wrapped) produces the same loglik and coefficients as
bare pdSymm(Q0+alpha~1) to optimizer tolerance.tmb-advanced-random-effects.Rmd covers all
random-effects structures beyond intercepts-only:
pdBlocked (Phase 3).get_subject_pars(fit, expanded = TRUE) and attr(re_q0_mat) /
re_alpha_mat access for power users.at argument and the
parameter-first marginalization convention.tmb-mixed-effects.Rmd stays as the intro tier
(intercept-only and basic 2-RE).calc_group_metrics.beezdemand_tmb() gains an at argument for
explicit conditioning on continuous covariates and factor levels.
When at = NULL (default), continuous covariates are evaluated at
their training mean and factors are marginalized across observed
levels (equal weights). The argument shape matches
get_demand_param_emms() / get_demand_comparisons():
at = list(age = 30, gender = "Male").cli::cli_warn() for "covariates held at 0" is retired
entirely. The new training-mean default is statistically defensible,
so the warning would only train users to ignore warnings. The
conditioned_on field in the return list still labels the actual
conditioning point (covariate values used; per-factor treatment, with
"marginal" for the default and the supplied level when at is
given)..tmb_build_emm_ref_grid() helper,
marginalized with equal weights, then Pmax/Omax/Qmax are derived from
the marginalized parameter values. This matches the parameter-level
marginalization convention used by get_demand_param_emms(). (Note:
this differs from "compute metrics per cell, then average"; the two
approaches give different answers for nonlinear transforms.)summary.beezdemand_tmb() now prints a single line under the
Population Demand Metrics block (Metrics conditioned at: <cov>=<X>, <factor>=marginal) that states the conditioning point so a printed
summary is self-describing.get_subject_pars() gains an opt-in expanded argument. When
expanded = TRUE and the fit's random-effects design varies within
id (e.g. M1-style multi-block fits with condition - 1 slopes), the
function returns a long-form table with one row per (subject,
factor-level) combination. Columns include the within-subject factor
names plus model-derived per-cell Q0, alpha, Pmax, and Omax.
Numeric within-id-varying RE-RHS terms are conditioned at the
subject's mean rather than expanded.expanded = FALSE returns the wide one-row-per-subject table
unchanged; predict(), ranef(), and other consumers that depend on
unique IDs continue to work bit-for-bit.Z_q0 / Z_alpha
(random-effects design) columns in addition to X_q0 / X_alpha. Prior
to this fix, M1-style fits where condition appeared only in
random_effects (not in factors) silently returned first-observed-
row Q0/alpha values without warning. Now those fits emit the
subject_pars validation warning and set affected subjects'
Q0/alpha/Pmax/Omax to NA in the default wide table, with a
pointer to expanded = TRUE for per-(subject, condition) values.subject_pars$Q0 / $alpha directly
(plot(fit, type = "individual") and
calculate_amplitude_persistence()) now abort with a targeted message
when those columns are NA, pointing users at
get_subject_pars(fit, expanded = TRUE). Native expanded-shape
support in those consumers is deferred to a follow-up release.get_subject_pars() and calc_group_metrics()
generics gain ... so the new expanded and at arguments dispatch
through UseMethod() correctly. All existing methods updated.formula(), model.matrix(), and update() methods for
beezdemand_tmb, plus formula() and model.matrix() for
beezdemand_hurdle (five new S3 methods total).formula(fit_tmb) returns list(Q0, alpha, random): one-sided
formulas for Q0 and alpha (reconstructed from
fit$formula_details$rhs_q0 / $rhs_alpha, so they reflect any
asymmetric collapse_levels) plus the original random_effects
spec preserved at fit time and round-trippable back to
fit_demand_tmb().model.matrix(fit_tmb) returns a named list of four matrices
(X_q0, X_alpha, Z_q0, Z_alpha); use what = ... to select
one. The list-rather-than-matrix return is intentional and
documented: the TMB tier has two fixed-effect linear predictors
(one per nonlinear parameter) rather than one. X_q0 / X_alpha are
zero-copy references to fit$formula_details; Z_q0 / Z_alpha
are recomputed via the internal .tmb_build_z_matrices() helper.
Degenerate Z requests (e.g., what = "Z_alpha" on a Q0-only fit)
return NULL with an informational message.update(fit_tmb, ...) re-fits with named arguments substituted
into the original call (e.g., update(fit, factors = NULL)).
Honors evaluate = FALSE per the stats::update.default
convention. Does not support formula-update syntax
(. - term); fit_demand_tmb() is argument-driven, not
formula-driven.formula(fit_hurdle) returns list(binary, consumption, random)
with both component formulas intercept-only today. Future support
for factor/covariate effects on hurdle components will enrich
these without changing the API.model.matrix(fit_hurdle) returns intercept-only design matrices
(X_binary, X_consumption) for parity with
model.matrix.beezdemand_tmb().beezdemand_tmb fits now store the original fit_demand_tmb()
call as fit$call so update() can rebuild it. fit_demand_hurdle()
already captured this slot.emmeans / effects / ggeffects
(need recover_data + emm_basis methods), drop1 / add1 /
MuMIn::dredge (need terms() + formula-update update() form),
or any tool that dispatches via coef() after a future
coef() default change. Those each remain follow-up tickets.vcov(), fitted(), and residuals() methods for both
beezdemand_tmb and beezdemand_hurdle (six new S3 methods total).
TMB methods default to the model's native likelihood scale (log scale
for "exponential", natural/LL4 scale for others), matching
broom::augment(fit)$.fitted / $.resid. scale = "natural" opts
into back-transformed values; type = "pearson" divides by the
residual SD on the model scale. Requesting type = "pearson" with
scale = "natural" falls back to type = "response" with an
informational message because a response-scale residual SD is not
identified for the exponential/zben variants without a separate
variance assumption.vcov.beezdemand_hurdle() returns the joint fixed-effect VCOV with
row/column names prefixed by component (zero_probability.<term>,
consumption.<term>, variance.<term>) so downstream tools can
index by component without an ad-hoc string match.fitted.beezdemand_hurdle(marginal = TRUE) (default) returns
marginal expected consumption P(y > 0) * E(y | y > 0) (the
.fitted column of predict(fit, type = "demand"));
marginal = FALSE returns the conditional-on-positive expectation.coef.beezdemand_tmb() gained a forward-compatible type = "internal"
alias that returns the optimizer's flat parameterization (current
default behavior). This preserves the numeric-vector escape hatch
consumed by car::deltaMethod, multcomp::glht, and similar tooling
across a future coef() default change to a per-subject tibble.augment.beezdemand_tmb() was refactored to share an internal
.tmb_fitted_resid() helper with the new fitted() and residuals()
methods, eliminating the duplicate predict() call and guaranteeing
the three accessors cannot drift apart on scale convention. Output
is bit-identical to the previous implementation.parameters::standard_error
(needs insight::get_parameters class-aware dispatch),
DHARMa::simulateResiduals (needs a simulate() method), or
performance::check_residuals (needs class-specific dispatch). The
explicit numeric-vector forms of car::deltaMethod and
multcomp::glht are unlocked when the caller pre-extracts
coef(fit, type = "internal").nobs.beezdemand_hurdle() for universal-accessor parity with
nobs.beezdemand_tmb() and the cross-price classes. broom::glance(fit)$nobs
and BIC(fit) were already correct via their own paths
(param_info$n_obs and the nobs attribute on logLik(), respectively);
this method closes the gap for any caller that consumes nobs() directly.anova.beezdemand_tmb(): joint Wald-chi-square tests on grouped
fixed-effect terms for a single fit, and sequential likelihood-ratio
tests for nested fits (TICKET-013).summary(fit_tmb)$variance_components now reports the Q0 and
alpha random-effect SDs on the log10 scale. TMB estimates these SDs on the
natural-log scale internally (src/MixedDemand.h evaluates
Q0 = exp(log_q0)); summary() previously reported them raw, off by a
factor of log(10) ~= 2.303 from nlme::VarCorr() on a structurally
matched fit_demand_mixed() fit using the default param_space = "log10".
The two backends' random-effect SDs are now directly comparable. The
residual SD (on the model's likelihood scale) and the random-effect
correlations (scale-invariant) are unchanged.summary(fit_tmb)$variance_components$Estimate change in value by a factor
of 1 / log(10). Analysis code that divided TMB RE SDs by ~2.303 to
compare them with nlme::VarCorr() should drop that manual conversion --
it is now applied internally.predict.beezdemand_tmb() gains a level argument for
type = "response". level = "subject" (default) preserves the previous
behavior: it conditions on each subject's random effects, requires the
model's ID column in newdata, and returns a .fitted column. level = "population" evaluates at the fixed-effect coefficients with all random
effects set to zero (the population-mean curve), does not require the ID
column, and returns a predict.fixed column. Passing
c("population", "subject") returns both predict.fixed and
predict.id columns in one call, matching the
nlme::predict.lme(level = 0:1) schema so nlme-based plotting code runs
unchanged.predict.beezdemand_nlme(), which accepts the nlme-style numeric
level (0 / 1), the TMB method takes the character form only; a
numeric level is rejected with a match.arg()-style error.fitted() and residuals() for beezdemand_tmb fits honor the same
level argument: level = "population" now returns population-mean
fitted values and the corresponding residuals. Previously this argument
was an unimplemented stub that returned subject-level values.VarCorr() now has a beezdemand_tmb method. VarCorr(fit_tmb) returns
the random-effect variance components in the matrix layout produced by
nlme::VarCorr() (a "VarCorr.lme"-class object with Variance,
StdDev, and (for pdSymm fits) Corr columns plus a final Residual
row), so users coming from nlme or lme4 can introspect a TMB fit with
a familiar accessor. The values match
summary(fit_tmb)$variance_components: the Q0/alpha random-effect SDs on
the log10 scale and the residual SD on the model's likelihood scale.check_demand_model() on a beezdemand_tmb fit now reports random-effect
scale on the log10 scale, consistent with
summary(fit_tmb)$variance_components (the TICKET-015 convention).
Previously these were raw natural-log-scale SDs, a factor of log(10)
larger. The raw internal SDs (still used for the near-zero degeneracy
check) are now exposed separately as $random_effects$sd_internal_log.
(Superseded later in this release: those log10-scale SDs are now
$random_effects$sd_log10, and $random_effects$variances holds their
squares -- see "Inference gates and diagnostic reporting" above.)The tidy() and glance() introspection methods now expose the same
column names, default arguments, and component labels on the
beezdemand_nlme and beezdemand_tmb backends, so backend-agnostic code
needs no dispatch glue.
glance(fit_tmb)$equation is renamed to
equation_form, matching glance(fit_nlme). There is no aliased
equation column. The fit_demand_tmb() API is new in 0.3.0 (this
rename happened during its development), so no released code depends on
the old name.tidy(fit_tmb) labels fixed-effect rows
component == "fixed" instead of "consumption", matching
tidy(fit_nlme) and the nlme / lme4 convention. Code filtering TMB
tidy() output on component == "consumption" will return zero rows.
Hurdle methods are unchanged. (summary(fit_tmb)$coefficients was
harmonized to "fixed" separately in the TICKET-031 follow-up below.)tidy(fit_tmb, effects = "ran_pars") reports the
random-effect variance components on the same scale as
summary(fit_tmb)$variance_components (Q0/alpha RE SDs on the log10
scale, residual SD on the likelihood scale) rather than the raw
internal logsigma optimizer coefficients. std.error is NA for
these rows, as it is for tidy(fit_nlme).tidy.beezdemand_tmb() gains an effects argument
(c("fixed", "ran_pars"), both by default) matching
tidy.beezdemand_nlme(): effects = "fixed" returns the fixed-effect
rows, effects = "ran_pars" returns the variance-component rows. An
invalid value is rejected with a match.arg()-style error.glance(fit_nlme) gains an n_random_effects column (the count of
random-effect terms), so the canonical glance() columns (
model_class, backend, equation_form, nobs, n_subjects,
n_random_effects, converged, logLik, AIC, BIC) are now
identical across both backends.tidy(fit_nlme, effects = "ran_pars")$estimate now reports random-effect
standard deviations (pulled from nlme::VarCorr(model)[, "StdDev"]),
not variances. This matches the tidy(fit_tmb) sibling (post-TICKET-015)
and the broom.mixed::tidy.lme upstream convention, closing the
cross-backend divergence on "ran_pars" rows. Migration: callers that
consumed the previous value as a variance should square the estimate
(estimate^2) or read nlme::VarCorr(fit$model)[, "Variance"] directly.
Hurdle and fixed tiers are unaffected.summary(fit_tmb)$coefficients$component now also emits "fixed" for
q0 / alpha / log_k rows, matching tidy(fit_tmb) (renamed in TICKET-017
above) and summary(fit_nlme)$coefficients. Code filtering
summary(fit_tmb)$coefficients on component == "consumption" will
return zero rows. summary(fit_tmb)$derived_metrics$component is
deliberately left as "consumption" because those rows describe derived
demand metrics (pmax, omax, q_at_pmax, elasticity_at_pmax), not fitted
coefficients, and a future ticket may rename them to "derived" or
"metric". Hurdle methods are unchanged.These sections capture the original 0.3.0 release scope (TMB mixed-effects modeling tier, hurdle improvements, bug fixes, quality / tooling). The TICKET-011 phases above were added under the same 0.3.0 development cycle.
fit_demand_tmb() is the new modern mixed-effects path for behavioral
economic demand models, alongside the existing NLME tier
(fit_demand_mixed()). Backed by Template Model Builder (TMB) with
automatic differentiation and Laplace approximation. Supports four
equations (exponential, exponentiated, simplified, zben),
1-RE (Q0 only) or 2-RE (correlated Q0 + alpha) random-effect structures,
estimated or fixed k, factor and continuous covariates with asymmetric
collapse_levels, and three data-adaptive starting-value strategies via
multi_start = TRUE (default).
fit_demand_tmb() exposes full optimizer controls via tmb_control:
optimizer: "nlminb" (default) or "L-BFGS-B" for recovering from
convergence failures (code 1 or 8).rel_tol: Convergence tolerance for nlminb (default 1e-10).lower / upper: Named numeric vectors for parameter bounds on the
optimizer scale; applied to all occurrences of repeated parameter names
(e.g., beta_q0).warm_start: Restart optimization from a previous fit$opt$par vector.
Automatically disables multi_start.trace: Optimizer trace output (default 0).The beezdemand_tmb class has a full S3 method suite:
print, summary, coef, fixef, ranef, logLik, AIC, BIC,
nobs, predict (response, parameters, demand), confint,
residuals, fitted, vcov, broom (tidy, glance, augment),
get_subject_pars, calc_group_metrics, get_demand_param_emms,
get_demand_comparisons, and visualization (plot, plot_qq,
plot_loss_surface, plot_loss_profile, plot_re_diagnostics,
plot_alpha_distribution, plot_elasticity, plot_expenditure,
plot_demand_overlay, and the cross-model forest plot
plot_model_comparison(), which compares parameter estimates and CIs
across any fitted demand models via their tidy() methods).
vignettes/tmb-mixed-effects.Rmd walks through the full TMB workflow
(equations, random-effect structures, diagnostics) with cache-aware
fast/full mode.
vignettes/convergence-guide.Rmd documents convergence troubleshooting
for both the TMB and NLME tiers.
predict.beezdemand_hurdle(type = "probability", marginal = TRUE)
computes population-averaged P(zero) by integrating over the random
intercept distribution. Methods: "kde" (default), "normal",
"empirical".
plot.beezdemand_hurdle(type = "probability") now shows the marginal
(population-averaged) P(zero) curve by default. marginal = FALSE
reverts to the old conditional (RE = 0) behavior.
calc_group_metrics.beezdemand_hurdle() now returns both conditional
and unconditional Pmax / Omax (TICKET-003). Conditional ($Pmax,
$Omax) keeps its long-standing Part-II-only meaning. New
$Pmax_unconditional / $Omax_unconditional come from optimizing
p * (1 - P0(p)) * Q(p) over the observed price domain. Subject-level
subject_pars likewise gains Pmax_unconditional /
Omax_unconditional columns. summary()'s derived_metrics reports
both sets, with component = "unconditional" on the new rows.
plot_expenditure.beezdemand_hurdle() gains
demand_type = c("unconditional", "conditional") (default
unconditional). The displayed expenditure curve and the Pmax/Omax
reference lines now come from the same metric set, so they always
align, which fixes the visible misalignment where the curve used
(1 - P0) * Q but the reference lines used the Part-II-only metrics.
New calculate_amplitude_persistence.beezdemand_tmb() method
(TICKET-004) lets users compute amplitude/persistence factors directly
from fit_demand_tmb() results (default persistence components:
c("Pmax", "Omax", "alpha"), because TMB subject_pars does not include
breakpoint).
New cross-price S3 methods (TICKET-005): print, augment, confint
(lm + lmer), and nobs for cp_model_nls, cp_model_lm, and
cp_model_lmer. Augment for lmer additionally includes a .fixed
column (population-level prediction with random effects = 0). All
methods return safe fallbacks (e.g., NA or an empty tibble) instead of
erroring when the underlying model is NULL.
Linearized marginal NLL surface for plot_loss_surface() on NLME
models (commit 2bedd29).
Visualization helpers added to TMB and NLME vignettes (commit
75a202a).
get_demand_param_emms.beezdemand_nlme() gains a param argument
("both", "Q0", "alpha") for API parity with the
beezdemand_tmb method (TICKET-012). Default "both" preserves the
historical return shape; "Q0" and "alpha" narrow the output to a
single parameter's columns for easier pivoting and plotting.
fit_demand_tmb() now drops rows with NA values in any modeling
column (id, price, response, factors, continuous covariates) before
entering the TMB pipeline, matching the fit_demand_mixed() behavior.
Previously a single NA in the response crashed .tmb_prepare_data()
with "missing value where TRUE/FALSE needed", and NA in factors or
covariates could propagate into model.matrix() and trigger a TMB
segfault during MakeADFun().
predict.beezdemand_tmb() now rebuilds the fixed-effect linear
predictor from newdata instead of reusing training-time
subject_pars$Q0 / alpha. Predictions for any model with factors or
continuous covariates now correctly reflect the values supplied in
newdata; an unknown subject id at level = "subject" is an error
(use level = "population" for the random-effects-at-zero prediction).
Previously the function silently used cached subject parameters for
known subjects and the reference-level intercepts for unknown ones,
producing systematically biased .fitted values.
augment.beezdemand_tmb() inherits the fix. Predict now also errors
clearly when newdata is missing a required modeling column or
contains factor levels not seen in training.
get_demand_param_emms.beezdemand_tmb() and
get_demand_comparisons.beezdemand_tmb() now include continuous
covariates in the reference grid, matching the dimensionality of the
fitted beta coefficients. Covariates default to their training-data
mean (matching emmeans::ref_grid) and can be overridden via
at = list(covname = value). Previously, TMB fits that mixed factors
and continuous covariates produced non-conformable arguments in the
Wald variance calculation or silently used the wrong model basis.
get_demand_comparisons() now restricts pairwise contrasts to observed
factor combinations. Previously, with unbalanced designs (e.g.,
different dose levels per drug), the function computed contrasts on the
full factorial grid, producing phantom comparisons for non-existent
factor combinations and identical estimates across contrast_by groups
in additive models.
summary.beezdemand_nlme() / tidy() keep the Wald statistic/p.value on
the estimation scale when report_space != internal_space, so natural-scale
inference is nlme's native containment-t (DF-aware) test rather than a
recomputed (and degenerate) natural-scale Wald test. Only estimate/std.error
are back-transformed (broom convention; supersedes the earlier TICKET-006
delta-method recompute).
summary.beezdemand_hurdle()$coefficients_matrix is now labelled
"z value" (was "t value") to match the pnorm-based p-value
computation. TMB-based hurdle models use Laplace approximation, so the
z-test is the correct asymptotic inference (matches the glmmTMB
convention); only the label was wrong (TICKET-006).
fit_demand_tmb() and fit_demand_hurdle() now expose hessian_pd on
the fit object and warn at fit time when the Hessian is not positive
definite (pdHess = FALSE). summary() adds a corresponding note;
tidy() sets a hessian_warning attribute on its output. Previously
these models silently reported unreliable standard errors / p-values /
Wald intervals when the Hessian was singular (TICKET-008). Backwards
compatible: legacy fit objects without the field still work.
check_demand_model() no longer crashes when called on
beezdemand_tmb objects (TICKET-002; the underlying name + residual
fields were already corrected in commit 719c0ed; this release adds
the regression coverage that pins the fix).
fit_demand_mixed() correctly handles NA values in the input data
and strips spurious names from start values (commit b39c24a).
Visualization improvements (smoothing, value clamping, and APA styling
for diagnostic plots; commit bf125f6).
Package audit fixes (boundary detection, data
validation, heuristic improvements; commit 60b13a2).
Bare stop() / warning() / message() calls in non-legacy R files
replaced with their cli equivalents (cli::cli_abort, cli::cli_warn,
cli::cli_inform) and the package's structured error helpers
(validation_error, fitting_error, missing_package_error;
TICKET-009). Errors now carry stable class tags (e.g.,
"beezdemand_validation_error") for programmatic catch-handling. Bare
calls in legacy R/analyze.R are intentionally left in place; that
file is slated for removal in v1.0. The error helpers themselves were
refactored to use cli::cli_abort internally so callers can pass cli
inline markup ({.arg}, {.field}, {.val}, {.fn}).
Test coverage added for four previously-untested exported functions
(TICKET-007): cp_posthoc_slopes(), cp_posthoc_intercepts(),
extract_coefficients(), and get_demand_param_trends(). 23 new
tests across three new test files using the etm and ko example
datasets.
New regression tests for check_demand_model.beezdemand_tmb() (1-RE
and 2-RE), calculate_amplitude_persistence.beezdemand_tmb(),
summary.beezdemand_nlme() p-value preservation,
summary.beezdemand_hurdle() coefficient-matrix labelling,
hessian_pd propagation on TMB and hurdle fits, the cross-price S3
methods, hurdle unconditional Pmax/Omax, and the structured error
class hierarchy.
FitCurves() is now superseded by fit_demand_fixed(). FitCurves() will
continue to work but emits a soft deprecation warning. The new function
provides a modern S3 interface with summary(), tidy(), glance(),
predict(), and plot() methods. See vignette("migration-guide") for
migration instructions.
FitMeanCurves() is now superseded by fit_demand_fixed(agg = "Mean") or
fit_demand_fixed(agg = "Pooled").
fit_demand_mixed() now supports the Koffarnus et al. (2015) exponentiated
equation via equation_form = "koff". This enables fitting demand curves
using the same equation form available in FitCurves() within the modern
hierarchical mixed-effects framework. The k parameter can be user-specified
or auto-calculated from data range.augment() methods for all model classes provide fitted values and
residuals in a tidy tibble:
augment.beezdemand_fixed(): Returns .fitted, .residaugment.beezdemand_hurdle(): Returns .fitted, .fitted_link,
.fitted_prob, .resid, .resid_responseaugment.beezdemand_nlme(): Returns .fitted, .resid, .fixedNew compare_models() function for unified model comparison across all
beezdemand model classes. Reports AIC, BIC, delta_AIC, delta_BIC, and
performs likelihood ratio tests when models are from the same backend
and nested.
New anova() S3 methods for comparing nested models:
anova.beezdemand_hurdle(): LRT for nested hurdle modelsanova.beezdemand_nlme(): Delegates to nlme::anova.lme()New check_demand_model() generic with methods for all model classes. Performs
comprehensive diagnostics including convergence checks, boundary condition
detection, random effect variance assessment, and residual outlier detection.
Returns structured diagnostics object with issues and recommendations. (Named
check_demand_model() to avoid conflict with performance::check_model().)
New plot_residuals() function creates diagnostic plots: residuals vs fitted,
histogram of residuals, and Q-Q plots. Works with all model classes via
the augment() infrastructure.
New plot_qq() function creates Q-Q plots for random effects to assess
normality assumptions in hurdle and NLME models.
alpha_star (normalized alpha, Strategy B;
Rzeszutek et al., 2025), which makes the elasticity parameter comparable
across different values of k. Available in FitCurves() output (columns
alpha_star and alpha_star_se), tidy() on beezdemand_fixed objects,
and tidy() on beezdemand_hurdle objects. Standard errors are obtained
via the delta method. See ?param-registry for details.New get_empirical_measures() as a modern replacement for GetEmpirical().
Returns a beezdemand_empirical S3 object; access the results via
$measures.
New get_descriptive_summary() as a modern replacement for
GetDescriptives(). Returns a beezdemand_descriptive S3 object; access
the results via $statistics.
New get_k() as a modern replacement for GetK(). Returns a single numeric
k value with optional verbose output.
New confint() methods for extracting confidence intervals from all model
classes: beezdemand_fixed, beezdemand_hurdle, beezdemand_nlme, and
cp_model_nls.
New migration guide vignette (vignette("migration-guide")) documenting the
transition from FitCurves() to fit_demand_fixed().
summary() methods for beezdemand_hurdle and beezdemand_nlme now return
structured summary objects instead of printing directly. Use
print(summary(fit)) for console output. Programmatic access is now
possible: s <- summary(fit); s$coefficients.
fit_demand_hurdle() now fits demand parameters in natural-log space
(log_q0, log_alpha, log_k) and reports back-transformed values; the
param_space argument has been removed.
fit_cp_nls() now uses log10-parameterized optimizer coefficients
(log10_qalone, I, log10_beta) across equation forms; the "exponential"
form fits on the log10(y) response scale and filters y <= 0 with a warning.
predict.cp_model_nls() now always returns y_pred on the natural y scale;
for "exponential" it additionally returns y_pred_log10 (and no longer returns
y_pred_natural).
New fit_demand_fixed() function provides a modern interface for individual
demand curve fitting. Returns a structured S3 object with summary(),
tidy(), and glance() methods. This wrapper offers the same functionality
as FitCurves() but with a standardized API.
New systematicity wrappers with unified output vocabulary:
check_systematic_demand() for purchase task data (wraps CheckUnsystematic())check_systematic_cp() for cross-price data (wraps check_unsystematic_cp())Both return beezdemand_systematicity objects with identical column
schemas (differing only in NA values for domain-specific fields).
First-class tidy() and glance() support is now guaranteed across all
beezdemand model classes. All methods return tibbles with standardized
columns including model_class and backend.
All summary objects now inherit from beezdemand_summary base class,
enabling shared fallback behavior and consistent field availability.
This release introduces Stability Contracts for all model classes:
summary() objects now return structured S3 objects with class
c("summary.<class>", "beezdemand_summary"). Required fields include:
call, model_class, backend, nobs, n_subjects, converged,
logLik, AIC, BIC, coefficients (tibble), notes.
tidy() methods return tibbles with columns: term, estimate,
std.error, statistic, p.value. Multi-part models include a
component column (e.g., "fixed", "variance", "derived").
glance() methods return 1-row tibbles with columns: model_class,
backend, nobs, n_subjects, converged, logLik, AIC, BIC.
fit_cp_nls() and fit_cp_linear() now accept x_var/y_var to map
non-standard column names to canonical ones ("x", "y"). fit_cp_linear()
additionally accepts id_var, group_var, and target_var. Default behavior
is unchanged when these arguments are omitted.
fit_cp_linear() gains explicit filter_target and target_level top-level
arguments (previously these were handled implicitly via validate_cp_data()).
Existing calls without these arguments are unaffected.
fit_cp_nls(start_vals=) is deprecated in favor of start_values=. The old
argument still works but emits a deprecation warning.
The following deprecations will take effect in version 0.2.0:
beezdemand::pull() is deprecated in favor of dplyr::pull(). The beezdemand
version was a legacy helper that predates the dplyr function.
The inverse_fun argument in summary.cp_model_nls(), plot.cp_model_nls(),
and predict.cp_model_nls() is deprecated in favor of inv_fun for consistency
with mixed-effects model methods.
Standardized argument names across cross-price model methods (inv_fun instead
of inverse_fun)
Cross-price plot methods now have consistent argument ordering across
plot.cp_model_nls(), plot.cp_model_lm(), and plot.cp_model_lmer()
Key user-facing functions now return tibbles for better compatibility with
tidyverse workflows: predict.cp_model_nls(), predict.cp_model_lm(),
tidy.cp_model_nls(), glance.cp_model_nls()
Added standardized error helpers (validation_error(), fitting_error(),
missing_package_error()) for consistent error messaging
check_unsystematic_cp() now returns an object of class cp_unsystematic
with proper summary() method dispatch (no longer overrides summary.tbl_df())
Added comprehensive hurdle model functionality using TMB (Template Model Builder):
fit_demand_hurdle(): Fit two-part hurdle models with 2 or 3 random effects
Part I models probability of zero consumption (logistic regression with random intercept)
Part II models log-consumption given positive response (nonlinear mixed effects)
S3 methods for beezdemand_hurdle objects:
print(), summary(), coef(), logLik(), AIC(), BIC()
predict(): Extract subject parameters or predict demand/probability
plot(): Visualize demand curves, zero probability, parameter distributions
Utility functions:
calc_omax_pmax(): Calculate Pmax and Omax from demand parameters
get_subject_pars(): Extract subject-specific parameter estimates
compare_hurdle_models(): Likelihood ratio test for model comparison
get_hurdle_param_summary(): Summary statistics for individual parameters
Simulation functions:
simulate_hurdle_data(): Generate synthetic hurdle model data
run_hurdle_monte_carlo(): Monte Carlo simulation studies
New vignette "Hurdle Demand Models" with comprehensive examples
New dataset apt_full: Full alcohol purchase task data with 1,100 subjects and demographic covariates
Added comprehensive cross-price demand model functionality:
check_unsystematic_cp(): Check for unsystematic data in cross-price models
fit_cp_nls(): Nonlinear model fitting for cross-price demand data
fit_cp_linear(): Linear model fitting for cross-price demand data with options for fixed effects and mixed-effects models
New utility functions for cross-price model objects:
summary(), plot(), glance(), and tidy() methodsextract_coefficients(): Extract model coefficients in tidy format
cp_posthoc_slopes() and cp_posthoc_intercepts(): Post-hoc comparisons for model parameters
validate_cp_data(): Validate and filter cross-price demand data
Added new vignette "How to Use Cross-Price Demand Model Functions" demonstrating:
Required data structure for cross-price analyses
Checking for unsystematic data
Both two-stage and pooled model fitting approaches
Linear and mixed-effects modeling options
Model visualization and coefficient extraction
Post-hoc comparisons
No longer relies on nlmrt and instead relies on nlsr
Fixes an issue where CheckUnsystematic may not flag certain cases when data are passed as tibble
Fixes deprecated arguments in ggplot2
Add ability to specify a start value for alpha in ExtraF() function
FitCurves(). These arguments are constrainq0, startq0, and startalpha. These arguments allow Q0 to be constrained so alpha is the only fitted parameter and allow for user-specified starting values.One major change that might affect previous scripts is that in output summary tables, the column formally named ID is now named id (lowercase)
Cleaned up a few things here and there. The package is close for submission to CRAN as it passes R CMD check with no errors, warnings, or notes
GetSharedK() updated to work better and faster at finding a reasonable value
Internal helper functions added to optimize GetSharedK()
ExtraF() now compares alpha and Q0
A number of functions now allows you to specify the column names
FitCurves() correctly pulls alpha and q0 standard errors when k is fitted as a free parameter. Also no longer accepts data transformations. Must be done prior to fitting using ChangeData().
FitCurves() now fits mean/pooled data based on method argument.
GetSharedK() no longer accepts data transformations.
New ChangeData() will soon serve as the replacement to
ReplaceZeros() and other arguments specified in FitCurves().
For the time being, FitCurves() will actually output a list
object. This may cause failures with old scripts. Try modifying
scripts to take the first element out of the list. This will soon be
taken care of.
New FitMeanCurves() will fit curve to averaged or pooled
data. Can also make plots.
FitCurves() can now make plots.
GetDescriptives() can make box and whisker plots.