I read that ANOVA and linear regression are the same thing. How can that be, considering that the output of ANOVA is some $F$ value and some $p$-value based on which you conclude if the sample means across the different samples are same or different.
But assuming the means are not equal (reject null hypothesis), ANOVA tells you nothing about the coefficients of the linear model. So how is linear regression the same as ANOVA?
ANOVA and linear regression are equivalent when the two models test against the same hypotheses and use an identical encoding. The models differ in their basic aim: ANOVA is mostly concerned to present differences between categories' means in the data while linear regression is mostly concern to estimate a sample mean response and an associated $\sigma^2$.
Somewhat aphoristically one can describe ANOVA as a regression with dummy variables. We can easily see that this is the case in the simple regression with categorical variables. A categorical variable will be encoded as a indicator matrix (a matrix of 0/1 depending on whether a subject is part of a given group or not) and then used directly for the solution of the linear system described by a linear regression.
Let's see an example with 5 groups. For the sake of argument I will assume that the mean of group1 equals 1, the mean of group2 equals 2, ... and the mean of group5 equals 5. (I use MATLAB, but the exact same thing is equivalent in R.)
rng(123); % Fix the seed
X = randi(5,100,1); % Generate 100 random integer U[1,5]
Y = X + randn(100,1); % Generate my response sample
Xcat = categorical(X); % Treat the integers are categories
% One-way ANOVA
[anovaPval,anovatab,stats] = anova1(Y,Xcat);
% Linear regression
fitObj = fitlm(Xcat,Y);
% Get the group means from the ANOVA
ANOVAgroupMeans = stats.means
% ANOVAgroupMeans =
% 1.0953 1.8421 2.7350 4.2321 5.0517
% Get the beta coefficients from the linear regression
LRbetas = [fitObj.Coefficients.Estimate']
% LRbetas =
% 1.0953 0.7468 1.6398 3.1368 3.9565
% Rescale the betas according the intercept
scaledLRbetas = [LRbetas(1) LRbetas(1)+LRbetas(2:5)]
% scaledLRbetas =
% 1.0953 1.8421 2.7350 4.2321 5.0517
% Check if the two results are numerically equivalent
abs(max( scaledLRbetas - ANOVAgroupMeans))
% ans =
% 2.6645e-15
As it can be seen in this scenario the results where exactly the same. The minute numerical difference is due to the design not being perfectly balanced as well as the underlaying estimation procedure; the ANOVA accumulates numerical errors a bit more aggressively. To that respect we fit an intercept, LRbetas(1); we could fit an intercept-free model but that would not be a "standard" linear regression. (The results would be even closer to ANOVA in that case though.)
The $F$-statistic (a ratio of the means) in the case of the ANOVA and in the case of linear regression will be also be the same for the above example:
abs( fitObj.anova.F(1) - anovatab{2,5} )
% ans =
% 2.9132e-13
This is because procedures test the same hypothesis but with different wordings: ANOVA will qualitatively check if "the ratio is high enough to suggest that no grouping is implausible" while linear regression will qualitatively check if "the ratio is high enough to suggest an intercept only model is possibly inadequate".
(This is a somewhat free interpretation of the "possibility to see a value equal or greater than the one observed under the null hypothesis" and it is not meant to be a text-book definition.)
Coming back to the final part of your question about "ANOVA tell(ing) you nothing about the coefficients of the linear model (assuming the means are not equal") I hope you can now see that the ANOVA, in the case that your design is simple/balanced enough, tells you everything that a linear model would. The confidence intervals for group means will be the same you have for your $\beta$, etc. Clearly when ones starts adding multiple covariate in his regression model, a simple one-way ANOVA does not have a direct equivalence. In that case one augments the information used to calculate the linear regression's mean response with information that are not directly available for a one way ANOVA. I believe that one can re-express things in ANOVA terms once more but it is mostly an academic exercise.
An interesting paper on the matter is Gelman's 2005 paper titled: Analysis of Variance - Why it is more important than ever. Some important points raised; I am not fully supportive of the paper (I think I personally align much more with McCullach's view) but it can be a constructive read.
As a final note: The plot thickens when you have mixed effects models. There you have different concepts about what can be considered a nuisance or actual information regarding the grouping of your data. These issues are outside the scope of this question but I think they are worthy of a nod.
Let me put some color into the idea that OLS with categorical (dummy-coded) regressors is equivalent to the factors in ANOVA. In both cases there are levels (or groups in the case of ANOVA).
In OLS regression it is most usual to have also continuous variables in the regressors. These logically modify the relationship in the fit model between the categorical variables and the dependent variable (D.C.). But not to the point of making the parallel unrecognizable.
Based on the mtcars data set we can first visualize the model lm(mpg ~ wt + as.factor(cyl), data = mtcars) as the slope determined by the continuous variable wt (weight), and the different intercepts projecting the effect of the categorical variable cylinder (four, six or eight cylinders). It is this last part that forms a parallel with a one-way ANOVA.
Let's see it graphically on the sub-plot to the right (the three sub-plots to the left are included for side-to-side comparison with the ANOVA model discussed immediately afterwards):
Each cylinder engine is color coded, and the distance between the fitted lines with different intercepts and the data cloud is the equivalent of within-group variation in an ANOVA. Notice that the intercepts in the OLS model with a continuous variable (weight) is not mathematically the same as the value of the different within-group means in ANOVA, due to the effect of weight and the different model matrices (see below): the mean mpg for 4-cylinder cars, for example, is mean(mtcars$mpg[mtcars$cyl==4]) #[1] 26.66364, whereas the OLS "baseline" intercept (reflecting by convention cyl==4 (lowest to highest numerals ordering in R)) is markedly different: summary(fit)$coef[1] #[1] 33.99079. The slope of the lines is the coefficient for the continuous variable weight.
If you try to suppress the effect of weight by mentally straightening these lines and returning them to the horizontal line, you'll end up with the ANOVA plot of the model aov(mtcars$mpg ~ as.factor(mtcars$cyl)) on the three sub-plots to the left. The weight regressor is now out, but the relationship from the points to the different intercepts is roughly preserved - we are simply rotating counter-clockwise and spreading out the previously overlapping plots for each different level (again, only as a visual device to "see" the connection; not as a mathematical equality, since we are comparing two different models!).
Each level in the factor cylinder is separate, and the vertical lines represent the residuals or within-group error: the distance from each point in the cloud and the mean for each level (color-coded horizontal line). The color gradient gives us an indication of how significant the levels are in validating the model: the more clustered the data points are around their group means, the more likely the ANOVA model will be statistically significant. The horizontal black line around $\small 20$ in all the plots is the mean for all the factors. The numbers in the $x$-axis are simply the placeholder number/identifier for each point within each level, and don't have any further purpose than to separate points along the horizontal line to allow a plotting display different to boxplots.
And it is through the sum of these vertical segments that we can manually calculate the residuals:
mu_mpg <- mean(mtcars$mpg) # Mean mpg in dataset
TSS <- sum((mtcars$mpg - mu_mpg)^2) # Total sum of squares
SumSq=sum((mtcars[mtcars$cyl==4,"mpg"]-mean(mtcars[mtcars$cyl=="4","mpg"]))^2)+
sum((mtcars[mtcars$cyl==6,"mpg"] - mean(mtcars[mtcars$cyl=="6","mpg"]))^2)+
sum((mtcars[mtcars$cyl==8,"mpg"] - mean(mtcars[mtcars$cyl=="8","mpg"]))^2)
The result: SumSq = 301.2626 and TSS - SumSq = 824.7846. Compare to:
Call:
aov(formula = mtcars$mpg ~ as.factor(mtcars$cyl))
Terms:
as.factor(mtcars$cyl) Residuals
Sum of Squares 824.7846 301.2626
Deg. of Freedom 2 29
Exactly the same result as testing with an ANOVA the linear model with only the categorical cylinder as regressor:
fit <- lm(mpg ~ as.factor(cyl), data = mtcars)
summary(fit)
anova(fit)
Analysis of Variance Table
Response: mpg
Df Sum Sq Mean Sq F value Pr(>F)
as.factor(cyl) 2 824.78 412.39 39.697 4.979e-09 ***
Residuals 29 301.26 10.39
What we see, then, is that the residuals - the part of the total variance not explained by the model - as well as the variance are the same whether you call an OLS of the type lm(DV ~ factors), or an ANOVA (aov(DV ~ factors)): when we strip the model of continuous variables we end up with an identical system. Similarly, when we evaluate the models globally or as an omnibus ANOVA (not level by level), we naturally get the same p-value F-statistic: 39.7 on 2 and 29 DF, p-value: 4.979e-09.
This is not to imply that the testing of individual levels is going to yield identical p-values. In the case of OLS, we can invoke summary(fit) and get:
lm(formula = mpg ~ as.factor(cyl), data = mtcars)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 26.6636 0.9718 27.437 < 2e-16 ***
as.factor(cyl)6 -6.9208 1.5583 -4.441 0.000119 ***
as.factor(cyl)8 -11.5636 1.2986 -8.905 8.57e-10 ***
This is not possible in ANOVA, which is more of an omnibus test. To get these types of $p$-value assessments we need to run a Tukey Honest Significant Difference test, which will try to reduce the possibility of a type I error as a result of performing multiple pairwise comparisons (hence, "p adjusted"), resulting in a completely different output:
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = mtcars$mpg ~ as.factor(mtcars$cyl))
$`as.factor(mtcars$cyl)`
diff lwr upr p adj
6-4 -6.920779 -10.769350 -3.0722086 0.0003424
8-4 -11.563636 -14.770779 -8.3564942 0.0000000
8-6 -4.642857 -8.327583 -0.9581313 0.0112287
Ultimately, nothing is more reassuring than taking a peek at the engine under the hood, which is none other than the model matrices, and the projections in the column space. These are actually quite simple in the case of an ANOVA:
$$\small\begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ \vdots \\\vdots\\\vdots\\.\\y_n \end{bmatrix} = \begin{bmatrix} \color{magenta} 1 & 0 & 0 \\ \color{magenta}1 & 0 & 0 \\ \vdots & \vdots & \vdots \\ \color{magenta} 0 & 1 & 0 \\ \color{magenta}0 & 1 & 0 \\ \vdots & \vdots & \vdots \\ .&.&.\\\color{magenta} 0 & 0 & 1 \\ \color{magenta}0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} \mu_1\\ \mu_2\\ \mu_3 \end{bmatrix} +\begin{bmatrix} \varepsilon_1 \\ \varepsilon_2\\ \varepsilon_3\\ \vdots\\ \vdots\\ \vdots\\ .\\ \varepsilon_n \end{bmatrix}\tag 1$$
This would be the one-way ANOVA model matrix with three levels (e.g. cyl 4, cyl 6, cyl 8), summarized as $\small y_{ij} = \mu_i + \epsilon_{ij}$, where $\mu_i$ is the mean at each level or group: when the error or residual for the observation $j$ of the group or level $i$ is added, we obtain the actual DV $y_{ij}$ observation.
On the other hand, the model matrix for an OLS regression is:
$$\small\begin{bmatrix}y_1 \\ y_2 \\ y_3 \\ y_4 \\ \vdots \\ y_n \end{bmatrix} = \begin{bmatrix} 1 & x_{12} & x_{13}\\ 1 & x_{22} & x_{23} \\ 1 & x_{32} & x_{33} \\ 1 & x_{42} & x_{43} \\ \vdots & \vdots & \vdots \\1 & x_{n2} & x_{n3} \end{bmatrix} \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \end{bmatrix} + \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \\ \varepsilon_4 \\ \vdots \\ \varepsilon_n \end{bmatrix}$$
This is of the form $ \small y_i = \beta_0 + \beta_1\, x_{i1} + \beta_2\, x_{i2} + \epsilon_i $ with a single intercept $\beta_0$ and two slopes ($\beta_1$ and $\beta_2$) each for a different continuous variables, say weight and displacement.
The trick now is to see how we can create different intercepts, as in the initial example, lm(mpg ~ wt + as.factor(cyl), data = mtcars) - so let's get rid of the second slope and stick to the original single continuous variable weight (in other words, one single column besides the column of ones in the model matrix; the intercept $\beta_0$ and the slope for weight, $\beta_1$). The column of $\color{brown}1$'s will by default correspond to the cyl 4 intercept. Again, its value is not identical to the ANOVA within-group mean for cyl 4, an observation that shouldn't be surprising comparing the column of $\color{brown}1$'s in the OLS model matrix (below) to the first column of $\color{magenta}1$'s in the ANOVA model matrix $(1),$ which only selects examples with 4-cylinders. The intercept will be shifted via dummy coding to explain the effect of cyl 6 and cyl 8as follows:
$$\small\begin{bmatrix}y_1 \\ y_2 \\ y_3 \\ y_4\\ y_5 \\ \vdots \\ y_n\end{bmatrix} = \begin{bmatrix} \color{brown}1 & x_1 \\ \color{brown}1 & x_2 \\\color{brown} 1 & x_3 \\ \color{brown}1 & x_4 \\ \color{brown}1 & x_5 \\ \vdots & \vdots \\\color{brown}1 & x_n \end{bmatrix} \begin{bmatrix} \beta_0 \\ \beta_1 \end{bmatrix}+ \begin{bmatrix}\color{red}1&0\\\color{red}1&0\\\color{red}1&0\\0&\color{blue}1\\0&\color{blue}1\\ \vdots & \vdots\\0&\color{blue}1\end{bmatrix} \begin{bmatrix} \tilde\mu_2 \\ \tilde\mu_3 \end{bmatrix} + \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \\ \varepsilon_4 \\ \varepsilon_5\\ \vdots \\ \varepsilon_n \end{bmatrix}$$
Now when the third column is $\color{red}1$ we'll be systematically shifting the intercept by $\tilde\mu_2.$ The $\tilde\cdot$ indicates that, as in the case of the "baseline" intercept in the OLS model not being identical to the group mean of 4-cylinder cars, but reflecting it, the differences between levels in the OLS model are not mathematically the between-groups differences in means:
fit <- lm(mpg ~ wt + as.factor(cyl), data = mtcars)
summary(fit)$coef[3] #[1] -4.255582 (difference between intercepts cyl==4 and cyl==6 in OLS)
fit <- lm(mpg ~ as.factor(cyl), data = mtcars)
summary(fit)$coef[2] #[1] -6.920779 (difference between group mean cyl==4 and cyl==6)
Likewise, when the fourth column is $\color{blue}1$, a fixed value $\tilde\mu_3$ will added to the intercept. The matrix equation, hence, will be $\small y_i = \beta_0 + \beta_1\, x_i + \tilde\mu_i + \epsilon_i $. Therefore, going with this model to the ANOVA model is just a matter of getting rid of the continuous variables, and understanding that the default intercept in OLS reflects the first level in ANOVA.
Antoni Parellada and usεr11852 had very good answer. I will address your question for coding perspective with R.
ANOVA tells you nothing about the coefficients of the linear model. So how is linear regression the same as ANOVA?
In fact, we can aov function in R can be used as same as lm. Here are some examples.
> lm_fit=lm(mpg~as.factor(cyl),mtcars)
> aov_fit=aov(mpg~as.factor(cyl),mtcars)
> coef(lm_fit)
(Intercept) as.factor(cyl)6 as.factor(cyl)8
26.663636 -6.920779 -11.563636
> coef(aov_fit)
(Intercept) as.factor(cyl)6 as.factor(cyl)8
26.663636 -6.920779 -11.563636
> all(predict(lm_fit,mtcars)==predict(aov_fit,mtcars))
[1] TRUE
As you can see, not only we can get coefficient from ANOVA model, but also we can use it for prediction, just like the linear model.
If we check the help file for aov function it says
This provides a wrapper to lm for fitting linear models to balanced or unbalanced experimental designs. The main difference from lm is in the way print, summary and so on handle the fit: this is expressed in the traditional language of the analysis of variance rather than that of linear models.
The analysis of variance (ANOVA) is a method that occurs within regression models. The technique gives a particular set of outputs for the model that analyse the estimated variance of different parts, and use this to make inferences about whether or not there are relationships between the explanatory variables and the response variable. Comparison of linear regression to ANOVA is an "apples and oranges" comparison, since the former is a model and the latter is a method of analysis that occurs within a model. Sometimes you will see references to ANOVA as if it were a model, in which case there must be some underlying model to which the method is applied.
The ANOVA method is fairly simple to understand when you look at it holistically in the context of a general regression model. The technique is based on the law of iterated variance. Suppose you are working in the context of some regression model:
$$Y_i = f(\mathbf{x}_i, \theta) + \varepsilon_i \quad \quad \quad \varepsilon_1,...,\varepsilon_n \sim \text{IID Dist}(0, \sigma^2).$$
Using the law of iterated variance we can write the marginal variance of $Y_i$ as:
$$\begin{equation} \begin{aligned} \mathbb{V}(Y_i) &= \mathbb{V}(\mathbb{E}(Y_i|\mathbf{X}_i)) + \mathbb{E}(\mathbb{V}(Y_i|\mathbf{X}_i)) \\[6pt] &= \mathbb{V}(f(\mathbf{X}_i, \theta)) + \mathbb{E}(\mathbb{V}(\varepsilon_i)) \\[6pt] &= \mathbb{V}(f(\mathbf{X}_i, \theta)) + \mathbb{E}(\sigma^2) \\[6pt] &= \mathbb{V}(f(\mathbf{X}_i, \theta)) + \sigma^2. \\[6pt] \end{aligned} \end{equation}$$
Now, if the explanatory vector does not have a relationship with the response variable then the regression function does not depend on the explanatory vector, and so $\mathbb{V}(f(\mathbf{X}_i, \theta)) = 0$, which implies $\mathbb{V}(Y_i) = \sigma^2$. On the other hand, if the explanatory vector does have a relationship with the response variable, then we will generally have $\mathbb{V}(f(\mathbf{X}_i, \theta)) > 0$, which implies $\mathbb{V}(Y_i) > \sigma^2$. Thus, generally speaking, a larger gap between the estimated variance of the response variable, and the estimated variance of the error term, constitutes evidence in favour of the hypothesis that there is a relationship between the explanatory vector and the response variable.
ANOVA can be taken further than this by breaking down the variance contributions of each of the explanatory variables in the model (or groups of these variables, etc.) to allow you to further test whether there is a plausible relationship between the response variable and an individual explanatory variable or group of explanatory variables. The lovely answer by Antoni Parellada shows you a colourful illustration of estimating the variance in three groups from a categorical explanatory variable.
The above decomposition from the law of iterated variance is the basic insight that underlies ANOVA. It is used to construct formal ANOVA tests to determine whether or not there is evidence of a relationship between the explanatory vector and the response variable. By conditioning on parts of the explanatory vector, this basic method can also be used to test for a relationship between particular subsets of explanatory variables and the response variable. In summary, ANOVA is a particular method that is used within the context of regression analysis to test for relationships between variables.
When is ANOVA "equivalent" to linear regression: The technique of ANOVA gives you a breakdown of the estimated variance of the components in the data, and this is often augmented with formal F-tests that use those estimated variance components. This is "equivalent" to performing the F-tests in a regression model, since that is what you are doing in ANOVA. (In simple linear regression there is only one explanatory variable so the F-test gives the same p-value as the T-test for this coefficient. In this case the ANOVA test is also equivalent to the individual T-test for the sole explanatory variable.)
You are correct that there are many aspects of the regression model that fall outside the scope of ANOVA (e.g., estimation of the coefficients for the individual explanatory variables). Again, that is because ANOVA is a method that occurs within the context of a regression model, not a model in its own right. When you see a whole model referred to as an "ANOVA model", it is more accurate to think of it as an ANOVA method applied to an underlying regression model.
If we take all data entries and arrange them into one single column Y, with the rest of the columns being indicator variables 1{ith data is element of the jth column in the original anova arrangement} then by taking a simple linear regression of Y on any of the other columns (say column B), you should obtain the same DF, SS, MS and F test statistic as in your ANOVA problem.
Thus ANOVA can be 'treated as' Linear Regression by writing the data with binary variables. Also note that the coefficient of regression for, say, a regression of Y on B should be the same as the avg. of the column B, computed with the original data.
