Code
using CairoMakie
using DataFrames
using LinearAlgebra
using MixedModels
using MixedModelsMakie
const progress = isinteractive()Convergence, singularity and all that
After working through this page you will be able to:
- explain what it means for a converged mixed model to be singular;
- recognize singularity from the estimated variance components and correlation parameters;
- relate singularity to model complexity and the likelihood criterion that balances fit against complexity.
NoteBefore you start
Prerequisites: Analysis of the sleepstudy data.
Datasets used: sleepstudy, plus small simulated examples (see the dataset catalog).
1 What does it mean for a converged model to be singular?
Add the packages to be used
Fit a model for reaction time in the sleepstudy example.
m01 = lmm(@formula(reaction ~ 1 + days + (1 + days|subj)), MixedModels.dataset(:sleepstudy); progress)print(m01)Linear mixed model fit by maximum likelihood
reaction ~ 1 + days + (1 + days | subj)
logLik -2 logLik AIC AICc BIC
-875.9697 1751.9393 1763.9393 1764.4249 1783.0971
Variance components:
Column Variance Std.Dev. Corr.
subj (Intercept) 565.52071 23.78068
days 32.68242 5.71685 +0.08
Residual 654.94016 25.59180
Number of obs: 180; levels of grouping factors: 18
Fixed-effects parameters:
──────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|)
──────────────────────────────────────────────────
(Intercept) 251.405 6.6323 37.91 <1e-99
days 10.4673 1.50224 6.97 <1e-11
──────────────────────────────────────────────────The covariance matrix for each subject’s random effects is evaluated from its “matrix square root”, called the Cholesky factor.
λ = only(m01.λ)2×2 LowerTriangular{Float64, Matrix{Float64}}:
0.92923 ⋅
0.0181644 0.222646The transpose of \(\lambda\), written \(\lambda'\), is an upper triangular matrix generated by “flipping” \(\lambda\) about the main diagonal.
λ'2×2 UpperTriangular{Float64, Adjoint{Float64, Matrix{Float64}}}:
0.92923 0.0181644
⋅ 0.222646The product \(\lambda * \lambda'\) will be symmetric. The covariance matrix of the random effects, \(\Sigma\), is this symmetric matrix scaled by \(\sigma^2\)
Σ = m01.σ^2 * λ * λ'2×2 Matrix{Float64}:
565.521 11.0547
11.0547 32.6824The estimated variances of the random effects, which are the diagonal elements of \(\Sigma\), correspond to the values shown in the table. To evaluate the covariance, isolate the correlation
# m01.σρs extracts the `σρs` property of the model.
# This property is a NamedTuple where the names
# correspond to grouping factors - in this case, `subj`.
# So `m01.σρs.subj.ρ` is the estimated correlation(s) for
# this grouping factor. Because there is only one such correlation
# we can extract it with `only()`, which also verifies that
# there is exactly one.
ρ = only(m01.σρs.subj.ρ)0.08131378907082809and multiply by the standard deviations
ρ * sqrt(first(Σ) * last(Σ))11.054663967907072 The factor is generated from a parameter vector
In practice we optimize the log-likelihood with respect to a parameter vector called \(\theta\) that generates \(\lambda\).
m01.θ3-element Vector{Float64}:
0.9292304825690973
0.018164375827454657
0.22264643709604015The elements of this parameter vector are subject to constraints. In particular, two of the three elements have a lower bound of zero.
lowerbd(m01)3-element Vector{Float64}:
0.0
-Inf
0.0That is, the first and third elements of \(\theta\), corresponding to diagonal elements of \(\lambda\), must be non-negative, whereas the second component is unconstrained (has a lower bound of \(-\infty\)).
3 Progress of the iterations
The optsum.fitlog property of the model is table that contains the value of the \(\theta\) vector and the value of the objective at that \(\theta\).
m01.optsum.fitlogTable with 2 columns and 72 rows:
θ objective
┌────────────────────────────────────────────
1 │ [1.0, 0.0, 1.0] 1784.64
2 │ [2.0, 0.0, 1.0] 1792.86
3 │ [1.0, 1.0, 1.0] 1799.0
4 │ [1.0, 0.0, 2.0] 1808.55
5 │ [0.0, 0.0, 1.0] 1806.24
6 │ [1.0, -1.0, 1.0] 1798.6
7 │ [1.0, 0.0, 0.0] 1794.77
8 │ [1.22435, -0.00691503, 0.797624] 1778.03
9 │ [1.46448, -0.0169209, 0.611336] 1771.86
10 │ [1.80499, -0.0631689, 0.109571] 1774.69
11 │ [1.46736, 0.101338, 0.607374] 1772.35
12 │ [1.56923, -0.151459, 0.422853] 1763.91
13 │ [1.66673, -0.627902, 0.274894] 1783.3
14 │ [1.65242, -0.14957, 0.478319] 1767.42
15 │ [1.49341, -0.0708007, 0.32143] 1757.94
16 │ [1.51395, 0.228361, 0.337085] 1763.58
17 │ [1.40516, -0.106856, 0.291229] 1756.52
⋮ │ ⋮ ⋮There were 72 evaluations of the objective before convergence was declared, according to rather stringent convergence criteria. We can see that the last 10 evaluations only produced changes in the fourth decimal place of the objective or even smaller. That is, effective convergence occurred a fair bit sooner.
4 A model that converges to a degenerate distribution
When the iterations converge to a value of \(\theta\) that is on the boundary, which means that one of the diagonal elements of \(\lambda\) is exactly zero, we say that \(\lambda\) is singular, in the sense that its inverse does not exist.
Another, more evocative characterization, is that the distribution of the random effects is a degenerate distribution, in the sense that the values of the random effects are constrained to a plane (technically, a linear subspace) in the space of all possible random effects.
This is shown in section 3.5.3 of Embrace Uncertainty, especially Figure 3.13.
5 Appendix
These are some appendices, mostly so that I can keep track of them
5.1 Evaluating the random effects correlation from \(\theta\)
There is a short-cut for evaluating the correlation which is to “normalize” the second row of \(\lambda\), in the sense that the row is scaled so that it has unit length.
normed = normalize!(λ[2, :])2-element Vector{Float64}:
0.08131378907082806
0.9966885510062535providing the correlation as
first(normed)0.081313789070828065.2 Optimizing with a fixed correlation
To profile the correlation we need optimize the objective while fixing a value of the correlation. The way we will do this is to determine \(\theta_2\) as a function of the fixed \(\rho\) and \(\theta_3\).
We need to solve \[ \rho = \frac{\theta_2}{\sqrt{\theta_2^2 + \theta_3^2}} \]
for \(\theta_2\) as a function of \(\rho\) and \(\theta_3\).
Notice that \(\theta_2\) and \(\rho\) have the same sign. Thus it is sufficient to determine the absolute value of \(\theta_2\) then transfer the sign from \(\rho\). \[ \theta_2^2=\rho^2(\theta_2^2 + \theta_3^2) \] which implies \[ \theta_2^2 = \frac{\rho^2}{1-\rho^2}\theta_3^2, \quad \theta_3\ge 0 \] and thus \[ \theta_2=\frac{\rho}{\sqrt{1-\rho^2}}\theta_3 \]
5.3 Exercises
- Provoke a singular fit. Fit the random-slope model to a small subset of
sleepstudy(say the first six subjects). Is the fit singular? Inspect the estimated correlation parameter.
NoteSolution
With little data the model often cannot estimate the intercept–slope correlation away from the boundary, so issingular returns true and the correlation sits at ±1.
df = DataFrame(dataset(:sleepstudy))
sub = filter(:subj => in(first(unique(df.subj), 3)), df)
m = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days | subj)), sub)
issingular(m)- Interpretation. A converged model reports an intercept–slope correlation of exactly +1. What does that tell you about the estimated random-effects covariance, and what is a reasonable next step?
NoteSolution
A correlation of ±1 means the estimated random-effects covariance is rank-deficient (degenerate) — the random intercept and slope are perfectly collinear, so the data do not support two independent random-effect dimensions. A reasonable response is to simplify the random-effects structure (e.g. a zero-correlation parametrization zerocorr(...) or dropping the random slope) rather than treating the correlation as a substantive estimate. We will discuss this more in subsequent lessons.
This page was rendered from git revision e563d55 using Quarto 1.9.38 and Julia 1.12.6.