51
$\begingroup$

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?

$\endgroup$
46
+50
$\begingroup$

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.

$\endgroup$
  • 6
    $\begingroup$ The accepted answer on this Cross Validated page also shows the relation between ANOVA and regression pretty nicely, via a mathematical approach that nicely complements the practical approach of this answer. $\endgroup$ – EdM Oct 2 '15 at 21:12
  • $\begingroup$ +1. Oh yeah, @MichaelHardy's answer is quite good in that thread. Thanks for mentioning it! $\endgroup$ – usεr11852 says Reinstate Monic Oct 2 '15 at 21:15
  • $\begingroup$ +1, in addition, I feel this figure in this answer is really helpful to bridge the gap between ANOVA and linear regression $\endgroup$ – Haitao Du May 5 '17 at 15:44
  • $\begingroup$ Would you agree that ANOVA is a Gaussian GLM with categorical predictors? $\endgroup$ – Digio Mar 8 at 16:43
  • $\begingroup$ @Digio: No, it would oversimplify the appropriateness of their use; I would keep GLM out of picture. $\endgroup$ – usεr11852 says Reinstate Monic Mar 8 at 22:38
56
+100
$\begingroup$

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):

enter image description here

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.

$\endgroup$
  • 6
    $\begingroup$ +1, I love your graphical illustration!! publication quality! $\endgroup$ – Haitao Du May 5 '17 at 15:48
  • $\begingroup$ @hxd1011 This is very nice of you. I appreciate it. $\endgroup$ – Antoni Parellada May 5 '17 at 15:51
6
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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