# 1-way ANOVA {#ANOVA}

Comparing a continuous variable across more than two levels of a categorical variable 

## Introduction to ANOVA

In linear regression we were concerned with fitting models relating a continuous response variable $Y$ to another continuous variable $X$, and we saw that we can use either the $t$--test or the $F$--test to find out if taking $X$ into account adds any information. But what if the $X$ variable is not intrinsically continuous?! 

If $X$ is a categorical variable with two levels, then the solution is to use the two--sample $t$--test described in chapter _sec:two_sample_t_tests_. In this chapter, we'd like to do the trick for $X$ variables that have more than two categories. The method of choice in this case is called ANOVA which is short for analysis of variance. We will see that ANOVA is just a special case of the $F$--test in linear regression. From a different angle, it can also be seen as a natural generalisation of the two--sample $t$--test -- so relax and enjoy, there's not much in this chapter that you don't know already!!

<!-- , fig.cap = "Scatterplot and OLS solution for Eye--sight score and Age in our Martian sample from site I."  -->


```{r scatterploteyesightage, fig.cap="XX", fig.width = large_width, echo = FALSE}
library(tidyverse)
martian <- read.delim("data/martian.txt")

lm_orig <- lm(Eye ~ Age, data = martian[martian$Site == "Site I",])
mars.lm <- summary(lm_orig)
mars.anova.eye <- anova(lm_orig)

martian %>% 
  filter(Site == "Site I") %>%
  ggplot(aes(x = Age, y = Eye)) +
    geom_point(col = DarkYellow, alpha = 0.6) +
    stat_smooth(method = "lm", se = FALSE, col = "#C42126") +
    scale_x_continuous("Age (years)", expand = c(0, 0)) +
    scale_y_continuous("Eye-sight (score)", expand = c(0, 0)) +
    coord_cartesian(xlim = c(0, 100), ylim = c(0, 100)) +
    theme(legend.position = c(0.8, 0.8),
      panel.background = element_blank(),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      axis.ticks.x = element_line(colour = "black"),
      axis.ticks.y = element_line(colour = "black"),
      axis.text.x = element_text(colour = "black"),
      axis.text.y = element_text(colour = "black"),
      axis.line = element_line(colour = "black")
)
```


Let's first emphasise the connection between linear regression and ANOVA. For this we'll use yet another linear regression example that we will then turn into a ANOVA--type problem. The question we'll be investigating here is: _Does Martian eye--sight depend on their age?_ For now, both variables (eye--sight and age) are continuous.

Figure \@ref(fig:scatterploteyesightage) shows a scatter plot of the two variables, together with the regression line we get from OLS. As we know now this has the form $f(X) = \hat \beta_0 + \hat \beta_1 X$, where $\hat \beta_0$ and  $\hat \beta_1$ are the OLS estimates (determined by our data) of the underlying parameters  $\beta_0$ and  $\beta_1$ that describe the linear relationship between $X$ and $Y$ in the population. As in the previous chapter, our main interest here is to determine whether or not the observed value $\hat \beta_1$ is big enough relative to the noise in the data so that we can reject the null hypothesis $H_0: \beta_1 = 0$, i.e. the scenario where looking at age does not help when trying to predict the eye-sight of a Martian.

We probably all agree the scatterplot shows a lovely linear relationship between the two variables and that a lot can be said about the eye-sight of a Martian once we know their age. The linear regression output shown in Table \@ref(tab:AgeEyeOLS) confirms this suspicion. The $p$--value for age is very small!



```{r AgeEyeOLS, echo = FALSE}

# broom::tidy(mars.lm) %>% 
broom::tidy(lm_orig) %>%   
  mutate(term = c("$\\beta_0$ (intercept)", "$\\beta_1$ (slope)")) %>% 
  knitr::kable(col.names = c("Term", "Estimate", "Standard Error", "Test statistic", "p-value"),
               caption = "Summary of the OLS model for Eye--sight as described by Age.")

```

Even before we formally carry out the $t$--test, we can guess that:

1. the estimated value for $\beta_1$ will be some negative number
2. it will be large compared to it's standard deviation
3. the $p$--value for $H_0: \beta_1 = 0$ will be small (this is a consequence of 1 and 2)

```{exercise, name="Understanding your data"}
Explain why it's easy to guess these things.
```


Table \@ref(tab:AgeEyeOLS) also tells us that our best guess for the linear relationship between age and eye-sight has the form

$$ \text{(predicted eye-sight score when age} = x) = f(x) = 126.9 - 1.4 x,$$

i.e. for each year of age, the expected score for eye-sight decreases on average by 1.4 points. 

<!-- % Also, from the associated $p$--value we know that the scatter around the OLS regression line is small enough to reject the hypothesis that our estimate $\hat \beta_1 = -1.4$ could just have occurred by chance. So we can confidently say that looking at age helps us to explain the variation which is going on in eye-sight. -->

As we learned in the previous chapter, the $F$--test provides another way to look at the problem. As opposed to the $t$--test, it doesn't provide us with a formula for (our best guess of) the linear relationship between age and eye-sight, but it does look at the question of whether or not the variable age adds information when we try to explain the variation in eye-sight.



Table table below provides us with the output of the $F$--test for our example. As expected, we see a giant value for the $F$--statistic. This means that the residuals of the OLS solution are significantly smaller than the residuals for the horizontal null model. In other words, the improvement we get from taking age into account is a lot bigger than we'd expect if age didn't play a role.


```{r echo = FALSE}

broom::tidy(mars.anova.eye) %>% 
  mutate(term = c("Model", "Residuals")) %>% 
  knitr::kable(col.names = c("Term", "Degrees of freedom", "Sum of the Squared...", "Mean of the Squared...", "Test statistic", "p-value"),
               caption = "Summary of the F--test applied to Eye-sight described by Age.", digits = 2)

```

Now imagine that age was an ordinal instead of a continuous variable that consists of the categories youth, adult and senior.^[It is always possible to turn a continuous variable into an ordinal variable by defining consecutive categories.] Figure \@ref(fig:tripleAgeEyePlot) shows the new situation where $X$ is a categorical variable with three levels. So far we have covered continuous $X$ and categorical $X$ with two levels in which cases you'd use linear regression and the two--sample $t$--test, respectively. Our new situation is sighly different and corresponds to a classical situation in which you'd use ANOVA!

```{r tripleAgeEyePlot, fig.cap = "Transitioning from a continuous X variable to a categorical X variable with three levels.", echo = FALSE}

# Calculate summaries
martian %>% 
  dplyr::filter(Site == "Site I") %>% 
  group_by(EyeIndex) %>% 
  dplyr::summarise(mean = mean(Age),
            stdev = sd(Age)) -> eye.mean

martian %>% 
  dplyr::filter(Site == "Site I") %>% 
  group_by(AgeIndex) %>% 
  dplyr::summarise(mean = mean(Eye),
            stdev = sd(Eye)) -> age.mean

martian %>% 
  filter(Site == "Site I") %>%
  ggplot(aes(x = AgeIndex, y = Eye)) +
  geom_jitter(position = position_jitter(width = .2), col = DarkYellow, alpha = 0.6) +
  scale_y_continuous("Eye-sight (score)", expand = c(0, 0)) +
  scale_x_discrete("Age Group", limits = c("Youth", "Adult", "Senior")) +
  coord_cartesian(ylim = c(0, 100)) +
  geom_point(data = age.mean,
             aes(x = AgeIndex, y = mean),
             col = "#C42126") +
  geom_errorbar(data = age.mean,
                aes(x = AgeIndex,
                    ymax = mean + stdev,
                    ymin = mean - stdev),
                width = 0.1, col = "#C42126", inherit.aes = FALSE) +
  theme(legend.position = c(0.8, 0.8),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(colour = "black"),
        axis.ticks.y = element_line(colour = "black"),
        axis.text.x = element_text(colour = "black"),
        axis.text.y = element_text(colour = "black"),
        axis.line = element_line(colour = "black")
  )

```


<!-- % section correlation_on_martian_data_set (end) -->

<!-- % Is there a formal way to assess whether or not age \emph{group} has a significant impact on eye-sight? We obviously wouldn't have started a new chapter if there wasn't, so the answer is yes! And as we said, the technique that is appropriate for this type of question is called ANOVA. -->

Actually, ANOVA is an umbrella term for a family of methods. Here we will cover the simplest case which is called one--way ANOVA. One--way ANOVA answers the question if taking into account the levels of one categorical variable helps to explain the variation of a continuous variable.^[Other ways of saying `if taking into account the levels of a categorical variable helps to explain the variation of a continuous variable' are `if the values of a continuous variable significantly differ between the groups defined by a categorical variable' or `if the values of a categorical variable help us to predict the value of a continuous variable'.] The statistical test that runs in the background when we perform ANOVA is the $F$--test that we already encountered. Moreover, we will show that the same analysis can also be carried out in a linear regression framework. 

<!-- % As we said in the beginning of this chapter, -->
<!-- % -->
<!-- % The good news is that the test which is used to answer this question is the $F$--test that we heard about before. And as we will see, ANOVA can even be performed by fitting a linear model in the same way as we did in the previous chapter! It's just for historical reasons that ANOVA is mainly associated with categorical explanatory variables and $F$--tests whereas linear models and the corresponding $t$--tests are mainly thought of in the context of continuous explanatory variables. In the end, both cases can be written as a linear model and both $F$--test and $t$--test can be used in either situation! -->

## A simple example

Before we start looking at the actual Martian example, let's first examine a simple dataset to observe the underlying concepts. The raw data are shown in the left--hand side of Figure \ref{simpleNolines}. We have two variables, $X$ and $Y$, where $X$ is a categorical variable with three groups, and $Y$ is a continuous variable. 

Remember, we can analyse the question whether or not the levels of $X$ can help us explain the variation in $Y$ either in the linear regression framework or with the $F$--test.

### Linear regression with categorical explanatory variables

We first show how the data can be examined in a linear regression framework! In linear regression we were interested if we can find a line through the scatterplot of $Y$ against $X$ that better explains the data than the horizontal line that goes through the observed mean value $\bar{y}$. Here, we are asking the same question except that it's a bit meaningless to look for a straight line because $X$ is not continuous and can only take three distinct values. So instead of looking for a line, we are asking whether or not we can find three group--specific $y$--values that are distinct from each other and (in terms of the residuals) fit the data better than just representing each group by the overall mean $\bar{y}$. 

As you can guess, the group--specific $y$--values that will fit the data better  than the overall mean $\bar{y}$ are the group mean values $\bar{y}_1$, $\bar{y}_2$ and $\bar{y}_3$. However, the question is if this improvement is _significantly_ better, i.e. if going the extra mile (reporting three instead of one value) is worth the effort.

```{r simpleNolines, fig.cap = "A fake dataset to help us explore the concepts in ANOVA. The data consist of a continuous variable described by three groups of a categorical variable.", fig.width = large_width, echo = FALSE}

data.frame(group = rep(1:3,each = 5),
           value = c(3,2,1,1,4,
                     5,2,4,2,3,
                     7,4,5,3,6),
           unique = rep(1:5, 3)) %>% 
  group_by(group) %>% 
  mutate(avg = mean(value),
            sd = sd(value)) -> playANOVA

# First just draw the individual points:
m <- ggplot(playANOVA, aes(x = group, y = value, group = factor(unique), col = factor(group))) +
  guides(colour = FALSE) +
  geom_point(position =  position_dodge(.3), size = 4) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous("", limits = c(0,8), expand = c(0,0)) +
  scale_x_continuous("", limits = c(0.7,3.2), labels = c("Group 1", "Group 2", "Group 3"), breaks = 1:3) +
      theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_line(colour = "black"),
        axis.line = element_line(colour = "black"),
#        axis.line.x = element_blank(),
        axis.text = element_text(colour = "black")
    )

col.y <- "#F5B319"
col.r <- "#C42126"

playANOVA.lm <- lm(value ~ factor(group), data = playANOVA)
playANOVA.lm.summary <- summary(playANOVA.lm)
playANOVA.lm.anova <- anova(playANOVA.lm)

# Simple
q1 <- m + ggtitle("Data points")

# Grand mean line:
q2 <- m + geom_hline(yintercept = mean(playANOVA$value), col = col.r) +
ggtitle("Null Hypothesis")

# Individual mean lines:
q3 <- m + geom_segment(aes(x = group - .15, y = avg, xend = group + .15, yend = avg), color = col.y) +
ggtitle("Model")

plot_grid(q1, q2, q3, nrow = 1, align = "h")
```



Figure \@ref(fig:simpleNolines) displays the raw data and the two models we seek to compare. The middle plot shows the null model which is equal to the horizontal line through $\bar{y}$, and the right hand side shows the plot with the three individual group means which takes the role of the OLS regression line in our situation. When we write down the two models as a formula, the null model is pretty simple:

\begin{eqnarray*}
\hat{y}_i = \bar{y} & \text{ for all individuals $i$}.
\end{eqnarray*}

The model we want to assess is a little bit more complicated as we need to distinguish between the three groups

\begin{eqnarray*}
\hat{y}_i = \begin{cases} \bar{y}_1 &  \text{if $i$th individual belongs to group 1} \\
\bar{y}_2 &  \text{if $i$th individual belongs to group 2} \\
\bar{y}_3 &  \text{if $i$th individual belongs to group 3} . \end{cases}
\end{eqnarray*}

The distinction between the groups is a little annoying as the formula doesn't look nice and concise. But luckily, with a bit of effort, we can rewrite it in a more compressed form as

$$ \hat{y}_i = \bar{y}_1 + (\bar{y}_2 - \bar{y}_1) \cdot x^{g2}_i + (\bar{y}_3 - \bar{y}_1) \cdot x^{g3}_i, $$

where the weird looking symbols  $x^{g2}_i$ and $x^{g3}_i$ are so called _indicator variables_ that take the values 0 and 1 to signify whether or not the $i$th individual belongs to group 2 and group 3, respectively.

For example if the 6$th$ individual belongs to group 2, then $x^{g2}_6 = 1$ and $x^{g3}_6 = 0$ so that the whole equation becomes

\begin{eqnarray*}
\hat{y}_6 & = & \bar{y}_1 + (\bar{y}_2 - \bar{y}_1) \cdot 1 + (\bar{y}_3 - \bar{y}_1) \cdot 0 \\
 & = &  \bar{y}_1 + (\bar{y}_2 - \bar{y}_1) \\
 & = & \bar{y}_2.
\end{eqnarray*}

Similarly, if the 11$th$ individual belongs to group 3, then  $x^{g2}_{11} = 0$ and $x^{g2}_{11} = 1$ so that the whole equation becomes

\begin{eqnarray*}
\hat{y}_{11} & = & \bar{y}_1 + (\bar{y}_2 - \bar{y}_1) \cdot 0 + (\bar{y}_3 - \bar{y}_1) \cdot 1 \\
 & = &  \bar{y}_1 + (\bar{y}_3 - \bar{y}_1) \\
 & = & \bar{y}_3.
 \end{eqnarray*}

```{exercise, name="Defining our model"}
Write down the equation for the 3$rd$ individual, if it's in group 1.
```

If you understand the trick with the indicator variables, you actually understand the connection between ANOVA and linear models, because the only thing that's left is a tiny substitution: if we replace $\bar{y}_1$ by $\hat{\beta}_0$, $(\bar{y}_2 - \bar{y}_1)$  by $\hat{\beta}_1$ and $(\bar{y}_3 - \bar{y}_1)$  by $\hat{\beta}_2$, our nice compressed equation above becomes
$$\hat{y}_i = \hat{\beta}_0 +   \hat{\beta}_1 x^{g2}_i + \hat{\beta}_2  x^{g3}_i,$$

which is just what we'd expect a linear model to look like! And we even know the values of the estimated coefficients already! In our little toy example $\bar {y}_1 = 2.2$, $\bar y_2 = 3.2$ and $\bar y_3 = 5.0$, so that $\hat{\beta}_0 = 2.2$ $\hat{\beta}_1 = 1.0$ and $\hat{\beta}_3 = 2.8$. With this, the entire regression equation becomes


$$\hat{y}_i = 2.2 +   1.0 \cdot x^{g2}_i + 2.8 \cdot  x^{g3}_i. $$

Therefore the only thing left to close the circle with the linear regression as you know, is the $t$--test.

In the linear regression we were interested in the null hypothesis $H_0: \beta_1 = 0$ because if we can reject this based on the data, it means that the OLS solution with the non-zero slope fits the data better than the null (horizontal) line. Here, we want to do the same thing, except that we are interested in two null hypotheses, namely  $H_0: \beta_1 = 0$ AND $H_0: \beta_2 = 0$. If either one of them can be rejected based on the data, then we know that the group means are significantly different from each other. Hence, there is a significant improvement in taking group membership into account! 

```{r results = "asis", echo = FALSE}

# XT <- xtable::xtable(LibidoDose.lm.summary, caption = "Linear regression output for the toy data shown in Figure 15.3.")
# 
# print(XT, type = "html")

```

<!-- , label = "ANOVALM1" -->

Table _ANOVALM1_ shows the computer output that we get when we apply linear regression to our data. The first value is the intercept which corresponds to $\hat \beta_0$. The fact that it has a small $p$--value tells us that the mean of the first group is significantly different from zero. This is not the important point in our investigation because we are interested in the _differences_ between the group means which are given by $\hat \beta_1$ and $\hat \beta_2$. The $p$--value for $\hat \beta_2$ is quite large, so there is no evidence in the data that group 1 and group 2 have different mean values. However, the $p$--value for $\hat \beta_2$ is very small. As $\hat \beta_2$ measures the mean difference between group 1 and group 3, we have evidence that the observed difference is bigger than what we'd expect by pure chance! Therefore, we can conclude that the more complex linear model provides a better fit to the data than the null model.

```{exercise, name="A non-significant result"}
Can you draw a picture similar to Figure _simpleNolines_ were the data don't provide evidence against the null model?
```

```{exercise, name="A linear model with two levels in the predictor"}
IMPORTANT!! Do you think it's possible to write down a similar linear model when we only want to compare two groups?! What would happen?
```

## Classical ANOVA

We said in the beginning of this chapter, we can carry out the same investigation with an $F$--test. Doing so would be the classical approach to assess mean differences for more than two groups. So when people talk about (one--way) ANOVA, this is what they probably mean. Let's have a look at what the $F$--test does in our simple situation. Just like in the case where we compared two continuous variables against each other, the $F$--test works with the residuals of the two models under comparison.

```{r simplePlusResiduals, fig.cap="Three plots of our fake ANOVA dataset depicting SSR, SST and SSM, respectively. SSR considers the distance between each data point and the within--group mean, SST considers the distance between each data point and the grand mean of all data points, and SSM considers the distance between the within--group means and the grand mean.", fig.width = large_width, echo = FALSE}

q1 <- m + geom_errorbar(aes(ymin=value, ymax=avg), width=0, position =  position_dodge(.3), col = "black") +
  geom_segment(aes(x = group - .15, y = avg, xend = group + .15, yend = avg), color = col.y) +
  ggtitle("SSR")


q2 <- m + geom_errorbar(aes(ymin=value, ymax=mean(playANOVA$value)), width=0, position =  position_dodge(.3), col = "black") +
  geom_hline(yintercept = mean(playANOVA$value), col = col.r) +
  ggtitle("SST")

q3 <- m + geom_errorbar(aes(ymin=mean(playANOVA$value), ymax = avg), width=0, position =  position_dodge(.3), col = "black") +
  geom_hline(yintercept = mean(playANOVA$value), col = col.r) +
  geom_segment(aes(x = group - .15, y = avg, xend = group + .15, yend = avg), color = col.y) +
  ggtitle("SSM")

plot_grid(q1, q2, q3, nrow = 1, align = "h")

```



```{r ANOVASummary, echo = FALSE}

# CleanData <- ddply(viagraData, c("dose"), transform,
#                    mean = mean(libido))
# 
# CleanData <- CleanData[-1]
# CleanData <- CleanData[c(2,1,4)]
# names(CleanData)[2] <- "value"
# CleanData$Grand <- round(mean(CleanData$value), 2)
# grandMean <- mean(playANOVA$value)
# 

data.frame(group = rep(1:3,each = 5),
           value = c(3,2,1,1,4,
                     5,2,4,2,3,
                     7,4,5,3,6),
           unique = rep(1:5, 3)) %>% 
  group_by(group) %>% 
  mutate(avg = mean(value),
            sd = sd(value)) %>%
  ungroup() %>% 
  mutate(grand = mean(value)) %>% 
  group_by(group) %>% 
  mutate(y_hat = mean(value)) %>% 
  ungroup() %>%
  mutate(SR = round((value - y_hat)^2, 2),
         ST = round((value - grand)^2, 2),
         SM = round((grand - y_hat)^2, 2)) -> playANOVA
  
SSR <- round(sum(playANOVA$SR),1)
SSM <- round(sum(playANOVA$SM),1)
SST <- round(sum(playANOVA$ST),1)

playANOVA %>% 
  select(group, value, y_hat, grand, SR, ST, SM) %>% 
  knitr::kable(col.names = c("x", "y", "group-wise mean", "global mean", "y - group sq", "y - grand sq", "bar - group sq"),
               caption = "Summary of the ANOVA model.")


```

FakeDataTable  The values used in our example ANOVA calculations.

Figure \@ref(fig:simplePlusResiduals) shows the residuals of both the fitted model (left--hand--side) and the null model (middle). And similar to Figure XX in the previous chapter, it also displays the differences between the fitted and the null model (right). If you square each of the values and sum them up in each of the three cases, you get our dear friends $SSR$, $SST$ and $SSM$, respectively. In our example, these quantities are:

<div class="fooNarrow">

| Term  | Value   |
|-------|---------|
| $SSR$ | `r SSR` |
| $SST$ | `r SST` |
| $SSM$ | `r SSM` |

</div>

For reference, Table \@ref(tab:FakeDataTable) shows all the values we are working with. Note that $SSR$, $SST$ and $SSM$ are the sums of the values in columns 5--7. ^[If you know what to look for you can see that things can even be calculated by hand!!]

Have you noticed anything with the three sums of squares? Like before, we have the relationship

$$ SST = SSM + SSR, $$

so the total variation in the data, $SST$, can be split into the variation that can be explained by fitting our model, $SSM$, and the variation that we cannot do anything about, $SSR$. As before, we could calcuate $R^2 = SSM/SST$ as the proportion of the total variance that can be explained by our model. However, this number just gives us an impression whether or not or model fits the data significantly better than the null model. If we want a statistical test, we need to employ the $F$--test which compares $SSM$ with $SSR$. The corresponding output is shown in Table \@ref(fig:FakeAnova). This actually is the type of output you'll usually encounter when ANOVA is performed.

<!-- <<echo = FALSE, results = 'asis'>>= -->
<!-- xtable(LibidoDose.lm.anova, caption = "One-way ANOVA output for the toy data shown in Figure 15.3.", label = "FakeAnova") -->
<!-- @ -->

In the table we see our familiar values $SSM$ and $SSR$, with the associated degrees of freedom $df_1$ and $df_2$. Like in the previous chapter, these numbers can be used to calculate the mean sum of squared for the model and the total, i.e.

$$MSM = \frac{SSM}{df_1} = \frac{ `r SSM`}{2}  \quad \quad \text{and} \quad \quad MSR = \frac{SSR}{df_2} = \frac{`r SSR`}{12}. $$

The ratio of these two numbers then gives the $F$--statistic. If it is large, we know that the variance explained by the model is large relative to the variance that remains after fitting it. In other words the group mean values are far apart from each other while the scatter of the original data points around the group mean values is small. Statisticians say "the variance between the groups is significantly larger than the variance within the groups". ^[We talked about how to derive the correct values for $df_1$ and $df_2$ already a little on page _explainDf_. In one--way ANOVA with $k$ groups we have $df_1 = k-1$ and $df_2 = n-k$. In our example we have $n = 15$ observations in total and $k = 3$ groups.]



You know the rest of the story already. To obtain the $p$--value, the observed value of the $F$--statistic is compared to the corresponding null distribution which is the $F$--distribution with $df_1$ and $df_2$ degrees of freedom. If it is large, there is not much area under the curve to the right of the observed $F$--statistic and the $p$--value is small. If it's smaller than you chosen signficance level $\alpha$, you can reject the null model with ($1-\alpha$)% confidence.

Here the $p$--value is quite small, so it seems safe the reject the null model. But if linear regression and ANOVA are supposed to be the same, why do we only have one $p$--value on the ANOVA--type output while we had three in the linear regression output?! The answer is that linear regression looks at the data more closely and tests the diffrences between the individual groups (specifically, it tests the differences between group 1 and the other groups in turn) whereas the $F$--test takes a more unifying approach and examines the question whether there is at least one significant difference. So as a rule of thumb, if at least one $t$--test in the linear regression framework is significant (APART FROM THE ONE FOR THE INTERCEPT), then the $F$--test will also be ^[Note that the $p$--value for the $F$--test is sort of a compromise between the two interesting $p$--values for the $t$--tests in linear regression].

## ANOVA on our Martian dataset

Let's return to our burning question whether Martian age group has an influence on eye--sight. Here, we'll do it the other way round and first use the $F$--test to take a global look at whether there are _any_ differences in eye--sight between the age groups. If we do find any differences, then it is worth looking at the problem in a more detailed way to see where exactly the differences are.

```{r echo = FALSE}
martian$AgeIndex <- factor(martian$AgeIndex, levels=c("Youth","Adult","Senior"))

AgeEyeIndex.lm <- lm(Eye ~ AgeIndex, data=martian)
AgeEyeIndex.lm.summary <- summary(AgeEyeIndex.lm)
AgeEyeIndex.anova <- anova(AgeEyeIndex.lm)

broom::tidy(AgeEyeIndex.anova) %>% 
  mutate(term = c("Model", "Residuals")) %>% 
  knitr::kable(col.names = c("Term", "Degrees of freedom", "Sum of the Squared...", "Mean of the Squared...", "Test statistic", "p-value"),
               caption = "Summary of the F--test applied to Eye-sight described by Age.", 
               digits = 2)


```


Table \@ref(tab:ANOVAeyeage) shows the ANOVA output that contains the relevant information for the $F$--test. As expected from Figure \@ref(fig:tripleAgeEyePlot), the model that fits different mean values for each age group offers a big improvement over the null model that fits all data points using the overall mean. The $MSM$ is much bigger than the $MSR$ so that we know that the variation of the eye--sight scores _between_ the age groups is much bigger than the variation _within_ the age groups. The observed $F$--statistic therefore takes an epic value of over $34.94$ which obviously results in a very small $p$--value. So there basically is no chance that the null model is true and we observe group differences that are as big (or bigger) than the ones in our data.


The $F$--test confirms our visual suspicion that the differences between the groups are significant. Let's put our linear regression glasses on and see _where_ the differences are. Table \@ref(tab:linregeyeage) shows the corresponding output. It contains three very small $p$--values. The first $p$--value tells us the (not so interesting) fact that the underlying population mean of the eye--sight scores of the young Martians is significantly different from zero. The other two $p$--values are more interesting, since they tell us that both the adult and the senior Martians reveal eye--scores that are significantly different from the scores of the young Martians.


```{r echo = FALSE}
broom::tidy(AgeEyeIndex.lm) %>% 
  mutate(term = c("$\\hat\\beta_0$ ($\\bar{y}_{youth}$)", 
                  "$\\hat\\beta_1$ ($\\bar{y}_{adult}$ - $\\bar{y}_{youth}$)", 
                  "$\\hat\\beta_2$ ($\\bar{y}_{senior}$ - $\\bar{y}_{youth}$)")) %>% 
  knitr::kable(col.names = c("Term", "Estimate", "Standard Error", "Test statistic", "p-value"),
               caption = "Summary of the ANOVA model for Eye--sight as described by Age, categorical.", 
               digits = 2)

```

To find out if there are significant differences between the adult and the senior Martians you can carry out a separate $t$--test. Or you fit a different linear model where either the adult or the senior group is the reference group, which then corresponds to the intercept. A word of caution though: if you do many tests on the same dataset, you run into problems of multiple testing. That is beyond the scope of this course!




```{r echo = FALSE}

AgeEyeIndex.aov <- aov(Eye ~ AgeIndex, data=martian)
AgeEyeIndex.aov.summary <- summary(AgeEyeIndex.aov)
AgeEyeIndex.aov.Tukey <- TukeyHSD(AgeEyeIndex.aov)

AgeEyeIndex.aov.Tukey$AgeIndex %>% 
  as.data.frame() %>% 
  mutate(Comparison = row.names(.)) %>% 
  select(Comparison, diff, lwr, upr, `p adj`) %>% 
  # mutate(term = c("$\\beta_0$ ($\\\bar{y}_youth$)", "$\\beta_1$ ($\\\bar{y}_aduld$ - $\\\bar{y}_youth$)", "$\\beta_2$ ($\\\bar{y}_senior$ - $\\\bar{y}_youth$)")) %>% 
  knitr::kable(col.names = c("Comparison", "Difference in means", "Lower 95", "Upper 95", "Adjusted p-value"),
               caption = "All pair-wise group comparisons via Tukey's HSD (Honestly Significant Differences) as a post-hoc test.")

```



