After working through this page you will be able to:
MixedModel;Prerequisites: Analysis of the sleepstudy data.
Datasets used: sleepstudy, kb07, dyestuff (see the dataset catalog).
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 in power simulation as a convenient way to simulate new data, potentially with different parameter values, for power analysis.
Attach the packages to be used
using AlgebraOfGraphics
using CairoMakie
using DataFrames
using MixedModels
using MixedModelsMakie
using Random
using MixedModelsDatasets: dataset
using SMLP2026: fit_or_restore
using AlgebraOfGraphics: AlgebraOfGraphics as AoG
const progress = isinteractive()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.
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 Int16Convert the table to a DataFrame for summary.
kb07 = DataFrame(kb07)
describe(kb07)| 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.6727 | 77.2965 | 28.22 | <1e-99 | 301.8071 | 362.1424 |
| spkr: old | 67.7486 | 18.2896 | 3.70 | 0.0002 | 42.9089 | 40.7136 |
| prec: maintain | -333.9211 | 47.1581 | -7.08 | <1e-11 | 62.0140 | 246.9255 |
| load: yes | 78.7704 | 19.5368 | 4.03 | <1e-04 | 65.1404 | 42.3626 |
| 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 | ||
| Residual | 668.4995 |
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) | 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 |
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)VarCorr(m1)| Column | Variance | Std.Dev | Corr. | |
|---|---|---|---|---|
| item | (Intercept) | 133027.438 | 364.729 | |
| prec: maintain | 63841.834 | 252.669 | -0.70 | |
| subj | (Intercept) | 88870.013 | 298.111 | |
| Residual | 460948.572 | 678.932 |
These two models are nested and can be compared with a likelihood-ratio test.
MixedModels.likelihoodratiotest(m0, m1)| model-dof | -2 logLik | χ² | χ²-dof | P(>χ²) | |
|---|---|---|---|---|---|
| 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 | |||
| rt_trunc ~ 1 + spkr + prec + load + spkr & prec + spkr & load + prec & load + spkr & prec & load + (1 | subj) + (1 + prec | item) | 13 | -28658 | 21 | -16 | 0.1649 |
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.
const RNG = MersenneTwister(42)
samp = parametricbootstrap(RNG, 5_000, m1; progress)
tbl = samp.tblTable with 18 columns and 5000 rows:
obj β1 β2 β3 β4 β5 β6 ⋯
┌────────────────────────────────────────────────────────────────────
1 │ 28666.1 2173.0 55.5993 -359.306 87.2052 -35.3047 33.0377 ⋯
2 │ 28710.1 1973.91 88.8248 -313.74 83.7709 -23.1799 -15.9286 ⋯
3 │ 28582.8 2167.34 48.2284 -297.197 79.0159 -22.0111 41.9171 ⋯
4 │ 28718.5 2259.43 66.0501 -386.658 82.6574 -35.3235 14.2615 ⋯
5 │ 28589.3 2182.94 44.9803 -410.528 109.299 -19.1385 42.9881 ⋯
6 │ 28517.0 2167.97 61.043 -360.569 75.7417 15.1462 61.9219 ⋯
7 │ 28708.8 2210.25 70.856 -349.679 99.5077 -32.944 16.7733 ⋯
8 │ 28731.5 2235.07 59.2516 -353.734 58.049 -35.2583 10.0729 ⋯
9 │ 28731.8 2043.11 72.2606 -300.714 71.0784 -14.1545 10.6162 ⋯
10 │ 28606.7 2011.05 76.431 -254.577 63.9764 -37.9284 3.92921 ⋯
11 │ 28643.4 2166.5 71.0285 -324.491 58.3968 -21.5751 -31.4075 ⋯
12 │ 28660.0 2158.28 64.38 -336.465 80.4072 -39.4399 1.98082 ⋯
13 │ 28641.5 2147.14 70.6071 -399.78 59.9985 -35.6204 28.5064 ⋯
14 │ 28682.0 2138.16 75.5171 -262.3 97.0663 -32.1177 26.8657 ⋯
15 │ 28668.1 2239.36 65.2487 -373.311 83.4486 -29.7478 5.37272 ⋯
16 │ 28605.0 2088.58 62.0174 -268.102 99.6468 -27.3314 49.5304 ⋯
17 │ 28531.1 2313.5 61.0588 -338.01 74.7445 -14.5737 48.7751 ⋯
⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱optsum_overrides: faster replicates at reduced precision
The keyword argument optsum_overrides passes additional arguments controlling the model fitting of each replicate, e.g. parametricbootstrap(RNG, 5_000, m1; optsum_overrides=(; ftol_rel=1e-8)). ftol_rel=1e-8 lowers the threshold for changes in the objective – more directly, the optimizer considers the model converged if the change in the deviance is less than \(10^{-8}\), which is a very small change, but larger than the default \(10^{-12}\). Because the majority of the optimization time is spent in final fine-tuning, changing this threshold can greatly speed up the fitting time at the cost of a small loss of quality in fit. For a stochastic process like the bootstrap, that change in quality just adds to the general noise, but that’s an acceptable tradeoff in order to get many more replicates. The MixedModels.jl documentation also includes a little note on the reduced precision bootstrap.
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.37 2332.46
β2 │ 35.0188 98.4562
β3 │ -430.963 -246.167
β4 │ 48.9126 112.362
β5 │ -51.5002 12.1719
β6 │ -13.6978 49.6548
β7 │ -27.2307 34.5997
β8 │ -6.39206 55.7104
ρ1 │ -0.898241 -0.458759
σ │ 654.495 700.684
σ1 │ 260.548 445.231
σ2 │ 172.313 313.725
σ3 │ 231.042 361.333draw(
data(samp.β) * mapping(:β; color=:coefname) * AoG.density();
figure=(; size=(800, 450)),
)
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)
Let’s consider the classic dyestuff dataset:
dyestuff = dataset(:dyestuff)
mdye = fit(MixedModel, @formula(yield ~ 1 + (1 | batch)), dyestuff)sampdye = parametricbootstrap(MersenneTwister(1234321), 10_000, mdye)
tbldye = sampdye.tblTable 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.52057
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.7733 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.08345 0.194589
16 │ 335.375 1530.32 57.4864 38.9066 0.676797
17 │ 318.347 1511.28 41.0926 39.1128 0.95182
⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮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)sampsleep = parametricbootstrap(MersenneTwister(666), 10_000, msleep)
tblsleep = sampsleep.tblTable with 10 columns and 10000 rows:
obj β1 β2 σ σ1 σ2 ρ1 ⋯
┌─────────────────────────────────────────────────────────────────────
1 │ 1694.07 258.109 10.7928 22.1713 13.639 5.42874 0.430476 ⋯
2 │ 1731.59 241.267 9.77437 23.5943 23.8822 7.07987 -0.277573 ⋯
3 │ 1706.51 254.403 10.6731 23.0742 21.4347 3.78776 0.493341 ⋯
4 │ 1737.33 257.833 7.61743 24.5027 27.2653 4.85669 -0.0191323 ⋯
5 │ 1764.5 244.154 13.0045 25.8322 30.9872 7.78257 -0.532291 ⋯
6 │ 1730.85 247.672 13.2977 24.3583 18.6133 5.87767 -0.0542254 ⋯
7 │ 1709.66 256.847 8.60196 22.4272 26.8869 6.62172 -0.637998 ⋯
8 │ 1776.49 252.028 9.56227 27.3707 27.7005 6.53338 -0.353799 ⋯
9 │ 1732.12 248.687 11.5173 23.7153 22.2867 6.69778 0.304229 ⋯
10 │ 1748.89 249.387 10.2581 25.3961 24.7506 5.44657 0.00470547 ⋯
11 │ 1735.34 247.16 8.57303 25.5931 19.2141 3.9918 1.0 ⋯
12 │ 1757.34 248.074 9.7696 25.9833 29.2858 4.74155 0.216293 ⋯
13 │ 1727.32 248.007 12.8942 23.9316 32.176 3.50424 0.152037 ⋯
14 │ 1753.38 252.684 9.88844 25.2124 21.3448 7.64865 -0.00704676 ⋯
15 │ 1731.05 262.218 7.70813 25.0737 18.1157 4.89353 -0.280415 ⋯
16 │ 1742.68 252.273 12.7376 24.4489 20.684 7.47363 0.0242163 ⋯
17 │ 1727.37 243.956 12.107 24.6337 26.502 3.89649 -0.411467 ⋯
⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱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.56524 13.4072
ρ1 │ -0.40997 1.0
σ │ 22.6674 28.5759
σ1 │ 10.5608 33.258
σ2 │ 3.08371 7.73195A 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)312Close examination of the histogram shows a few values of -1.
count(tblsleep.ρ1 .≈ -1)0Furthermore 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)1There 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))313sleepstudy random-slope model and report a 95% confidence interval for the days effect.Use a fixed RNG for reproducibility and shortestcovint (or the confint method) on the bootstrap object.
using Random
df = DataFrame(dataset(:sleepstudy))
m = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days | subj)), df)
boot = parametricbootstrap(MersenneTwister(42), 1000, m)
DataFrame(shortestcovint(boot))count(issingular, boot) (or inspecting the relevant column) gives the number of singular replicates. A large fraction of singular fits means the correlation parameter is poorly identified, so its bootstrap distribution piles up at ±1 and its confidence interval should be read with caution.
This page was rendered from git revision 29e8d33 using Quarto 1.10.18 and Julia 1.12.7.