Code
using AlgebraOfGraphics
using CairoMakie
using DataFrames
using MixedModels
using MixedModelsMakie
using Random
using SMLP2025: dataset
using AlgebraOfGraphics: AlgebraOfGraphics as AoG
const progress=false
Phillip Alday, Douglas Bates, and Reinhold Kliegl
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
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
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.
The kb07
data (Kronmüller & Barr, 2007) are one of the datasets provided by the MixedModels
package.
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.
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.
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.
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 |
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.
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
.
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.
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:
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
Often the intercept will be on a different scale and potentially less interesting, so we can stop it from being included in the plot:
Let’s consider the classic dysetuff dataset:
Est. | SE | z | p | σ_batch | |
(Intercept) | 1527.5000 | 17.6946 | 86.33 | <1e-99 | 37.2603 |
Residual | 49.5101 |
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 |
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.
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,
Close examination of the histogram shows a few values of -1
.
Furthermore there are even a few cases where the estimate of the standard deviation of the random effect for the intercept is zero.
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.
This page was rendered from git revision 0e4cb20 using Quarto 1.7.29.