Parametric bootstrap for mixed-effects models

Author

Phillip Alday, Douglas Bates, and Reinhold Kliegl

Published

2025-05-12

The speed of MixedModels.jl relative to its predecessors makes the parametric bootstrap much more computationally tractable. This is valuable because the parametric bootstrap can be used to produce more accurate confidence intervals than methods based on standard errors or profiling of the likelihood surface.

This page is adapted from the MixedModels.jl docs

1 The parametric bootstrap

Bootstrapping is a family of procedures for generating sample values of a statistic, allowing for visualization of the distribution of the statistic or for inference from this sample of values. Bootstrapping also belongs to a larger family of procedures called resampling, which are based on creating new samples of data from an existing one, then computing statistics on the new samples, in order to examine the distribution of the relevant statistics.

A parametric bootstrap is used with a parametric model, m, that has been fit to data. The procedure is to simulate n response vectors from m using the estimated parameter values and refit m to these responses in turn, accumulating the statistics of interest at each iteration.

The parameters of a LinearMixedModel object are the fixed-effects parameters (β), the standard deviation (σ), of the per-observation noise, and the covariance parameter (θ), that defines the variance-covariance matrices of the random effects. A technical description of the covariance parameter can be found in the MixedModels.jl docs. Lisa Schwetlick and Daniel Backhaus have provided a more beginner-friendly description of the covariance parameter in the documentation for MixedModelsSim.jl. For today’s purposes – looking at the uncertainty in the estimates from a fitted model – we can simply use values from the fitted model, but we will revisit the parametric bootstrap as a convenient way to simulate new data, potentially with different parameter values, for power analysis.

Attach the packages to be used

Code
using AlgebraOfGraphics
using CairoMakie
using DataFrames
using MixedModels
using MixedModelsMakie
using Random
using SMLP2025: dataset

using AlgebraOfGraphics: AlgebraOfGraphics as AoG
const progress=false

Note that the precise stream of random numbers generated for a given seed can change between Julia versions. For exact reproducibility, you either need to have the exact same Julia version or use the StableRNGs package.

2 A model of moderate complexity

The kb07 data (Kronmüller & Barr, 2007) are one of the datasets provided by the MixedModels package.

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

Convert the table to a DataFrame for summary.

kb07 = DataFrame(kb07)
describe(kb07)
7×7 DataFrame
Row variable mean min median max nmissing eltype
Symbol Union… Any Union… Any Int64 DataType
1 subj S030 S103 0 String
2 item I01 I32 0 String
3 spkr new old 0 String
4 prec break maintain 0 String
5 load no yes 0 String
6 rt_trunc 2182.2 579 1940.0 5171 0 Int16
7 rt_raw 2226.24 579 1940.0 15923 0 Int16

The experimental factors; spkr, prec, and load, are two-level factors.

contrasts = Dict(:spkr => EffectsCoding(),
                 :prec => EffectsCoding(),
                 :load => EffectsCoding(),
)

The EffectsCoding contrast is used with these to create a ±1 encoding.

We can look at an initial fit of moderate complexity:

form = @formula(rt_trunc ~ 1 + spkr * prec * load +
                          (1 + spkr + prec + load | subj) +
                          (1 + spkr + prec + load | item))
m0 = fit(MixedModel, form, kb07; contrasts, progress)
Est. SE z p σ_subj σ_item
(Intercept) 2181.6729 77.3131 28.22 <1e-99 301.6736 362.3182
spkr: old 67.7490 18.2747 3.71 0.0002 42.5441 40.7073
prec: maintain -333.9206 47.1388 -7.08 <1e-11 61.9978 246.8078
load: yes 78.7702 19.5387 4.03 <1e-04 64.9874 42.5139
spkr: old & prec: maintain -21.9655 15.8070 -1.39 0.1646
spkr: old & load: yes 18.3838 15.8070 1.16 0.2448
prec: maintain & load: yes 4.5333 15.8070 0.29 0.7743
spkr: old & prec: maintain & load: yes 23.6072 15.8070 1.49 0.1353
Residual 668.5371

The default display in Quarto uses the pretty MIME show method for the model and omits the estimated correlations of the random effects.

The VarCorr extractor displays these.

VarCorr(m0)
Column Variance Std.Dev Corr.
subj (Intercept) 91006.9887 301.6736
spkr: old 1810.0012 42.5441 +0.79
prec: maintain 3843.7268 61.9978 -0.59 +0.02
load: yes 4223.3647 64.9874 +0.36 +0.85 +0.54
item (Intercept) 131274.4558 362.3182
spkr: old 1657.0830 40.7073 +0.44
prec: maintain 60914.1141 246.8078 -0.69 +0.35
load: yes 1807.4303 42.5139 +0.32 +0.15 -0.14
Residual 446941.7901 668.5371

None of the two-factor or three-factor interaction terms in the fixed-effects are significant. In the random-effects terms only the scalar random effects and the prec random effect for item appear to be warranted, leading to the reduced formula

# formula f4 from https://doi.org/10.33016/nextjournal.100002
form = @formula(rt_trunc ~ 1 + spkr * prec * load + (1 | subj) + (1 + prec | item))

m1 = fit(MixedModel, form, kb07; contrasts, progress)
Est. SE z p σ_item σ_subj
(Intercept) 2181.7582 77.4709 28.16 <1e-99 364.7286 298.1109
spkr: old 67.8114 16.0526 4.22 <1e-04
prec: maintain -333.8582 47.4629 -7.03 <1e-11 252.6687
load: yes 78.6849 16.0525 4.90 <1e-06
spkr: old & prec: maintain -21.8802 16.0525 -1.36 0.1729
spkr: old & load: yes 18.3214 16.0526 1.14 0.2537
prec: maintain & load: yes 4.4710 16.0526 0.28 0.7806
spkr: old & prec: maintain & load: yes 23.5219 16.0525 1.47 0.1428
Residual 678.9318
VarCorr(m1)
Column Variance Std.Dev Corr.
item (Intercept) 133026.918 364.729
prec: maintain 63841.496 252.669 -0.70
subj (Intercept) 88870.080 298.111
Residual 460948.432 678.932

These two models are nested and can be compared with a likelihood-ratio test.

MixedModels.likelihoodratiotest(m0, m1)
model-dof deviance χ² χ²-dof P(>χ²)
rt_trunc ~ 1 + spkr + prec + load + spkr & prec + spkr & load + prec & load + spkr & prec & load + (1 | subj) + (1 + prec | item) 13 28658
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) 29 28637 21 16 0.1650

The p-value of approximately 20% leads us to prefer the simpler model, m1, to the more complex, m0.

3 Bootstrap basics

To bootstrap the model parameters, first initialize a random number generator then create a bootstrap sample and extract the table of parameter estimates from it.

const RNG = MersenneTwister(42)
samp = parametricbootstrap(RNG, 5_000, m1)
tbl = samp.tbl
Table with 18 columns and 5000 rows:
      obj      β1       β2       β3        β4       β5        β6        ⋯
    ┌────────────────────────────────────────────────────────────────────
 1  │ 28666.6  2173.04  55.6593  -359.246  87.1657  -35.2652  32.9777   ⋯
 2  │ 28709.1  1973.88  88.8108  -313.755  83.8083  -23.2173  -15.9147  ⋯
 3  │ 28575.8  2167.28  48.2509  -297.175  79.0748  -22.07    41.8946   ⋯
 4  │ 28721.3  2259.36  65.9821  -386.726  82.7248  -35.3908  14.3295   ⋯
 5  │ 28590.9  2182.95  45.0315  -410.478  109.293  -19.132   42.9369   ⋯
 6  │ 28518.9  2167.98  61.0523  -360.56   75.7268  15.1612   61.9125   ⋯
 7  │ 28702.9  2210.25  70.8323  -349.703  99.4985  -32.9348  16.797    ⋯
 8  │ 28726.1  2235.04  59.248   -353.737  58.0711  -35.2804  10.0765   ⋯
 9  │ 28734.7  2042.97  72.1467  -300.829  71.2183  -14.2944  10.7301   ⋯
 10 │ 28609.2  2011.1   76.404   -254.606  63.9238  -37.8758  3.95621   ⋯
 11 │ 28636.7  2166.35  71.1184  -324.401  58.5484  -21.7266  -31.4973  ⋯
 12 │ 28659.2  2158.36  64.3944  -336.452  80.3281  -39.3608  1.96649   ⋯
 13 │ 28637.3  2147.09  70.6536  -399.734  60.0463  -35.6682  28.4599   ⋯
 14 │ 28672.7  2138.23  75.4572  -262.36   96.9987  -32.0501  26.9257   ⋯
 15 │ 28662.5  2239.34  65.1886  -373.371  83.4639  -29.7632  5.43283   ⋯
 16 │ 28604.9  2088.54  62.0142  -268.106  99.684   -27.3686  49.5337   ⋯
 17 │ 28539.6  2313.47  61.1166  -337.951  74.769   -14.5982  48.7173   ⋯
 ⋮  │    ⋮        ⋮        ⋮        ⋮         ⋮        ⋮         ⋮      ⋱

An empirical density plot of the estimates of the residual standard deviation is obtained as

plt = data(tbl) * mapping(:σ) * AoG.density()
draw(plt; axis=(;title="Parametric bootstrap estimates of σ"))

A density plot of the estimates of the standard deviation of the random effects is obtained as

plt = data(tbl) * mapping(
  [:σ1, :σ2, :σ3] .=> "Bootstrap replicates of standard deviations";
  color=dims(1) => renamer(["Item intercept", "Item speaker", "Subj"])
) * AoG.density()
draw(plt; figure=(;supertitle="Parametric bootstrap estimates of variance components"))

The bootstrap sample can be used to generate intervals that cover a certain percentage of the bootstrapped values. We refer to these as “coverage intervals”, similar to a confidence interval. The shortest such intervals, obtained with the confint extractor, correspond to a highest posterior density interval in Bayesian inference.

We generate these for all random and fixed effects:

confint(samp; method=:shortest)
DictTable with 2 columns and 13 rows:
 par   lower      upper
 ────┬─────────────────────
 β1  │ 2033.36    2332.44
 β2  │ 34.9886    98.421
 β3  │ -430.861   -246.169
 β4  │ 48.9002    112.424
 β5  │ -51.4927   12.1684
 β6  │ -13.6705   49.681
 β7  │ -27.2138   34.6244
 β8  │ -6.36896   55.7056
 ρ1  │ -0.899979  -0.460383
 σ   │ 654.51     700.692
 σ1  │ 264.114    448.733
 σ2  │ 174.044    316.132
 σ3  │ 228.265    359.429
draw(
  data(samp.β) * mapping(:β; color=:coefname) * AoG.density();
  figure=(; resolution=(800, 450)),
)
Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
@ Makie ~/.julia/packages/Makie/KcEWO/src/scenes.jl:238

For the fixed effects, MixedModelsMakie provides a convenience interface to plot the combined coverage intervals and density plots

ridgeplot(samp)

Often the intercept will be on a different scale and potentially less interesting, so we can stop it from being included in the plot:

ridgeplot(samp; show_intercept=false)

4 Singularity

Let’s consider the classic dysetuff dataset:

dyestuff = dataset(:dyestuff)
mdye = fit(MixedModel, @formula(yield ~ 1 + (1 | batch)), dyestuff)
Est. SE z p σ_batch
(Intercept) 1527.5000 17.6946 86.33 <1e-99 37.2603
Residual 49.5101
sampdye = parametricbootstrap(MersenneTwister(1234321), 10_000, mdye)
tbldye = sampdye.tbl
Table with 5 columns and 10000 rows:
      obj      β1       σ        σ1       θ1
    ┌─────────────────────────────────────────────
 1  │ 316.738  1528.14  43.5735  22.6691  0.520249
 2  │ 336.101  1552.01  58.1535  39.5166  0.679521
 3  │ 322.046  1501.55  51.857   0.0      0.0
 4  │ 326.893  1525.24  43.6512  66.3744  1.52056
 5  │ 331.544  1522.05  50.7415  51.0074  1.00524
 6  │ 326.892  1550.23  53.5174  19.0928  0.356759
 7  │ 313.786  1504.16  45.1878  0.0      0.0
 8  │ 322.636  1494.05  49.8556  17.7732  0.356494
 9  │ 323.582  1529.74  46.4678  35.2094  0.757716
 10 │ 321.639  1532.6   45.3573  32.4904  0.71632
 11 │ 331.565  1548.25  55.8562  28.7512  0.514736
 12 │ 310.076  1481.8   34.9418  38.4381  1.10006
 13 │ 314.492  1541.95  40.3751  28.3694  0.702646
 14 │ 340.876  1530.48  63.006   42.6409  0.676776
 15 │ 309.778  1496.17  41.5411  8.08346  0.19459
 16 │ 335.375  1530.32  57.4864  38.9066  0.676797
 17 │ 318.347  1511.28  41.0926  39.1128  0.951821
 ⋮  │    ⋮        ⋮        ⋮        ⋮        ⋮
plt = data(tbldye) * mapping(:σ1) * AoG.density()
draw(plt; axis=(;title="Parametric bootstrap estimates of σ_batch"))

Notice that this density plot has a spike, or mode, at zero. Although this mode appears to be diffuse, this is an artifact of the way that density plots are created. In fact, it is a pulse, as can be seen from a histogram.

plt = data(tbldye) * mapping(:σ1) * AoG.histogram(;bins=100)
draw(plt; axis=(;title="Parametric bootstrap estimates of σ_batch"))

A value of zero for the standard deviation of the random effects is an example of a singular covariance. It is easy to detect the singularity in the case of a scalar random-effects term. However, it is not as straightforward to detect singularity in vector-valued random-effects terms.

For example, if we bootstrap a model fit to the sleepstudy data

sleepstudy = dataset(:sleepstudy)
msleep = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days | subj)),
             sleepstudy)
Est. SE z p σ_subj
(Intercept) 251.4051 6.6323 37.91 <1e-99 23.7805
days 10.4673 1.5022 6.97 <1e-11 5.7168
Residual 25.5918
sampsleep = parametricbootstrap(MersenneTwister(666), 10_000, msleep)
tblsleep = sampsleep.tbl
Table with 10 columns and 10000 rows:
      obj      β1       β2       σ        σ1       σ2       ρ1           ⋯
    ┌─────────────────────────────────────────────────────────────────────
 1  │ 1694.07  258.109  10.7928  22.1713  13.6388  5.42872  0.430499     ⋯
 2  │ 1731.59  241.267  9.77434  23.5943  23.882   7.07979  -0.277561    ⋯
 3  │ 1706.51  254.403  10.6731  23.0742  21.4346  3.78782  0.493345     ⋯
 4  │ 1737.33  257.833  7.61745  24.5028  27.2647  4.8567   -0.0190962   ⋯
 5  │ 1764.5   244.154  13.0045  25.8322  30.9873  7.78254  -0.532286    ⋯
 6  │ 1730.85  247.672  13.2977  24.3584  18.6105  5.87738  -0.0535883   ⋯
 7  │ 1709.66  256.847  8.60197  22.4272  26.8884  6.62183  -0.638036    ⋯
 8  │ 1776.49  252.028  9.56225  27.3708  27.7002  6.53333  -0.353783    ⋯
 9  │ 1732.12  248.687  11.5173  23.7153  22.2865  6.69781  0.304236     ⋯
 10 │ 1748.89  249.387  10.2581  25.3961  24.7505  5.44654  0.00471327   ⋯
 11 │ 1735.34  247.16   8.57302  25.593   19.2151  3.99188  1.0          ⋯
 12 │ 1757.34  248.074  9.76959  25.9834  29.2859  4.74129  0.216426     ⋯
 13 │ 1727.32  248.007  12.8942  23.9314  32.1774  3.50451  0.152022     ⋯
 14 │ 1753.38  252.684  9.88845  25.2124  21.3446  7.64867  -0.00704067  ⋯
 15 │ 1731.05  262.218  7.70815  25.0737  18.1155  4.89347  -0.280397    ⋯
 16 │ 1742.68  252.273  12.7376  24.4489  20.6843  7.47369  0.0242807    ⋯
 17 │ 1727.37  243.956  12.107   24.6337  26.502   3.89643  -0.411438    ⋯
 ⋮  │    ⋮        ⋮        ⋮        ⋮        ⋮        ⋮          ⋮       ⋱

the singularity can be exhibited as a standard deviation of zero or as a correlation of ±1.

confint(sampsleep)
DictTable with 2 columns and 6 rows:
 par   lower      upper
 ────┬───────────────────
 β1  │ 238.265    264.28
 β2  │ 7.56525    13.4072
 ρ1  │ -0.409958  1.0
 σ   │ 22.6673    28.5759
 σ1  │ 10.5606    33.2568
 σ2  │ 3.08374    7.73233

A histogram of the estimated correlations from the bootstrap sample has a spike at +1.

plt = data(tblsleep) * mapping(:ρ1) * AoG.histogram(;bins=100)
draw(plt; axis=(;title="Parametric bootstrap samples of correlation of random effects"))

or, as a count,

count(tblsleep.ρ1  .≈ 1)
309

Close examination of the histogram shows a few values of -1.

count(tblsleep.ρ1 .≈ -1)
0

Furthermore there are even a few cases where the estimate of the standard deviation of the random effect for the intercept is zero.

count(tblsleep.σ1 .≈ 0)
4

There is a general condition to check for singularity of an estimated covariance matrix or matrices in a bootstrap sample. The parameter optimized in the estimation is θ, the relative covariance parameter. Some of the elements of this parameter vector must be non-negative and, when one of these components is approximately zero, one of the covariance matrices will be singular.

The issingular method for a MixedModel object that tests if a parameter vector θ corresponds to a boundary or singular fit.

This operation is encapsulated in a method for the issingular function that works on MixedModelBootstrap objects.

count(issingular(sampsleep))
313

5 References

Kronmüller, E., & Barr, D. J. (2007). Perspective-free pragmatics: Broken precedents and the recovery-from-preemption hypothesis. Journal of Memory and Language, 56(3), 436–455. https://doi.org/10.1016/j.jml.2006.05.002

This page was rendered from git revision 0e4cb20 using Quarto 1.7.29.

Back to top