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.6208167 0.6208167 1.55695 0.2586
Residuals  6 2.3924333 0.3987389               

Error: Within
          Df       Sum Sq      Mean Sq  F value     Pr(>F)    
C          2 0.0000750000 0.0000375000  0.96429 0.40892911    
A:C        2 0.0010583333 0.0005291667 13.60714 0.00082115 ***
Residuals 12 0.0004666667 0.0000388889                        
---
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.5569503852 0.2586013071796       2.059970247e-01
3      C   2  12  0.9642857143 0.4089291051673       3.134173988e-05
4    A:C   2  12 13.6071428571 0.0008211483483     * 4.420851101e-04

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

$`Sphericity Corrections`
  Effect          GGe          p[GG] p[GG]<.05          HFe          p[HF] p[HF]<.05
3      C 0.5517241379 0.370943601025           0.5845070423 0.374919141101          
4    A:C 0.5517241379 0.007837640541         * 0.5845070423 0.006627509443         *

#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.6208167 0.6208167 13.93528 0.0097023 **
Residuals  6 0.2673000 0.0445500                      
---
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.0000750 0.00003750 0.00021 0.99979
A:C        2 0.0010583 0.00052917 0.00299 0.99702
Residuals 12 2.1256000 0.17713333                

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