Understanding Convergence in Demand Models

Optimization in Nonlinear Models

Unlike linear regression, which has a closed-form solution, nonlinear demand models must be fitted iteratively. The optimizer starts from initial parameter values and takes steps through parameter space, trying to minimize an objective function:

• TMB models: minimize negative log-likelihood (NLL)
• NLME models: minimize restricted/maximum likelihood (REML/ML)
• NLS models: minimize residual sum of squares (RSS)

At each iteration, the optimizer evaluates the objective function, computes gradients (the direction of steepest descent), and takes a step. It stops when one of several stopping rules fires:

1. Gradient near zero: the objective function is flat, suggesting a minimum
2. Parameter change below tolerance: successive iterations produce nearly identical estimates
3. Iteration limit reached: the budget of iterations is exhausted
4. Step size collapses: the optimizer can’t make progress

The convergence status you see in model output tells you which stopping rule fired. It does not tell you whether the solution is the global minimum, whether the parameters are scientifically plausible, or whether the model fits the data well.

TMB Convergence Codes (nlminb)

TMB models fitted with fit_demand_tmb() and fit_demand_hurdle() use the nlminb optimizer by default. It reports the following convergence codes:

Code 0 is what you want but doesn’t guarantee the solution is correct — the optimizer may have converged to a local minimum or a boundary.

Code 1 (iteration limit) is often benign. If the objective function was still decreasing, increasing iter_max in tmb_control may resolve it. If the objective had plateaued, the current estimates are likely at or very near the minimum.

Code 8 (false convergence) sounds alarming but is common in well-specified demand models. TMB provides exact gradients via automatic differentiation, so code 8 is rarer than with finite-difference methods. When it does occur, it typically means the solution is at a feasibility boundary (e.g., a random effect variance near zero) rather than that the optimizer is stuck.

NLME Convergence Warnings

NLME models fitted with fit_demand_mixed() use a fundamentally different optimization architecture than TMB. Rather than a single optimizer minimizing a single objective, NLME alternates between two procedures:

1. EM algorithm: estimates variance components (random effect variances and correlations)
2. PNLS (penalized nonlinear least squares): estimates fixed effects and predicts random effects, conditional on the variance components

These two steps alternate until the overall likelihood converges. After convergence, the approximate variance-covariance matrix (apVar) is computed from the Hessian at the solution to provide standard errors for the variance components.

NLME signals problems through R warnings rather than numeric codes. The package captures these during fitting and stores them in the fit_warnings field of the returned object. The following table summarizes the most common warnings:

# After fitting an NLME model
fit <- fit_demand_mixed(data, y_var = "y_ll4", equation_form = "zben")

# apVar is a matrix when successful, a character string describing the error
# when it fails
is_apvar_ok <- is.matrix(fit$model$apVar)

NLS Convergence (Fixed-Effect Models)

fit_demand_fixed() handles convergence automatically through a multi-backend fallback strategy. For each subject, the package tries:

1. nls.multstart() with random starting values
2. If that fails, minpack.lm::nlsLM() (Levenberg-Marquardt)
3. If that fails, nlsr::wrapnlsr()

Models either converge or fail outright; partial convergence doesn’t occur. The returned object tracks per-subject success/failure in the Notes column. Use check_demand_model() to see the failure rate:

fit <- fit_demand_fixed(apt, equation = "hs")
diagnostics <- check_demand_model(fit)
print(diagnostics)

A high failure rate (>20%) suggests data quality problems — insufficient price range, too many zero-consumption observations, or extreme outliers — rather than modeling issues.

The Key Insight

Convergence status reflects the optimizer’s stopping rule, not the quality of the solution. A non-converged fit with plausible parameters and valid standard errors can be more trustworthy than a converged fit at a boundary or local minimum. This principle applies equally to TMB numeric codes and NLME warning messages.

The next section provides a systematic checklist for evaluating fit quality regardless of convergence status.

1. Parameter Plausibility

The most important check is whether the estimated parameters make scientific sense. For demand models: - Q0 (intensity) should be positive and in the range of observed consumption at low prices - alpha (elasticity) should be positive and not extreme - k (scaling constant) typically falls between 1 and 4 for purchase task data

fit_tmb <- fit_demand_tmb(data, y_var = "y", x_var = "x", id_var = "id")

# Population-level (fixed) parameters
coef(fit_tmb)

# Full summary with standard errors
summary(fit_tmb)

# Subject-level parameters (fixed + random effects)
subject_pars <- get_subject_pars(fit_tmb)
head(subject_pars)

# Check ranges
summary(subject_pars$Q0)
summary(subject_pars$alpha)

fit_nlme <- fit_demand_mixed(data, y_var = "y_ll4", equation_form = "zben")

# Population-level fixed effects
fixef(fit_nlme$model)

# Subject-level deviations from population mean
ranef(fit_nlme$model)

# Subject-level parameters (fixed + random)
subject_pars <- get_subject_pars(fit_nlme)
head(subject_pars)

fit_hurdle <- fit_demand_hurdle(data, y_var = "y", x_var = "x", id_var = "id")

# Population parameters
coef(fit_hurdle)

# Subject-level parameters
subject_pars <- get_subject_pars(fit_hurdle)
head(subject_pars)

2. Standard Errors and Variance Estimates

Finite, reasonable standard errors indicate that the likelihood surface is well-behaved near the solution.

# Check if standard errors were successfully computed
fit_tmb$se_available

# Hessian positive-definiteness gate (since 0.3.0): TRUE means the inverse
# Hessian is well-conditioned, so SEs / p-values / Wald CIs are trustworthy.
# FALSE means they are not — and `fit_demand_tmb()` / `fit_demand_hurdle()`
# also emit a cli warning at fit time when verbose >= 1.
fit_tmb$hessian_pd

# View SEs in summary output (notes section flags pdHess = FALSE)
summary(fit_tmb)

# Check apVar status
is_apvar_ok <- is.matrix(fit_nlme$model$apVar)
cat("apVar OK:", is_apvar_ok, "\n")

# If apVar is OK, get confidence intervals on variance components
if (is_apvar_ok) {
  intervals(fit_nlme$model, which = "var-cov")
}

# Fixed effect standard errors are always available (from the PNLS step)
summary(fit_nlme$model)$tTable

3. Gradient Norm (TMB-Specific)

For TMB models, you can check the gradient at the solution directly. A small gradient confirms the optimizer stopped at (or very near) a stationary point.

# Compute gradient at the solution
grad <- fit_tmb$tmb_obj$gr(fit_tmb$opt$par)
max_grad <- max(abs(grad))
cat("Maximum absolute gradient:", max_grad, "\n")

# Rule of thumb:
# < 0.001: excellent — optimizer found a clean minimum
# 0.001-0.01: good — likely near a minimum
# > 0.01: concerning — may not be at a minimum

NLME does not expose gradients directly, so this check is TMB-only. For NLME models, rely on the other five checks.

4. Hessian / apVar Status

The Hessian matrix (second derivatives of the objective function) tells you about the curvature at the solution. A positive-definite Hessian confirms the solution is at a local minimum (not a saddle point or maximum).

# Preferred (since 0.3.0): the dedicated `hessian_pd` accessor on the fit
# object. TRUE = positive definite; FALSE = not; NA = sdreport failed.
fit_tmb$hessian_pd

# Equivalent low-level path through the sdreport object:
if (!is.null(fit_tmb$sdr)) {
  isTRUE(fit_tmb$sdr$pdHess)
}

out <- tidy(fit_tmb)
attr(out, "hessian_warning")

# Check apVar (NLME's version of the Hessian check)
cat("apVar OK:", is.matrix(fit_nlme$model$apVar), "\n")

# Check random effects correlation for near-singularity
vc <- nlme::VarCorr(fit_nlme$model)
print(vc)
# Look at the Corr column — values > |0.99| indicate near-singularity

5. Comparison with Nested Converged Models

If a non-converged complex model beats a simpler converged model by a large AIC margin, the parameters are in a good region of the likelihood surface. AIC and BIC depend on the log-likelihood at the parameter values, not on convergence status.

# Fit a simpler model (e.g., fewer random effects)
fit_simple <- fit_demand_tmb(data, y_var = "y", x_var = "x", id_var = "id",
                             random_effects = "q0")
fit_complex <- fit_demand_tmb(data, y_var = "y", x_var = "x", id_var = "id",
                              random_effects = c("q0", "alpha"))

# Compare
AIC(fit_simple, fit_complex)

# Fit with random effect on Q0 only vs. Q0 + alpha
fit_q0 <- fit_demand_mixed(data, y_var = "y_ll4", equation_form = "zben",
                           random_effects = Q0 ~ 1)
fit_both <- fit_demand_mixed(data, y_var = "y_ll4", equation_form = "zben",
                             random_effects = Q0 + alpha ~ 1)

# Compare AIC
AIC(fit_q0$model, fit_both$model)

# Or use anova for a formal likelihood ratio test
anova(fit_q0$model, fit_both$model)

6. Stability Across Starting Values

If different starting values lead to the same (or very similar) parameter estimates, you can be confident the solution is not a local minimum.

# Fit with default multi-start
fit1 <- fit_demand_tmb(data, y_var = "y", x_var = "x", id_var = "id")

# Fit with warm start from the first fit
fit2 <- fit_demand_tmb(data, y_var = "y", x_var = "x", id_var = "id",
                       tmb_control = list(warm_start = fit1$opt$par))

# Compare NLL values
cat("Fit 1 NLL:", fit1$opt$objective, "\n")
cat("Fit 2 NLL:", fit2$opt$objective, "\n")

# Fit with pooled NLS starting values (data-driven)
fit_pooled <- fit_demand_mixed(data, y_var = "y_ll4", equation_form = "zben",
                               start_values = "pooled_nls")

# Fit with different manual starting values
fit_manual <- fit_demand_mixed(data, y_var = "y_ll4", equation_form = "zben",
                               start_values = list(
                                 fixed = c(Q0 = 5, alpha = 0.001)
                               ))

# Compare fixed effects
fixef(fit_pooled$model)
fixef(fit_manual$model)

TMB Decision Table

Use this table after running the diagnostic checklist on a TMB model (fit_demand_tmb() or fit_demand_hurdle()):

NLME Decision Table

Use this table after running the diagnostic checklist on an NLME model (fit_demand_mixed()):

NLME Remedies

NLME convergence issues typically fall into a few categories. The table below maps symptoms to solutions:

fit <- fit_demand_mixed(
  data,
  y_var = "y_ll4",
  equation_form = "zben",
  nlme_control = list(
    msMaxIter = 500,  # Max iterations for the outer optimization (default: 200)
    maxIter = 500,    # Max iterations for the PNLS step (default: 200)
    niterEM = 200     # Number of EM iterations before switching to Newton (default: 100)
  )
)

fit <- fit_demand_mixed(
  data,
  y_var = "y_ll4",
  equation_form = "zben",
  nlme_control = list(
    pnlsTol = 1e-2,     # Loosen PNLS tolerance (default: 1e-3)
    tolerance = 1e-5,    # Slightly loosen overall tolerance (default: 1e-6)
    niterEM = 150        # More EM iterations can help with difficult surfaces
  )
)

# Default uses pdDiag (diagonal); try pdSymm for unstructured if you have
# enough data, or stay with pdDiag if convergence is problematic
fit_diag <- fit_demand_mixed(
  data,
  y_var = "y_ll4",
  equation_form = "zben",
  covariance_structure = "pdDiag"  # No correlation between random effects
)

# Compare with unstructured (pdSymm) if convergence allows
fit_symm <- fit_demand_mixed(
  data,
  y_var = "y_ll4",
  equation_form = "zben",
  covariance_structure = "pdSymm"  # Estimate correlation between random effects
)

# Random effects on both Q0 and alpha (more complex)
fit_both <- fit_demand_mixed(
  data,
  y_var = "y_ll4",
  equation_form = "zben",
  random_effects = Q0 + alpha ~ 1
)

# Random effect on Q0 only (simpler, often sufficient)
fit_q0_only <- fit_demand_mixed(
  data,
  y_var = "y_ll4",
  equation_form = "zben",
  random_effects = Q0 ~ 1
)

# Compare fit
AIC(fit_both$model, fit_q0_only$model)

TMB Remedies

TMB convergence issues have their own set of solutions:

fit <- fit_demand_tmb(
  data,
  y_var = "y", x_var = "x", id_var = "id",
  tmb_control = list(
    iter_max = 2000,      # Increase iteration limit (default: 1000)
    eval_max = 4000,      # Increase function evaluation limit (default: 2000)
    rel_tol = 1e-12       # Tighten relative tolerance for more precision
  )
)

fit_lbfgsb <- fit_demand_tmb(
  data,
  y_var = "y", x_var = "x", id_var = "id",
  tmb_control = list(optimizer = "L-BFGS-B")
)

# First fit (may not fully converge)
fit1 <- fit_demand_tmb(data, y_var = "y", x_var = "x", id_var = "id")

# Use as warm start for a second attempt
fit2 <- fit_demand_tmb(data, y_var = "y", x_var = "x", id_var = "id",
                       tmb_control = list(warm_start = fit1$opt$par))

# Reduce random effects: 2-RE instead of 3-RE
fit_2re <- fit_demand_hurdle(
  data,
  y_var = "y", x_var = "x", id_var = "id",
  random_effects = c("zeros", "q0")  # Drop alpha random effect
)

# Try a different Part II equation
fit_snd <- fit_demand_hurdle(
  data,
  y_var = "y", x_var = "x", id_var = "id",
  part2 = "simplified_exponential"  # Fewer parameters (no k)
)

Fixed-Effect NLS Remedies

fit_demand_fixed() handles convergence automatically via multi-backend fallback. When many subjects fail to converge, the issue is almost always data quality rather than model specification:

• Too few price points: demand curves need at least 4-5 prices to fit reliably
• No variation in consumption: if consumption is constant across prices, the model has nothing to fit
• Extreme outliers: a single wildly atypical data point can prevent convergence
• Too many zeros: consider fit_demand_hurdle() instead, which handles zeros explicitly

# Check data quality first
systematic <- CheckUnsystematic(data, deltaq = 0.025, bounce = 0.1, reversals = 0)
table(systematic$TotalPass)

# Try a different equation if "hs" fails for many subjects
fit_koff <- fit_demand_fixed(data, equation = "koff")

TMB Example

fit_hurdle <- fit_demand_hurdle(apt, y_var = "y", x_var = "x", id_var = "id")
diagnostics <- check_demand_model(fit_hurdle)
print(diagnostics)

The output shows:

• Convergence: status and any optimizer messages
• Random Effects: variance estimates and whether any are near zero
• Residuals: summary statistics and outlier count
• Issues/Recommendations: actionable suggestions

NLME Example

fit_nlme <- fit_demand_mixed(data, y_var = "y_ll4", equation_form = "zben")
diagnostics <- check_demand_model(fit_nlme)
print(diagnostics)

For NLME models, the diagnostics additionally check for:

• Stored fit warnings (e.g., “false convergence”)
• apVar status (Hessian positive-definiteness)
• Random effect correlation near-singularity (|corr| > 0.99)

Live Example with Fixed-Effect Model

Here is a live example using the built-in apt dataset:

data(apt)
fit <- fit_demand_fixed(apt, equation = "hs")
diagnostics <- check_demand_model(fit)
print(diagnostics)
#> 
#> Model Diagnostics
#> ================================================== 
#> Model class: beezdemand_fixed 
#> 
#> Convergence:
#>   Status: Converged
#> 
#> Residuals:
#>   Mean: 0.0284
#>   SD: 0.5306
#>   Range: [-1.458, 2.228]
#>   Outliers: 3 observations
#> 
#> --------------------------------------------------
#> Issues Detected (1):
#>   1. Detected 3 potential outliers across subjects

Accessing Fit Warnings Directly

For NLME models, you can also inspect the captured warnings directly:

# Warnings captured during NLME fitting
fit_nlme$fit_warnings

# These are also surfaced by glance()
glance_result <- broom::glance(fit_nlme)
glance_result$converged

Cross-References

• vignette("hurdle-demand-models") — Fitting TMB hurdle models
• vignette("mixed-demand") — Fitting NLME mixed-effects models
• vignette("mixed-demand-advanced") — Advanced NLME topics including control tuning
• vignette("fixed-demand") — Fitting fixed-effect NLS models
• vignette("model-selection") — Choosing between model tiers