kb07.qmd

---
title: "Bootstrapping a fitted model"
jupyter: julia-1.9
---

Begin by loading the packages to be used.

```{julia}
#| code-fold: true
#| output: false
using AlgebraOfGraphics
using CairoMakie
using DataFrames
using MixedModels
using ProgressMeter
using Random
using SMLP2023: dataset

CairoMakie.activate!(; type="svg")

ProgressMeter.ijulia_behavior(:clear)
```

# 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.
Furthermore, `Grouping` contrasts are assigned to the `subj` and `item` factors.
This is not a contrast per-se but an indication that these factors will be used as grouping factors for random effects and, therefore, there is no need to create a contrast matrix.
For large numbers of levels in a grouping factor, an attempt to create a contrast matrix may cause memory overflow.

It is not important in these cases but a good practice in any case.

```{julia}
#| output: false
contrasts = merge(
  Dict(nm => EffectsCoding() for nm in (:spkr, :prec, :load)),
  Dict(nm => Grouping() for nm in (:subj, :item)),
)
```

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)
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)
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}
:::