Shrinkage plots

Authors

Douglas Bates

Phillip Alday

Published

2026-06-21

After working through this page you will be able to:

read a shrinkage plot and explain what its two sets of random-effects coefficients 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 SMLP2026: dataset
using MixedModelsMakie
using Random
using ProgressMeter

const progress = isinteractive()

false

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)	91064.7403	301.7693
	spkr: old	1850.5250	43.0177	+0.78
	prec: maintain	3843.8825	61.9991	-0.59	+0.03
	load: yes	4246.0532	65.1617	+0.36	+0.83	+0.53
item	(Intercept)	131173.0381	362.1782
	spkr: old	1656.7129	40.7027	+0.44
	prec: maintain	60946.6098	246.8737	-0.69	+0.35
	load: yes	1795.0255	42.3677	+0.32	+0.16	-0.14
Residual		446886.0133	668.4953

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.5613  28637.1227  28695.1227  28696.1119  28854.3156

Variance components:
             Column       Variance  Std.Dev.   Corr.
subj     (Intercept)      91064.7403 301.7693
         spkr: old         1850.5250  43.0177 +0.78
         prec: maintain    3843.8825  61.9991 -0.59 +0.03
         load: yes         4246.0532  65.1617 +0.36 +0.83 +0.53
item     (Intercept)     131173.0381 362.1782
         spkr: old         1656.7129  40.7027 +0.44
         prec: maintain   60946.6098 246.8737 -0.69 +0.35
         load: yes         1795.0255  42.3677 +0.32 +0.16 -0.14
Residual                 446886.0133 668.4953
 Number of obs: 1789; levels of grouping factors: 56, 32

  Fixed-effects parameters:
───────────────────────────────────────────────────────────────────────────────
                                             Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────────────────────────────────
(Intercept)                             2181.67        77.299   28.22    <1e-99
spkr: old                                 67.7485      18.2934   3.70    0.0002
prec: maintain                          -333.921       47.1493  -7.08    <1e-11
load: yes                                 78.7703      19.5383   4.03    <1e-04
spkr: old & prec: maintain               -21.9657      15.806   -1.39    0.1646
spkr: old & load: yes                     18.3843      15.806    1.16    0.2448
prec: maintain & load: yes                 4.53389     15.806    0.29    0.7742
spkr: old & prec: maintain & load: yes    23.6074      15.806    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.451416    ⋅          ⋅          ⋅ 
  0.0502725  0.0401697   ⋅          ⋅ 
 -0.0550321  0.0726579  0.01714     ⋅ 
  0.0351643  0.0850281  0.0321733  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.541781    ⋅           ⋅          ⋅ 
  0.0268651  0.0546398    ⋅          ⋅ 
 -0.253039   0.268006    0.0229063   ⋅ 
  0.0200358  0.00129292  0.0601135  0.0

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}:
  91064.7   10141.5     -11101.7     7093.74
  10141.5    1850.52        67.9472  2316.37
 -11101.7      67.9472    3843.88    2142.48
   7093.74   2316.37      2142.48    4246.05

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

4-element Vector{Float64}:
 91064.74027228545
  1850.5249616692
  3843.882496563724
  4246.0532111504635

sqrt.(diag(Σ₁))

4-element Vector{Float64}:
 301.76934945796836
  43.017728457802136
  61.9990523844012
  65.16174653238251

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.438	364.729
	prec: maintain	63841.835	252.669	-0.70
subj	(Intercept)	88870.014	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.905	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)	2022.91	2334.0
2	β	missing	prec: maintain	-430.239	-239.802
3	β	missing	spkr: old	34.0592	96.72
4	β	missing	load: yes	46.5349	109.526
5	σ	item	(Intercept)	270.065	451.979
6	σ	item	prec: maintain	181.739	325.127
7	ρ	item	(Intercept), prec: maintain	-0.907255	-0.490094
8	σ	subj	(Intercept)	233.834	364.504
9	σ	residual	missing	657.341	702.655

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.15	2335.45
2	β	missing	prec: maintain	-427.769	-248.337
3	β	missing	spkr: old	35.6938	97.8081
4	β	missing	load: yes	45.3095	107.368
5	σ	item	(Intercept)	260.498	451.649
6	σ	item	prec: maintain	179.495	315.414
7	ρ	item	(Intercept), prec: maintain	-0.901831	-0.471356
8	σ	subj	(Intercept)	236.684	360.363
9	σ	residual	missing	657.178	702.739

VarCorr(m4)

	Column	Variance	Std.Dev	Corr.
item	(Intercept)	133014.905	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.599	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 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 e563d55 using Quarto 1.9.38 and Julia 1.12.6.