ranovaposthoc

Planned contrasts using ezANOVA output in R


I've been looking into using planned contrasts as opposed to post-hoc t-tests. I typically use ezANOVA (Type III ANOVA) but it seems that conducting planned contrasts using ezANOVA is not currently catered for.

aov() on the other hand is a Type I ANOVA (I don't want to get into a debate about which type is best for which type of design). It is straight forward to conduct planned contrasts using aov() (for between group designs) but I want to conduct a Type III ANOVA in a repeated measures and to be frank ezANOVA has a much more user friendly output.

Bearing in mind ezANOVA has the option to include return_aov = TRUE does anyone know of a way to use the information provided by ezANOVA to conduct planned contrasts?

Note: return_aov = TRUE allows access to the output of aov via something along these lines:

summary.lm(ModelName$aov$'Participant:IndependentVariable1')

Participant above is an example variable added to wid in ezANOVA:

wid = .(Participant)

summary.lm() is typically used when presenting the results of planned contrasts in aov, granted between groups ANOVAs as opposed to repeated measures.

I'm particularly interested in using the output to conduct planned contrasts for repeated measures ANOVA.

BOUNTY GOALS

Goals I would like to achieve from this bounty:

1) Use the output of ezANOVA to conduct planned contrasts in a repeated measures ANOVA.

1A) Use the output of ezANOVA to conduct planned contrasts on between subjects ANOVA (this one should be relatively easy & therefore not a requisite to claim the bounty.)

Any dummy data should suffice but here is a reminder of the format for ezANOVA repeated measures ANOVA:

ModelName <- ezANOVA(
data = DataSet,
dv = .(DependentVariable), 
wid = .(Participant), 
within = .(IndependentVariable1, IndependentVariable2), 
type=3, 
detailed = TRUE, 
return_aov = TRUE)

This is a Related Question with reproducible data and code which can be used in relation to this problem.

You can find a PDF here giving some background on planned contrasts and what they do.


Solution

  • The emmeans package provides the appropriate functionality to calculate custom contrasts/arbitrary linear functions of the estimates marginal means (EMMs) for aov and aovlist objects (see here for a full list of supported models).

    In the following I use the ANT data set which comes with the ez package.

    First we set up a mixed factorial ANOVA using ezANOVA. Note that one needs to set orthogonal contrasts in order to get meaningful Type-III tests (see, e.g., John Fox' answer here).

    library("ez")
    library("emmeans")
    
    # set orthogonal contrasts
    options(contrasts = c("contr.sum", "contr.poly"))
    
    data(ANT)
    rt_anova <- ezANOVA(data = ANT[ANT$error == 0, ], 
                        dv = rt,
                        wid = subnum, 
                        within = .(cue, flank),
                        between = group,
                        type = 3,
                        return_aov = TRUE)
    

    We can then calculate the EMMs for, say all group-flank combinations.

    emm <- emmeans(rt_anova$aov, ~ group * flank)
    emm
    ## group     flank         emmean       SE    df lower.CL upper.CL
    ## Control   Neutral     381.5546 1.735392 53.97 378.0753 385.0339
    ## Treatment Neutral     379.9286 1.735392 53.97 376.4493 383.4079
    ## Control   Congruent   381.6363 1.735392 53.97 378.1570 385.1155
    ## Treatment Congruent   379.7520 1.735392 53.97 376.2727 383.2313
    ## Control   Incongruent 466.6770 1.735392 53.97 463.1977 470.1563
    ## Treatment Incongruent 452.2352 1.735392 53.97 448.7559 455.7145
    

    Now it is easy to calculate all pairwise comparisons or any desired contrast on these EMMs.
    See also this book chapter and my answer here if you need a bit more detail on how one can derive the contrasts weights from the hypotheses.

    # all pairwise comparisons 
    pairs(emm, adjust = "Holm")
    ## contrast                                        estimate       SE    df t.ratio p.value
    ## Control,Neutral - Treatment,Neutral           1.62594836 2.454215 53.97   0.663  1.0000
    ## Control,Neutral - Control,Congruent          -0.08167403 2.473955 36.00  -0.033  1.0000
    ## Control,Neutral - Treatment,Congruent         1.80259257 2.454215 53.97   0.734  1.0000
    ## Control,Neutral - Control,Incongruent       -85.12239797 2.473955 36.00 -34.407  <.0001
    ## Control,Neutral - Treatment,Incongruent     -70.68062093 2.454215 53.97 -28.800  <.0001
    ## Treatment,Neutral - Control,Congruent        -1.70762239 2.454215 53.97  -0.696  1.0000
    ## Treatment,Neutral - Treatment,Congruent       0.17664421 2.473955 36.00   0.071  1.0000
    ## Treatment,Neutral - Control,Incongruent     -86.74834633 2.454215 53.97 -35.347  <.0001
    ## Treatment,Neutral - Treatment,Incongruent   -72.30656929 2.473955 36.00 -29.227  <.0001
    ## Control,Congruent - Treatment,Congruent       1.88426660 2.454215 53.97   0.768  1.0000
    ## Control,Congruent - Control,Incongruent     -85.04072394 2.473955 36.00 -34.374  <.0001
    ## Control,Congruent - Treatment,Incongruent   -70.59894690 2.454215 53.97 -28.766  <.0001
    ## Treatment,Congruent - Control,Incongruent   -86.92499054 2.454215 53.97 -35.419  <.0001
    ## Treatment,Congruent - Treatment,Incongruent -72.48321351 2.473955 36.00 -29.299  <.0001
    ## Control,Incongruent - Treatment,Incongruent  14.44177704 2.454215 53.97   5.884  <.0001
    ## 
    ## Results are averaged over the levels of: cue 
    ## P value adjustment: holm method for 15 tests 
    
    # custom contrasts
    contrast(
      emm, 
      list(c1 = c(1, -1, 0, 0, 0, 0), # reproduces first pairwise comparison
           # emmean of row 1 - (emmean of row 1 + emmean of row 2) / 2; see EMMs table
           # 381.5546 - (379.9286 + 381.6363) / 2
           c2 = c(1, -0.5, -0.5, 0, 0, 0))
     )
     ## contrast  estimate       SE    df t.ratio p.value
     ## c1       1.6259484 2.454215 53.97   0.663  0.5105
     ## c2       0.7721372 2.136825 43.84   0.361  0.7196
    

    The same applies for purely within-subjects ANOVAs or between-subjects ANOVAs.

    # within-subjects ANOVA
    rt_anova_wi <- ezANOVA(data = ANT[ANT$error == 0, ], 
                        dv = rt,
                        wid = subnum, 
                        within = .(cue, flank),
                        type = 3,
                        return_aov = TRUE)
    
    emm <- emmeans(rt_anova_wi$aov, ~ cue * flank)
    contrast(
      emm, 
      list(c1 = c(1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 
           c2 = c(1, -0.5, -0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0))
    ) 
    ## contrast estimate       SE     df t.ratio p.value
    ## c1       47.31005 3.802857 170.34  12.441  <.0001
    ## c2       50.35320 3.293371 170.34  15.289  <.0001
    
    # between-subjects ANOVA
    rt_anova_bw <- ezANOVA(data = ANT[ANT$error == 0, ], 
                           dv = rt,
                           wid = subnum, 
                           within_full = .(cue, flank), 
                           between = group,
                           type = 3,
                           return_aov = TRUE)
    
    emm_bw <- emmeans(rt_anova_bw$aov, ~ group)
    # custom linear function
    contrast(
      emm_bw, 
      list(c1 = c(2/3, 1/2)) 
    )
    ## contrast estimate        SE df t.ratio p.value
    ## c1       475.2899 0.8213448 18 578.673  <.0001