ANOVA

Analysis of variance (ANOVA) is a powerful statistical technique. First-year graduate students in statistics are taught ANOVA almost immediately because of its importance, both theoretical and practical. We are often amazed, however, at the extent to which people outside the field are unaware of its purpose and value.

Regression creates a model, and ANOVA is one method of evaluating such models. The mathematics of ANOVA are intertwined with the mathematics of regression, so statisticians usually present them together; we follow that tradition here.

ANOVA is actually a family of techniques that are connected by a common mathematical analysis. This chapter mentions several applications:

One-way ANOVA

This is the simplest application of ANOVA. Suppose you have data samples from several populations and are wondering whether the populations have different means. One-way ANOVA answers that question. If the populations have normal distributions, use the oneway.test function (“Performing One-Way ANOVA”); otherwise, use the nonparametric version, the kruskal.test function (“Performing Robust ANOVA (Kruskal–Wallis Test)”).

Model comparison

When you add or delete a predictor variable from a linear regression, you want to know whether that change did or did not improve the model. The anova function compares two regression models and reports whether they are significantly different (“Comparing Models by Using ANOVA”).

ANOVA table

The anova function can also construct the ANOVA table of a linear regression model, which includes the F statistic needed to gauge the model’s statistical significance (“Getting Regression Statistics”). This important table is discussed in nearly every textbook on regression.

The See Also section below contain more about the mathematics of ANOVA.

Example Data

In many of the examples in this chapter, we start with creating example data using R’s pseudo random number generation capabilities. So at the beginning of each recipe you may see something like the following:

set.seed(42)
x <- rnorm(100)
e <- rnorm(100, mean=0, sd=5)
y <- 5 + 15 * x + e

We use set.seed to set the random number generation seed so that if you run the example code on your machine you will get the same answer. In the above example, x is a vector of 100 draws from a standard normal (mean=0, sd=1) distribution. Then we create a little random noise called e from a normal distribution with mean=0 and sd=5. y is then calculated as 5 + 15 * x + e. The idea behind creating example data rather than using “real world” data is that with simulated “toy” data you can change the coefficients and parameters in the example data and see how the change impacts the resulting model. For example, you could increase the standard deviation of e in the example data and see what impact that has on the R^2 of your model.

See Also

There are many good texts on linear regression. One of our favorites is Applied Linear Regression Models (4th ed.) by Kutner, Nachtsheim, and Neter (McGraw-Hill/Irwin). We generally follow their terminology and conventions in this chapter.

We also like Linear Models with R by Julian Faraway (Chapman & Hall), because it illustrates regression using R and is quite readable. Earlier versions of Faraday’s work are available free on-line, too (e.g., http://cran.r-project.org/doc/contrib/Faraway-PRA.pdf).

Problem

You have two vectors, x and y, that hold paired observations: (x₁, y₁), (x₂, y₂), …, (x_n, y_n). You believe there is a linear relationship between x and y, and you want to create a regression model of the relationship.

Solution

The lm function performs a linear regression and reports the coefficients:

set.seed(42)
x <- rnorm(100)
e <- rnorm(100, mean = 0, sd = 5)
y <- 5 + 15 * x + e

lm(y ~ x)
#>
#> Call:
#> lm(formula = y ~ x)
#>
#> Coefficients:
#> (Intercept)            x
#>        4.56        15.14

Discussion

Simple linear regression involves two variables: a predictor (or independent) variable, often called x; and a response (or dependent) variable, often called y. The regression uses the ordinary least-squares (OLS) algorithm to fit the linear model:

• y_i = β₀ + β₁x_i + ε_i

where β₀ and β₁ are the regression coefficients and the ε_i are the error terms.

The lm function can perform linear regression. The main argument is a model formula, such as y ~ x. The formula has the response variable on the left of the tilde character (~) and the predictor variable on the right. The function estimates the regression coefficients, β₀ and β₁, and reports them as the intercept and the coefficient of x, respectively:

Coefficients:
(Intercept)            x
      4.558       15.136

In this case, the regression equation is:

• y_i = 17.72 + 3.25x_i + ε_i

It is quite common for data to be captured inside a data frame, in which case you want to perform a regression between two data frame columns. Here, x and y are columns of a data frame dfrm:

df <- data.frame(x, y)
head(df)
#>        x     y
#> 1  1.371 31.57
#> 2 -0.565  1.75
#> 3  0.363  5.43
#> 4  0.633 23.74
#> 5  0.404  7.73
#> 6 -0.106  3.94

The lm function lets you specify a data frame by using the data parameter. If you do, the function will take the variables from the data frame and not from your workspace:

lm(y ~ x, data = df)          # Take x and y from df
#>
#> Call:
#> lm(formula = y ~ x, data = df)
#>
#> Coefficients:
#> (Intercept)            x
#>        4.56        15.14

Problem

You have several predictor variables (e.g., u, v, and w) and a response variable y. You believe there is a linear relationship between the predictors and the response, and you want to perform a linear regression on the data.

Solution

Use the lm function. Specify the multiple predictors on the righthand side of the formula, separated by plus signs (+):

lm(y ~ u + v + w)

Discussion

Multiple linear regression is the obvious generalization of simple linear regression. It allows multiple predictor variables instead of one predictor variable and still uses OLS to compute the coefficients of a linear equation. The three-variable regression just given corresponds to this linear model:

• y_i = β₀ + β₁u_i + β₂v_i + β₃w_i + ε_i

R uses the lm function for both simple and multiple linear regression. You simply add more variables to the righthand side of the model formula. The output then shows the coefficients of the fitted model:

set.seed(42)
u <- rnorm(100)
v <- rnorm(100, mean = 3,  sd = 2)
w <- rnorm(100, mean = -3, sd = 1)
e <- rnorm(100, mean = 0,  sd = 3)

y <- 5 + 4 * u + 3 * v + 2 * w + e

lm(y ~ u + v + w)
#>
#> Call:
#> lm(formula = y ~ u + v + w)
#>
#> Coefficients:
#> (Intercept)            u            v            w
#>        4.77         4.17         3.01         1.91

The data parameter of lm is especially valuable when the number of variables increases, since it’s much easier to keep your data in one data frame than in many separate variables. Suppose your data is captured in a data frame, such as the df variable shown here:

df <- data.frame(y, u, v, w)
head(df)
#>       y      u     v     w
#> 1 16.67  1.371 5.402 -5.00
#> 2 14.96 -0.565 5.090 -2.67
#> 3  5.89  0.363 0.994 -1.83
#> 4 27.95  0.633 6.697 -0.94
#> 5  2.42  0.404 1.666 -4.38
#> 6  5.73 -0.106 3.211 -4.15

When we supply df to the data parameter of lm, R looks for the regression variables in the columns of the data frame:

lm(y ~ u + v + w, data = df)
#>
#> Call:
#> lm(formula = y ~ u + v + w, data = df)
#>
#> Coefficients:
#> (Intercept)            u            v            w
#>        4.77         4.17         3.01         1.91

See Also

See “Performing Simple Linear Regression” for simple linear regression.

Problem

You want the critical statistics and information regarding your regression, such as R², the F statistic, confidence intervals for the coefficients, residuals, the ANOVA table, and so forth.

Solution

Save the regression model in a variable, say m:

m <- lm(y ~ u + v + w)

Then use functions to extract regression statistics and information from the model:

anova(m)

ANOVA table

coefficients(m)

Model coefficients

coef(m)

Same as coefficients(m)

confint(m)

Confidence intervals for the regression coefficients

deviance(m)

Residual sum of squares

effects(m)

Vector of orthogonal effects

fitted(m)

Vector of fitted y values

residuals(m)

Model residuals

resid(m)

Same as residuals(m)

summary(m)

Key statistics, such as R², the F statistic, and the residual standard error (σ)

vcov(m)

Variance–covariance matrix of the main parameters

Discussion

When we started using R, the documentation said use the lm function to perform linear regression. So we did something like this, getting the output shown in “Performing Multiple Linear Regression”:

lm(y ~ u + v + w)
#>
#> Call:
#> lm(formula = y ~ u + v + w)
#>
#> Coefficients:
#> (Intercept)            u            v            w
#>        4.77         4.17         3.01         1.91

How disappointing! The output was nothing compared to other statistics packages such as SAS. Where is R²? Where are the confidence intervals for the coefficients? Where is the F statistic, its p-value, and the ANOVA table?

Of course, all that information is available—you just have to ask for it. Other statistics systems dump everything and let you wade through it. R is more minimalist. It prints a bare-bones output and lets you request what more you want.

The lm function returns a model object that you can assign to a variable:

m <- lm(y ~ u + v + w)

From the model object, you can extract important information using specialized functions. The most important function is summary:

summary(m)
#>
#> Call:
#> lm(formula = y ~ u + v + w)
#>
#> Residuals:
#>    Min     1Q Median     3Q    Max
#> -5.383 -1.760 -0.312  1.856  6.984
#>
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)
#> (Intercept)    4.770      0.969    4.92  3.5e-06 ***
#> u              4.173      0.260   16.07  < 2e-16 ***
#> v              3.013      0.148   20.31  < 2e-16 ***
#> w              1.905      0.266    7.15  1.7e-10 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 2.66 on 96 degrees of freedom
#> Multiple R-squared:  0.885,  Adjusted R-squared:  0.882
#> F-statistic:  247 on 3 and 96 DF,  p-value: <2e-16

The summary shows the estimated coefficients. It shows the critical statistics, such as R² and the F statistic. It shows an estimate of σ, the standard error of the residuals. The summary is so important that there is an entire recipe devoted to understanding it (“Understanding the Regression Summary”).

There are specialized extractor functions for other important information:

Model coefficients (point estimates)

    coef(m)
#> (Intercept)           u           v           w
#>        4.77        4.17        3.01        1.91

Confidence intervals for model coefficients

    confint(m)
#>             2.5 % 97.5 %
#> (Intercept)  2.85   6.69
#> u            3.66   4.69
#> v            2.72   3.31
#> w            1.38   2.43

Model residuals

    resid(m)
#>       1       2       3       4       5       6       7       8       9
#> -0.5675  2.2880  0.0972  2.1474 -0.7169 -0.3617  1.0350  2.8040 -4.2496
#>      10      11      12      13      14      15      16      17      18
#> -0.2048 -0.6467 -2.5772 -2.9339 -1.9330  1.7800 -1.4400 -2.3989  0.9245
#>      19      20      21      22      23      24      25      26      27
#> -3.3663  2.6890 -1.4190  0.7871  0.0355 -0.3806  5.0459 -2.5011  3.4516
#>      28      29      30      31      32      33      34      35      36
#>  0.3371 -2.7099 -0.0761  2.0261 -1.3902 -2.7041  0.3953  2.7201 -0.0254
#>      37      38      39      40      41      42      43      44      45
#> -3.9887 -3.9011 -1.9458 -1.7701 -0.2614  2.0977 -1.3986 -3.1910  1.8439
#>      46      47      48      49      50      51      52      53      54
#>  0.8218  3.6273 -5.3832  0.2905  3.7878  1.9194 -2.4106  1.6855 -2.7964
#>      55      56      57      58      59      60      61      62      63
#> -1.3348  3.3549 -1.1525  2.4012 -0.5320 -4.9434 -2.4899 -3.2718 -1.6161
#>      64      65      66      67      68      69      70      71      72
#> -1.5119 -0.4493 -0.9869  5.6273 -4.4626 -1.7568  0.8099  5.0320  0.1689
#>      73      74      75      76      77      78      79      80      81
#>  3.5761 -4.8668  4.2781 -2.1386 -0.9739 -3.6380  0.5788  5.5664  6.9840
#>      82      83      84      85      86      87      88      89      90
#> -3.5119  1.2842  4.1445 -0.4630 -0.7867 -0.7565  1.6384  3.7578  1.8942
#>      91      92      93      94      95      96      97      98      99
#>  0.5542 -0.8662  1.2041 -1.7401 -0.7261  3.2701  1.4012  0.9476 -0.9140
#>     100
#>  2.4278

Residual sum of squares

    deviance(m)
#> [1] 679

ANOVA table

    anova(m)
#> Analysis of Variance Table
#>
#> Response: y
#>           Df Sum Sq Mean Sq F value  Pr(>F)
#> u          1   1776    1776   251.0 < 2e-16 ***
#> v          1   3097    3097   437.7 < 2e-16 ***
#> w          1    362     362    51.1 1.7e-10 ***
#> Residuals 96    679       7
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

If you find it annoying to save the model in a variable, you are welcome to use one-liners such as this:

summary(lm(y ~ u + v + w))

Or you can use Magritr pipes:

lm(y ~ u + v + w) %>%
  summary

See Also

See “Understanding the Regression Summary”. See “Identifying Influential Observations” for regression statistics specific to model diagnostics.

Problem

You created a linear regression model, m. However, you are confused by the output from summary(m).

Discussion

The model summary is important because it links you to the most critical regression statistics. Here is the model summary from “Getting Regression Statistics”:

summary(m)
#>
#> Call:
#> lm(formula = y ~ u + v + w)
#>
#> Residuals:
#>    Min     1Q Median     3Q    Max
#> -5.383 -1.760 -0.312  1.856  6.984
#>
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)
#> (Intercept)    4.770      0.969    4.92  3.5e-06 ***
#> u              4.173      0.260   16.07  < 2e-16 ***
#> v              3.013      0.148   20.31  < 2e-16 ***
#> w              1.905      0.266    7.15  1.7e-10 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 2.66 on 96 degrees of freedom
#> Multiple R-squared:  0.885,  Adjusted R-squared:  0.882
#> F-statistic:  247 on 3 and 96 DF,  p-value: <2e-16

Let’s dissect this summary by section. We’ll read it from top to bottom—even though the most important statistic, the F statistic, appears at the end:

Call

    summary(m)$call

This shows how lm was called when it created the model, which is important for putting this summary into the proper context.

Residuals statistics

    # Residuals:
    #     Min      1Q  Median      3Q     Max
    # -5.3832 -1.7601 -0.3115  1.8565  6.9840

Ideally, the regression residuals would have a perfect, normal distribution. These statistics help you identify possible deviations from normality. The OLS algorithm is mathematically guaranteed to produce residuals with a mean of zero.[‸1] Hence the sign of the median indicates the skew’s direction, and the magnitude of the median indicates the extent. In this case the median is negative, which suggests some skew to the left.

If the residuals have a nice, bell-shaped distribution, then the first quartile (1Q) and third quartile (3Q) should have about the same magnitude. In this example, the larger magnitude of 3Q versus 1Q (1.3730 versus 0.9472) indicates a slight skew to the right in our data, although the negative median makes the situation less clear-cut.

The Min and Max residuals offer a quick way to detect extreme outliers in the data, since extreme outliers (in the response variable) produce large residuals.

Coefficients

summary(m)$coefficients
#>             Estimate Std. Error t value Pr(>|t|)
#> (Intercept)     4.77      0.969    4.92 3.55e-06
#> u               4.17      0.260   16.07 5.76e-29
#> v               3.01      0.148   20.31 1.58e-36
#> w               1.91      0.266    7.15 1.71e-10

The column labeled Estimate contains the estimated regression coefficients as calculated by ordinary least squares.

Theoretically, if a variable’s coefficient is zero then the variable is worthless; it adds nothing to the model. Yet the coefficients shown here are only estimates, and they will never be exactly zero. We therefore ask: Statistically speaking, how likely is it that the true coefficient is zero? That is the purpose of the t statistics and the p-values, which in the summary are labeled (respectively) t value and Pr(>|t|).

The p-value is a probability. It gauges the likelihood that the coefficient is not significant, so smaller is better. Big is bad because it indicates a high likelihood of insignificance. In this example, the p-value for the u coefficient is a mere 0.00106, so u is likely significant. The p-value for w, however, is 0.05744; this is just over our conventional limit of 0.05, which suggests that w is likely insignificant.[^2] Variables with large p-values are candidates for elimination.

A handy feature is that R flags the significant variables for quick identification. Do you notice the extreme righthand column containing double asterisks (**), a single asterisk (*), and a period(.)? That column highlights the significant variables. The line labeled "Signif. codes" at the bottom gives a cryptic guide to the flags’ meanings:

+

--------- ---------------------------------- *** p-value between 0 and 0.001 ** p-value between 0.001 and 0.01 * p-value between 0.01 and 0.05 . p-value between 0.05 and 0.1 (blank) p-value between 0.1 and 1.0 --------- ----------------------------------

+

The column labeled Std. Error is the standard error of the estimated coefficient. The column labeled t value is the t statistic from which the p-value was calculated.

Residual standard error::
+
[source, r]

# Residual standard error: 2.66 on 96 degrees of freedom

+
-------------------------------------------------------------------
This reports the standard error of the residuals (*σ*)—that is, the
sample standard deviation of *ε*.
-------------------------------------------------------------------

_R_^2^ (coefficient of determination)::
+
[source, r]

# Multiple R-squared:  0.8851,  Adjusted R-squared:  0.8815

+
-------------------------------------------------------------------
*R*^2^ is a measure of the model’s quality. Bigger is better.
Mathematically, it is the fraction of the variance of *y* that is
explained by the regression model. The remaining variance is not
explained by the model, so it must be due to other factors (i.e.,
unknown variables or sampling variability). In this case, the model
explains 0.4981 (49.81%) of the variance of *y*, and the remaining
0.5019 (50.19%) is unexplained.

That being said, we strongly suggest using the adjusted rather than
the basic *R*^2^. The adjusted value accounts for the number of
variables in your model and so is a more realistic assessment of
its effectiveness. In this case, then, we would use 0.8815,
not 0.8851s
-------------------------------------------------------------------

_F_ statistic::
+
[source, r]

# F-statistic: 246.6 on 3 and 96 DF,  p-value: < 2.2e-16

+
--------------------------------------------------------------------
The *F* statistic tells you whether the model is significant
or insignificant. The model is significant if any of the
coefficients are nonzero (i.e., if *β*~*i*~ ≠ 0 for some *i*). It is
insignificant if all coefficients are zero (*β*~1~ = *β*~2~ = … =
*β*~*n*~ = 0).

Conventionally, a *p*-value of less than 0.05 indicates that the
model is likely significant (one or more *β*~*i*~ are nonzero)
whereas values exceeding 0.05 indicate that the model is likely
not significant. Here, the probability is only 0.000391 that our
model is insignificant. That’s good.

Most people look at the *R*^2^ statistic first. The statistician
wisely starts with the *F* statistic, for if the model is not
significant then nothing else matters.
--------------------------------------------------------------------

[[see_also-id240]]
==== See Also

See <<recipe-id231>> for more on extracting statistics and information from the
model object.

[[recipe-id205]]
=== Performing Linear Regression Without an Intercept

[[problem-id205]]
==== Problem

You want to perform a linear regression, but you want to force the
intercept to be zero.

[[solution-id205]]
==== Solution

Add "`+` `0`" to the righthand side of your regression formula. That
will force `lm` to fit the model with a zero intercept:

[source, r]

lm(y ~ x + 0)

The corresponding regression equation is:

++++
<ul class="simplelist">
  <li><em>y</em><sub><em>i</em></sub> = <em>βx</em><sub><em>i</em></sub> + <em>ε</em><sub><em>i</em></sub></li>
</ul>
++++

[[discussion-id205]]
==== Discussion

Linear regression ordinarily includes an intercept term, so that is the
default in R. In rare cases, however, you may want to fit the data while
assuming that the intercept is zero. In this you make a modeling
assumption: when _x_ is zero, _y_ should be zero.

When you force a zero intercept, the `lm` output includes a coefficient
for _x_ but no intercept for _y_, as shown here:

[source, r]

lm(y _{x + 0) #> #> Call: #> lm(formula = y} x + 0) #> #> Coefficients: #> x #> 4.3

We strongly suggest you check that modeling assumption before
proceeding. Perform a regression with an intercept; then see if the
intercept could plausibly be zero. Check the intercept’s confidence
interval. In this example, the confidence interval is (6.26, 8.84):

[source, r]

confint(lm(y ~ x)) #> 2.5 % 97.5 % #> (Intercept) 6.26 8.84 #> x 2.82 5.31

Because the confidence interval does not contain zero, it is NOT
statistically plausible that the intercept could be zero. So in this
case, it is not reasonable to rerun the regression while forcing a zero
intercept.

[[title-highcor]]
=== Regressing Only Variables that Highly Correlate with your Dependent Variable

[[problem-highcor]]
==== Problem

You have a data frame with many variables and you want to build a
multiple linear regression using only the variables that are highly
correlated to your response (dependent) variable.

[[solution-highcor]]
==== Solution

If `df` is our data frame containing both our response (dependent) and
all our predictor (independent) variables and `dep_var` is our response
variable, we can figure out our best predictors and then use them in a
linear regression. If we want the top 4 predictor variables, we can use
this recipe:

[source, r]

best_pred ← df %>% select(-dep_var) %>% map_dbl(cor, y = df$dep_var) %>% sort(decreasing = TRUE) %>% .[1:4] %>% names %>% df[.]

mod ← lm(df$dep_var ~ as.matrix(best_pred))

This recipe is a combination of many differnt pieces of logic used
elsewhere in this book. We will describe each step here then walk
through it in the discussion using some example data.

First we drop the response variable out of our pipe chain so that we
have only our predictor variables in our data flow:

[source, r]

df %>% select(-dep_var)

Then we use `map_dbl` from `purrr` to perform a pairwise correlation on
each column relative to the response variable.

[source, r]

map_dbl(cor, y = df$dep_var) %>%

We then take the resulting correlations and sort them in decreasing
order:

[source, r]

sort(decreasing = TRUE) %>%

We want only the top 4 correlated variables so we select the top 4
records in the resulting vector:

[source, r]

.[1:4] %>%

And we don't need the correlation values, only the names of the rows
which are the variable names from our original data frame `df`:

[source, r]

names %>%

Then we can pass those names into our subsetting brackets to select only
the columns with names matching the ones we want:

[source, r]