kkl15.qmd

---
title: 'RePsychLing Kliegl, Kuschela, & Laubrock (2015)- Reduction of Model Complexity'
engine: julia
julia:
  exeflags: ["--project", "--threads=auto"]
---

# Background

@Kliegl2015 is a follow-up to @Kliegl2011 (see also script `kwdyz11.qmd`) from an experiment looking at a variety of effects of visual cueing under four different cue-target relations (CTRs). In this experiment two rectangles are displayed (1) in horizontal orientation , (2) in vertical orientation, (3) in left diagonal orientation, or in (4) right diagonal orientation relative to a central fixation point. Subjects react to the onset of a small or a large visual target occurring at one of the four ends of the two rectangles. The target is cued validly on 70% of trials by a brief flash of the corner of the rectangle at which it appears; it is cued invalidly at the three other locations 10% of the trials each. This implies a latent imbalance in design that is not visible in the repeated-measures ANOVA, but we will show its effect in the random-effect structure and conditional modes.

There are a couple of differences between the first and this follow-up experiment, rendering it more a conceptual than a direct replication. First, the original experiment was carried out at Peking University and this follow-up at Potsdam University. Second, diagonal orientations of rectangles and large target sizes were not part of the design of @Kliegl2011.

We specify three contrasts for the four-level factor CTR that are derived from spatial, object-based, and attractor-like features of attention. They map onto sequential differences between appropriately ordered factor levels. Replicating @Kliegl2011, the attraction effect was not significant as a fixed effect, but yielded a highly reliable variance component (VC; i.e., reliable individual differences in positive and negative attraction effects cancel the fixed effect). Moreover, these individual differences in the attraction effect were negatively correlated with those in the spatial effect.

This comparison is of interest because a few years after the publication of @Kliegl2011, the theoretically critical correlation parameter (CP) between the spatial effect and the attraction effect was determined as the source of a non-singular LMM in that paper. The present study served the purpose to estimate this parameter with a larger sample and a wider variety of experimental conditions.

Here we also include two additional experimental manipulations of target size and orientation of cue rectangle. A similar analysis was reported in the parsimonious mixed-model paper [@Bates2015]; it was also used in a paper of GAMEMs [@Baayen2017]. Data and R scripts of those analyses are also available in [R-package RePsychLing](https://github.com/dmbates/RePsychLing/tree/master/data/).

The analysis is based on reaction times `rt` to maintain compatibility with @Kliegl2011.

In this vignette we focus on the reduction of model complexity. And we start with a quote:

“Neither the [maximal] nor the [minimal] linear mixed models are appropriate for most repeated measures analysis. Using the [maximal] model is generally wasteful and costly in terms of statiscal power for testing hypotheses. On the other hand, the [minimal] model fails to account for nontrivial correlation among repeated measurements. This results in inflated [T]ype I error rates when non-negligible correlation does in fact exist. We can usually find middle ground, a covariance model that adequately accounts for correlation but is more parsimonious than the [maximal] model. Doing so allows us full control over [T]ype I error rates without needlessly sacrificing power.”

Stroup, W. W. (2012, p. 185). _Generalized linear mixed models: Modern concepts, methods and applica?ons._ CRC Press, Boca Raton.

# Packages

```{julia}
#| code-fold: true
#| output: false
using AlgebraOfGraphics
using AlgebraOfGraphics: density
using BoxCox
using CairoMakie
using CategoricalArrays
using Chain
using DataFrameMacros
using DataFrames
using MixedModels
using MixedModelsMakie
using Random
using SMLP2024: dataset
using StatsBase
```

# Read data, compute and plot means

```{julia}
dat = DataFrame(dataset(:kkl15))
describe(dat)
```

```{julia}
dat_subj = combine(
  groupby(dat, [:Subj, :CTR]),
  nrow => :n,
  :rt => mean => :rt_m,
  :rt => (c -> mean(log, c)) => :lrt_m,
)
dat_subj.CTR = categorical(dat_subj.CTR, levels=levels(dat.CTR))
describe(dat_subj)
```

```{julia}
#| code-fold: true
#| fig-cap: Comparative boxplots of mean response time by subject under different conditions
#| label: fig-bxpltsubjcond
boxplot(
  dat_subj.CTR.refs,
  dat_subj.lrt_m;
  orientation=:horizontal,
  show_notch=true,
  axis=(;
    yticks=(
      1:4,
      [
        "valid cue",
        "same obj/diff pos",
        "diff obj/same pos",
        "diff obj/diff pos",
      ]
    )
  ),
  figure=(; resolution=(800, 300)),
)
```

Mean of reaction times for four cue-target relations. Targets appeared at (a) the cued position (valid) in a rectangle, (b) in the same rectangle cue, but at its other end, (c) on the second rectangle, but at a corresponding horizontal/vertical physical distance, or (d) at the other end of the second rectangle, that is $\sqrt{2}$ of horizontal/vertical distance diagonally across from the cue, that is also at larger physical distance compared to (c).

# Contrasts

```{julia}
contrasts = Dict(
  :CTR => SeqDiffCoding(; levels=["val", "sod", "dos", "dod"]),
  :cardinal => EffectsCoding(; levels=["cardinal", "diagonal"]),
  :size => EffectsCoding(; levels=["big", "small"])
)
```

# Maximum LMM

This is the maximum LMM for the design; `size` is a between-subject factor, 
ignoring other information such as trial number, age and gender of subjects.

```{julia}
m_max = let
  form = @formula rt ~ 1 + CTR * cardinal * size +
                           (1 + CTR * cardinal | Subj)
  fit(MixedModel, form, dat; contrasts)
end;
```

```{julia}
issingular(m_max)
```

```{julia}
only(MixedModels.PCA(m_max))
```

```{julia}
VarCorr(m_max)
```

The LMM `m_max` is overparameterized but it is not immediately apparent why. 

# Reduction strategy 1

##  Zero-correlation parameter LMM (1)

Force CPs to zero. _Reduction strategy 1_ is more suited for reducing model w/o theoretical expectations about CPs. The better reduction strategy for the present experiment with an _a priori_ interest in CPs is described as _Reduction strategy 2_.

```{julia}
m_zcp1 = let
  form = @formula rt ~ 1 + CTR * cardinal * size +
                   zerocorr(1 + CTR * cardinal | Subj)
  fit(MixedModel, form, dat; contrasts)
end;
```

```{julia}
issingular(m_zcp1)
```

```{julia}
only(MixedModels.PCA(m_zcp1))
```

```{julia}
VarCorr(m_zcp1)
```

The LMM `m_zcp1` is also overparameterized, but now there is clear evidence for absence of evidence for the VC of one of the interactions and the other two interaction-based VCs are also very small. 

## Reduced zcp LMM

Take out VCs for interactions.

```{julia}
m_zcp1_rdc = let
  form = @formula rt ~ 1 + CTR * cardinal * size +
                   zerocorr(1 + CTR + cardinal | Subj)
  fit(MixedModel, form, dat; contrasts)
end;
```

```{julia}
issingular(m_zcp1_rdc)
```

```{julia}
only(MixedModels.PCA(m_zcp1_rdc))
```

```{julia}
VarCorr(m_zcp1_rdc)
```

LMM `m_zcp_rdc` is ok . We add in CPs. 

## Parsimonious LMM (1)

Extend zcp-reduced LMM with CPs

```{julia}
m_prm1 = let
  form = @formula rt ~ 1 + CTR * cardinal * size +
                           (1 + CTR + cardinal | Subj)
  fit(MixedModel, form, dat; contrasts)
end;
```

```{julia}
issingular(m_prm1)
```

```{julia}
only(MixedModels.PCA(m_prm1))
```

LMM `m_zcp_rdc` is ok . We add in CPs. 

```{julia}
VarCorr(m_prm1)
```

We note that the critical correlation parameter between spatial (`sod`) and attraction (`dod`) is now estimated at .54 -- not that close to the 1.0 boundary that caused singularity in @Kliegl2011.

## Model comparison 1

```{julia}
gof_summary = let
  nms = [:m_zcp1_rdc, :m_prm1, :m_max]
  mods = eval.(nms)
  lrt = MixedModels.likelihoodratiotest(m_zcp1_rdc, m_prm1, m_max)
  DataFrame(;
    name = nms,
    dof=dof.(mods),
    deviance=round.(deviance.(mods), digits=0),
    AIC=round.(aic.(mods),digits=0),
    AICc=round.(aicc.(mods),digits=0),
     BIC=round.(bic.(mods),digits=0),
    χ²=vcat(:., round.(lrt.tests.deviancediff, digits=0)),
    χ²_dof=vcat(:., round.(lrt.tests.dofdiff, digits=0)),
    pvalue=vcat(:., round.(lrt.tests.pvalues, digits=3))
  )
end
```

AIC prefers LMM `m_prm1` over `m_zcp1_rdc`; BIC LMM `m_zcp1_rdc`. 
As the CPs were one reason for conducting this experiment, AIC is the criterion of choice.

# Reduction strategy 2

## Complex LMM
Relative to LMM `m_max`, first we take out interaction VCs and associated CPs, because these VCs are very small. This is the same as LMM `m_prm1` above. 

```{julia}
m_cpx = let
  form = @formula rt ~ 1 + CTR * cardinal * size +
                      (1 + CTR + cardinal | Subj)
  fit(MixedModel, form, dat; contrasts)
end;
```

## Zero-correlation parameter LMM (2)
Now we check the significance of ensemble of CPs.

```{julia}
m_zcp2 = let
  form = @formula rt ~ 1 + CTR * cardinal * size  +
              zerocorr(1 + CTR + cardinal | Subj)
  fit(MixedModel, form, dat; contrasts)
end;
```


```{julia}
VarCorr(m_zcp2)
```

## Parsimonious LMM (2)

The cardinal-related CPs are quite small. Do we need them?

```{julia}
m_prm2 = let
  form = @formula(rt ~ 1 + CTR * cardinal * size  +
                      (1 + CTR | Subj) + (0 + cardinal | Subj))
  fit(MixedModel, form, dat; contrasts)
end;
```

```{julia}
VarCorr(m_prm2)
```

## Model comparison 2

```{julia}
gof_summary = let
  nms = [:m_zcp2, :m_prm2, :m_cpx, :m_max]
  mods = eval.(nms)
  lrt = MixedModels.likelihoodratiotest(m_zcp2, m_prm2, m_cpx, m_max)
  DataFrame(;
    name = nms,
    dof=dof.(mods),
    deviance=round.(deviance.(mods), digits=0),
    AIC=round.(aic.(mods),digits=0),
    AICc=round.(aicc.(mods),digits=0),
     BIC=round.(bic.(mods),digits=0),
    χ²=vcat(:., round.(lrt.tests.deviancediff, digits=0)),
    χ²_dof=vcat(:., round.(lrt.tests.dofdiff, digits=0)),
    pvalue=vcat(:., round.(lrt.tests.pvalues, digits=3))
  )
end
```

The cardinal-related CPs could be removed w/o loss of goodness of fit. 
However, there is no harm in keeping them in the LMM. 
The data support both LMM `m_prm2` and  `m_cpx` (same as: `m_prm1`).
We keep the slightly more complex LMM `m_cpx` (`m_prm1`).

# Diagnostic plots of LMM residuals

Do model residuals meet LMM assumptions? Classic plots are

  - Residual over fitted
  - Quantiles of model residuals over theoretical quantiles of normal distribution

## Residual-over-fitted plot

The slant in residuals show a lower and upper boundary of reaction times, that is we have have too few short and too few long residuals. Not ideal, but at least width of the residual band looks similar across the fitted values, that is there is no evidence for heteroskedasticity.

```{julia}
#| code-fold: true
#| label: fig-m1fittedresid
#| fig-cap: Residuals versus fitted values for model m1
CairoMakie.activate!(; type="png")
scatter(fitted(m_prm1), residuals(m_prm1); alpha=0.3)
```

With many observations the scatterplot is not that informative. Contour plots or heatmaps may be an alternative.

```{julia}
#| code-fold: true
#| label: fig-m1fittedresid2
#| fig-cap: Heatmap of residuals versus fitted values for model m1
set_aog_theme!()
draw(
  data((; f=fitted(m_prm1), r=residuals(m_prm1))) *
  mapping(
    :f => "Fitted values from m1", :r => "Residuals from m1"
  ) *
  density();
)
```

## Q-Q plot

The plot of quantiles of model residuals over corresponding quantiles of the normal distribution should yield a straight line along the main diagonal.

```{julia}
#| code-fold: true
#| label: fig-qqnormm1
#| fig-cap: Quantile-quantile plot of the residuals for model m1 versus a standard normal
CairoMakie.activate!(; type="png")
qqnorm(
  residuals(m_prm1);
  qqline=:none,
  axis=(;
    xlabel="Standard normal quantiles",
    ylabel="Quantiles of the residuals from model m1",
  ),
)
```

# Conditional modes

## Caterpillar plot

```{julia}
#| code-fold: true
#| label: fig-caterpillarm1
#| fig-cap: Prediction intervals of the subject random effects in model m1
cm1 = only(ranefinfo(m_prm1))
caterpillar!(Figure(; resolution=(800, 1200)), cm1; orderby=2)
```

## Shrinkage plot


```{julia}
#| code-fold: true
#| label: fig-caterpillarm1L
#| fig-cap: Shrinkage plots of the subject random effects in model m1L
shrinkageplot!(Figure(; resolution=(1000, 1200)), m_prm1)
```

# Parametric bootstrap

Here we

  - generate a bootstrap sample
  - compute shortest covergage intervals for the LMM parameters
  - plot densities of bootstrapped parameter estimates for residual, fixed effects, variance components, and correlation parameters

## Generate a bootstrap sample

We generate 2500 samples for the 15 model parameters (4 fixed effect, 7 VCs, 15 CPs, and 1 residual).

```{julia}
samp = parametricbootstrap(MersenneTwister(1234321), 2500, m_prm1;
                           optsum_overrides=(; ftol_rel=1e-8));
```

```{julia}
tbl = samp.tbl
```

## Shortest coverage interval

```{julia}
confint(samp)
```

We can also visualize the shortest coverage intervals for fixed effects with the `ridgeplot()` command:

```{julia}
#| code-fold: true
#| label: fig-bsridgem1
#| fig-cap: Ridge plot of fixed-effects bootstrap samples from model m1L
ridgeplot(samp; show_intercept=false)
```

## Comparative density plots of bootstrapped parameter estimates

### Residual

```{julia}
#| code-fold: true
#| label: fig-sigmadensitym1
#| fig-cap: '  Kernel density estimate from bootstrap samples of the residual standard  deviation for model m_prm1'
draw(
  data(tbl) *
  mapping(:σ => "Residual") *
  density();
  figure=(; resolution=(800, 400)),
)
```

### Fixed effects and associated variance components (w/o GM)

The shortest coverage interval for the `GM` ranges from x to x ms and the associate variance component from .x to .x. To keep the plot range small we do not include their densities here.

```{julia}
#| code-fold: true
#| label: fig-betadensitym1
#| fig-cap: '  Kernel density estimate from bootstrap samples of the fixed effects for model  m_prm1'
rn = renamer([
  "(Intercept)" => "GM",
  "CTR: sod" => "spatial effect",
  "CTR: dos" => "object effect",
  "CTR: dod" => "attraction effect",
  "(Intercept), CTR: sod" => "GM, spatial",
  "(Intercept), CTR: dos" => "GM, object",
  "CTR: sod, CTR: dos" => "spatial, object",
  "(Intercept), CTR: dod" => "GM, attraction",
  "CTR: sod, CTR: dod" => "spatial, attraction",
  "CTR: dos, CTR: dod" => "object, attraction",
])
draw(
  data(tbl) *
  mapping(
    [:β02, :β03, :β04] .=> "Experimental effect size [ms]";
    color=dims(1) =>
    renamer(["spatial effect", "object effect", "attraction effect"]) =>
    "Experimental effects",
  ) *
  density();
  figure=(; resolution=(800, 350)),
)
```

The densitiies correspond nicely with the shortest coverage intervals.

```{julia}
#| code-fold: true
#| label: fig-sigmasdensitym1
#| fig-cap: '  Kernel density estimate from bootstrap samples of the standard deviations for  model m1L (excluding Grand Mean)'
draw(
  data(tbl) *
  mapping(
    [:σ2, :σ3, :σ4] .=> "Standard deviations [ms]";
    color=dims(1) =>
    renamer(["spatial effect", "object effect", "attraction effect"]) =>
    "Variance components",
  ) *
  density();
  figure=(; resolution=(800, 350)),
)
```

The VC are all very nicely defined.

### Correlation parameters (CPs)

```{julia}
#| code-fold: true
#| label: fig-corrdensitym1
#| fig-cap: '  Kernel density estimate from bootstrap samples of the standard deviations for  model m1L'
draw(
  data(tbl) *
  mapping(
    [:ρ01, :ρ02, :ρ03, :ρ04, :ρ05, :ρ06] .=> "Correlation";
    color=dims(1) =>
    renamer(["GM, spatial", "GM, object", "spatial, object",
    "GM, attraction", "spatial, attraction", "object, attraction"]) =>
    "Correlation parameters",
  ) *
  density();
  figure=(; resolution=(800, 350)),
)
```

Three CPs stand out positively, the correlation between GM and the spatial effect, GM and attraction effect, and the correlation between spatial and attraction effects.
The second CP was positive, but not significant in the first study.
The third CP replicates a CP that was judged questionable in script `kwdyz11.jl`.
The three remaining CPs are not well defined for reaction times.

# References

::: {#refs}
:::