One-Way ANOVA with Random Effects

Sven Rieger

December 1, 2023

Preface: Software I

  • The following packages are used:
mlmRaPkg <- c("merTools",
              "lme4",
              "flextable",
              "ggplot2")
  • Install packages when not already installed:
lapply(mlmRaPkg,
        function(x) 
          if(!x %in% rownames(installed.packages())) {
            install.packages(x)
            }
      )
  • Load (a subset of) the required package(s) into the R session.
library(lme4)
library(ggplot2)
library(flextable)

Preface: Software II

Print list of packages and cite them via Pandoc citation.

Show/hide fenced code
```{r}
#| label: write-mlmPkg
#| code-fold: true
#| code-summary: "Show/hide fenced code"
#| output-location: fragment
#| output: asis

for (i in 1:length(mlmRaPkg)) {
  
  cat(paste0(i, ". ",
             mlmRaPkg[i],
             " [", "v", utils::packageVersion(mlmRaPkg[i]),", @R-", mlmRaPkg[i],
             "]\n"))
}
```
  1. merTools (v0.6.1, Knowles & Frederick, 2024)
  2. lme4 (v1.1.34, Bates et al., 2024)
  3. flextable (v0.9.2, Gohel & Skintzos, 2024)
  4. ggplot2 (v3.4.3, Wickham et al., 2023)

Preface: Dataset

Recall, we use the HSB dataset from the merTools package (Knowles & Frederick, 2024). For a more detailed description see :

dat <- merTools::hsb
head(dat, 6)
  schid minority female    ses mathach size schtype meanses
1  1224        0      1 -1.528   5.876  842       0  -0.428
2  1224        0      1 -0.588  19.708  842       0  -0.428
3  1224        0      0 -0.528  20.349  842       0  -0.428
4  1224        0      0 -0.668   8.781  842       0  -0.428
5  1224        0      0 -0.158  17.898  842       0  -0.428
6  1224        0      0  0.022   4.583  842       0  -0.428

One-Way ANOVA with Random Effects I

The one-Way ANOVA with Random Effects Model is a model without predictors (often called the empty model [ger: “Nullmodell”]); the variance of a variable is decomposed into level 1 (within-) and level 2 (between-) variance components and it is the first step in (every?) multilevel analysis.

Level-specific equations

\[ \text{Level 1: } Y_{ij} = \beta_{0j} + r_{ij} \qquad(1)\]

\[ \text{Level 2: } \beta_{0j} = \gamma_{00} + u_{0j} \qquad(2)\]

Combined equation

Substituting Equation 2 into Equation 1 yields the combined model equation (see Equation 3):

\[ Y_{ij} = \gamma_{00} + u_{0j} + r_{ij} \qquad(3)\]

One-Way ANOVA with Random Effects II

Variance decomposition of \(Y_{ij}\)

\[ \begin{aligned} VAR(Y_{ij}) &= VAR(u_{0j}) + VAR(r_{ij}) \\ &= \tau_{00} + \sigma^2 \end{aligned} \qquad(4)\]


Intraclass correlation coefficient (ICC):

\[ ICC(\rho) = \frac{ \tau_{00} } {\tau_{00} + \sigma^2 } \qquad(5)\]

→ It measures the proportion of the variance in the outcome that is between level 2 units; if \(ICC = 0\) there is no need for multilevel modeling!

One-Way ANOVA with Random Effects in R I

raMod <- dat |>
  lmer(formula = mathach ~ 1 + 
                  (1 | schid),
       REML = FALSE)

summary(raMod)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: mathach ~ 1 + (1 | schid)
   Data: dat

     AIC      BIC   logLik deviance df.resid 
 47121.8  47142.4 -23557.9  47115.8     7182 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-3.06262 -0.75365  0.02676  0.76070  2.74184 

Random effects:
 Groups   Name        Variance Std.Dev.
 schid    (Intercept)  8.553   2.925   
 Residual             39.148   6.257   
Number of obs: 7185, groups:  schid, 160

Fixed effects:
            Estimate Std. Error t value
(Intercept)  12.6371     0.2436   51.87
dat |>
  lmer(formula = mathach ~ 1 + 
        (1 | schid),
       REML = FALSE) |>
  as_flextable()
Table 1: Random ANOVA Table generated by the flextable package (Gohel & Skintzos, 2024)

group

Estimate

Standard Error

statistic

Fixed effects

(Intercept)

12.637

0.244

51.873

Random effects

schid

sd__(Intercept)

2.925

Residual

sd__Observation

6.257

Signif. codes: 0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05

square root of the estimated residual variance: 6.3

data's log-likelihood under the model: -23,557.9

Akaike Information Criterion: 47,121.8

Bayesian Information Criterion: 47,142.4

 
raFixedEff <- summary(raMod) |>
  coefficients()

raFixedEff
            Estimate Std. Error  t value
(Intercept) 12.63707  0.2436171 51.87266
raRandomEff <- raMod |>
  VarCorr() |>
  as.data.frame() |>
  subset(select = "vcov")
  
raRandomEff
       vcov
1  8.553464
2 39.148400

Parameters

Recall the combined equation of the model (Equation 3): \(Y_{ij} = \gamma_{00} + u_{0j} + r_{ij}\)

Fixed effects:

\(\gamma_{00}\) = 12.637

Random effects:

\(VAR(u_{0j}) = \tau_{00}\) = 8.553

\(VAR(r_{ij}) = \sigma^2\) = 39.148

TBC …

References

Bates, D., Maechler, M., Bolker, B., & Walker, S. (2024). lme4: Linear mixed-effects models using eigen and S4. https://github.com/lme4/lme4/
Gohel, D., & Skintzos, P. (2024). Flextable: Functions for tabular reporting. https://ardata-fr.github.io/flextable-book/
Knowles, J. E., & Frederick, C. (2024). merTools: Tools for analyzing mixed effect regression models. https://CRAN.R-project.org/package=merTools
Wickham, H., Chang, W., Henry, L., Pedersen, T. L., Takahashi, K., Wilke, C., Woo, K., Yutani, H., & Dunnington, D. (2023). ggplot2: Create elegant data visualisations using the grammar of graphics. https://ggplot2.tidyverse.org