DRM.jl

Appearance

Profile-likelihood intervals ​

Status — Stable

Mirrors drmTMB's Profile-likelihood intervals. In DRM.jl today: profile-likelihood confidence intervals via confint(fit; method = :profile), the auditable profile_result object behind them, and the profile_curve data for a likelihood-ratio diagnostic plot.

A profile-likelihood interval asks a different question from a Wald interval. A Wald interval uses the local curvature at the fitted estimate — fast, but only as good as the quadratic approximation there. A profile interval fixes one target at a sequence of values, re-optimises every other parameter at each one, and keeps the values whose likelihood-ratio distance is still compatible with the fitted model. When a parameter sits near a boundary or its likelihood is skewed, the two intervals can differ noticeably, and a successfully computed profile can reveal asymmetry that Wald misses.

Fit a small model ​

The example asks whether a continuous response changes with temperature while the residual standard deviation is constant:

julia
using DRM, Random
Random.seed!(20260608)

n = 70
temperature = 2.4 .* rand(n) .- 1.2                 # Uniform(-1.2, 1.2)
growth = 0.5 .+ 0.8 .* temperature .+ 0.55 .* randn(n)
dat = (; growth, temperature)

fit = drm(bf(@formula(growth ~ temperature), @formula(sigma ~ 1)),
          Gaussian(); data = dat)
coef(fit, :mu)
2-element Vector{Float64}:
 0.5310004510619495
 0.9114840433792056

Wald first, then profile ​

confint returns one row per coefficient — (param, coef, estimate, lower, upper) — on each parameter's working scale (μ on the response scale, σ on log σ). The default is the Wald interval:

julia
confint(fit)                          # method = :wald (the default)
3-element Vector{@NamedTuple{param::Symbol, coef::String, estimate::Float64, lower::Float64, upper::Float64}}:
 (param = :mu, coef = "(Intercept)", estimate = 0.5310004510619495, lower = 0.4001712447241781, upper = 0.6618296573997209)
 (param = :mu, coef = "temperature", estimate = 0.9114840433792056, lower = 0.7184994464556558, upper = 1.1044686403027555)
 (param = :sigma, coef = "(Intercept)", estimate = -0.5861537908239811, lower = -0.7518009808888497, upper = -0.4205066007591125)

Switch to the profile interval with method = :profile. Each endpoint is found by re-optimising the nuisance parameters along a likelihood-ratio root-search, so it costs more than Wald but does not assume a quadratic log-likelihood:

julia
confint(fit; method = :profile)
3-element Vector{@NamedTuple{param::Symbol, coef::String, estimate::Float64, lower::Float64, upper::Float64}}:
 (param = :mu, coef = "(Intercept)", estimate = 0.5310004510619495, lower = 0.3983556429750908, upper = 0.6636452591488087)
 (param = :mu, coef = "temperature", estimate = 0.9114840433792056, lower = 0.715821274038889, upper = 1.1071468127195219)
 (param = :sigma, coef = "(Intercept)", estimate = -0.5861537908239811, lower = -0.7431378591127786, upper = -0.41083252372690837)

For a single target, pass parm to profile just that parameter block — here the constant residual scale, on the log σ scale:

julia
confint(fit; method = :profile, parm = :sigma)
1-element Vector{@NamedTuple{param::Symbol, coef::String, estimate::Float64, lower::Float64, upper::Float64}}:
 (param = :sigma, coef = "(Intercept)", estimate = -0.5861537908239811, lower = -0.7431378591127786, upper = -0.41083252372690837)

The auditable profile object ​

confint(...; method = :profile) is a thin wrapper over profile_result, which returns the full object so a profile can be audited, not just trusted. Alongside ci it reports the per-endpoint work counts, the wall-clock time, and which nuisance-gradient backend was used:

julia
res = profile_result(fit; parm = :sigma)
res.ci                                # same rows as confint(...; method = :profile)
1-element Vector{@NamedTuple{param::Symbol, coef::String, estimate::Float64, lower::Float64, upper::Float64}}:
 (param = :sigma, coef = "(Intercept)", estimate = -0.5861537908239811, lower = -0.7431378591127786, upper = -0.41083252372690837)
julia
(attempted = res.attempted, used = res.used, failed = res.failed,
 autodiff = res.autodiff, level = res.level)
(attempted = 1, used = 1, failed = 0, autodiff = :forward, level = 0.95)

The stats field carries one row per coefficient, with diagnostics for both endpoints: likelihood evaluations, bracket expansions, endpoint failures and searched-range limits. The generic profiler also records each arm's last nuisance method, fallback use and acceptance reason:

julia
res.stats
1-element Vector{@NamedTuple{param::Symbol, coef::String, evaluations::Int64, gradient_evaluations::Int64, bracket_expansions::Int64, root_iterations::Int64, lower_unbounded::Bool, upper_unbounded::Bool, nonmonotone::Bool, lower_endpoint_failed::Bool, upper_endpoint_failed::Bool, lower_nuisance_method::Symbol, upper_nuisance_method::Symbol, lower_nuisance_fallback::Bool, upper_nuisance_fallback::Bool, lower_nuisance_reason::Symbol, upper_nuisance_reason::Symbol}}:
 (param = :sigma, coef = "(Intercept)", evaluations = 14, gradient_evaluations = 14, bracket_expansions = 4, root_iterations = 8, lower_unbounded = 0, upper_unbounded = 0, nonmonotone = 0, lower_endpoint_failed = 0, upper_endpoint_failed = 0, lower_nuisance_method = :lbfgs_forward, upper_nuisance_method = :lbfgs_forward, lower_nuisance_fallback = 0, upper_nuisance_fallback = 0, lower_nuisance_reason = :accepted, upper_nuisance_reason = :accepted)

Check lower_endpoint_failed and upper_endpoint_failed first. A failed arm has no valid endpoint: its signed infinity is a failure placeholder, not evidence for an unbounded statistical interval. confint warns when an endpoint fails; use profile_result to retain the diagnostics. The generic path rejects non-finite or unsuccessfully terminated nuisance solves. Optimizer termination alone does not prove a global optimum or a sufficiently small score.

If lower_unbounded or upper_unbounded is true without endpoint failure, the threshold was not crossed within the searched range. This is a useful diagnostic, not proof that the mathematical interval is infinite. Specialized location-only and location-scale profilers retain their separate solver contracts; the generic nuisance diagnostics do not certify those paths.

Inspect the likelihood-ratio curve ​

To see the profile rather than only its endpoints, profile_curve returns the data for a 1-D likelihood-ratio plot of one coefficient (by its global index k). At each grid value of θ[k] the nuisance parameters are re-optimised. Failed generic grid solves raise an error identifying the grid location. A profiled objective materially below the reference fit also raises an error instead of being clamped to zero. It returns the grid x, the profile deviance (2(ℓ̂ − ℓ_profile)), the fitted estimate, and the χ²₁ cutoff the interval uses.

We profile the log σ coefficient. Its global index is the one whose coefficient block is :sigma:

julia
# global coefficient index of the (single) sigma coefficient
kσ = only(only(r for (p, r) in fit.blocks if p === :sigma))

curve = profile_curve(fit, kσ; npoints = 21)
(estimate = curve.estimate, cutoff = curve.cutoff, param = curve.param, coef = curve.coef)
(estimate = -0.5861537908239811, cutoff = 3.8414588206941245, param = :sigma, coef = "(Intercept)")

The interval endpoints are where the deviance curve crosses cutoff. The deviance is zero at the MLE and rises on both sides; the values of x where it equals cutoff are exactly the profile endpoints reported above:

julia
[(x = xi, deviance = di) for (xi, di) in zip(curve.x, curve.deviance)][1:5]
5-element Vector{@NamedTuple{x::Float64, deviance::Float64}}:
 (x = -0.8397000672425361, deviance = 10.735474938492786)
 (x = -0.8143454396006806, deviance = 8.53803376818675)
 (x = -0.788990811958825, deviance = 6.6247574200493204)
 (x = -0.7636361843169696, deviance = 4.9815953642220165)
 (x = -0.7382815566751141, deviance = 3.5951917993064626)

Crossings on both sides are useful checks, but do not by themselves certify optimizer accuracy, interval coverage or the absence of other optima. Inspect solver diagnostics and the curve. A nonmonotone curve requires investigation; the numerical search does not guarantee finding its first crossing.

When to reach for which ​

  • Wald (method = :wald, the default) — fast, fine when the log-likelihood is near-quadratic around the estimate and the parameter is well inside its range.

  • Profile (method = :profile) — slower but curvature-honest; prefer it for scale and correlation parameters, near boundaries, or whenever a Wald interval runs past a parameter's natural range.

  • Bootstrap (bootstrap_ci) — when even the profile's asymptotics are in doubt; see the figure gallery and the post-fit tour.

See also ​