2
$\begingroup$

In a mixed effects model, I'm trying to understand the distinction between the following two models:

y ~ group + (1+group|subject)

and

y ~ group + (1|subject) + (1|subject:group)

where group is a fixed factor containing (say) two levels, subject is a random factor containing n levels, and each level of group appears multiple times within each subject. In the first model, there is a random intercept for each subject and a random slope for group within each subject. In the second model, there is a random intercept for each subject and there is a random intercept for each combination of subject and group. Conceptually, these two models appear to be doing the same thing -- but they produce different results using lmer in the lme4 R package.

I'm struggling to understand why using n random slopes (where each slope is across the two group levels within each subject) isn't essentially the same as using 2n random intercepts (one for each combination of subject and group).

Can anybody explain why these two models are not equivalent, and under which circumstances it would be appropriate to use each one?


Sample R code

# Load required libraries
library(lme4)
library(lmerTest)

# Generate dummy data
# Group 1:
n1 <- c(15,15,20,70,25,35)
m1 <- c(-1.9,-1.3,-1.4,-1.6,-1.1,-1.7)
s1 <- c(0.6,0.5,0.8,1,0.7,0.9)
# Group 2:
n2 <- c(100,200,250,480,130,270)
m2 <- c(-2,-1.9,-1.7,-2.9,-2.2,-2.2)
s2 <- c(0.8,1.1,1,1,1.2,1.1)
# Set random number seed and populate data frame:
set.seed(123)
dat <- data.frame(y=NULL,subject=NULL,group=NULL)
for(i in 1:6){
  dat <- rbind(dat,
               cbind(y=rnorm(n1[i],m1[i],s1[i]),
                     subject=i,
                     group=1))
  dat <- rbind(dat,
               cbind(y=rnorm(n2[i],m2[i],s2[i]),
                     subject=i,
                     group=2))
}
# Subject and group are factors
dat$subject <- factor(dat$subject)
dat$group <- factor(dat$group)

# Model 1: random slopes
out.slope <- lmer(y ~ group + (1+group|subject), data=dat)
summary(out.slope)
# Linear mixed model fit by REML t-tests use Satterthwaite approximations to degrees
#   of freedom [lmerMod]
# Formula: y ~ group + (1 + group | subject)
#    Data: dat
# 
# REML criterion at convergence: 4638.1
# 
# Scaled residuals: 
#     Min      1Q  Median      3Q     Max 
# -2.9831 -0.6452  0.0085  0.6617  3.4625 
# 
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr
#  subject  (Intercept) 0.02132  0.1460       
#           group2      0.14267  0.3777   0.03
#  Residual             1.02692  1.0134       
# Number of obs: 1610, groups:  subject, 6
# 
# Fixed effects:
#             Estimate Std. Error       df t value Pr(>|t|)    
# (Intercept) -1.55090    0.09988  3.55800 -15.528 0.000214 ***
# group2      -0.57190    0.17632  5.14500  -3.244 0.021938 *  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Correlation of Fixed Effects:
#        (Intr)
# group2 -0.349

# Model 2: random intercept for each combination of subject and group
out.interact <- lmer(y ~ group + (1|subject) + (1|subject:group), data=dat)
summary(out.interact)
# Linear mixed model fit by REML t-tests use Satterthwaite approximations to degrees
#   of freedom [lmerMod]
# Formula: y ~ group + (1 | subject) + (1 | subject:group)
#    Data: dat
# 
# REML criterion at convergence: 4640.1
# 
# Scaled residuals: 
#     Min      1Q  Median      3Q     Max 
# -2.9852 -0.6382  0.0111  0.6650  3.4675 
# 
# Random effects:
#  Groups        Name        Variance Std.Dev.
#  subject:group (Intercept) 0.09235  0.3039  
#  subject       (Intercept) 0.02319  0.1523  
#  Residual                  1.02633  1.0131  
# Number of obs: 1610, groups:  subject:group, 12; subject, 6
# 
# Fixed effects:
#             Estimate Std. Error      df t value Pr(>|t|)    
# (Intercept)  -1.5444     0.1625 14.9030  -9.502 1.03e-07 ***
# group2       -0.5796     0.1971  5.4940  -2.941   0.0288 *  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Correlation of Fixed Effects:
#        (Intr)
# group2 -0.704
$\endgroup$
  • $\begingroup$ The best way is to write two math models, and check the which one is more suitable to the variance-covariance structure of your data. If you go from math model to software model, you are the master of software; if you go directly to software model, you are the slave of the software. $\endgroup$ – user158565 Nov 29 '18 at 4:36

Your Answer

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

Browse other questions tagged or ask your own question.