6.7.3: Using R

R - Mixed Effects Models

• Load the schools data.
• Obtain the ANOVA for the mixed effects model.
• Obtain estimators and CIs for means for each combination of region and school type.
• Obtain a means plot for each combination of region and school type.
• Obtain Tukey’s multiple comparisons CIs.

1. Load the schools data by using the following commands:

setwd("~/path-to-folder/")
schools_data <- read.table("schools_data.txt",header=T)
attach(schools_data)

2. Obtain the ANOVA for the mixed effects model by using the following commands:

library(lmerTest)
library(lme4)
mixed_schools<-lmer(SR_score ~ region + school_type + region:school_type + (1 | teacher : (region:school_type)) , schools_data)
summary(mixed_schools) # Partial output
#Random effects:
# Groups                       Name        Variance Std.Dev.
# (region:school_type):teacher (Intercept) 9.375    3.062   
# Residual                                 4.687    2.165   
# Number of obs: 16, groups:  (region:school_type):teacher, 8
anova(mixed_schools)
#Type III Analysis of Variance Table with Satterthwaites method
#                    Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
#region             112.812 112.812     1     4 24.0667 0.008011 **
#school_type         15.312  15.312     1     4  3.2667 0.144986   
#region:school_type  52.812  52.812     1     4 11.2667 0.028395 *
#---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Note that the command lmer() gives the ANOVA table only for the fixed effects. Therefore, in this example, since there are fixed effects, we get the ANOVA table with their \(F\) values and \(p\)-values.

In the "Random effects" section of the output, under the column variance, we get the estimates for \(\sigma_{\gamma}^{2}\) and \(\sigma^{2}\) which are equal to 9.375 and 4.687 respectively.

Alternatively, we can use the command aov() which gives a partial ANOVA table.

mixed_schools1<-aov(SR_score ~ region + school_type + region*school_type + Error((region*school_type)/teacher),schools_data)
summary(mixed_schools1)
#Error: region
#       Df Sum Sq Mean Sq
#region  1  564.1   564.1
#Error: school_type
#            Df Sum Sq Mean Sq
#school_type  1  76.56   76.56
#Error: region:school_type
#                   Df Sum Sq Mean Sq
#region:school_type  1  264.1   264.1
#Error: region:school_type:teacher
#          Df Sum Sq Mean Sq F value Pr(>F)
#Residuals  4  93.75   23.44               
#Error: Within
#          Df Sum Sq Mean Sq F value Pr(>F)
#Residuals  8   37.5   4.688               

3. Obtain estimators, CIs , and multiple comparisons CIs for means for each combination of region and school type by using the following commands:

library(emmeans)
pairwise_conf_intervals<-emmeans(mixed_schools,list(pairwise~region:school_type),adjust="Tukey")
CI<-confint(pairwise_conf_intervals)
$`emmeans of region, school_type`
# region school_type emmean   SE df lower.CL upper.CL
# EastUS Private       85.8 2.42  4     79.0     92.5
# WestUS Private       89.5 2.42  4     82.8     96.2
# EastUS Public        73.2 2.42  4     66.5     80.0
# WestUS Public        93.2 2.42  4     86.5    100.0
#Degrees-of-freedom method: kenward-roger
#Confidence level used: 0.95
$`pairwise differences of region, school_type`
# 1                               estimate   SE df lower.CL upper.CL
# EastUS Private - WestUS Private    -3.75 3.42  4   -17.69    10.19
# EastUS Private - EastUS Public     12.50 3.42  4    -1.44    26.44
# EastUS Private - WestUS Public     -7.50 3.42  4   -21.44     6.44
# WestUS Private - EastUS Public     16.25 3.42  4     2.31    30.19
# WestUS Private - WestUS Public     -3.75 3.42  4   -17.69    10.19
# EastUS Public - WestUS Public     -20.00 3.42  4   -33.94    -6.06
#Degrees-of-freedom method: kenward-roger
#Confidence level used: 0.95
#Conf-level adjustment: tukey method for comparing a family of 4 estimates

4. Obtain means plot for each combination of region and school type by using the following commands:

library(plotrix)
region_means<-as.data.frame(CI$`emmeans of region, school_type`)
region<-region_means$region
school_type<-region_means$school_type
region_school_type<-paste(region,school_type)
plotCI(x=region_means$emmean,y=NULL,li=region_means$lower.CL,ui=region_means$upper.CL,xaxt="n",xlab="Region*SchoolType",ylab="SR scores")
axis(1,at=1:4,labels=region_school_type)

Figure \(\PageIndex{1}\): SR scores mean plot for Region*SchoolType.

5. Obtain Tukey’s multiple comparisons plot by using the following commands:

diff_comp<-as.data.frame(CI$`pairwise differences of region, school_type`)
diff_reg_sch<-diff_comp[,1]
plotCI(x=diff_comp$estimate,y=NULL,li=diff_comp$lower.CL,ui=diff_comp$upper.CL,xaxt="n",xlab="",ylab="Differences of means")
abline(h=0)
axis(1,at=1:6,labels=diff_reg_sch,las=1,cex.axis=0.6)
detach(schools_data)

Figure \(\PageIndex{2}\): Tukey comparisons differences of means plot.