# Everything You Always Wanted to Know About ANOVA*

**R on I Should Be Writing**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

The goals of this post are to (1) examine what ANOVAs are and are not, (2) demonstrate what the various types of ANOVAs are, and (3) familiarize you with how R does ANOVA.

Here are some assumptions I make about you, the reader, in this post:

- You’re familiar with the ideas of multi-factor ANOVAs (what a main effect is, what interactions are…).
- You know some
`R`

– how to fit a linear model, how to wrangle some data. ~~You are IID and normally distributed.~~

# What *is* ANOVA?

ANOVA tables are a way of summarizing a model – any model – by presenting the results grouped by the model’s terms / effects.
These tables contain a test statistic (often *F*), which represents the combined significance of all the parameters associated with a term.
This test is sometimes called the omnibus test, and is often accompanied by a measure of effect size, informing us how much of the variation in the outcome can be attributed to this term.
In other words, these answer the question: *how much of the variation in \(Y\) is related to the variation in \(X\)?*

This type of model summary is often paired with **factorial designs** – designs in which all the predictors (independent variables) are categorical, and all possible combinations between the levels of the different categories are represented in the data (often in a balanced way).

In fact, ANOVAs and factorial designs are paired so often that the two have (wrongfully) become synonymous and ANOVAs are commonly seen as “a type of model” that is applicable to factorial data. A consequence of this is that it is common knowledge^{1} that ANOVAs are:

- A type of model, that is a special case of linear regression (i.e. with a conditionally normal outcome variable), where
- All predictors are categorical, and
- The model is maximal – it has all possible main effects and interactions.

These 3 points are in fact ** false**.

Let’s look at some examples…

# ANOVAs in R

##
*Toy Data*

d <- tibble::tribble( ~id, ~group, ~X, ~Z, ~Rx, ~condition, ~Y, 1L, "Gb", 102, 1L, "Placebo", "Ca", 584.07, 2L, "Ga", 52, 1L, "Placebo", "Cb", 790.29, 3L, "Gb", 134, 2L, "Dose100", "Ca", 875.76, 4L, "Gb", 128, 3L, "Dose100", "Cb", 848.37, 5L, "Ga", 78, 1L, "Dose250", "Ca", 270.42, 6L, "Gb", 150, 2L, "Dose250", "Cb", 999.87, 7L, "Ga", 73, 1L, "Placebo", "Ca", 364.1, 8L, "Ga", 87, 7L, "Placebo", "Cb", 420.84, 9L, "Gb", 115, 6L, "Dose100", "Ca", 335.78, 10L, "Gb", 113, 4L, "Dose100", "Cb", 627, 11L, "Gc", 148, 3L, "Dose250", "Ca", 607.79, 12L, "Gc", 82, 3L, "Dose250", "Cb", 329.32, 13L, "Ga", 139, 1L, "Placebo", "Ca", 335.56, 14L, "Ga", 65, 2L, "Placebo", "Cb", 669.04, 15L, "Gb", 139, 1L, "Dose100", "Ca", 405.04, 16L, "Gc", 96, 1L, "Dose100", "Cb", 367.15, 17L, "Gb", 50, 5L, "Dose250", "Ca", 27.37, 18L, "Gc", 90, 2L, "Dose250", "Cb", 468.69, 19L, "Ga", 90, 2L, "Placebo", "Ca", 584.67, 20L, "Ga", 116, 2L, "Placebo", "Cb", 277.71, 21L, "Gb", 78, 2L, "Dose100", "Ca", 266.01, 22L, "Gb", 60, 1L, "Dose100", "Cb", 0.04, 23L, "Gc", 112, 4L, "Dose250", "Ca", 593.25, 24L, "Ga", 63, 4L, "Dose250", "Cb", 512.26, 25L, "Ga", 89, 1L, "Placebo", "Ca", 635.57, 26L, "Ga", 97, 2L, "Placebo", "Cb", 468.69, 27L, "Gc", 76, 3L, "Dose100", "Ca", 514.66, 28L, "Gb", 83, 1L, "Dose100", "Cb", 264.87, 29L, "Gc", 84, 4L, "Dose250", "Ca", 220.34, 30L, "Gb", 88, 1L, "Dose250", "Cb", 216.54 ) d$id <- factor(d$id) d$group <- factor(d$group) d$Rx <- factor(d$Rx, levels = c("Placebo", "Dose100", "Dose250")) d$condition <- factor(d$condition)

dplyr::glimpse(d) #> Rows: 30 #> Columns: 7 #> $ id <fct> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1~ #> $ group <fct> Gb, Ga, Gb, Gb, Ga, Gb, Ga, Ga, Gb, Gb, Gc, Gc, Ga, Ga, Gb, ~ #> $ X <dbl> 102, 52, 134, 128, 78, 150, 73, 87, 115, 113, 148, 82, 139, ~ #> $ Z <int> 1, 1, 2, 3, 1, 2, 1, 7, 6, 4, 3, 3, 1, 2, 1, 1, 5, 2, 2, 2, ~ #> $ Rx <fct> Placebo, Placebo, Dose100, Dose100, Dose250, Dose250, Placeb~ #> $ condition <fct> Ca, Cb, Ca, Cb, Ca, Cb, Ca, Cb, Ca, Cb, Ca, Cb, Ca, Cb, Ca, ~ #> $ Y <dbl> 584.07, 790.29, 875.76, 848.37, 270.42, 999.87, 364.10, 420.~ m <- lm(Y ~ group + X, data = d)

This is a multiple regression model with a covariate `X`

and a 3-level factor `group`

. We can summarize the results in a coefficient table (aka a “regression” table):

summary(m) #> #> Call: #> lm(formula = Y ~ group + X, data = d) #> #> Residuals: #> Min 1Q Median 3Q Max #> -372.51 -166.47 -53.67 134.57 451.16 #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 118.608 145.557 0.815 0.42256 #> groupGb -102.591 94.853 -1.082 0.28937 #> groupGc -92.384 107.302 -0.861 0.39712 #> X 4.241 1.504 2.820 0.00908 ** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> Residual standard error: 218.8 on 26 degrees of freedom #> Multiple R-squared: 0.2383, Adjusted R-squared: 0.1504 #> F-statistic: 2.711 on 3 and 26 DF, p-value: 0.06556

Or we can look at an ANOVA table:

anova(m) #> Analysis of Variance Table #> #> Response: Y #> Df Sum Sq Mean Sq F value Pr(>F) #> group 2 8783 4391 0.0918 0.912617 #> X 1 380471 380471 7.9503 0.009077 ** #> Residuals 26 1244265 47856 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

While the `group`

term had 2 parameters in the coefficient table, it now has a single test with a `Df`

of 2 - the test represents the *total* significance of the two parameters combined! That is, these *F*-tests are **omnibus tests** for all the parameters of a given term.

**Note** that the model being summarized here is neither completely factorial (`X`

is a continuous covariate), nor maximal (the `group:X`

interaction is missing)!

By default `R`

calculates *type 1* sums of squares (SS) - these are also called *sequential SS*, because each term’s SS is calculated for its addition on top of the PREVIOUS terms - sequentially!
In our example, the effect of `X`

represents only what `X`

explains on top of what `group`

was already able to explain / predict - it is the *unique* contribution of `X`

, controlling for `group`

; the effect of `group`

however does *not* represent its unique contribution, but instead its *total* contribution.

This means that because the following models have the same terms *in a different order* they will produce different *type 1* ANOVA tables!

anova(lm(Y ~ group + X, data = d)) #> Analysis of Variance Table #> #> Response: Y #> Df Sum Sq Mean Sq F value Pr(>F) #> group 2 8783 4391 0.0918 0.912617 #> X 1 380471 380471 7.9503 0.009077 ** #> Residuals 26 1244265 47856 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 anova(lm(Y ~ X + group, data = d)) #> Analysis of Variance Table #> #> Response: Y #> Df Sum Sq Mean Sq F value Pr(>F) #> X 1 325745 325745 6.8067 0.01486 * #> group 2 63509 31754 0.6635 0.52353 #> Residuals 26 1244265 47856 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We can recreate the ANOVA table above by building a sequence of models, and comparing them:

m0 <- lm(Y ~ 1, data = d) m1 <- lm(Y ~ group, data = d) anova(m0, m1, m) #> Analysis of Variance Table #> #> Model 1: Y ~ 1 #> Model 2: Y ~ group #> Model 3: Y ~ group + X #> Res.Df RSS Df Sum of Sq F Pr(>F) #> 1 29 1633519 #> 2 27 1624737 2 8783 0.0918 0.912617 #> 3 26 1244265 1 380471 7.9503 0.009077 ** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 anova(m) # same SS values #> Analysis of Variance Table #> #> Response: Y #> Df Sum Sq Mean Sq F value Pr(>F) #> group 2 8783 4391 0.0918 0.912617 #> X 1 380471 380471 7.9503 0.009077 ** #> Residuals 26 1244265 47856 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Simultaneous Sum of Squares

There are also *type 2 SS* - also called *simultaneous SS*, because the SS of each term is calculated for its addition over all other terms - its unique contribution while controlling for the other terms.
Type 2 SS can be obtained with the `Anova()`

function from the {`car`

} package:

car::Anova(m, type = 2) #> Anova Table (Type II tests) #> #> Response: Y #> Sum Sq Df F value Pr(>F) #> group 63509 2 0.6635 0.523533 #> X 380471 1 7.9503 0.009077 ** #> Residuals 1244265 26 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We can recreate the above ANOVA table by building two sequences of models, and comparing them:

m_sans_X <- lm(Y ~ group, data = d) m_sans_group <- lm(Y ~ X, data = d) anova(m_sans_group, m) # Same SS as the type 2 test for group #> Analysis of Variance Table #> #> Model 1: Y ~ X #> Model 2: Y ~ group + X #> Res.Df RSS Df Sum of Sq F Pr(>F) #> 1 28 1307774 #> 2 26 1244265 2 63509 0.6635 0.5235 anova(m_sans_X, m) # Same SS as the type 2 test for X #> Analysis of Variance Table #> #> Model 1: Y ~ group #> Model 2: Y ~ group + X #> Res.Df RSS Df Sum of Sq F Pr(>F) #> 1 27 1624737 #> 2 26 1244265 1 380471 7.9503 0.009077 ** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Because the order of terms is usually of little importance, type 1 tests are rarely used in the analysis of factorial designs…

Unfortunately, they are R’s default…

## Adding Interactions

Things get a bit more complicated when interactions are involved, as type 2 SS treat interactions differently than main effects:

m_int <- lm(Y ~ group * X, data = d)

The SS of each main effect is calculated over all other main effects (simultaneously, as we’ve already seen) but *without* accounting for the interaction, while the SS of the interaction term is calculated on top of the underlying main effects (sequentially). So we get the unique contribution of each main effect when controlling only for the other main effect, and the unique contribution of the interaction controlling for the already-included contribution of the main effects.

car::Anova(m_int, type = 2) #> Anova Table (Type II tests) #> #> Response: Y #> Sum Sq Df F value Pr(>F) #> group 63509 2 1.3555 0.2768607 #> X 380471 1 16.2410 0.0004884 *** #> group:X 682026 2 14.5566 7.246e-05 *** #> Residuals 562240 24 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Is the same as^{2}:

anova(m_sans_group, m) # Same SS as the type 2 test for group #> Analysis of Variance Table #> #> Model 1: Y ~ X #> Model 2: Y ~ group + X #> Res.Df RSS Df Sum of Sq F Pr(>F) #> 1 28 1307774 #> 2 26 1244265 2 63509 0.6635 0.5235 anova(m_sans_X, m) # Same SS as the type 2 test for X #> Analysis of Variance Table #> #> Model 1: Y ~ group #> Model 2: Y ~ group + X #> Res.Df RSS Df Sum of Sq F Pr(>F) #> 1 27 1624737 #> 2 26 1244265 1 380471 7.9503 0.009077 ** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 anova(m, m_int) # Same SS as the type 2 test for group:X #> Analysis of Variance Table #> #> Model 1: Y ~ group + X #> Model 2: Y ~ group * X #> Res.Df RSS Df Sum of Sq F Pr(>F) #> 1 26 1244265 #> 2 24 562240 2 682026 14.557 7.246e-05 *** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In designs of higher order, each “order” is tested in a similar *simultaneous-sequential* manner. E.g., in a 3-way design, all main effects (1st order) are tested *simultaneously* (accounting for one another), **then** all 2-way interactions (2nd order) are tested *simultaneously* (accounting for the main effects *and* one another), and **then** the 3-way interaction is tested (accounting for all main effects and 2-way interactions).

### Type 2 vs Type 3

There is another type of simultaneous SS - the type 3 test, which does not treat interactions differently than main effects: the SS for each main effect / interaction is calculated as its contribution on top of all other main effects AND interactions. So the effect of `X`

is its unique contribution while controlling both for `group`

*and* for `group:X`

!

However, remember that we previously saw that these methods actually produce omnibus tests for the combined effect of the *parameters of each term*. But in the `m_int`

model the parameters labeled `X`

, `groupGb`

, and `groupGc`

are no longer parameters of the main effects - instead they are parameters of simple (i.e., conditional) effects!

summary(m_int) #> #> Call: #> lm(formula = Y ~ group * X, data = d) #> #> Residuals: #> Min 1Q Median 3Q Max #> -370.18 -67.15 33.68 111.35 172.14 #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 833.130 173.749 4.795 6.99e-05 *** #> groupGb -1308.898 233.218 -5.612 8.90e-06 *** #> groupGc -752.718 307.747 -2.446 0.0222 * #> X -4.041 1.942 -2.081 0.0482 * #> groupGb:X 13.041 2.419 5.390 1.55e-05 *** #> groupGc:X 7.731 3.178 2.432 0.0228 * #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> Residual standard error: 153.1 on 24 degrees of freedom #> Multiple R-squared: 0.6558, Adjusted R-squared: 0.5841 #> F-statistic: 9.146 on 5 and 24 DF, p-value: 5.652e-05

`X`

is the slope of`X`

(*when*`group=Ga`

`Ga`

is the reference level)

`groupGb`

is the difference between`group=Ga`

and`group=Gb`

*when*`X=0`

`groupGc`

is the difference between`group=Ga`

and`group=Gc`

*when*`X=0`

(Pay attention to the “when” - this is what makes them conditional.)

We can see that changing the reference group changes the test for `X`

:

d$group <- relevel(d$group, ref = "Gb") m_int2 <- lm(Y ~ group * X, data = d) car::Anova(m_int, type = 3) #> Anova Table (Type III tests) #> #> Response: Y #> Sum Sq Df F value Pr(>F) #> (Intercept) 538630 1 22.9922 6.994e-05 *** #> group 738108 2 15.7536 4.269e-05 *** #> X 101495 1 4.3325 0.04823 * #> group:X 682026 2 14.5566 7.246e-05 *** #> Residuals 562240 24 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 car::Anova(m_int2, type = 3) #> Anova Table (Type III tests) #> #> Response: Y #> Sum Sq Df F value Pr(>F) #> (Intercept) 219106 1 9.3528 0.005402 ** #> group 738108 2 15.7536 4.269e-05 *** #> X 910646 1 38.8722 1.918e-06 *** #> group:X 682026 2 14.5566 7.246e-05 *** #> Residuals 562240 24 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

How can we resolve this?

### Centering

By centering our predictors! Centering is transforming our data in such a way that 0 represents the overall mean. When this is done, conditioning on 0 is the same as conditioning on the overall mean = looking at the main effect!

For covariables this is easy enough:

d$X_c <- d$X - mean(d$X) # or scale(d$X, center = TRUE, scale = FALSE)

But how do we center a factor??

The answer is - use some type of orthogonal coding, for example `contr.sum()`

(effects coding). This makes the coefficients harder to interpret ^{3}, but we’re not looking at those anyway!

contrasts(d$group) <- contr.sum

Now when looking at type 3 tests, the main effects terms actually are main effects!

m_int3 <- lm(Y ~ group*X_c, data = d) car::Anova(m_int3, type = 3) #> Anova Table (Type III tests) #> #> Response: Y #> Sum Sq Df F value Pr(>F) #> (Intercept) 4743668 1 202.4902 3.401e-13 *** #> group 19640 2 0.4192 0.66231 #> X_c 143772 1 6.1371 0.02067 * #> group:X_c 682026 2 14.5566 7.246e-05 *** #> Residuals 562240 24 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Remember: ** ABC - Always Be Centering** (your predictors) - type 3 ANOVA tables make little sense without centering!

##
*Recreate the Type 3 ANOVA*

Unfortunately, we can’t just build a model without an interaction term and use it to recreate the type 3 ANOVA. Instead, we need to actually build the model matrix (i.e., the design matrix), and drop the correct columns:

mm <- model.matrix(m_int2) head(mm) #> (Intercept) groupGa groupGc X groupGa:X groupGc:X #> 1 1 0 0 102 0 0 #> 2 1 1 0 52 52 0 #> 3 1 0 0 134 0 0 #> 4 1 0 0 128 0 0 #> 5 1 1 0 78 78 0 #> 6 1 0 0 150 0 0

A type 3 test for a term, is equal to the comparison between a model without the parameters associated with that term and the **full model**:

m_sans_group <- lm(Y ~ mm[,-(2:3)], data = d) m_sans_X <- lm(Y ~ mm[,-4], data = d) m_sans_int <- lm(Y ~ mm[,-(5:6)], data = d) anova(m_sans_group, m_int2) # Same SS as type 3 for group #> Analysis of Variance Table #> #> Model 1: Y ~ mm[, -(2:3)] #> Model 2: Y ~ group * X #> Res.Df RSS Df Sum of Sq F Pr(>F) #> 1 26 1300348 #> 2 24 562240 2 738108 15.754 4.269e-05 *** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 anova(m_sans_X, m_int2) # Same SS as type 3 for X #> Analysis of Variance Table #> #> Model 1: Y ~ mm[, -4] #> Model 2: Y ~ group * X #> Res.Df RSS Df Sum of Sq F Pr(>F) #> 1 25 1472886 #> 2 24 562240 1 910646 38.872 1.918e-06 *** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 anova(m_sans_int, m_int2) # Same SS as type 3 for group:X #> Analysis of Variance Table #> #> Model 1: Y ~ mm[, -(5:6)] #> Model 2: Y ~ group * X #> Res.Df RSS Df Sum of Sq F Pr(>F) #> 1 26 1244265 #> 2 24 562240 2 682026 14.557 7.246e-05 *** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Compare to:

car::Anova(m_int2, type = 3) #> Anova Table (Type III tests) #> #> Response: Y #> Sum Sq Df F value Pr(>F) #> (Intercept) 219106 1 9.3528 0.005402 ** #> group 738108 2 15.7536 4.269e-05 *** #> X 910646 1 38.8722 1.918e-06 *** #> group:X 682026 2 14.5566 7.246e-05 *** #> Residuals 562240 24 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Okay, that’s some ugly stuff. Let us never look at that again.

## Balanced vs. Unbalanced Data

The distinction between types 1, 2 and 3 SS is only relevant when there is some dependency between predictors (aka some collinearity). In our example, we can see that `group`

and `X`

are somewhat co-linear (VIF / tolerance are not strictly 1):

performance::check_collinearity(m) #> # Check for Multicollinearity #> #> Low Correlation #> #> Term VIF Increased SE Tolerance #> group 1.08 1.04 0.92 #> X 1.08 1.04 0.92

In a factorial design, we might call this dependence / collinearity among our predictors an “unbalanced design” (the number of observations differs between cells), and when the predictors are completely independent we would call this a “balanced design” (equal number of observations in all cells).

Let’s look at two examples:

### Balanced data

We can see that `Rx`

and `condition`

are balanced:

table(d$Rx, d$condition) #> #> Ca Cb #> Placebo 5 5 #> Dose100 5 5 #> Dose250 5 5 chisq.test(d$Rx, d$condition)$statistic # Chisq is exactly 0 #> X-squared #> 0

And so type 2 and type 3 ANOVAs are identical:

contrasts(d$Rx) <- contr.sum contrasts(d$condition) <- contr.sum m_balanced <- lm(Y ~ condition * Rx, data = d) anova(m_balanced) #> Analysis of Variance Table #> #> Response: Y #> Df Sum Sq Mean Sq F value Pr(>F) #> condition 1 13666 13666 0.2162 0.6461 #> Rx 2 41379 20690 0.3273 0.7240 #> condition:Rx 2 61444 30722 0.4860 0.6210 #> Residuals 24 1517030 63210 car::Anova(m_balanced, type = 2) #> Anova Table (Type II tests) #> #> Response: Y #> Sum Sq Df F value Pr(>F) #> condition 13666 1 0.2162 0.6461 #> Rx 41379 2 0.3273 0.7240 #> condition:Rx 61444 2 0.4860 0.6210 #> Residuals 1517030 24 car::Anova(m_balanced, type = 3) #> Anova Table (Type III tests) #> #> Response: Y #> Sum Sq Df F value Pr(>F) #> (Intercept) 6422803 1 101.6112 4.204e-10 *** #> condition 13666 1 0.2162 0.6461 #> Rx 41379 2 0.3273 0.7240 #> condition:Rx 61444 2 0.4860 0.6210 #> Residuals 1517030 24 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

### Unbalanced data

However, `condition`

and `group`

are NOT balanced:

table(d$group, d$condition) #> #> Ca Cb #> Gb 6 6 #> Ga 5 6 #> Gc 4 3 chisq.test(d$group, d$condition)$statistic # Chisq is NOT 0 #> Warning in chisq.test(d$group, d$condition): Chi-squared approximation may be #> incorrect #> X-squared #> 0.2337662

And so type 2 and type 3 ANOVAs are NOT identical:

contrasts(d$group) <- contr.sum contrasts(d$condition) <- contr.sum m_unbalanced <- lm(Y ~ condition * group, data = d) anova(m_unbalanced) #> Analysis of Variance Table #> #> Response: Y #> Df Sum Sq Mean Sq F value Pr(>F) #> condition 1 13666 13666 0.2087 0.6519 #> group 2 7165 3582 0.0547 0.9469 #> condition:group 2 41204 20602 0.3146 0.7330 #> Residuals 24 1571485 65479 car::Anova(m_unbalanced, type = 2) #> Anova Table (Type II tests) #> #> Response: Y #> Sum Sq Df F value Pr(>F) #> condition 12048 1 0.1840 0.6718 #> group 7165 2 0.0547 0.9469 #> condition:group 41204 2 0.3146 0.7330 #> Residuals 1571485 24 car::Anova(m_unbalanced, type = 3) #> Anova Table (Type III tests) #> #> Response: Y #> Sum Sq Df F value Pr(>F) #> (Intercept) 5858845 1 89.4773 1.439e-09 *** #> condition 3452 1 0.0527 0.8203 #> group 8913 2 0.0681 0.9344 #> condition:group 41204 2 0.3146 0.7330 #> Residuals 1571485 24 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

As the deviation from perfect balance (i.e. independence) is larger, so will the differences between the types increase.

### Why Does This Happen?

For type 1 SS, each effect is estimated sequentially - so in the model `Y ~ A * B`

, the effect of B is its unique contribution after controlling for and partialling out the effect of A.
But if A and B are independent, then there is nothing from A to partial out of B, and vice versa - so the order does not matter!

The distinction between types 2 and 3 comes from how they estimate main effects in the presence of interactions:
Type 2 SS looks at the SS between the marginal means of `group`

, across the levels of `condition`

.
So the mean of the first group is:

Type 3 SS however looks at the SS between the means of `group`

, weighted by `condition`

. So the mean of group `a`

is

This makes type 3 SS invariant to the cell frequencies! (However, both equations would give the same result when the data are balanced.)

A lot has been said about type 2 vs type 3. I will not go into the weeds here, but it is important to note that

- Most statistical softwares (SAS, Stata, SPSS, …) default to type 3 SS with orthogonal factor coding (but covariables are
*not*mean-centered in most cases by default) (see Langsrud, 2003). This makes`R`

inconsistent as we’ve seen it defaults to type 1 ANOVA and treatment coding.

- Often in factorial designs, any imbalance in the design is incidental, so it is often beneficial to have a method that is invariant to such imbalances. (Though this might not be true if the data is observational.)

- Coefficient tables give results that are analogous to type 3 SS when all terms are covariables:

m_covs <- lm(Y ~ X * Z, data = d) car::Anova(m_covs, type = 2) #> Anova Table (Type II tests) #> #> Response: Y #> Sum Sq Df F value Pr(>F) #> X 324676 1 6.7657 0.01513 * #> Z 2430 1 0.0506 0.82373 #> X:Z 57640 1 1.2011 0.28315 #> Residuals 1247704 26 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 car::Anova(m_covs, type = 3) #> Anova Table (Type III tests) #> #> Response: Y #> Sum Sq Df F value Pr(>F) #> (Intercept) 83130 1 1.7323 0.1996 #> X 5866 1 0.1222 0.7294 #> Z 59917 1 1.2486 0.2740 #> X:Z 57640 1 1.2011 0.2831 #> Residuals 1247704 26 summary(m_covs) # same p-values as type 3 #> #> Call: #> lm(formula = Y ~ X * Z, data = d) #> #> Residuals: #> Min 1Q Median 3Q Max #> -384.97 -153.72 -23.98 141.28 422.79 #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 366.352 278.348 1.316 0.200 #> X 1.014 2.900 0.350 0.729 #> Z -112.681 100.843 -1.117 0.274 #> X:Z 1.175 1.072 1.096 0.283 #> #> Residual standard error: 219.1 on 26 degrees of freedom #> Multiple R-squared: 0.2362, Adjusted R-squared: 0.1481 #> F-statistic: 2.68 on 3 and 26 DF, p-value: 0.06772

**Note** once again that the model being summarized here is not factorial at all - both `Z`

and `X`

are continuous covariables!

## ANOVA Made Easy

Making sure that our factors are orthogonally coded is a pain in the @$$.

Thankfully, we have the `afex`

package which turns this:

d$group <- factor(d$group) d$condition <- factor(d$condition) contrasts(d$group) <- contr.sum contrasts(d$condition) <- contr.sum m_lm <- lm(Y ~ group * condition, d) car::Anova(m_lm, type = 3) #> Anova Table (Type III tests) #> #> Response: Y #> Sum Sq Df F value Pr(>F) #> (Intercept) 5858845 1 89.4773 1.439e-09 *** #> group 8913 2 0.0681 0.9344 #> condition 3452 1 0.0527 0.8203 #> group:condition 41204 2 0.3146 0.7330 #> Residuals 1571485 24 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Into this:

afex::aov_car(Y ~ group * condition + Error(id), data = d) #> Anova Table (Type 3 tests) #> #> Response: Y #> Effect df MSE F ges p.value #> 1 group 2, 24 65478.53 0.07 .006 .934 #> 2 condition 1, 24 65478.53 0.05 .002 .820 #> 3 group:condition 2, 24 65478.53 0.31 .026 .733 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Much easier!

# Other Types of Models

So far we’ve seen how ANOVAs are applied to linear OLS models. However the idea of omnibus tests per term can be extended to many other types of models. For models where SS cannot be calculated, analogous methods based on deviance or likelihood are used instead (read more in the `car::Anova()`

docs). Here are some examples:

## GLMs

### Logistic

m_logistic <- glm(condition ~ group * X_c, data = d, family = binomial()) car::Anova(m_logistic, type = 2) #> Analysis of Deviance Table (Type II tests) #> #> Response: condition #> LR Chisq Df Pr(>Chisq) #> group 0.15679 2 0.9246 #> X_c 0.75134 1 0.3861 #> group:X_c 1.10225 2 0.5763

### Poisson

m_poisson <- glm(Z ~ group * X_c, data = d, family = poisson()) car::Anova(m_poisson, type = 3) #> Analysis of Deviance Table (Type III tests) #> #> Response: Z #> LR Chisq Df Pr(>Chisq) #> group 0.88074 2 0.6438 #> X_c 0.04370 1 0.8344 #> group:X_c 0.11357 2 0.9448

### Ordinal

m_ordinal <- ordinal::clm(group ~ X_c * condition, data = d) anova(m_ordinal) #> Type I Analysis of Deviance Table with Wald chi-square tests #> #> Df Chisq Pr(>Chisq) #> X_c 1 0.4469 0.5038 #> condition 1 0.1584 0.6907 #> X_c:condition 1 0.6129 0.4337

**Note** that all of the models summarized here with an ANOVA table do not have a continuous (conditionally normal) outcome, and not *all* of their predictors are categorical.
We may have been inclined to summarize these models with a coefficient table, but it is equally valid to present an ANOVA table!

## (G)LMMs

m_mixed <- lmerTest::lmer(Y ~ X_c * group + (group | Z), data = d) anova(m_mixed, type = 2) #> Type II Analysis of Variance Table with Satterthwaite's method #> Sum Sq Mean Sq NumDF DenDF F value Pr(>F) #> X_c 290878 290878 1 23.645 14.4210 0.0008951 *** #> group 7198 3599 2 9.520 0.1784 0.8393243 #> X_c:group 582071 291036 2 21.938 14.4288 0.0001001 *** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Also for GLMMs:

m_mixed2 <- lme4::glmer(condition ~ X_c * group + (1 | Z), data = d, family = binomial()) car::Anova(m_mixed2, type = 3) #> Analysis of Deviance Table (Type III Wald chisquare tests) #> #> Response: condition #> Chisq Df Pr(>Chisq) #> (Intercept) 0.0886 1 0.7660 #> X_c 1.1917 1 0.2750 #> group 0.0829 2 0.9594 #> X_c:group 0.9972 2 0.6074

# Concluding Remarks

Although we might be more inclined to summarize our model with an ANOVA table when our model contains categorical predictors, especially when interactions are involved, we’ve seen that (contrary to what we may have learned) ANOVA tables do not require any special data (factorial, balanced, normal outcome), and can be used as an alternative to a coefficient table.

Regardless of how we summarize our model - with a coefficient table or with an ANOVA table, using type 1, 2 or 3 SS, with orthogonal or treatment (dummy) coding, with centered or uncentered covariables - our underlying model is equivalent - and will produce the same estimated simple effects, marginal means and contrasts. That is, the method we use to summarize our model will not have any bearing on whatever follow-up analysis we may wish to carry out (using `emmeans`

of course! Check out the materials from my R course: **Analysis of Factorial Designs**).^{4}

Hopefully now you know what ANOVAs ARE - just ~~ANOVA~~ another way of summarizing a statistical model, term by term…

FIN

**leave a comment**for the author, please follow the link and comment on their blog:

**R on I Should Be Writing**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.