support for GLMMs with a dispersion parameter - #914
Draft
palday wants to merge 18 commits into
Draft
Conversation
palday
force-pushed
the
pa/dispersion-again
branch
from
August 12, 2026 14:35
52893df to
348b884
Compare
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #914 +/- ##
==========================================
+ Coverage 94.88% 95.02% +0.14%
==========================================
Files 38 38
Lines 3946 4004 +58
==========================================
+ Hits 3744 3805 +61
+ Misses 202 199 -3
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
The Laplace and AGQ objectives carried the penalty as sum(u²) when it should be sum(u²)/ϕ. θ is relative to the scale parameter in this package -- Var(b) = ϕΛΛ', which is the σ factor σs and VarCorr were already applying -- so u ~ N(0, ϕI) and the penalty belongs on the ϕ scale along with the data term. The check that pins this down: for a Normal family pwrss/n is the conditional MLE of ϕ, so profiling the objective over ϕ at fixed β and θ must return pwrss/n. It does with uss/ϕ; with uss it misses by 16%. Equivalently, the corrected expression specialises to objective(::LinearMixedModel) for a Normal response and identity link, which the previous one did not. That identity is now a test. For AGQ the same rescaling applies in three coupled places: the integrand becomes (u² + Σdev)/ϕ, the quadrature nodes are spaced by √ϕ/L.diag because the conditional sd of u carries the same ϕ as the prior, and Cϕ picks up a compensating nᵤ·log ϕ for the resulting Jacobian term. _agq_deviance(m, 1) == _laplace_deviance(m) still holds exactly, which is what fixes those constants. This changes θ̂ and ϕ̂ for every dispersion-family fit; the references in test/pirls.jl are re-recorded. On the sleepstudy Gamma fit θ̂ goes from [0.944, -0.061, 0.143] to [1.248, -0.006, 0.217] and σ̂ from 0.0885 to 0.0808. The new values deliberately do not match lme4: glmer optimises the ϕ ≡ 1 penalty while reporting VarCorr on the relative scale, and that inconsistency is what issue #206 is about. Families without a dispersion parameter are untouched, bit for bit. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
palday
force-pushed
the
pa/dispersion-again
branch
2 times, most recently
from
August 12, 2026 17:19
2212770 to
3fc2457
Compare
ϕ is now estimated one of two ways, selected by `fast`:
fast=true plugs in the Pearson moment estimator pwrss/n, matching
lme4's sigma(). This is what every dispersion fit did
previously, and it is unchanged.
fast=false carries log ϕ as a parameter of the outer optimisation
alongside β and θ, so it converges to the conditional MLE.
The two are the same estimator only for Normal, where pwrss/n provably
is the conditional MLE -- the outer optimisation just rediscovers it,
and agreeing to 7e-5 on the sleepstudy fit is a useful check on the
whole arrangement. For Gamma and InverseGaussian the MLE solves a
digamma equation instead of a moment condition, and the two answers
separate (0.18% apart on the sleepstudy Gamma fit). Solving that
equation numerically through the outer optimiser avoids having to
special-case the score for each family.
log ϕ rather than ϕ keeps the parameter unbounded and comparably scaled
to θ (ϕ̂ ≈ 0.0065 there, three orders below the θ entries). It also
keeps ϕ out of the zero-snapping pass in fit!, which is keyed on a lower
bound of zero and would otherwise be free to send ϕ to 1.
Note that ϕ still does not enter PIRLS. The penalised objective is
(Σwᵢrᵢ² + ‖u‖²)/ϕ, so the 1/ϕ factors out of the minimisation over u:
the conditional modes are the same in both regimes and only the deviance
differs. _dispersion is the single place either estimator is read from,
so the objective, loglikelihood and dispersion cannot drift apart.
Supporting changes: GeneralizedLinearMixedModel gains a ϕ field, empty
when ϕ is plugged in and one element when it is free, which is also the
flag distinguishing the regimes; setβθ! reads the trailing log ϕ when
present; unfit! clears it so refit! can move between regimes;
restoreoptsum! recognises the third parameter-vector length; and the
PRIMA scale vector is sized for the extra parameter.
Tests cover both regimes, that Normal reproduces pwrss/n while Gamma
does not, and that the regime survives refit! and a saveoptsum /
restoreoptsum! round trip.
ϕ̂ from the joint fit is only reproducible to ~0.15% -- between runs, and
between a fit and its own refit! -- because β and θ move and ϕ̂ follows.
That is too loose to lock ϕ̂ against a literal value, and close enough to
the 0.18% moment-vs-MLE gap to be worth distinguishing. Conditionally on
the fitted β and θ, though, the outer optimiser matches the true
conditional MLE to 0.0009%, and the gap to the moment estimator is worth
3e-4 in deviance against an ftol that resolves 2e-9. The tests therefore
compare ϕ̂ to a golden-section conditional MLE within a single fit, and
assert the moment estimator is distinguishable from it, rather than
pinning either to a constant.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
palday
force-pushed
the
pa/dispersion-again
branch
from
August 13, 2026 06:37
3fc2457 to
bb8ebb2
Compare
How ϕ was estimated is a property of the fitted model, so `show` is a
better place for it than an `@info` in `fit!`: it stays available for as
long as the model does, rather than scrolling past once, and it no longer
prints one copy per replicate when `parametricbootstrap` refits in a loop.
The line is gated on `dispersion_parameter` and reads the same flag the
objective does, so it cannot disagree with the ϕ actually used:
Distribution: Gamma{Float64}
Link: LogLink()
Dispersion parameter ϕ: estimated jointly with β and θ
The `@suppress` wrappers added to the dispersion tests are removed again --
they were working around this noise, and suppressing stdout in tests also
hides genuine warnings. `@test_logs` on the fit now asserts that nothing
is emitted at all, and a new testset covers the `show` output for both
estimators and for a family without a dispersion parameter.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
`nothing` versus a value says what the two ϕ regimes are far more directly
than an empty versus a one-element vector, and it matches how `sigma` is
spelled in `OptSummary`.
The vector was never a positive choice: `GeneralizedLinearMixedModel` is an
immutable struct, and ϕ has to be written on every objective evaluation by
`setβθ!` as well as switched between regimes by `fit!`, `unfit!`, `refit!`
and `_simulate!`. A bare `Union{T,Nothing}` field cannot be reassigned, so
the vector was standing in for a mutable cell. A `Ref` is that cell, and
says so. The union is read once per objective evaluation, which is nothing
next to the PIRLS iteration around it.
isempty(m.ϕ) -> isnothing(m.ϕ[])
first(m.ϕ) -> m.ϕ[]
push!(m.ϕ, x) -> m.ϕ[] = x
empty!(m.ϕ) -> m.ϕ[] = nothing
length(m.ϕ) -> !isnothing(m.ϕ[]) (PRIMA scale vector)
No behavioural change.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
`_dispersion` gains a leading branch: a non-nothing `optsum.sigma` wins over both estimators and returns σ² directly. `fit!` then leaves ϕ out of the parameter vector entirely, so a fixed-ϕ fit optimises over β and θ alone. On the sleepstudy Gamma fit with σ = 0.09 that is 5 parameters rather than 6, ϕ̂ is 0.0081 exactly, and the deviance is 1735.099 against 1732.002 for the unconstrained fit -- worse, as a constrained fit must be, which is the check that ϕ is really pinned and not quietly re-estimated.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
docs/NEWS-update.jlto update the cross-references.closes #206, closes #291, closes #787