---
title: "TMB Mixed-Effects Demand Models"
author: "Brent Kaplan"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{TMB Mixed-Effects Demand Models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  dev = "png",
  dpi = 144,
  fig.width = 7,
  fig.height = 5,
  warning = FALSE,
  message = FALSE,
  cache = TRUE,
  cache.path = "tmb-mixed-effects_cache/",
  autodep = TRUE
)

library(beezdemand)
library(dplyr)
library(ggplot2)

data("apt_full", package = "beezdemand", envir = environment())

# Set `BEEZDEMAND_VIGNETTE_MODE=full` to run the slower examples.
render_fast <- tolower(Sys.getenv("BEEZDEMAND_VIGNETTE_MODE", "fast")) != "full"

cache_key_object <- function(x) {
  tmp <- tempfile(fileext = ".rds")
  saveRDS(x, tmp)
  on.exit(unlink(tmp), add = TRUE)
  unname(tools::md5sum(tmp))
}

knitr::opts_chunk$set(
  cache.extra = list(
    beezdemand_version = as.character(utils::packageVersion("beezdemand")),
    apt_full = cache_key_object(apt_full)
  )
)
```

```{r data-setup, include = FALSE}
# Consistent subset used throughout the vignette
set.seed(42)
ids_sub <- sample(unique(apt_full$id), 100)
dat <- apt_full[apt_full$id %in% ids_sub, ]
```

## Introduction

`fit_demand_tmb()` fits continuous mixed-effects demand models using
[Template Model Builder](https://github.com/kaskr/adcomp) (TMB). It is the
modern alternative to `fit_demand_mixed()` (which uses `nlme`) and provides
several advantages:

- **Automatic differentiation** -- exact gradients via compiled C++, replacing
  the numerical finite-difference approximations used by `nlme`
- **Laplace approximation** -- integrates over random effects analytically
  rather than relying on iterative linearization
- **Multi-start optimization** -- automatically tries multiple starting-value
  sets and keeps the best fit
- **Four equation forms** -- exponential, exponentiated, simplified, and
  zero-bounded exponential (zben)
- **Factor and covariate support** -- design matrices for group comparisons
  with estimated marginal means (EMMs) and pairwise contrasts

This vignette covers `fit_demand_tmb()` for continuous consumption data. For
**two-part hurdle models** that explicitly model zero consumption, see
`vignette("hurdle-demand-models")`. For **individual NLS curves**, see
`vignette("fixed-demand")`.

## Quick Start

```{r quick-fit}
fit <- fit_demand_tmb(
  dat,
  y_var = "y", x_var = "x", id_var = "id",
  equation = "exponential",
  random_effects = c("q0", "alpha"),
  verbose = 0
)
fit
```

```{r quick-summary}
summary(fit)
```

```{r quick-plot, fig.cap = "Population demand curve from the exponential equation."}
plot(fit, type = "demand")
```

The `exponential` equation (Hursh & Silberberg, 2008) models log-transformed
consumption and is the most reliable choice for 2-random-effect models. It
automatically drops zero-consumption observations.

## Choosing an Equation

`fit_demand_tmb()` supports four demand equations:

| Equation | Response | Zeros | k | Best for |
|----------|----------|-------|---|----------|
| `"exponential"` | log(Q) | Dropped | Estimated or fixed | Most datasets; robust 2-RE convergence |
| `"exponentiated"` | Raw Q | Allowed | Estimated or fixed | Data with few zeros; 1-RE models |
| `"simplified"` | Raw Q | Allowed | None | Simpler model without k |
| `"zben"` | LL4(Q) | Allowed (via transform) | None | Wide dynamic range with LL4 compression |

### Mathematical Specifications

**Exponential** (Hursh & Silberberg, 2008):
$$\log(Q_{ij}) = \log(Q_{0i}) + k \left(e^{-\alpha_i \cdot Q_{0i} \cdot C_j} - 1\right) + \varepsilon_{ij}$$

**Exponentiated** (Koffarnus et al., 2015):
$$Q_{ij} = Q_{0i} \cdot 10^{k \left(e^{-\alpha_i \cdot Q_{0i} \cdot C_j} - 1\right)} + \varepsilon_{ij}$$

**Simplified**:
$$Q_{ij} = Q_{0i} \cdot e^{-\alpha_i \cdot Q_{0i} \cdot C_j} + \varepsilon_{ij}$$

**Zero-bounded exponential (zben)** (see `src/MixedDemand.h` for full
specification):
$$\log_{10}(Q_{0i}) \cdot e^{-\frac{\alpha_i}{\log_{10}(Q_{0i})} \cdot Q_{0i} \cdot C_j} + \varepsilon_{ij}$$

where $Q_{0i} = \exp(\mathbf{x}_i^\top \boldsymbol{\beta}_{Q_0} + b_i)$,
$\alpha_i = \exp(\mathbf{x}_i^\top \boldsymbol{\beta}_\alpha + c_i)$, and
$(b_i, c_i)$ are correlated random effects.

```{r eq-exponential}
fit_exp <- fit_demand_tmb(dat, equation = "exponential",
                          random_effects = c("q0", "alpha"), verbose = 0)
```

```{r eq-exponentiated}
# exponentiated converges reliably with 1-RE
fit_expon <- fit_demand_tmb(dat, equation = "exponentiated",
                            random_effects = "q0", verbose = 0)
```

```{r eq-simplified}
fit_simp <- fit_demand_tmb(dat, equation = "simplified",
                           random_effects = "q0", verbose = 0)
```

```{r eq-zben}
dat$y_ll4 <- ll4(dat$y)
fit_zben <- fit_demand_tmb(dat, y_var = "y_ll4", equation = "zben",
                           random_effects = "q0", verbose = 0)
```

```{r eq-compare}
data.frame(
  equation = c("exponential", "exponentiated", "simplified", "zben"),
  random_effects = c("2-RE", "1-RE", "1-RE", "1-RE"),
  converged = c(fit_exp$converged, fit_expon$converged,
                fit_simp$converged, fit_zben$converged),
  AIC = round(c(AIC(fit_exp), AIC(fit_expon),
                AIC(fit_simp), AIC(fit_zben)), 1)
)
```

**Note:** AIC values are not directly comparable across equations because they
model different response scales (log Q vs raw Q vs LL4(Q)).

**Convergence tip:** The `exponential` equation is the most reliable for
2-random-effect models. The `exponentiated` equation works well with 1-RE but
can struggle to converge with 2-RE, especially with smaller samples.

## Random Effects Structure

`fit_demand_tmb()` supports two configurations:

- `random_effects = "q0"` -- **1-RE**: random intercept on $Q_0$ only;
  $\alpha$ is constant across subjects
- `random_effects = c("q0", "alpha")` -- **2-RE**: random effects on both
  $Q_0$ and $\alpha$ with an estimated correlation

```{r re-compare}
fit_1re <- fit_demand_tmb(dat, equation = "exponential",
                          random_effects = "q0", verbose = 0)
fit_2re <- fit  # reuse from Quick Start (exponential, 2-RE)
```

```{r re-compare-output}
compare_models(fit_1re, fit_2re)
```

A significant LRT p-value and substantially lower AIC/BIC for the 2-RE model
indicate that subjects differ meaningfully in both intensity ($Q_0$) and
elasticity ($\alpha$).

```{r re-ranef}
head(nlme::ranef(fit_2re))
```

The `b_i` column is the random deviation on log($Q_0$) and `c_i` is the
random deviation on log($\alpha$) for each subject.

## Estimating vs. Fixing k

For the `exponential` and `exponentiated` equations, $k$ scales the range of
the demand curve. By default it is estimated (`estimate_k = TRUE`). You can
fix it at the conventional value of 2 (Hursh & Silberberg, 2008):

```{r k-compare}
fit_k_free <- fit_demand_tmb(dat, equation = "exponential",
                              random_effects = c("q0", "alpha"),
                              estimate_k = TRUE, verbose = 0)
fit_k_fixed <- fit_demand_tmb(dat, equation = "exponential",
                               random_effects = c("q0", "alpha"),
                               estimate_k = FALSE, k = 2, verbose = 0)

data.frame(
  k = c("estimated", "fixed at 2"),
  converged = c(fit_k_free$converged, fit_k_fixed$converged),
  AIC = round(c(AIC(fit_k_free), AIC(fit_k_fixed)), 1),
  k_value = round(c(exp(coef(fit_k_free)[["log_k"]]), 2), 3)
)
```

The `simplified` and `zben` equations do not use a $k$ parameter.

## Examining Results

### Summary and Coefficients

```{r results-summary}
summary(fit_2re)
```

The fixed-effect coefficient table reports `estimate` and `std.error` on the
scale set by `report_space`: the default, `"natural"`, back-transforms $Q_0$,
$\alpha$, and $k$ off the log estimation scale. The `statistic` and `p.value`
columns are always computed on the estimation (log) scale, where the Wald test
is well defined; they are not recomputed from the back-transformed `estimate`
and `std.error` (that recompute is degenerate for a strictly positive
parameter -- it would test "ratio = 0" rather than "ratio = 1"). One
consequence: `estimate / std.error` from the default table does not reproduce
the reported `statistic`. This matches the `broom` and
`emmeans` convention; pass `report_space = "internal"` (or `"log10"`) to read
the estimate and SE on the same scale as the test.

The `variance_components` block reports the $Q_0$ and $\alpha$ random-effect
standard deviations on the **log10 scale**. TMB estimates these SDs on the
natural-log scale internally; `summary()` divides them by $\log(10)$ so they
are directly comparable with `nlme::VarCorr()` on a `fit_demand_mixed()` fit
using the default `param_space = "log10"`. The residual SD is reported on the
model's likelihood scale, which is equation-dependent (log-consumption for
`exponential`, the LL4-transformed `y_var` for `zben`); the random-effect
correlations are scale-invariant.

For users coming from `nlme` or `lme4`, the `VarCorr()` accessor returns these
same variance components in the familiar `nlme::VarCorr()` matrix layout -- a
`Variance` / `StdDev` matrix (with a `Corr` column for `pdSymm` fits) and a
final `Residual` row:

```{r results-varcorr}
VarCorr(fit_2re)
```

Note that `tidy()` (shown next) reports the raw internal optimizer parameters
instead -- the `logsigma` rows are the natural log of each RE SD, not the
log10-scale SDs from `summary()$variance_components`.

```{r results-tidy}
tidy(fit_2re)
```

```{r results-glance}
glance(fit_2re)
```

### Subject-Specific Parameters

Each subject gets empirical Bayes estimates of $Q_0$, $\alpha$, $P_{max}$, and
$O_{max}$:

```{r results-subject-pars}
spars <- get_subject_pars(fit_2re)
head(spars)
```

```{r results-subject-summary}
spars |>
  summarise(
    across(c(Q0, alpha, Pmax, Omax),
           list(median = median, mean = mean, sd = sd),
           .names = "{.col}_{.fn}")
  ) |>
  tidyr::pivot_longer(everything(), names_to = c("parameter", "stat"),
                      names_sep = "_") |>
  tidyr::pivot_wider(names_from = stat, values_from = value) |>
  mutate(across(where(is.numeric), \(x) round(x, 4)))
```

### Amplitude–Persistence Decomposition

`calculate_amplitude_persistence()` collapses the per-subject parameters into
two latent factors used in behavioral economic indices: **Amplitude**
(intensity of demand, primarily $Q_0$) and **Persistence** (sensitivity to
price, drawn from $P_{max}$, $O_{max}$, and $1/\alpha$). The TMB method
extracts subject parameters from `fit$subject_pars` and delegates to the
default Z-score implementation, so the result is comparable across model
tiers.

```{r results-amplitude-persistence}
ap <- calculate_amplitude_persistence(fit_2re)
head(ap)
```

```{r results-amplitude-persistence-summary}
ap |>
  summarise(
    Amplitude_mean = mean(Amplitude, na.rm = TRUE),
    Persistence_mean = mean(Persistence, na.rm = TRUE),
    Amplitude_sd = sd(Amplitude, na.rm = TRUE),
    Persistence_sd = sd(Persistence, na.rm = TRUE)
  )
```

By construction, both factors are sample-standardized (mean 0, SD 1 within the
fitted dataset). For cross-sample comparisons, supply external `basis_means`
and `basis_sds` so the standardization uses a fixed reference.

### Confidence Intervals

```{r results-confint}
confint(fit_2re)
```

By default, confidence intervals are on the internal (log) scale. Use
`report_space = "natural"` for back-transformed intervals:

```{r results-confint-natural}
confint(fit_2re, report_space = "natural")
```

For a quick diagnostic on the Gaussian (Wald) approximation, request Monte
Carlo intervals with `method = "simulate"`. These draw `R` samples from the
joint asymptotic posterior `N(coef, vcov)` and report empirical quantiles.
They are asymptotically Wald-equivalent, so a large discrepancy between the
two flags a fit where the Gaussian approximation is suspect. This is a
diagnostic comparison, **not** an accuracy improvement: the `simulate`
method does not capture non-Gaussian posterior shape and carries no
positivity guarantee on the internal scale. Set `seed` for reproducibility.

```{r results-confint-simulate}
ci_wald <- confint(fit_2re)
ci_sim <- confint(fit_2re, method = "simulate", R = 2000, seed = 42)
data.frame(
  term = ci_wald$term,
  width_wald = ci_wald$conf.high - ci_wald$conf.low,
  width_sim = ci_sim$conf.high - ci_sim$conf.low
)
```

### Variance-Covariance and the Delta Method

`vcov()` returns the fixed-effect variance-covariance matrix from the TMB
sdreport — the inverse of the negative Hessian at the MLE, restricted to
fixed effects after Laplace-marginalizing the random effects. Combined
with the optimizer's internal parameter vector (`coef(fit, type = "internal")`),
it lets you apply the delta method to any nonlinear function of the
parameters via `car::deltaMethod`. Pass the parameter vector explicitly
so the call is stable across the planned `coef()` default change in a
future release.

```{r results-vcov-delta, eval = requireNamespace("car", quietly = TRUE)}
beta <- coef(fit_2re, type = "internal")
car::deltaMethod(
  beta,
  paste0(names(beta)[1], " * 1"),
  vcov. = vcov(fit_2re)
)
```

`fitted(fit)` and `residuals(fit)` are also exposed as direct accessors,
matching the model-scale convention used by `broom::augment()` (log scale
for `"exponential"`, natural/LL4 scale for others). Pass
`scale = "natural"` to back-transform.

### Predictions

Three prediction types are available:

```{r predict-response}
# Fitted values for observed data
pred_resp <- predict(fit_2re, type = "response")
head(pred_resp)
```

```{r predict-demand}
# Population demand curve at specific prices
predict(fit_2re, type = "demand", prices = c(0, 0.5, 1, 2, 5, 10, 20))
```

```{r predict-parameters}
# Subject-level parameter estimates (same as get_subject_pars)
pred_pars <- predict(fit_2re, type = "parameters")
head(pred_pars)
```

For `type = "response"`, the `level` argument controls whether predictions
condition on subject random effects. `level = "subject"` (the default)
conditions on each subject's random effects and needs an `id` column in
`newdata`; `level = "population"` sets the random effects to zero and needs
no `id`. Requesting both at once returns `predict.fixed` (population) and
`predict.id` (subject) side by side, matching the
`nlme::predict.lme(level = 0:1)` layout so `nlme`-based plotting code runs
unchanged.

```{r predict-levels}
px <- exp(seq(log(0.05), log(20), length.out = 60))
ids_show <- unique(dat$id)[1:6]

pred_levels <- predict(
  fit_2re,
  newdata = expand.grid(x = px, id = ids_show),
  level = c("population", "subject"),
  scale = "natural"
)
head(pred_levels)
```

The population-mean curve is identical for every subject, so a single
`predict.fixed` column overlays the per-subject `predict.id` curves.

```{r predict-levels-plot, fig.cap = "Subject-conditional demand curves (grey) around the population-mean curve (dark)."}
ggplot(pred_levels, aes(x = x)) +
  geom_line(aes(y = predict.id, group = id), colour = "grey65") +
  geom_line(
    data = subset(pred_levels, id == ids_show[1]),
    aes(y = predict.fixed), colour = "#2c3e50", linewidth = 1.2
  ) +
  scale_x_log10() +
  labs(x = "Price", y = "Predicted consumption") +
  theme_minimal()
```

### Population Metrics

```{r results-group-metrics}
calc_group_metrics(fit_2re)
```

### Visualization

```{r plot-demand, fig.cap = "Population-level demand curve with confidence band."}
plot(fit_2re, type = "demand")
```

```{r plot-individual, fig.cap = "Individual demand curves for a random sample of subjects."}
plot(fit_2re, type = "individual")
```

```{r plot-parameters, fig.cap = "Distribution of subject-specific demand parameters."}
plot(fit_2re, type = "parameters")
```

## Diagnostics

After fitting, assess model health with the built-in diagnostic tools.

The fit object exposes a `hessian_pd` field that flags whether the Hessian
returned by `TMB::sdreport()` is positive definite. When `FALSE`, standard
errors, p-values, and Wald confidence intervals derived from the Hessian
should be treated as unreliable; the warning is also surfaced as a note in
`summary()` output and as a `hessian_warning` attribute on `tidy()` output.

```{r diagnostics-hessian}
fit_2re$hessian_pd
```

```{r diagnostics-check}
# Model health check — convergence, variance components, residual stats,
# and (since 0.3.0) Hessian positive-definiteness.
check_demand_model(fit_2re)
```

```{r diagnostics-augment}
# Augment with fitted values and residuals
aug <- augment(fit_2re)
head(aug[, c("id", "x", "y", ".fitted", ".resid", ".std_resid")])
```

```{r diagnostics-qq, fig.cap = "Q-Q plots of subject-level random effects."}
# Q-Q plot of random effects
plot_qq(fit_2re)
```

```{r diagnostics-re, fig.cap = "Random effects diagnostic panels."}
# Random effects diagnostic panels
plot_re_diagnostics(fit_2re)
```

```{r diagnostics-resid, fig.cap = "Residuals vs fitted values for the 2-RE exponential model."}
# Residual plot — standard in every modeling workflow
plot_residuals(fit_2re, type = "fitted")
```

These diagnostics help identify:

- **Q-Q plots**: Non-normality of random effects (heavy tails, outliers)
- **`check_demand_model()`**: Convergence issues, boundary estimates, residual patterns
- **Augmented data**: Observation-level residuals for identifying poorly fitting subjects
- **Residual plots**: Heteroscedasticity, non-linearity, or outliers in the fitted vs residual pattern

## Advanced Visualization

`beezdemand` provides several specialized plots beyond the standard `plot()` method.
These work with any `beezdemand_tmb` object.

### Expenditure and Elasticity

The expenditure curve shows total spending ($P \times Q$) as a function of price,
with vertical and horizontal reference lines at $P_{max}$ and $O_{max}$:

```{r viz-expenditure, fig.cap = "Expenditure curve with Pmax and Omax reference lines."}
plot_expenditure(fit_2re)
```

The elasticity curve shows how responsive demand is to price changes. The dashed
line at $-1$ marks unit elasticity — prices above this threshold produce elastic
demand:

```{r viz-elasticity, fig.cap = "Own-price elasticity curve with unit elasticity reference."}
plot_elasticity(fit_2re)
```

### Loss Surface and Profile

The loss surface visualizes the sum-of-squared-residuals landscape over a grid
of $(Q_0, \alpha)$ values, with the MLE marked. This helps assess
identifiability — a sharp, well-defined minimum indicates good identification:

```{r viz-loss-surface, fig.cap = "2D loss surface (SSR) over the Q0-alpha parameter space."}
plot_loss_surface(fit_2re)
```

Profile plots show 1D slices through the loss surface, fixing one parameter at
its MLE and varying the other:

```{r viz-loss-profile, fig.cap = "1D loss profiles for Q0 and alpha."}
plot_loss_profile(fit_2re)
```

### Subject Heterogeneity

Visualize the distribution of subject-level $\alpha$ estimates to assess
heterogeneity in price sensitivity:

```{r viz-alpha-dist, fig.cap = "Distribution of subject-level alpha (elasticity) estimates."}
plot_alpha_distribution(fit_2re)
```

### Multi-Model Comparison

Overlay demand curves from multiple models to visualize how different
specifications affect the predicted demand function:

```{r viz-overlay, fig.cap = "Demand curves from 1-RE and 2-RE models overlaid."}
plot_demand_overlay(fit_1re, fit_2re, labels = c("1-RE", "2-RE"))
```

Compare parameter distributions side by side:

```{r viz-model-comparison, fig.cap = "Side-by-side parameter estimates from 1-RE and 2-RE models."}
plot_model_comparison(fit_1re, fit_2re, labels = c("1-RE", "2-RE"))
```

## Group Comparisons

To test whether demand parameters differ by group, pass factor variables via
the `factors` argument. This adds fixed effects to the design matrices for both
$Q_0$ and $\alpha$.

```{r factor-setup}
# Filter to Male/Female for a clean two-level comparison
dat_mf <- dat |> filter(gender %in% c("Male", "Female"))

fit_gender <- fit_demand_tmb(
  dat_mf,
  y_var = "y", x_var = "x", id_var = "id",
  equation = "exponential",
  factors = "gender",
  random_effects = c("q0", "alpha"),
  verbose = 0
)
fit_gender
```

### Joint Tests

`anova(fit)` reports a joint Wald-χ² test for each parameter × factor block —
here, whether `gender` shifts $Q_0$ and $\alpha$. Pass additional fits for a
nested likelihood-ratio test.

```{r results-anova}
anova(fit_gender)
```

### Estimated Marginal Means

```{r emms-q0}
get_demand_param_emms(fit_gender, param = "Q0")
```

```{r emms-alpha}
get_demand_param_emms(fit_gender, param = "alpha")
```

EMMs are reported on both the natural scale (`estimate`) and log scale
(`estimate_log`). The natural-scale estimates represent the population-average
$Q_0$ or $\alpha$ for each group.

### Pairwise Comparisons

```{r comparisons-q0}
get_demand_comparisons(fit_gender, param = "Q0")
```

```{r comparisons-alpha}
get_demand_comparisons(fit_gender, param = "alpha")
```

The `estimate_ratio` column gives the multiplicative ratio between groups on
the natural scale (e.g., a ratio of 1.34 means Group A's $Q_0$ is 34% higher
than Group B's).

```{r factor-interaction, eval = !render_fast, include = !render_fast}
# Factor interaction: gender x age group
dat_mf$age_group <- ifelse(
  dat_mf$age < median(dat_mf$age, na.rm = TRUE), "younger", "older"
)

fit_interact <- fit_demand_tmb(
  dat_mf, y_var = "y", x_var = "x", id_var = "id",
  equation = "exponential",
  factors = c("gender", "age_group"),
  factor_interaction = TRUE,
  random_effects = c("q0", "alpha"),
  verbose = 0
)
get_demand_param_emms(fit_interact, param = "Q0")
```

```{r asis-interaction-note, eval = render_fast, results = 'asis', echo = FALSE}
cat("*Factor interaction and continuous covariate sections are skipped in fast",
    "render mode. Set `BEEZDEMAND_VIGNETTE_MODE=full` to include them.*\n")
```

```{r continuous-covariates, eval = !render_fast, include = !render_fast}
# Continuous covariates enter the design matrices for Q0 and alpha
fit_cov <- fit_demand_tmb(
  dat_mf, y_var = "y", x_var = "x", id_var = "id",
  equation = "exponential",
  continuous_covariates = "age",
  random_effects = c("q0", "alpha"),
  verbose = 0
)
tidy(fit_cov)
```

```{r collapse-levels, eval = !render_fast, include = !render_fast}
# Asymmetric factor collapsing: different groupings for Q0 vs alpha
dat_mf$age_group <- cut(
  dat_mf$age,
  breaks = c(0, 25, 40, Inf),
  labels = c("young", "mid", "older")
)

fit_collapsed <- fit_demand_tmb(
  dat_mf, y_var = "y", x_var = "x", id_var = "id",
  equation = "exponential",
  factors = "age_group",
  # Collapse "young" and "mid" into "young_mid" for Q0,

  # but collapse "mid" with "older" for alpha
  collapse_levels = list(
    Q0 = list(age_group = c("young" = "young_mid", "mid" = "young_mid",
                            "older" = "older")),
    alpha = list(age_group = c("young" = "young", "mid" = "older",
                               "older" = "older"))
  ),
  random_effects = "q0",
  verbose = 0
)
tidy(fit_collapsed)
```

## Model Comparison

Use `compare_models()` to compare nested TMB models via likelihood ratio test,
AIC, and BIC:

```{r model-comparison}
compare_models(fit_1re, fit_2re)
```

**Valid comparisons** require models fit on the same data with the same
equation and response scale. Models with different equations (e.g., exponential
vs exponentiated) model different responses and cannot be compared via AIC.

### Building Nested Models with `update()`

`update(fit, ...)` re-fits with named arguments substituted into the
original `fit_demand_tmb()` call. Use it to build a candidate fit and
drop or add factors / covariates / random effects without typing the
full fit call again:

```{r update-nested}
# Drop the gender factor to test its joint significance:
fit_gender_null <- update(fit_gender, factors = NULL)
compare_models(fit_gender_null, fit_gender)
```

`evaluate = FALSE` returns the unevaluated call (for inspection) instead
of re-fitting, matching the convention of `stats::update.default`.

`formula(fit)` and `model.matrix(fit)` round out the introspection API.
`formula(fit)` returns a named list of one-sided formulas for `Q0` and
`alpha` plus the random-effect spec; `model.matrix(fit)` returns the
four design matrices the TMB template consumed (`X_q0`, `X_alpha`,
`Z_q0`, `Z_alpha`) — a named list rather than the single matrix `lm`
or `lme4` return, because the TMB tier truly has two fixed-effect
linear predictors.

## Convergence Tips

If a model fails to converge, try these strategies in order:

| Strategy | How | When |
|----------|-----|------|
| Reduce random effects | `random_effects = "q0"` | 2-RE models struggling |
| Use exponential equation | `equation = "exponential"` | Exponentiated/simplified not converging |
| Fix k | `estimate_k = FALSE, k = 2` | k estimate drifting to extreme values |
| Increase iterations | `tmb_control = list(iter_max = 2000)` | "false convergence" messages |
| Try L-BFGS-B | `tmb_control = list(optimizer = "L-BFGS-B")` | nlminb not making progress |
| Set parameter bounds | `tmb_control = list(lower = ..., upper = ...)` | Estimates at boundaries |
| Warm start | `tmb_control = list(warm_start = prev_fit$opt$par)` | Refining a near-converged fit |

```{r convergence-examples, eval = FALSE}
# Switch optimizer
fit <- fit_demand_tmb(
  dat, equation = "exponential",
  tmb_control = list(optimizer = "L-BFGS-B"),
  verbose = 0
)

# Warm-start from a previous fit
fit2 <- fit_demand_tmb(
  dat, equation = "exponential",
  tmb_control = list(warm_start = fit$opt$par),
  verbose = 0
)

# Apply parameter bounds
fit3 <- fit_demand_tmb(
  dat, equation = "exponential",
  tmb_control = list(
    lower = c(log_k = -2),
    upper = c(log_k = 4)
  ),
  verbose = 0
)

# Disable multi-start for faster iteration during exploration
fit4 <- fit_demand_tmb(
  dat, equation = "exponential",
  multi_start = FALSE,
  verbose = 2
)
```

### tmb_control Options

| Field | Default | Description |
|-------|---------|-------------|
| `optimizer` | `"nlminb"` | `"nlminb"` or `"L-BFGS-B"` |
| `iter_max` | 1000 | Maximum iterations |
| `eval_max` | 2000 | Maximum evaluations (nlminb only) |
| `rel_tol` | 1e-10 | Relative tolerance (nlminb only) |
| `lower` / `upper` | NULL | Named numeric bounds on log-scale parameters |
| `warm_start` | NULL | Starting values from a previous `fit$opt$par` |
| `trace` | 0 | Optimizer trace level |

## Choosing Between fit_demand_tmb() and fit_demand_mixed()

| Feature | `fit_demand_tmb()` | `fit_demand_mixed()` |
|---------|-------------------|---------------------|
| Backend | TMB (C++, automatic differentiation) | nlme (R, numerical gradients) |
| Equations | exponential, exponentiated, simplified, zben | zben, simplified |
| k parameter | Estimated or fixed | Not available |
| Random effects | 1 or 2 (Q0, alpha) | Configurable via nlme |
| Convergence | Robust (AD + Laplace + multi-start) | Can struggle with nonlinear equations |
| Speed | Fast (compiled C++) | Variable |
| Post-hoc EMMs | `get_demand_param_emms()` | `get_demand_param_emms()` (via emmeans) |
| Factors/covariates | Design matrices | Formula-based |

Prefer `fit_demand_tmb()` for new work. `fit_demand_mixed()` remains useful
when you need specific `nlme` features like custom correlation structures or
when working with existing pipelines.

## References

Hursh, S. R., & Silberberg, A. (2008). Economic demand and essential value.
*Psychological Review*, 115(1), 186--198.

Koffarnus, M. N., Franck, C. T., Stein, J. S., & Bickel, W. K. (2015). A
modified exponential behavioral economic demand model to better describe
consumption data. *Experimental and Clinical Psychopharmacology*, 23(6),
504--512.

Kristensen, K., Nielsen, A., Berg, C. W., Skaug, H., & Bell, B. M. (2016).
TMB: Automatic differentiation and Laplace approximation. *Journal of
Statistical Software*, 70(5), 1--21.

## See Also

- `vignette("fixed-demand")` -- Individual NLS demand curves
- `vignette("mixed-demand")` -- NLME-based mixed-effects models
- `vignette("mixed-demand-advanced")` -- Advanced topics: multi-factor designs, collapse_levels
- `vignette("hurdle-demand-models")` -- Two-part hurdle models for zero-heavy data
- `vignette("model-selection")` -- Choosing the right demand model
