-
Notifications
You must be signed in to change notification settings - Fork 3
/
kb07.qmd
152 lines (120 loc) · 3.61 KB
/
kb07.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
---
title: "Bootstrapping a fitted model"
engine: julia
julia:
exeflags: ["--project", "--threads=auto"]
---
Begin by loading the packages to be used.
```{julia}
#| code-fold: true
#| output: false
using AlgebraOfGraphics
using CairoMakie
using DataFrames
using MixedModels
using Random
using SMLP2024: dataset
const progress=false
```
# Data set and model
The `kb07` data [@Kronmueller2007] are one of the datasets provided by the `MixedModels` package.
```{julia}
kb07 = dataset(:kb07)
```
Convert the table to a DataFrame for summary.
```{julia}
describe(DataFrame(kb07))
```
The experimental factors; `spkr`, `prec`, and `load`, are two-level factors.
The `EffectsCoding` contrast is used with these to create a $\pm1$ encoding.
```{julia}
#| output: false
contrasts = Dict{Symbol,Any}(nm => EffectsCoding() for nm in (:spkr, :prec, :load))
```
The display of an initial model fit
```{julia}
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
```
does not include the estimated correlations of the random effects.
The `VarCorr` extractor displays these.
```{julia}
VarCorr(kbm01)
```
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
```{julia}
kbm02 = let
form = @formula(
rt_trunc ~
1 + spkr + prec + load + (1 | subj) + (1 + prec | item)
)
fit(MixedModel, form, kb07; contrasts, progress)
end
```
```{julia}
VarCorr(kbm02)
```
These two models are nested and can be compared with a likelihood-ratio test.
```{julia}
MixedModels.likelihoodratiotest(kbm02, kbm01)
```
The p-value of approximately 17% leads us to prefer the simpler model, `kbm02`, to the more complex, `kbm01`.
# 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.
```{julia}
Random.seed!(1234321)
kbm02samp = parametricbootstrap(2000, kbm02)
kbm02tbl = kbm02samp.tbl
```
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.
```{julia}
confint(kbm02samp)
```
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, @fig-kbm02fedens.
```{julia}
#| code-fold: true
#| fig-cap: "Comparative densities of the fixed-effects coefficients in kbm02samp"
#| label: fig-kbm02fedens
draw(
data(kbm02samp.β) * mapping(:β; color=:coefname) * AlgebraOfGraphics.density();
figure=(; resolution=(800, 450)),
)
```
```{julia}
#| code-fold: true
#| fig-cap: "Density plot of bootstrap samples standard deviation of random effects"
#| label: fig-kbm02sampsigmadens
let pars = ["σ1", "σ2", "σ3"]
draw(
data(kbm02tbl) *
mapping(pars .=> "σ"; color=dims(1) => renamer(pars)) *
AlgebraOfGraphics.density();
figure=(; resolution=(800, 450)),
)
end
```
```{julia}
#| code-fold: true
#| fig-cap: "Density plot of correlation parameters in bootstrap sample from model kbm02"
#| label: fig-kbm02sampcorrdens
draw(
data(kbm02tbl) *
mapping(:ρ1 => "Correlation") *
AlgebraOfGraphics.density();
figure=(; resolution=(800, 450)),
)
```
# References
::: {#refs}
:::