I thought it could be helpful to have a thread on ANOVA in R. As a statistical consultant, this is the most frequent FAQ I get from clients - how to run a linear model on their data, conduct hypothesis tests, extract predicted means and perform contrasts.
To put is plainly - running anova in R is a hot mess. It is a patchwork of functions from different packages and different authors with differing philosophies. 
The first thing to understand is formula notation, which generally takes the form: `formula(y ~ x)`. The formula itself is only a string of characters with the class "formula" attached to it - it is not (by itself) linked to any data.
Main effects: x1 + x2
interactions: x1|x2
both: x1*x2
And like in other languages, check the variable type to make sure factors are factors or characters and integers are integers - so year 2020 is not read as 2,020.
interactions: x1|x2
both: x1*x2
And like in other languages, check the variable type to make sure factors are factors or characters and integers are integers - so year 2020 is not read as 2,020.
I don't use aov() at all because the output is limited and it does not handle missing or unbalanced data well. Instead I run a linear model (or nonlinear or generalised model) followed by anova.
The function `Anova()` from the car package is preferred because it can handle type 1, 2 or 3 sums of squares easily.
There is a tremendous number of functions associated with an object of class "lm" - use the `methods()` to investigate this.
To evaluate if a reduced model is statistically different from a fuller model, use `anova()`.
For extracting treatment least squares means, the package "emmeans" provides an excellent set of options for obtaining the means, confidence intervals and performing all types of contrasts.
"emmeans" has many types of pre-built contrasts to enable all sorts of common contrasts - pairwise, comparison to a control, comparison of linear v polynomial contrasts...many useful options! Search `?contrast.emmGrid` and `?`emmc-functions`` to learn more.
A few notes on custom contrasts:
For mixed models, there are two main options: nlme and lme4. Both are useful and provide substantially equivalent results, provided the models are specified correctly. Here is an lme4 example. It calls an lmer version of `anova()`
Many of these functions with common names like "anova" are actually doing different things, depending on the object it acts upon:
objects from lme4 are "merMod" objects with a large number of accessory function that can act on them:
much of lme4 functionality is the same as "lm" objects. Likelihood ratio tests and the least squares of fixed effects can be accessed the same:
Accessing the variance components can be done with `VarCorr()`
What I have not talked about is assessing diagnostics. This is surprisingly easy - call `plot()` on lm, nlme or merMod objects to see some diagnostic plots. The function `bartlett.test()` can be used for testing for heteroscedasticity.
One of the most important things not discussed at all is having a discussion with the client about their research objects and what precisely they want to know - this guides my process more than anything. The next thing I do is plot the data prior to running any analyses.
This is the general framework I follow when helping a client conduct an analysis of variance. While this process may take different paths depending on the data - like a generalised linear model - the overall process is the same.
