hypothesis-tests.Rmd

---
title: "Hypothesis Testing"
author: "Math 271"
date: "Spring 2022"
output: 
  html_document:
    css: lab.css
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache=TRUE)
library(tidyverse)
library(magrittr)
library(moderndive)
set.seed(42)
```

## [Chapter 9: Hypothesis Testing](https://moderndive.com/9-hypothesis-testing.html){target="_blank"}

Confidence intervals : estimation of an unknown quantity

Hypothesis testing: One framework for making binary decisions based on data. (There are others: see Game Theory/Decision Theory for a broader view of the topic.)

full name: Null Hypothesis Significance Testing

Most appropriate for Man vs Nature type scenario.

Want to decide between two competing explanations for how some data was generated.

More precisely, usually we have one special explanation (called a __null hypothesis__)and want to see if the data "disproves" that explanation.

Decide between:

- Data is consistent with my hypothesis  (__fail to reject the null hypothesis__)
- Data is inconsistent with my hypothesis (__reject the null hypothesis as not reasonable__)

Using the null hypothesis we must be able to make some specific predictions about what I expect will happen when I observe a data set. (A probability model of some sort.)

Then we compare our actual observed data to the outcomes predicted by the hypothesis and decide if the theory is contradicted or not.



```{r}
promotions

promotions %>% count(gender,decision) %>% group_by(gender) %>% mutate(p=proportions(n))

promotions %>% group_by(gender) %>% summarize(promoted=mean(decision=="promoted"))

# This code computes the difference between the two groups (gender) based on their status (promoted/not)

promotions %>% 
  group_by(gender) %>% 
  summarize(promoted=mean(decision=="promoted")) %>% 
  summarize(p_diff=diff(promoted)) # does female - male


```

p_diff = -0.29 meaning that women were 30% less likely to be promoted

In a real life situation we establish innocence and seek to disprove that innocence (develop a hypothesis):

Hypothesize: There is no gender discrimination, males and females are equally likely to be hired/promoted. _the null hypothesis_ (aka: `gender` and the `promoted` variables have no relationship to each other (call this independence technically.))

Next step is: make some predictions about the behavior of `p_diff` if the null hypothesis is true.

`sample()` - column based version of `slice_sample()` 

The code chunk below takes in the decision column and shuffles (permutes) it which in turn breaks the relationship between gender and promotion (independence).
```{r}
promotions %>% mutate(decision=sample(decision))

get_p_diff <- . %>% group_by(gender) %>% 
  summarize(promoted=mean(decision=="promoted")) %>% 
  summarize(p_diff=diff(promoted))

# below we are running the shuffled data through the code that computes the difference between the two groups (gender) based on their status (promoted/not):

promotions %>% get_p_diff
promotions %>% mutate(decision=sample(decision)) %>% get_p_diff

# we now want to get an understanding of how 'diff_p' interacts with a large data set, we do this by rerunning it 2000 times:
permute_p_diffs <- rerun(2000, promotions %>% mutate(decision=sample(decision)) %>% get_p_diff) %>% bind_rows()

# developing a visualization of how p_diff responds to sampling when columns are independent 
# sampling distribution for diff_p _assuming the null is true_.
permute_p_diffs
ggplot(permute_p_diffs) + aes(x=p_diff) + geom_histogram(binwidth=1/12) +
  geom_vline(xintercept=-0.291666, col="red")
```

The next step for us to take is to bring in our actual data set and assess it's relationship to the resampled (permuted) data. Our results (red line found on histogram) indicate that our actual data isn't totally unreasonable. 

If we're trying to prove that discrimination against women exists within the workplace, we would be looking toward the negative values (left side) of the histogram. For males, we look toward the positive values (right side).

## Acceptance and Rejection Regions

- Which parts of the sampling distribution support the alternative hypothesis?
- Use `quantile` to identify a cutoff(s) to separate the histogram into two regions
    + __Rejection region__ parts of the histogram which are more compatible with the alternative hypothesis than the null hypothesis.
    + __Acceptance region__ everything else 

Alternative hypothesis: Women are discriminated against. -> Negative values of p_diff.

Left side of the histogram is evidence for the alternative hypothesis, everything else is the acceptance region; where you continue to accept (fail to reject?) the null hypothesis. 
Confidence level -> (1 - Significance level).
The usual value for confidence levels are 95%. The corresponding "thing" in a hypothesis test would be a 5% or 0.05 significance level. It is the reverse of the confidence level. 

```{r}
permute_p_diffs

permute_p_diffs %>% summarise(crit_value=quantile(p_diff, 0.05)) #Looking off the cutoff for the bottom 5%.

permute_p_diffs <- permute_p_diffs %>% mutate(reject = p_diff < -0.208) 


ggplot(permute_p_diffs) + aes(x=p_diff, fill=reject) + geom_histogram(binwidth=1/12) + geom_vline(xintercept=-0.291666, col="red")
```

The previous table, FALSE means you are on the right side of the significance level, consequently, TRUE means that you are on the left side. In the histogram, this is shown in Red as the acceptance region & Blue as the rejection region. 

If the null hypothesis is true we would consider this data set to be unusual, so we would reject the null hypothesis, and accept the alternative hypothesis. Based on this we can conclude that it is more likely that discrimination occurred in the creation of these data. 

## What if I wanted to identify any sort of discrimination (against men or women)

Very positive values would be discrimination against men while very negative values would be discrimination against women. So there should be rejection regions in both tails of the sampling distribution. We'll use the same approach as confidence intervals, to look up critical cutoff values on both sides of the distribution.

```{r}
ci_probs <- . %>% multiply_by(c(-1,1)) %>% add(1) %>% divide_by(2)


crit_vals <- permute_p_diffs %>% summarise(p=ci_probs(0.95), crit_value=quantile(p_diff, p)) %>% print


permute_p_diffs <- permute_p_diffs %>% 
  mutate(reject = ! between(p_diff, crit_vals$crit_value[1], crit_vals$crit_value[2]) ) 

ggplot(permute_p_diffs) + aes(x=p_diff, fill=reject) + geom_histogram(binwidth=1/12,center=0) + geom_vline(xintercept=-0.291666, col="red") + scale_x_continuous(breaks = seq(-0.55,0.55,by=0.1)) + geom_rug()
```



## P-Values


One way to decide if this is unusual or not is to compute the proportion of the histogram which is _"more unusual" than our actual observation_.

```{r}
permute_p_diffs %>% summarize(more_unusual=mean(p_diff <= -0.2916)) ## one-sided disc against women test

## to equivalent methods of getting a two-sided p-value
permute_p_diffs %>% summarize(p_value= mean(!between(p_diff, -0.2916, 0.2916) ))
permute_p_diffs %>% summarize(p_value= mean(abs(p_diff) > abs(-0.2916)))

```

The value 0.006 indicates that the value we saw in the actual dataset (-.29) is something that occurs less that 1% of the time when no relationship is assumed.

This positions us to make a judgment based on the findings. Do the above findings provide enough compelling evidence to suggest that gender discrimination does exist within the bank's workplace? That is our decision to make. Randomness is a potential explanation for the findings, but there is very small chance of that.

## Types of Error


- Type I Error: Rejection of a true null hypothesis. (False positive. Convict innocent defendant.)
- Type II Error: Acceptance of a false null hypothesis. (False negative. Acquit guilty defendant.)

The __significance level__ $\alpha$ for a hypothesis test represents an upper limit on the chance that Type I Error will occur. Note that this assumes the null hypothesis is true, and is based only on this assumption. $\alpha$ is the chance that the sampling process will produce an "unusual result". 

The chance of committing a Type II Error is usually denoted $\beta$. To compute this number one needs to make an assumption not only about the null hypothesis, but also a specific alternative hypothesis. It represents the chance that a sample drawn from the alternative hypothesis will produce a non-rejecting test statistic. The converse $1-\beta$ is the chance of correctly rejecting, and is called the __statistical power__ of the test to detect the alternative hypothesis.



## Testing a mean

This time we will compare IMDB ratings for movie genres, Action and Romance. The null hypothesis will be that there is no relationship between genre and the mean rating.

- Null hypothesis: mean romance rating = mean action rating $$\mu_r - \mu_a = 0$$
- Alternative: they are not the same $$\mu_r - \mu_a \neq 0$$

```{r}
movies_sample
ggplot(movies_sample) + aes(y=rating, x=genre) + geom_boxplot()

## observed statisic value
movies_sample %>% group_by(genre) %>% summarize(xbar=mean(rating)) %>% summarize(xbar_diff=diff(xbar))
```


```{r}
## statistic computing pipeline (input data set, output xbar_diff value)
xbar_diff_stat <- . %>% 
  group_by(genre) %>% 
  summarize(xbar=mean(rating)) %>%
  summarize(xbar_diff=diff(xbar))

observed_xbar_diff <- movies_sample %>% xbar_diff_stat

## generating a statistic under null hypothesis (no relationship between rating and genre)
movies_sample %>% mutate(rating=sample(rating)) %>% xbar_diff_stat

## Permuting removes any possible relationship between rating and genere.
## redo this many times to estimate the distribution under the null

movies_permuted <- rerun(2000, movies_sample %>% mutate(rating=sample(rating)) %>% xbar_diff_stat) %>% bind_rows

ggplot(movies_permuted) + aes(x=xbar_diff) + geom_histogram() + geom_vline(xintercept=1.047222, col="red")

## cutoff values for the acceptance/rejection regions of a two-sided test with alpha=0.05
movies_permuted %>% summarize(p=ci_probs(0.95), crit_vals=quantile(xbar_diff, p))


## two sided p-value is the proportion of permuted samples where xbar_diff is more extreme than the 
## original data
movies_permuted %>% summarize(p_value=mean(abs(xbar_diff) >= abs(1.047222)))
```

## Inference for regression problems

Consider the simple linear regression problem:

> Given a collection of $(x,y)$ values, find the coefficient pair $(a,b)$ that produces the "best fit" prediction equation $\hat y = a + bx$.

```{r}
evals_ch5 <- evals %>%
  select(ID, score, bty_avg, age)
glimpse(evals_ch5)

ggplot(evals_ch5) + aes(x=bty_avg, y=score) + geom_jitter(alpha=0.3) + geom_smooth(method="lm", fullrange=TRUE) + xlim(0,9)

lm(score~bty_avg, data=evals_ch5)
```

In this case, we are using the beauty score (bty_avg) to compute teaching score (score).

What is the idea of simple linear regression?
The basis of linear regression is to find a line of best fit for our data. From a mathematical point of view, we want to find the slope and intercept value that gives you the equation for the best fitting line. The line of best fit should minimize residuals.



Our overall plan will be to:

- Compute the coefficients for the data set we have
- Use bootstrapping/permuting to investigate the effect of sampling variability on the coefficients
- Visualize and interpret the sampling distributions to make inferences about the true values of the coefficients.

We first produce a pipeline that will compute the coefficients for us

```{r}
lm(score~bty_avg, data=evals_ch5)

evals_ch5 %>% lm(score~bty_avg, .) %>% coef() %>% t %>% as_tibble %>% rename(intercept="(Intercept)")

lm_coef_stat <- . %>% lm(score~bty_avg, .) %>% 
  coef() %>% t %>% as_tibble %>% 
  rename(intercept="(Intercept)")

evals_ch5 %>% lm_coef_stat()
```

First, we turn the linear model function into a pipeline which will compute the coefficients for us. You cannot directly pipe a formula into the `lm()` function. `lm()` doesn't work like the rest of the tidyverse functions, it likes the data to be the second argument in the function. By using the placeholder ".", we can tell the pipe where to place the incoming data. Then we can use `coef()` or `coefficients()` (depending on how much you like abbreviation) to pull out only the coefficient values from the linear model object generated by `lm()`. Finally, to turn our coefficients into an object which we can use with other tidyverse functions using either `enframe()` or a combination of `t()` (Which transposes data) and `as_tibble()`. The result is a two column data frame with the column name defining which coefficient is in the row.   

### Confidence Intervals

```{r}
evals_ch5 %>% slice_sample(prop=1, replace = TRUE) %>% lm_coef_stat()

evals_bootstrap <- rerun(2000, evals_ch5 %>% slice_sample(prop=1, replace = TRUE) %>% lm_coef_stat()) %>% bind_rows()

evals_bootstrap

ggplot(evals_bootstrap) + aes(x=intercept) + geom_histogram()
ggplot(evals_bootstrap) + aes(x=bty_avg) + geom_histogram()

```

We bootstrap the sample, which means we get a new set of coefficients for each bootstrapped sample. Then we rerun the bootstrapping a bunch of times to determine how the coefficients' variability change with the sampling process taken into account. Comparing the values we get from bootstrapping the coefficients to the plot we created in line 229, the confidence band surrounding our line of best fit matches up with the range of intercepts. While our best estimate is shown as the blue line, we can be 95% confident that the line of best fit falls within this confidence band.   

```{r}
#Calculating the critical values for our data
(evals_cis <- evals_bootstrap %>% summarize(p=ci_probs(0.95), 
                              crit_vals_icept=quantile(intercept, p), 
                              crit_vals_slope=quantile(bty_avg,p)))

lm(score~bty_avg, data=evals_ch5) %>% confint
```

Formally, we can interpret this result as:
"We can be 95% that the intercept of the real regression line is between 3.72 and 4.02"

The `confint` function will take an `lm` object and apply probability theory math to compute the same thing we computed when we bootstrapped our sample a bunch. In practice, most people will use the theory based method but being able to use two different methods and return the same result, you can be more confident in your answer.

```{r}
ggplot(evals_bootstrap) + aes(x=intercept, y=bty_avg) + geom_point() + 
  geom_vline(xintercept=evals_cis$crit_vals_icept, col="red") + 
  geom_hline(yintercept=evals_cis$crit_vals_slope, col="blue") + 
  stat_ellipse(col="darkorange1")
```

Together, how does our intercept and slope vary as computed by one data set.
The confidence interval, just by looking at the intercept, can be seen here as the red lines, and the CI, computed just by looking at the beta coefficient can be seen as the blue lines. By using both, you could compute a confidence oval, which combines the CI computed both variables to create a more complex confidence object.

```{r}
ggplot(evals_ch5) + aes(x=bty_avg, y=score) + geom_point() + 
  geom_abline(aes(intercept=intercept, slope=bty_avg), 
              data=evals_bootstrap %>% slice_sample(n=100), 
              alpha=0.1) + 
  geom_abline(col="red", intercept=3.729830, slope=0.03359138)
  
```

A large intercept tends to have a small slope value, and vice versa. Visually, if our intercept is higher on the y-axis it will tend to be flatter to go through all the points. 

```{r}
evals_ch5 %>% slice_sample(prop=1, replace = TRUE) %>% lm_coef_stat()

observed_model <- lm(score~bty_avg, data=evals_ch5)

## resample the residuals only (not the whole row)
evals_ch5 %>% mutate(yhat=predict(observed_model), 
                     resid=residuals(observed_model),
                     score=yhat + sample(resid)) %>% lm_coef_stat()

boot_residuals <- . %>% mutate(yhat=predict(observed_model), 
                     resid=residuals(observed_model),
                     score=yhat + sample(resid, replace = TRUE))

boot_resid_coefs <- rerun(2000, evals_ch5 %>% boot_residuals %>% lm_coef_stat) %>%
  bind_rows
ggplot(evals_ch5) + aes(x=bty_avg, y=score) + geom_point() + 
  geom_abline(aes(intercept=intercept, slope=bty_avg), 
              data=boot_resid_coefs %>% slice_sample(n=100), 
              alpha=0.1) 

boot_resid_coefs %>% summarize(p=ci_probs(0.95), 
                              crit_vals_icept=quantile(intercept, p), 
                              crit_vals_slope=quantile(bty_avg,p))
observed_model %>% confint()
```


### Hypothesis Tests

- Null: No relationship between variables: slope is zero 
- Alternative: there is a relationship: slope is not zero (provided the relationship is linear)

How can we simulate the null hypothesis. What can we do to the data to make that hypothesis true?

```{r}
evals_ch5
lm(score~bty_avg, evals_ch5)

evals_ch5 %>% mutate(score=sample(score))

ggplot(evals_ch5 %>% mutate(score=sample(score))) + aes(bty_avg, score) + geom_point() + geom_smooth(method="lm")

evals_ch5 %>% mutate(score=sample(score)) %>% lm_coef_stat

null_permute <- rerun(2000, evals_ch5 %>% mutate(score=sample(score)) %>% lm_coef_stat) %>%
  bind_rows()

#null_permute

ggplot(null_permute) + aes(bty_avg) + geom_histogram()

null_permute %>% summarize(p=ci_probs(0.95), crit_slope=quantile(bty_avg, p))

null_permute <- null_permute %>% mutate(reject= !between(bty_avg, 
                                         quantile(bty_avg,0.025), 
                                         quantile(bty_avg, 0.975)))

ggplot(null_permute) + aes(bty_avg, fill=reject) + geom_histogram() +
  geom_vline(xintercept=0.0666, color="red")

```


Our observed statistic is in the rejection region: conclude that the true slope value is non-zero.

To make a p-value: "What proportion of the sampling distribution is more "extreme" than the observed statistic?"
```{r}
null_permute %>% mutate(more_extreme= abs(bty_avg) > abs(0.06664)) %>%
  summarize(p_val=mean(more_extreme))

```
In our simulations this never occurred in 2000 attempts. It's not really zero, but I can say it's probably not bigger than say 1/2000 or so.

## Conditions for inference on regression

- __L__inear relationship: $\hat y = a + b\cdot x + c\cdot x^2$ or more generally $\hat y = \sum_{i=1}^k b_i\cdot x_i$.
- __I__independence of observations (technically, the residuals)
- __N__ormality of residuals
- __E__qual variance of all residuals

1. To check linearity, investigate the __predictions vs residuals__. For non-linear relationships there will often be a curved pattern.
2. Checking independence requires __thought about the data collection__ processes.
3. To check normality, make a __histogram of the residuals__.
4. To check equal variance (homoskedasticity if we're feeling fancy), look again at the __predictions vs resdiuals__. This time we're looking to see if there is a pattern in the variability of the dots.


```{r}
observed_model <- lm(score~bty_avg, evals_ch5)
fortified_evals <- evals_ch5 %>% 
  mutate(yhat=predict(observed_model), resid= resid(observed_model))
ggplot(fortified_evals) +  aes(x=yhat, y=resid) + geom_point() + geom_smooth()

ggplot(fortified_evals) + aes(x=resid) + geom_histogram() ## not really normal (left-skew)

```

## Testing equal means for multiple groups

In the _Testing a mean_ section above, we used the statistic $\bar x_1 - \bar x_2$ (difference of sample means) to explore whether two groups had the same average value for a response variable. This statistic seems to be the obvious choice when comparing two groups, but what if there are more than two groups? 

As in the two group case, the research question is often: __Does the grouping variable have any relationship to the response variable?__ With more than two groups, this becomes:

- null hypothesis: There is no relationship between the variables, the group means are __all equal to some common value__ ($\mu_1 = \mu_2 = \mu_3 =\ldots$)
- alternative hypothesis: There is some relationship between the two variables: __at least one pair of group means are different__.

The usual method used with this question is called ANOVA (ANalysis Of VAriance), the idea is based on a partitioning of quantities known as _Sums of Squares_.

### Sums of Squares, Variance, Standard deviation

Drawings for this section: https://geogebra.org/notes/m82na63q

The _Sum of Squares_ of a list of numbers is exactly what it sounds like. 

```{r}
iris %$% sum(Sepal.Length^2)
```

```{r}
squares_plot <- ggplot(iris %>% slice_sample(n=10)) + aes(x=Sepal.Length, y=Sepal.Length) + 
  geom_point()
squares_plot + 
  geom_rect(aes(xmin=0,ymin=0, xmax=Sepal.Length, ymax=Sepal.Length), alpha=0.1,fill=1:10) + 
  theme(aspect.ratio = 1) + 
  labs(title="Squares for Sepal.Length")
```


The _degrees of freedom_ of a list of numbers is, _informally_, tells you how many of the numbers in the list are free to change. For a list of independent data, the degrees of freedom is the same as the sample size.

```{r}
length(iris$Sepal.Length)
```

Conversely, suppose we run a very simple model, where we include only the intercept column in the model. The predictions are then all the same value, equal to the sample mean. This predictions vector has only one degree of freedom, matching the one column in the model matrix.

```{r}
iris %$% mean(Sepal.Length)
icept_only <- lm(Sepal.Length~1, iris)
icept_only %>% predict %>% unname
```

The residuals from the intercept-only model have $n-1$ degrees of freedom. A basic property of these residuals is that they must sum to zero, therefore if $n-1$ of the values are known, the last one is no longer free to change.
```{r}
icept_only %>% residuals
icept_only %>% residuals %>% sum ## not exactly zero due to small rounding errors
```

The core idea of ANOVA, is that the Total Sum of Squares for a response variable can be partitioned into a Sum of Squares from the predictions, plus the Sum of Squares from the residuals. 

```{r}
## The sum of squares partitioned into two parts: predictions + residuals
sum(predict(icept_only)^2) 
sum(residuals(icept_only)^2)
sum(predict(icept_only)^2) + sum(residuals(icept_only)^2)

iris %$% sum(Sepal.Length^2) ## Total Sum of Squares
```

```{r}
iris %>% summarize(mean(Sepal.Length))

squares_plot + geom_vline(xintercept=5.843) + geom_hline(yintercept=5.843) + 
  geom_rect(aes(xmin=5.843,ymin=5.843, 
                xmax=0, ymax=0), alpha=0.1,fill=1:10) + 
  theme(aspect.ratio = 1) +
  labs(title="Squares for predictions", subtitle="there is one square for each point") 


squares_plot + geom_vline(xintercept=5.843) + geom_hline(yintercept=5.843) +
  geom_rect(aes(xmin=5.843,ymin=5.843, 
                xmax=Sepal.Length, ymax=Sepal.Length), alpha=0.1,fill=1:10) + 
  theme(aspect.ratio = 1) + xlim(0,NA) + ylim(0,NA) +
  labs(title="Squares for residuals") 
```


Similarly, the total degrees of freedom in the response can be partitioned into predictions and residual freedom: $149 + 1 = 150$

#### Mean squares 

A _mean square_ is a sum of squares divided by its associated degree of freedom.

You may be familiar already with one particular mean square: the mean square of residual distances from the overall sample average. This quantity has a special name: The _variance_.

```{r}
## all these compute the same value
iris %$% sum((Sepal.Length - mean(Sepal.Length))^2) %>% divide_by(149)
sum(residuals(icept_only)^2) / (nrow(iris) - 1)
var(iris$Sepal.Length)
```

Note how the quantities we have just discussed all appear in the `anova()` table for the intercept-only linear model.

```{r}
anova(icept_only)
```

#### Standard Deviation

We won't use this, but note that the square root of the variance is called the _standard deviation_.

```{r}
sqrt(var(iris$Sepal.Length)) == sd(iris$Sepal.Length)
```

### The F-Statsitic

The basic statistic we wish to calculate is the _mean square for the predictions_ divided by the _mean square for the residuals_, which we call an _F-statistic_.

```{r}
msp <- sum(predict(icept_only)^2) / (1)
msr <- sum(residuals(icept_only)^2) / (nrow(iris) - 1)
msp/msr
```

If you are familiar with the $t$-test, this $F$-statistic is closely related to the $t$-value, in this case it is exactly the value $t^2$ when testing the hypothesis that the mean is zero.
```{r}
iris %$% t.test(Sepal.Length)

t.test(iris$Sepal.Length)$statistic^2

```

Note that I can get all of the numbers we've discussed by fitting a model with no predictors at all, and passing that and the intercept-only model to the `anova` function.

```{r}

iris %>% {anova(lm(Sepal.Length~0, .), 
                lm(Sepal.Length~1, .))}
```
- Model 1 has no predictors at all, and all of the SS and Df are in the residuals, showing the total sum of squares.
- Model 2 has only the intercept, this puts 1 Df into the prediction SS, leaving 149 for the residuals. We see the familiar prediction and residual sums of squres for this model, and R has also computed the F ratio.

#### Comparing group means

Now, we turn to our real problem: Comparing the means of 3 groups. We will add a third model to our `anova()` call, one that includes the `Species` categorical predictor (3 categories). 


```{r}

iris %>% lm(Sepal.Length ~ Species, .) %>% predict %>% unname
## this prediction vector has three degrees of freedom
iris %>% lm(Sepal.Length ~ Species, .) 

iris %>% lm(Sepal.Length ~ Species, .) %>% model.matrix

iris %>% lm(Sepal.Length ~ Species, .) %>% residuals %>% unname
## has 150 - 3 = 147 degrees of freedom in residuals
## 
```



```{r}
sum((iris %>% lm(Sepal.Length ~ 1 + Species, .) %>% predict)^2)
sum((iris %>% lm(Sepal.Length ~ 1 + Species, .) %>% residuals)^2)
iris %>% {anova(lm(Sepal.Length~0, .), 
                lm(Sepal.Length~1, .),
                lm(Sepal.Length~1+Species, .))} %>%
  print(digits=6)

(63.21 / 2) / (38.96 / 147)
```

Lines for models 1 and 2 have the same Df and SS values. Focus on line 3. We see that in our model with the three categories, the residual SS has decreased from 102.2 to 39.0: The difference in these Residual SS values, 63.2, is what appears in the SS column for model 3.

So in this ANOVA table, the `Sum of Squares` column is telling you: _This much of the total sum of squares is transferred from the residuals to the predictions when adding this term to the model._ Aka the amount of the prediction SS that can be attributed to that model term.

Of the total sum of squares of the response variable, `5223.9` we can partition this into 

- 5121.7 that comes from allowing an overall mean term in the model (one intercept column)
- 63.2 that comes from allowing the three groups to have different means (two species indicator columns)
- 39.0 remaining SS in the residuals (with 150 - 1 - 2 = 147 degrees of freedom remaining)

```{r}
species_model_predictions <- lm(Sepal.Length~1+Species, iris) %>% predict(newdata=squares_plot$data)

squares_plot + aes(color=Species) +
  geom_vline(xintercept=5.843) + 
  geom_hline(yintercept=5.843) +
  geom_rect(aes(xmin=5.843,ymin=5.843, 
                xmax=species_model_predictions,
                ymax=species_model_predictions), 
            alpha=0.1, fill=1:10) + 
  theme(aspect.ratio = 1) + xlim(0,NA) + ylim(0,NA) +
  labs(title="Squares for group means") 

squares_plot + aes(color=Species) + 
  geom_vline(xintercept=species_model_predictions) + 
  geom_hline(yintercept=species_model_predictions) +
  geom_rect(aes(xmin=Sepal.Length,ymin=Sepal.Length, 
                xmax=species_model_predictions,
                ymax=species_model_predictions), 
            alpha=0.1, fill=1:10) + 
  theme(aspect.ratio = 1) + xlim(0,NA) + ylim(0,NA) +
  labs(title="Squares for residuals", subtitle="from group means") 

```

The F values are now comparing the Mean Square for each term to the Mean Square for residuals in the bottom model.

```{r}
(5121.68/1) / (38.956 / 147) ## intercept row F value
(63.2121/2) / (38.956 / 147) ## Species row F value
```
If adding Species to the model improves the predictions greatly, then it will transfer a large amount of the total SS into the model and away from the residuals. So a large F ratio is evidence that the means are all different, while a small F ratio is evidence that the means are all the same.

Let's permute the species and come up with an idea of the sampling distribution of F under the null hypothesis.


```{r}
iris %>% lm(Sepal.Length~Species, .) %>% anova()

iris %>% lm(Sepal.Length~Species, .) %>% {anova(.)[1,4]} %>% tibble(F=.)


iris %>% lm(Sepal.Length~Species, .) %>% 
  {anova(.)["Species", "F value"]} %>%
  tibble(F=.)

species_F_stat <- . %>% lm(Sepal.Length~Species, .) %>% 
  {anova(.)["Species", "F value"]} %>% 
  tibble(F=.)

iris %>% species_F_stat() ## observed statistic

iris %>% mutate(Species=sample(Species)) %>% species_F_stat()

iris_species_permute <- rerun(2000, iris %>% mutate(Species=sample(Species)) %>% species_F_stat()) %>% bind_rows()

iris_species_permute
```

```{r}
ggplot(iris_species_permute) + aes(x=F) + 
  geom_histogram(aes(y=..density..), boundary=0) + 
  geom_vline(xintercept = quantile(iris_species_permute$F, 0.95), col="red") + 
  geom_function(fun=. %>% df(.,2,147)) + # the F-distribution with parameters 2 and 147
  geom_vline(xintercept=qf(0.95, 2, 147), col="blue", linetype=2) # We can get the critical value from probability theory 
  
```

The p-values in the anova table are computed by looking up a theory-based probability from the F-distribution with `Df` and `Res.Df` parameters.

### Regression with numerical + categorical predictors


```{r}
intercept_only <- lm(Sepal.Length ~ 1, iris)
no_species <- update(intercept_only, ~ . + Sepal.Width)
parallel_slopes <- update(no_species, ~ . + Species)
with_interactions <- update(parallel_slopes, ~ . + Species:Sepal.Width)

iris_plot <- ggplot(iris) + aes(x=Sepal.Width, y=Sepal.Length, color=Species) + geom_point() 
iris_plot + geom_line(aes(y=predict(intercept_only), color=NULL))
iris_plot + geom_line(aes(y=predict(no_species), color=NULL))
iris_plot + geom_line(aes(y=predict(parallel_slopes)))
iris_plot + geom_line(aes(y=predict(with_interactions)))

anova(with_interactions)
```


#### Permuting residuals

To test if a more complicated model is a significant improvement over a simpler model, our resampling scheme will be the following

- Add a reshuffled version of the residuals from the simpler model to the predictions from the simpler model to create a "new" data set that comes from the simpler model.
- Fit the more complicated model to the "new" permuted data set. Compute the F statistic for the improvement to model fit.
- Do the previous step many times to get a estimate of the null distribution. Compare your observed F stat for the real data set to the observed one.

```{r}
iris %>% mutate(Sepal.Length = predict(parallel_slopes) + residuals(parallel_slopes)) %>% lm(Sepal.Length~Species*Sepal.Width, .) %>% {anova(.)[3,"F value"]} %>% tibble(F=.)

iris %>% mutate(Sepal.Length = predict(parallel_slopes) + sample(residuals(parallel_slopes))) %>% lm(Sepal.Length~Species*Sepal.Width, .) %>% {anova(.)[3,"F value"]} %>% tibble(F=.)

interaction_f_null <- rerun(2000, iris %>% mutate(Sepal.Length = predict(parallel_slopes) + sample(residuals(parallel_slopes))) %>% lm(Sepal.Length~Species*Sepal.Width, .) %>% {anova(.)[3,"F value"]} %>% tibble(F=.)) %>% bind_rows()

ggplot(interaction_f_null) + aes(x=F) + geom_histogram(boundary=0) + geom_vline(xintercept=0.406, col="red")

interaction_f_null %>% summarize(p_value=mean(F>0.406))

```