One-Way ANOVA with Random Effects

Author

Sven Rieger

Last updated on

December 1, 2023

Abstract

This part introduces the One-Way ANOVA with Random Effects (or so-called empty model). The material relies on Chapter 2 of the book Hierarchical Linear Models: Applications and Data Analysis Methods (raudenbush2002?).

Revealjs Presentation

If you want to see the presentation in full screen go to Other Formats on the right.

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"))
}
```

merTools (v0.6.1, Knowles & Frederick, 2024)
lme4 (v1.1.34, Bates et al., 2024)
flextable (v0.9.2, Gohel & Skintzos, 2024)
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