Shrinkage plots

Authors

Douglas Bates

Phillip Alday

Published

2026-08-23

After working through this page you will be able to:

read a shrinkage plot and explain what its two sets of points represent;
connect the amount of shrinkage to the strength of a variance component;
use shrinkage plots to judge whether a random-effects term is contributing or “inert”.

Before you start

Prerequisites: Analysis of the sleepstudy data; Convergence, singularity and all that is helpful.

Datasets used: kb07 (see the dataset catalog).

I have stated that the likelihood criterion used to fit linear mixed-effects can be considered as balancing fidelity to the data (i.e. fits the observed data well) versus model complexity.
This is similar to some of the criterion used in Machine Learning (ML), except that the criterion for LMMs has a rigorous mathematical basis.
In the shrinkage plot we consider the values of the random-effects coefficients for the fitted values of the model versus those from a model in which there is no penalty for model complexity.
If there is strong subject-to-subject variation then the model fit will tend to values of the random effects similar to those without a penalty on complexity.
If the random effects term is not contributing much (i.e. it is “inert”) then the random effects will be shrunk considerably towards zero in some directions.

Code

using CairoMakie
using DataFrames
using LinearAlgebra
using MixedModels
using MixedModelsDatasets: dataset
using MixedModelsMakie
using Random
using ProgressMeter

const progress = isinteractive()

Load the kb07 data set.

kb07 = dataset(:kb07)

Arrow.Table with 1789 rows, 7 columns, and schema:
 :subj      String
 :item      String
 :spkr      String
 :prec      String
 :load      String
 :rt_trunc  Int16
 :rt_raw    Int16

contrasts = Dict(
  :spkr => HelmertCoding(),
  :prec => HelmertCoding(),
  :load => HelmertCoding(),
)
m1 = let
  form = @formula(
    rt_trunc ~
      1 +
      spkr * prec * load +
      (1 + spkr + prec + load | subj) +
      (1 + spkr + prec + load | item)
  )
  fit(MixedModel, form, kb07; contrasts, progress)
end

VarCorr(m1)

	Column	Variance	Std.Dev	Corr.
subj	(Intercept)	91087.5446	301.8071
	spkr: old	1841.1736	42.9089	+0.78
	prec: maintain	3845.7309	62.0140	-0.59	+0.02
	load: yes	4243.2715	65.1404	+0.36	+0.83	+0.53
item	(Intercept)	131147.0985	362.1424
	spkr: old	1657.6012	40.7136	+0.44
	prec: maintain	60972.1927	246.9255	-0.69	+0.35
	load: yes	1794.5887	42.3626	+0.32	+0.16	-0.14
Residual		446891.5597	668.4995

issingular(m1)

true

print(m1)

Linear mixed model fit by maximum likelihood
 rt_trunc ~ 1 + spkr + prec + load + spkr & prec + spkr & load + prec & load + spkr & prec & load + (1 + spkr + prec + load | subj) + (1 + spkr + prec + load | item)
    logLik   -2 logLik      AIC         AICc        BIC     
 -14318.5614  28637.1227  28695.1227  28696.1119  28854.3157

Variance components:
             Column       Variance  Std.Dev.   Corr.
subj     (Intercept)      91087.5446 301.8071
         spkr: old         1841.1736  42.9089 +0.78
         prec: maintain    3845.7309  62.0140 -0.59 +0.02
         load: yes         4243.2715  65.1404 +0.36 +0.83 +0.53
item     (Intercept)     131147.0985 362.1424
         spkr: old         1657.6012  40.7136 +0.44
         prec: maintain   60972.1927 246.9255 -0.69 +0.35
         load: yes         1794.5887  42.3626 +0.32 +0.16 -0.14
Residual                 446891.5597 668.4995
 Number of obs: 1789; levels of grouping factors: 56, 32

  Fixed-effects parameters:
──────────────────────────────────────────────────────────────────────────────
                                            Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────────────────────────────────
(Intercept)                             2181.67       77.2965  28.22    <1e-99
spkr: old                                 67.7486     18.2896   3.70    0.0002
prec: maintain                          -333.921      47.1581  -7.08    <1e-11
load: yes                                 78.7704     19.5368   4.03    <1e-04
spkr: old & prec: maintain               -21.9657     15.8061  -1.39    0.1646
spkr: old & load: yes                     18.3843     15.8061   1.16    0.2448
prec: maintain & load: yes                 4.5338     15.8061   0.29    0.7742
spkr: old & prec: maintain & load: yes    23.6075     15.8061   1.49    0.1353
──────────────────────────────────────────────────────────────────────────────

1 Expressing the covariance of random effects

Earlier today we mentioned that the parameters being optimized are from a “matrix square root” of the covariance matrix for the random effects. There is one such lower triangular matrix for each grouping factor.

l1 = first(m1.λ)   # Cholesky factor of relative covariance for subj

4×4 LowerTriangular{Float64, Matrix{Float64}}:
  0.45147     ⋅          ⋅          ⋅ 
  0.0502257  0.0399667   ⋅          ⋅ 
 -0.0550789  0.0728357  0.0163335   ⋅ 
  0.0351887  0.0855018  0.0307616  0.0

Notice the zero on the diagonal. A triangular matrix with zeros on the diagonal is singular.

l2 = last(m1.λ)    # this one is also singular

4×4 LowerTriangular{Float64, Matrix{Float64}}:
  0.541724    ⋅           ⋅          ⋅ 
  0.0269111  0.0546349    ⋅          ⋅ 
 -0.253071   0.268066    0.0230701   ⋅ 
  0.0200495  0.00126769  0.0600886  0.00121966

To regenerate the covariance matrix we need to know that the covariance is not the square of l1, it is l1 * l1' (so that the result is symmetric) and multiplied by σ̂²

Σ₁ = varest(m1) .* (l1 * l1')

4×4 Matrix{Float64}:
  91087.5   10133.4     -11112.6     7099.61
  10133.4    1841.17        64.6309  2316.95
 -11112.6      64.6309    3845.73    2141.45
   7099.61   2316.95      2141.45    4243.27

diag(Σ₁)  # compare to the variance column in the VarCorr output

4-element Vector{Float64}:
 91087.54462745202
  1841.1736312244882
  3845.730881202395
  4243.271473369348

sqrt.(diag(Σ₁))

4-element Vector{Float64}:
 301.807131505291
  42.90889920779241
  62.013957148390354
  65.14039816710786

2 Shrinkage plots

shrinkageplot!(Figure(; size=(500, 500)), m1)

Figure 1: Shrinkage plot of model m1 for the `:subj` grouping term

The upper left panel shows the perfect negative correlation for those two components of the random effects.

shrinkageplot!(Figure(; size=(500, 500)), m1, :item)

Figure 2: Shrinkage plot of model m1 for the `:item` grouping term

X1 = Int.(m1.X')

8×1789 Matrix{Int64}:
  1   1   1   1   1  1   1   1   1   1  …   1   1   1   1   1   1   1  1   1
 -1   1   1  -1  -1  1   1  -1  -1   1      1  -1  -1   1   1  -1  -1  1   1
 -1   1  -1   1  -1  1  -1   1  -1   1     -1   1  -1   1  -1   1  -1  1  -1
  1  -1  -1  -1  -1  1   1   1   1  -1      1   1   1  -1  -1  -1  -1  1   1
  1   1  -1  -1   1  1  -1  -1   1   1     -1  -1   1   1  -1  -1   1  1  -1
 -1  -1  -1   1   1  1   1  -1  -1  -1  …   1  -1  -1  -1  -1   1   1  1   1
 -1  -1   1  -1   1  1  -1   1  -1  -1     -1   1  -1  -1   1  -1   1  1  -1
  1  -1   1   1  -1  1  -1  -1   1  -1     -1  -1   1  -1   1   1  -1  1  -1

X1 * X1'

8×8 Matrix{Int64}:
 1789    -1    -1     3    -3     1     1     3
   -1  1789    -3     1    -1     3     3     1
   -1    -3  1789     1    -1     3     3     1
    3     1     1  1789     3    -1    -1    -3
   -3    -1    -1     3  1789     1     1     3
    1     3     3    -1     1  1789    -3    -1
    1     3     3    -1     1    -3  1789    -1
    3     1     1    -3     3    -1    -1  1789

3 How to interpret a shrinkage plot

Extreme shrinkage (shrunk to a line or to a point) is easy to interpret – the term is not providing benefit and can be removed.
When the range of the blue dots (shrunk values) is comparable to those of the red dots (unshrunk) it indicates that the term after shrinkage is about as strong as without shrinkage.
By itself, this doesn’t mean that the term is important. In some ways you need to get a feeling for the absolute magnitude of the random effects in addition to the relative magnitude.
Small magnitude and small relative magnitude indicate you can drop that term

4 Conclusions from these plots

Only the intercept for the subj appears to be contributing explanatory power
For the item both the intercept and the spkr appear to be contributing

m2 = let
  form = @formula(
    rt_trunc ~
      1 + prec * spkr * load + (1 | subj) + (1 + prec | item)
  )
  fit(MixedModel, form, kb07; contrasts, progress)
end

VarCorr(m2)

	Column	Variance	Std.Dev	Corr.
item	(Intercept)	133027.436	364.729
	prec: maintain	63841.834	252.669	-0.70
subj	(Intercept)	88870.013	298.111
Residual		460948.573	678.932

Code

shrinkageplot!(Figure(; size=(250, 250)), m2)

m3 = let
  form = @formula(
    rt_trunc ~
      1 + prec + spkr + load + (1 | subj) + (1 + prec | item)
  )
  fit(MixedModel, form, kb07; contrasts, progress)
end

VarCorr(m3)

	Column	Variance	Std.Dev	Corr.
item	(Intercept)	133014.904	364.712
	prec: maintain	63768.184	252.524	-0.70
subj	(Intercept)	88819.302	298.026
Residual		462443.261	680.032

rng = Random.seed!(1234321);

m3btstrp = parametricbootstrap(rng, 2000, m3);

DataFrame(shortestcovint(m3btstrp))

9×5 DataFrame

Row	type	group	names	lower	upper
	String	String?	String?	Float64	Float64
1	β	missing	(Intercept)	2023.0	2334.0
2	β	missing	prec: maintain	-430.22	-239.746
3	β	missing	spkr: old	33.9812	96.6883
4	β	missing	load: yes	46.4854	109.517
5	σ	item	(Intercept)	270.075	452.086
6	σ	item	prec: maintain	182.165	325.334
7	ρ	item	(Intercept), prec: maintain	-0.907612	-0.489103
8	σ	subj	(Intercept)	227.994	355.589
9	σ	residual	missing	657.347	702.618

ridgeplot(m3btstrp)

Figure 4: Ridge plot of the fixed-effects coefficients from the bootstrap sample

ridgeplot(m3btstrp; show_intercept=false)

Figure 5: Ridge plot of the fixed-effects coefficients from the bootstrap sample (without the intercept)

m4 = let
  form = @formula(
    rt_trunc ~
      1 + prec + spkr + load + (1 + prec | item) + (1 | subj)
  )
  fit(MixedModel, form, kb07; contrasts, progress)
end

m4bstrp = parametricbootstrap(rng, 2000, m4);

ridgeplot(m4bstrp; show_intercept=false)

DataFrame(shortestcovint(m4bstrp))

9×5 DataFrame

Row	type	group	names	lower	upper
	String	String?	String?	Float64	Float64
1	β	missing	(Intercept)	2034.24	2335.44
2	β	missing	prec: maintain	-428.711	-249.325
3	β	missing	spkr: old	35.803	97.9178
4	β	missing	load: yes	45.2571	107.383
5	σ	item	(Intercept)	260.507	451.799
6	σ	item	prec: maintain	178.443	315.113
7	ρ	item	(Intercept), prec: maintain	-0.900336	-0.47248
8	σ	subj	(Intercept)	232.089	357.732
9	σ	residual	missing	657.18	702.755

VarCorr(m4)

	Column	Variance	Std.Dev	Corr.
item	(Intercept)	133014.904	364.712
	prec: maintain	63768.184	252.524	-0.70
subj	(Intercept)	88819.302	298.026
Residual		462443.261	680.032

Code

let mods = [m1, m2, m4]
  DataFrame(;
    geomdof=(sum ∘ leverage).(mods),
    npar=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    BIC=bic.(mods),
    AICc=aicc.(mods),
  )
end

3×6 DataFrame

Row	geomdof	npar	deviance	AIC	BIC	AICc
	Float64	Int64	Float64	Float64	Float64	Float64
1	131.565	29	28637.1	28695.1	28854.3	28696.1
2	107.543	13	28658.5	28684.5	28755.8	28684.7
3	103.478	9	28663.9	28681.9	28731.3	28682.0

fig = Figure(; size=(400, 400))
ax = fig[1, 1] = Axis(fig)
scatter!(ax, fitted(m4), residuals(m4); alpha=0.5)
fig

Figure 6: Residuals versus fitted values for model m4

5 See Also

The introductory and longitudinal-data chapters of Embrace Uncertainty — the theory of conditional modes and further worked examples of shrinkage.

6 Exercises

Shrinkage in a well-supported model. Produce a shrinkage plot for the sleepstudy random-slope model. How much do the points move compared with the kb07 plots above, and what does that say about the strength of the random effects?

Solution

The sleepstudy random effects are well supported by the data, so the conditional modes shrink only modestly toward the origin — the points barely move. Strong subject-to-subject variation means the penalty for complexity has little to pull against.

df = DataFrame(dataset(:sleepstudy))
m = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days | subj)), df)
shrinkageplot!(Figure(; size=(500, 500)), m)

Reading shrinkage. In a different model, the points for one random-effects term collapse almost entirely onto a single line. What is that telling you about the term?

Solution

Heavy shrinkage toward a lower-dimensional subspace (a line, or the origin) means that random-effects term is contributing little — it is close to “inert”. The likelihood gains little from letting it vary, so the estimates are pulled strongly toward zero. This is a visual cue to consider simplifying the random-effects structure.

This page was rendered from git revision 29e8d33 using Quarto 1.10.18 and Julia 1.12.7.