ez
ez copied to clipboard
mike-lawrence
Type3 SS in ezANOVA not working in 4.2.0+
I am calculating a mixed-ANOVA with within and between subject variables. My groups are not completely equal in Ns, so I am trying to use Type 3 SS.
Starting from version 4.2.0, type=1 2 or 3 produces the same results. Whereas if I switch back to 4.1.1, it gives a more stringent results with using type 3.
Does anyone else experience the same issue?
I know it's been a long time since this first message, however I am still experiencing this in version 4.4.0. I was wondering if this was planning on being addressed.
I wasn't able to reproduce this when it was originally reported. I should have posted a message saying so. A user just emailed me separately with a reproducible example in which they were ignoring the standard result from ezANOVA(..., type=3) and instead looking at the result in ezANOVA(..., type=3, return_aov=TRUE)$aov, which I now see does indeed have an incorrect result. @mychan24 obviously you may not recall at all, but @mdefende do you recall if this latter usage matches what you were observing?
Yes, it looks like the reported sums of squares are not the same between the ANOVA and aov outputs for the ezANOVA(..., type = 3, return_aov=TRUE) call. One of the people I was working with determined the $ANOVA output was incorrect due to contrasts not being set properly for Type III sums of squares by default. Default is to use treatments contrasts (contr.treatment) which are not orthogonal, but Type III sums of squares have to have orthogonal contrasts to work correctly. He said using sums contrasts by calling options(contrasts=c(‘contr.sum', ‘contr.poly’)) will set the correct contrasts. This does indeed cause the reported sums of squares to be the same between $aov and $ANOVA outputs.
Hm, I can't replicate your colleague's observation that setting the contrasts option makes the SS's of $ANOVA and $aov match:
library(tidyverse)
library(ez)
#> Registered S3 methods overwritten by 'lme4':
#> method from
#> cooks.distance.influence.merMod car
#> influence.merMod car
#> dfbeta.influence.merMod car
#> dfbetas.influence.merMod car
#use the ANT dataset for demonstration
data(ANT)
#add an imbalanced second between-Ss variable
ANT$handedness = factor(ifelse(
as.numeric(ANT$subnum)%%3==1
, 'l'
, 'r'
))
#show the imbalanced design
ezDesign(
data = ANT
, y = subnum
, x = group
, col = handedness
)
#set treatment contrasts and run ezANOVA
options(contrasts = c('contr.treatment','contr.poly'))
out_treatment = ezANOVA(
data = ANT
, wid = subnum
, dv = rt
, between = .(group,handedness)
, type = 3
, detailed = T
, return_aov = TRUE
)
#> Warning: Data is unbalanced (unequal N per group). Make sure you specified
#> a well-considered value for the type argument to ezANOVA().
#> Warning: Collapsing data to cell means. *IF* the requested effects are a
#> subset of the full design, you must use the "within_full" argument, else
#> results may be inaccurate.
#> Coefficient covariances computed by hccm()
#combine, compare, and store for looking at later
tibble(
effect = out_treatment$ANOVA$Effect
, ANOVA_ssn = out_treatment$ANOVA$SSn
) %>%
dplyr::full_join(
tibble(
effect = str_trim(dimnames(summary(out_treatment$aov)[[1]])[[1]])
, aov_ssn = summary(out_treatment$aov)[[1]]$Sum
)
, by = 'effect'
) %>%
dplyr::mutate(
nearly_equal = near(ANOVA_ssn,aov_ssn)
) ->
treatment_compared
#set sum contrasts and run ezANOVA
options(contrasts = c('contr.sum','contr.poly'))
out_sum = ezANOVA(
data = ANT
, wid = subnum
, dv = rt
, between = .(group,handedness)
, type = 3
, detailed = T
, return_aov = TRUE
)
#> Warning: Data is unbalanced (unequal N per group). Make sure you specified
#> a well-considered value for the type argument to ezANOVA().
#> Warning: Collapsing data to cell means. *IF* the requested effects are a
#> subset of the full design, you must use the "within_full" argument, else
#> results may be inaccurate.
#> Coefficient covariances computed by hccm()
#combine, compare, and store for looking at later
tibble(
effect = out_sum$ANOVA$Effect
, ANOVA_ssn = out_sum$ANOVA$SSn
) %>%
dplyr::full_join(
tibble(
effect = str_trim(dimnames(summary(out_sum$aov)[[1]])[[1]])
, aov_ssn = summary(out_sum$aov)[[1]]$Sum
)
, by = 'effect'
) %>%
dplyr::mutate(
nearly_equal = near(ANOVA_ssn,aov_ssn)
) ->
sum_compared
#check if contrast type affects whether $ANOVA and $aov match
print(treatment_compared)
#> # A tibble: 5 x 4
#> effect ANOVA_ssn aov_ssn nearly_equal
#> <chr> <dbl> <dbl> <lgl>
#> 1 (Intercept) 2965387. NA NA
#> 2 group 98.8 141. FALSE
#> 3 handedness 0.0959 0.000933 FALSE
#> 4 group:handedness 17.5 17.5 TRUE
#> 5 Residuals NA 139. NA
print(sum_compared)
#> # A tibble: 5 x 4
#> effect ANOVA_ssn aov_ssn nearly_equal
#> <chr> <dbl> <dbl> <lgl>
#> 1 (Intercept) 2965387. NA NA
#> 2 group 98.8 141. FALSE
#> 3 handedness 0.0959 0.000933 FALSE
#> 4 group:handedness 17.5 17.5 TRUE
#> 5 Residuals NA 139. NA
#try setting contrasts explicitly ----
#set treatment contrasts and run ezANOVA
contrasts(ANT$group) = 'contr.treatment'
contrasts(ANT$handedness) = 'contr.treatment'
out_treatment = ezANOVA(
data = ANT
, wid = subnum
, dv = rt
, between = .(group,handedness)
, type = 3
, detailed = T
, return_aov = TRUE
)
#> Warning: Data is unbalanced (unequal N per group). Make sure you specified
#> a well-considered value for the type argument to ezANOVA().
#> Warning: Collapsing data to cell means. *IF* the requested effects are a
#> subset of the full design, you must use the "within_full" argument, else
#> results may be inaccurate.
#> Coefficient covariances computed by hccm()
#combine, compare, and store for looking at later
tibble(
effect = out_treatment$ANOVA$Effect
, ANOVA_ssn = out_treatment$ANOVA$SSn
) %>%
dplyr::full_join(
tibble(
effect = str_trim(dimnames(summary(out_treatment$aov)[[1]])[[1]])
, aov_ssn = summary(out_treatment$aov)[[1]]$Sum
)
, by = 'effect'
) %>%
dplyr::mutate(
nearly_equal = near(ANOVA_ssn,aov_ssn)
) ->
treatment_compared
#set sum contrasts and run ezANOVA
contrasts(ANT$group) = 'contr.sum'
contrasts(ANT$handedness) = 'contr.sum'
out_sum = ezANOVA(
data = ANT
, wid = subnum
, dv = rt
, between = .(group,handedness)
, type = 3
, detailed = T
, return_aov = TRUE
)
#> Warning: Data is unbalanced (unequal N per group). Make sure you specified
#> a well-considered value for the type argument to ezANOVA().
#> Warning: Collapsing data to cell means. *IF* the requested effects are a
#> subset of the full design, you must use the "within_full" argument, else
#> results may be inaccurate.
#> Coefficient covariances computed by hccm()
#combine, compare, and store for looking at later
tibble(
effect = out_sum$ANOVA$Effect
, ANOVA_ssn = out_sum$ANOVA$SSn
) %>%
dplyr::full_join(
tibble(
effect = str_trim(dimnames(summary(out_sum$aov)[[1]])[[1]])
, aov_ssn = summary(out_sum$aov)[[1]]$Sum
)
, by = 'effect'
) %>%
dplyr::mutate(
nearly_equal = near(ANOVA_ssn,aov_ssn)
) ->
sum_compared
#check if contrast type affects whether $ANOVA and $aov match
print(treatment_compared)
#> # A tibble: 5 x 4
#> effect ANOVA_ssn aov_ssn nearly_equal
#> <chr> <dbl> <dbl> <lgl>
#> 1 (Intercept) 499580. NA NA
#> 2 group 12.7 141. FALSE
#> 3 handedness 9.49 0.000933 FALSE
#> 4 group:handedness 17.5 17.5 TRUE
#> 5 Residuals NA 139. NA
print(sum_compared)
#> # A tibble: 5 x 4
#> effect ANOVA_ssn aov_ssn nearly_equal
#> <chr> <dbl> <dbl> <lgl>
#> 1 (Intercept) 2965387. NA NA
#> 2 group 98.8 141. FALSE
#> 3 handedness 0.0959 0.000933 FALSE
#> 4 group:handedness 17.5 17.5 TRUE
#> 5 Residuals NA 139. NA
Created on 2020-02-13 by the reprex package (v0.3.0)
Though that does remind me that ezANOVA internally sets options(contrasts=c('contr.sum','contr.poly') and I should really warn users about that (I'd even forgotten about it and was confused for a bit why setting the contrasts via options myself was yielding identical results).
Hey @mike-lawrence, not entirely recalling what I did, but pretty sure I was calling the $aov and extracting numbers from that instead. I think maybe just adding these info to the README would be sufficient for those who are digging around to understand why the output is not what they thought.
I have mostly switched to lmer and afex, but appreciate the development you put into this :)