Bootstrapping a fitted model

Author

Phillip Alday, Douglas Bates, and Reinhold Kliegl

Published

2024-09-09

Begin by loading the packages to be used.

Code

using AlgebraOfGraphics
using CairoMakie
using DataFrames
using MixedModels
using Random
using SMLP2024: dataset

const progress=false

1 Data set and model

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.

describe(DataFrame(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. The EffectsCoding contrast is used with these to create a \(\pm1\) encoding.

contrasts = Dict{Symbol,Any}(nm => EffectsCoding() for nm in (:spkr, :prec, :load))

The display of an initial model fit

kbm01 = 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

	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

does not include the estimated correlations of the random effects.

The VarCorr extractor displays these.

VarCorr(kbm01)

	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

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

	Est.	SE	z	p	σ_item	σ_subj
(Intercept)	2181.8526	77.4681	28.16	<1e-99	364.7125	298.0259
spkr: old	67.8790	16.0785	4.22	<1e-04
prec: maintain	-333.7906	47.4472	-7.03	<1e-11	252.5212
load: yes	78.5904	16.0785	4.89	<1e-05
Residual	680.0319

VarCorr(kbm02)

	Column	Variance	Std.Dev	Corr.
item	(Intercept)	133015.242	364.713
	prec: maintain	63766.937	252.521	-0.70
subj	(Intercept)	88819.438	298.026
Residual		462443.388	680.032

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

MixedModels.likelihoodratiotest(kbm02, kbm01)

	model-dof	deviance	χ²	χ²-dof	P(>χ²)
rt_trunc ~ 1 + spkr + prec + load + (1 \| subj) + (1 + prec \| item)	9	28664
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	27	20	0.1431

The p-value of approximately 17% leads us to prefer the simpler model, kbm02, to the more complex, kbm01.

2 A bootstrap sample

Create a bootstrap sample of a few thousand parameter estimates from the reduced model. The pseudo-random number generator is initialized to a fixed value for reproducibility.

Random.seed!(1234321)
kbm02samp = parametricbootstrap(2000, kbm02)
kbm02tbl = kbm02samp.tbl

Table with 14 columns and 2000 rows:
      obj      β1       β2       β3        β4       σ        σ1       σ2       ⋯
    ┌───────────────────────────────────────────────────────────────────────────
 1  │ 28769.9  2152.69  60.3921  -359.134  54.1231  702.317  421.683  212.915  ⋯
 2  │ 28775.3  2221.81  61.7318  -349.417  71.1992  697.028  442.809  278.934  ⋯
 3  │ 28573.6  2325.31  53.9668  -404.225  81.3407  657.746  364.354  274.392  ⋯
 4  │ 28682.9  2274.65  52.9438  -384.66   69.6788  682.801  398.9    319.726  ⋯
 5  │ 28616.1  2191.34  54.7785  -357.919  111.829  673.878  253.248  235.913  ⋯
 6  │ 28634.1  2176.73  95.9793  -378.383  82.1559  671.685  510.996  309.876  ⋯
 7  │ 28620.4  2128.07  57.6928  -267.747  60.4621  666.986  396.878  295.009  ⋯
 8  │ 28665.9  2152.59  55.8208  -274.382  76.5128  684.746  369.091  274.446  ⋯
 9  │ 28651.2  2088.3   76.4636  -296.14   95.6124  680.282  265.684  186.041  ⋯
 10 │ 28679.9  2375.23  72.2422  -423.922  94.0833  680.729  425.72   287.686  ⋯
 11 │ 28686.8  2034.11  59.9149  -224.859  59.5448  682.433  308.675  272.181  ⋯
 12 │ 28702.2  2130.98  83.3953  -271.61   79.974   689.679  401.974  264.255  ⋯
 13 │ 28625.0  2055.35  55.3526  -363.937  67.7128  676.817  309.154  265.796  ⋯
 14 │ 28570.3  2160.73  81.5575  -321.913  108.979  661.524  426.869  259.825  ⋯
 15 │ 28613.6  2237.73  36.9996  -362.443  82.1812  673.873  312.861  241.595  ⋯
 16 │ 28681.4  2306.64  43.8801  -429.983  103.082  686.437  360.536  204.514  ⋯
 17 │ 28577.2  2143.78  66.2959  -373.598  93.157   668.15   321.47   183.119  ⋯
 ⋮  │    ⋮        ⋮        ⋮        ⋮         ⋮        ⋮        ⋮        ⋮     ⋱

One of the uses of such a sample is to form “confidence intervals” on the parameters by obtaining the shortest interval that covers a given proportion (95%, by default) of the sample.

confint(kbm02samp)

DictTable with 2 columns and 9 rows:
 par   lower      upper
 ────┬─────────────────────
 β1  │ 2028.01    2337.92
 β2  │ 38.431     99.5944
 β3  │ -439.321   -245.864
 β4  │ 46.0262    107.511
 ρ1  │ -0.897984  -0.445598
 σ   │ 655.249    701.497
 σ1  │ 261.196    448.512
 σ2  │ 175.489    312.049
 σ3  │ 228.099    357.789

A sample like this can be used for more than just creating an interval because it approximates the distribution of the estimator. For the fixed-effects parameters the estimators are close to being normally distributed, Figure 1.

Code

draw(
  data(kbm02samp.β) * mapping(:β; color=:coefname) * AlgebraOfGraphics.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/WgbrE/src/scenes.jl:227

Figure 1: Comparative densities of the fixed-effects coefficients in kbm02samp

Code

let pars = ["σ1", "σ2", "σ3"]
  draw(
    data(kbm02tbl) *
    mapping(pars .=> "σ"; color=dims(1) => renamer(pars)) *
    AlgebraOfGraphics.density();
    figure=(; resolution=(800, 450)),
  )
end

┌ 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/WgbrE/src/scenes.jl:227

Figure 2: Density plot of bootstrap samples standard deviation of random effects

Code

draw(
  data(kbm02tbl) *
  mapping(:ρ1 => "Correlation") *
  AlgebraOfGraphics.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/WgbrE/src/scenes.jl:227

Figure 3: Density plot of correlation parameters in bootstrap sample from model kbm02

3 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