2
$\begingroup$

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.

Any advice?

$\endgroup$
1
$\begingroup$

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")
#> Loading required package: Matrix

# 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
$\endgroup$
  • $\begingroup$ How would this change if you left the correlation between random slope and intercept in the model? $\endgroup$ – dmacfour Sep 12 at 17:26
  • $\begingroup$ I'm also curious how you would calculate the ICC from the variances given in the model with the random slope. $\endgroup$ – dmacfour Sep 12 at 17:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.