estimating reliability from varying intercepts

reliability_multilevel.R

library(tidyverse)
# simulate unobserved trait and items with error
N <- 100
dat <- tibble(
  id = 1:N,
  trait = rnorm(N),
  item1 =  trait + rnorm(N),
  item2 =  trait + rnorm(N),
  item3 =  trait + rnorm(N),
  item4 =  trait + rnorm(N),
  item5 =  trait + rnorm(N)
)

# classic approach
psych::alpha(dat %>% select(item1:item5))
# cronbach's alpha = 0.82

# use a multilevel model with varying intercepts instead
library(brms)
datlong <- dat %>% pivot_longer(item1:item5)
m1 <- brm(value ~ (1|id), datlong)


rels <- coef(m1)$id[, , "Intercept"] %>% # extract varying intercept per ID
  as.data.frame() %>% 
  rownames_to_column("id") %>% 
  # divide std.error of the varying intercepts by the SD of the r. intercepts
  # square it and subtract from 1
  mutate(rel = 1 - (Est.Error/VarCorr(m1)$id$sd["Intercept", "Estimate"])^2)

# aggregate rel across IDs
rels %>% 
  pull(rel) %>% 
  mean()
# mean rel = 0.82