# Is there a better way to compare Intraclass correlations using mixed models?

The situation is this:

The lab I work for is building an intercept-only mixed model, in part to estimate the intraclass correlation related to a particular random effect. They want to build two of these models - one for a treatment condition and one for a control condition - and compare differences in ICC.

This seems odd to me. Is there a better way to do this? My first thought was to build a single model add a random slope for the treatment/control variable, but I can't find an straightforward way to compute the two ICC's of interest in this case.

You could indeed fit a mixed model for both groups simultaneously and allow for different fixed effects and random intercepts per group. This could also give you the means to statistically compare if there is a difference in the ICCs in the two groups. The R code below shows an example on how this could be done using the lme4 package and simulated data.

# simulate some data
# two groups A and B, with different variances for the random effects,
# and different average outcomes

set.seed(2019)
n <- 100
id <- rep(1:20, each = 5)
b1 <- rnorm(n, sd = 2)
b2 <- rnorm(n, sd = 4)

y1 <- 10 + b1[id] + rnorm(n, sd = 1)
y2 <- 20 + b2[id] + rnorm(n, sd = 1)

# we put all data in a data.frame in the long format
DF <- data.frame(y = c(y1, y2),
group = gl(2, 100, labels = c("A", "B")),
id = c(id, id + 200))

#############################################3

library("lme4")

# we fit the two models separately for two groups
fm_A <- lmer(y ~ 1 + (1 | id), data = DF, subset = group == "A")
fm_B <- lmer(y ~ 1 + (1 | id), data = DF, subset = group == "B")

# we fit the model in all data with different fixed and random effects per group
fm_both <- lmer(y ~ group + (0 + dummy(group, "A") + dummy(group, "B") || id),
data = DF)

# compare variance components from the separate and joint fits
VarCorr(fm_A)
#>  Groups   Name        Std.Dev.
#>  id       (Intercept) 2.2201
#>  Residual             1.0240
VarCorr(fm_B)
#>  Groups   Name        Std.Dev.
#>  id       (Intercept) 4.03766
#>  Residual             0.97875

VarCorr(fm_both)
#>  Groups   Name              Std.Dev.
#>  id       dummy(group, "A") 2.2222
#>  id.1     dummy(group, "B") 4.0365
#>  Residual                   1.0016

# You could also formally test if the ICCs are the same in the two groups
gm <- lmer(y ~ group + (1 | id), data = DF)

anova(gm , fm_both)
#> refitting model(s) with ML (instead of REML)
#> Data: DF
#> Models:
#> gm: y ~ group + (1 | id)
#> fm_both: y ~ group + ((0 + dummy(group, "A") | id) + (0 + dummy(group,
#> fm_both:     "B") | id))
#>         Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
#> gm       4 733.66 746.85 -362.83   725.66
#> fm_both  5 729.21 745.71 -359.61   719.21 6.4465      1    0.01112 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• How would this change if you left the correlation between random slope and intercept in the model? – dmacfour Sep 12 at 17:26
• I'm also curious how you would calculate the ICC from the variances given in the model with the random slope. – dmacfour Sep 12 at 17:51