6
$\begingroup$

While testing type I error rate using R, I found that I am getting higher than 5% false positives with a 2-way ANOVA.

However, I thought one of the reasons for using an ANOVA (as opposed to multiple tests) was to control the type I error. So, what am I doing wrong?

Here is the R code:

generate.dat <- function(N=30)
  # Generate data with 2 factors.
  # N is number of subjects per group.
{
  # id is the subject id, we use N*4 because there are four groups.
  # group is the full name of the group, A1, A2, B1, B2
  # value is the measured/observed result in our experiment, again N*4.
  # the mean value is 0 for all groups, i.e. there are no effects of group.
  dat <- data.frame(id=factor(1:(N*4)),
                    group=factor(c(
                          rep(c("A1", "A2"), each=N),
                          rep(c("B1", "B2"), each=N))),
                    value=rnorm(N*4))

  # split group into separate factors
  dat$factor1 <- factor(substring(dat$group, 1, 1))
  dat$factor2 <- factor(substring(dat$group, 2, 2))

  dat
}

# e.g.
# dat <- generate.dat(3)
# print(dat)
#   id group       value factor1 factor2
# 1   1    A1  0.42385602       A       1
# 2   2    A1 -0.34829466       A       1
# 3   3    A1 -1.40946883       A       1
# 4   4    A2 -0.09177423       A       2
# 5   5    A2 -1.26614034       A       2
# 6   6    A2 -1.18024188       A       2
# 7   7    B1 -0.86129559       B       1
# 8   8    B1 -1.30517594       B       1
# 9   9    B1  0.50221849       B       1
# 10 10    B2  0.42755864       B       2
# 11 11    B2 -0.03262990       B       2
# 12 12    B2  0.85593000       B       2

do.anova <- function(dat)
  # Perform a 2 way between-subjects anova on dat.
  # Return the lowest p value obtained.
  # (We are interested whether there is at least one significant effect
  # in this particular data set).
{
  results <- aov(value ~ factor1 * factor2, data=dat)
  min( summary(results)[[1]]$`Pr(>F)`[1:3] )
}


# Peform 1000 replications:
# each time generate some data and get the lowest p value
# from the resulting anova.
out <- replicate(1000, do.anova(generate.dat()))

# How many times out of the 1000 was there a significant effect?
sum(out < .05) / 1000 # about 14%

I doubt there is really a 14% chance of finding an effect, so what am I doing wrong?

Many thanks.

$\endgroup$
6
  • 2
    $\begingroup$ Note that $1-(1-0.05)^3 = 0.14$. $\endgroup$
    – amoeba
    Aug 11, 2014 at 10:34
  • 1
    $\begingroup$ @amoeba So if I did 3 separate tests, I would expect the type I error to be 14%. But I thought one of the reasons for using an ANOVA was to control the family-wise type I error rate, as an alternative to, say, multiple t-tests? $\endgroup$
    – trev
    Aug 11, 2014 at 11:42
  • 1
    $\begingroup$ What you are saying refers to a situation with one factor and many levels: e.g. if you do a single one-way ANOVA with 10 groups instead of doing 45 pairwise t-tests. Instead, here you have a situation with three factors (two simple and one interaction). I am not sure if ANOVA controls family-wise error rate when "family" refers to factors! Somebody more experienced in ANOVA testing will surely answer soon. $\endgroup$
    – amoeba
    Aug 11, 2014 at 12:53
  • 1
    $\begingroup$ @amoeba is quite right. Not the type-I error rate of each factor is controlled but the type-I error rate of the global F-test (i.e. the F-test that compares the full model with the empty model.) $\endgroup$
    – Michael M
    Aug 11, 2014 at 14:36
  • 1
    $\begingroup$ Really? Wow. But that means that a vast amount of studies have around a 14% chance of a type I error. Or more if there are more factors. I did not realize that. $\endgroup$
    – trev
    Aug 11, 2014 at 15:43

2 Answers 2

5
$\begingroup$

As we have already discussed in the comments above, ANOVA controls family-wise type I error rate across the "family" of levels, not across the "family" of factors. For example, one-way ANOVA with 10 groups (levels) controls the error rate, as opposed to performing all 45 pairwise t-tests that obviously runs into the multiple comparison problem. However, two-way ANOVA tests two factors and an interaction, so essentially performs three separate F tests. As is nicely demonstrated by your simulation, this also suffers from a multiple comparison problem, which is not controlled by ANOVA. Note that $1-(1-0.05)^3=0.14$, so your simulation provided a very precise estimate of the theoretically expected error rate.

I am no expert on that, but after I wrote my comments above I made a google search and stumbled upon Cramer et al. Hidden Multiplicity in Multiway ANOVA: Prevalence, Consequences, and Remedies, apparently a preprint from 2013. Let me quote their abstract verbatim:

Many empirical researchers do not realize that the common multiway analysis of variance (ANOVA) harbors a multiple comparison problem. In the case of two factors, three separate null hypotheses are subject to test (i.e., two main effects and one interaction). Consequently, the probability of a Type I error is 14% rather than 5%. We explain the multiple comparison problem and demonstrate that researchers almost never correct for it. We describe one of several correction procedures (i.e., sequential Bonferroni), and show that its application alters at least one of the substantive conclusions in 45 out of 60 articles considered. An alternative method to combat the multiplicity in multiway ANOVA is preregistration of hypotheses.

$\endgroup$
2
  • 1
    $\begingroup$ ANOVA can correct for multiple factors, or across the "family" of factors, but this usually requires an additional step that can be overlooked resulting in the conditions that you describe. $\endgroup$
    – Greg Snow
    Aug 11, 2014 at 17:33
  • $\begingroup$ +1 for a great find on the paper! There are also more references in there to the issue. $\endgroup$
    – trev
    Aug 12, 2014 at 9:14
3
$\begingroup$

The ANOVA procedures can control the family wise error across the entire set of factors and interactions, but you need to do it correctly, not look at the minimum of multiple p-values. One way to see the overall F test is to run summary.lm on your result object rather than just summary then look at the bottom of the print out. For simulation it can be easier to compare to a null model (intercept only) using the anova function. Here is a new version of your simulation function that does the overall test:

do.anova2 <- function(dat) {
  results <- aov(value ~ factor1 * factor2, data=dat)
  results0 <- update(results, .~1)
  anova(results0, results)[2,6]
}

When we run this version:

> out2 <- replicate(1000, do.anova2(generate.dat()))
> mean(out2 < 0.05)
[1] 0.043
> out2 <- replicate(1000, do.anova2(generate.dat()))
> mean(out2 < 0.05)
[1] 0.049
> out2 <- replicate(1000, do.anova2(generate.dat()))
> mean(out2 < 0.05)
[1] 0.051
>

We see the about 5% that you were expecting.

$\endgroup$
5
  • $\begingroup$ Is this the same as converting two-way ANOVA to one-way by "collapsing" all factors into one and treating all cells as levels of a single "combined" factor, and then looking at the p-value of this "flattened" ANOVA? $\endgroup$
    – amoeba
    Aug 11, 2014 at 19:45
  • $\begingroup$ @amoeba He is referring to comparing the full model to the null model with intercept only. This is one of the solutions in the Cramer paper you mentioned (p. 10). The problem is, this is not what researchers typically do in practise (based on the psychology journals they reviewed). $\endgroup$
    – trev
    Aug 12, 2014 at 9:16
  • $\begingroup$ I find it easier to explain this concept by writing the models in full. Not sure how clear this will be in the comments box, but: summary(aov(value ~ factor1*factor2, dat)) is the same as model1 <- lm(value ~ 1, dat); model2 <- lm(value ~ 1 + factor1, dat); model3 <- lm(value ~ 1 + factor1 + factor2, dat); model4 <- lm(value ~ 1 + factor1 + factor2 + factor1:factor2, dat); anova(model1, model2, model3, model4). What @Greg is doing is anova(model1, model4) $\endgroup$
    – trev
    Aug 12, 2014 at 9:19
  • $\begingroup$ @trev: Yes, but isn't the p-value of the two-way ANOVA compared to the model with intercept only (what Greg and you describe, and the authors of the paper above call an "omnibus" test) going to be the same as the p-value of the one-way ANOVA with all cells treated as levels of one factor? $\endgroup$
    – amoeba
    Aug 12, 2014 at 9:31
  • 1
    $\begingroup$ @amoeba Ah yes, I see what you mean. The resulting model will have the same number of terms (dummy variables, x1, x2 etc) and so the same p value, you are right. The interpretation of the individual model terms (x1, x2...) will be different of course, e.g. none of them represent the interaction between factors. $\endgroup$
    – trev
    Aug 12, 2014 at 10:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.