measure ICC

morning_call.R

library(lavaan)

X <- HolzingerSwineford1939[, paste0("x", 1:9)]
colMeans(X)
apply(X, 2, sd)

# min-max scaling version
X.p <- as.data.frame(apply(X, 2, function (x) (x - min(x)) / (max(x) - min(x))))
colMeans(X.p)
apply(X.p, 2, sd)

# create normal scores version
X.s <- as.data.frame(apply(X, 2, function (x) qnorm(rank(x) / (length(x) + 1))))
colMeans(X.s)  # about 0
apply(X.s, 2, sd)  # about 1

colnames(X.p) <- colnames(X.s) <- paste0("x.", 1:9)

# long form data
X.l.p <- reshape(X.p, direction = "long", varying = 1:9, timevar = "measure")
X.l.s <- reshape(X.s, direction = "long", varying = 1:9, timevar = "measure")

library(lme4)
library(nlme) # nlme provides a bit more flexibility

# ~ Current approach with min-max measures
(fit.orig <- lmer(x ~ (1 | measure) + (1 | id), X.l.p))

as.numeric(VarCorr(fit.orig)$measure)  # measure component
as.numeric(VarCorr(fit.orig)$id)  # person component
sigma(fit.orig) ^ 2  # random error variance component

# 16% of variance is due to "measure"
as.numeric(VarCorr(fit.orig)$measure) /
  (as.numeric(VarCorr(fit.orig)$measure) +
     as.numeric(VarCorr(fit.orig)$id) +
     sigma(fit.orig) ^ 2)

# Proposal with rankits
# only using ML rather than REML to better match psych results
(fit.prop <- lme(
  x ~ 1, random = ~ 1 | id, X.l.s, method = "ML"))

# ICC suggests that "measures" correlate about .26 on average:
as.numeric(getVarCov(fit.prop)) /
  (as.numeric(getVarCov(fit.prop)) + sigma(fit.prop) ^ 2)
# Average r using psych package:
psych::alpha(X.s[, 1:9])$total$average_r

# coefficient alpha using LME:
as.numeric(getVarCov(fit.prop)) /
  (as.numeric(getVarCov(fit.prop)) + sigma(fit.prop) ^ 2 / 9)
# coefficient alpha using psych:
psych::alpha(X.s[, 1:9])$total$raw_alpha

# rho parameter is average r:
gls(x ~ 1, correlation = corCompSymm(form = ~ 1 | id), X.l.s, method = "ML")