QK Box 17.4

We illustrate how “default” error bars produced by most software packages can be misleading when summarizing mixed models. We generated two simple data sets for a three-factor mixed model design. Factor A has two levels and is fixed. Factor B is a random effect with four levels and is nested in A. Factor C is crossed with A & C and is fixed with three levels. In the first data set (response variable Y1), there are strong effects of the random factor, and effects of factor C are consistent across that factor. The second data set (Y2) has identical observations for each combination of A & C, but data were permuted so the effects of C vary across levels of B.

Preliminaries

First, load the required packages (afex, car, lattice, lme4, lmerTest, nlme, VCA, ez, emmeans) and, for convenience, apaTables

Import errbar data file (errbar.csv)

errbar <- read.csv("../data/errbar.csv")
errbar

set contrasts from afex make individual a factor, make sure species a factor too

set_sum_contrasts()

errbar$B <- factor(errbar$B)

Fit full model with Y1

errbar1.aov <- aov(Y1~A*C+Error(B),errbar)
summary(errbar1.aov)

Error: B
          Df Sum Sq Mean Sq F value Pr(>F)
A          1 0.6208  0.6208   1.557  0.259
Residuals  6 2.3924  0.3987               

Error: Within
          Df    Sum Sq   Mean Sq F value   Pr(>F)    
C          2 0.0000750 0.0000375   0.964 0.408929    
A:C        2 0.0010583 0.0005292  13.607 0.000821 ***
Residuals 12 0.0004667 0.0000389                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

use ez for comparison with type 3 SS - same result as design is balanced

ezerrbar <- ezANOVA(data=errbar, dv=Y1, wid=B, within=C, between=A, type=3)

Warning: Converting "C" to factor for ANOVA.Warning: Converting "A" to factor for ANOVA.

print(ezerrbar)

$ANOVA
  Effect DFn DFd          F            p p<.05          ges
2      A   1   6  1.5569504 0.2586013072       2.059970e-01
3      C   2  12  0.9642857 0.4089291052       3.134174e-05
4    A:C   2  12 13.6071429 0.0008211483     * 4.420851e-04

$`Mauchly's Test for Sphericity`
  Effect      W         p p<.05
3      C 0.1875 0.0152231     *
4    A:C 0.1875 0.0152231     *

$`Sphericity Corrections`
  Effect       GGe       p[GG] p[GG]<.05      HFe       p[HF] p[HF]<.05
3      C 0.5517241 0.370943601           0.584507 0.374919141          
4    A:C 0.5517241 0.007837641         * 0.584507 0.006627509         *

#GG issue, so print clean table with corrections
apa.ezANOVA.table(ezerrbar)

ANOVA results
 

Note. df_num indicates degrees of freedom numerator. df_den indicates degrees of freedom denominator. 
Epsilon indicates Greenhouse-Geisser multiplier for degrees of freedom, 
p-values and degrees of freedom in the table incorporate this correction.
ges indicates generalized eta-squared.
 

Repeat for Y2

errbar2.aov <- aov(Y2~A*C+Error(B),errbar)
summary(errbar2.aov)

Error: B
          Df Sum Sq Mean Sq F value Pr(>F)   
A          1 0.6208  0.6208   13.94 0.0097 **
Residuals  6 0.2673  0.0446                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Within
          Df Sum Sq Mean Sq F value Pr(>F)
C          2 0.0001 0.00004   0.000  1.000
A:C        2 0.0011 0.00053   0.003  0.997
Residuals 12 2.1256 0.17713               

Graphics

Plot raw data for Y1 and Y2

graphY1<-ggplot(data=errbar, aes(x=C, y=Y1, color=A, group=B))+
    geom_point(size=2,show.legend = FALSE)+
  geom_line(show.legend=FALSE)+
  theme_qk()+scale_color_viridis(discrete=TRUE, option="H", begin=0.25)
graphY2<-ggplot(data=errbar, aes(x=C, y=Y2, color=A, group=B))+
    geom_point(size=2,show.legend = FALSE)+
  geom_line(show.legend=FALSE)+
  theme_qk()+scale_color_viridis(discrete=TRUE, option="H", begin=0.25)
emm1<-emmeans(errbar1.aov, ~A|C)      #extract means
emm2<-as.data.frame(emm1)                    #convert data table to frame

Plot means and standard errors The plot uses Y1, but could as easily be done with Y2

pd=position_dodge(width=0.1)
graphmean<-ggplot(emm2,aes(x=C,y=emmean,shape=A, group=A, color=A))+
  geom_point(size=3,show.legend = FALSE, position=pd)+    #size=3 for larger symbols, as means shown (so fewer points on fig)
  geom_errorbar(aes(ymin = emmean-SE, ymax = emmean+SE), position=pd, width=0, show.legend = FALSE)+
  geom_line(show.legend=FALSE, position=pd)+
  labs(y="Mean")+
   theme_qk()+scale_color_viridis(discrete=TRUE, option="H", begin=0.25)

Combined graph

graphY1+graphmean+graphY2