---
title: "Hypothesis testing details and practice (Notes)"
subtitle: "STAT 155"
author: "Your Name"
format:
  html:
    toc: true
    toc-depth: 2
    embed-resources: true
---


```{r setup}
#| include: false
knitr::opts_chunk$set(
  collapse = TRUE, 
  warning = FALSE,
  message = FALSE,
  error = TRUE,
  fig.height = 2.75, 
  fig.width = 4.25,
  fig.env = 'figure',
  fig.pos = 'h',
  fig.align = 'center')
```



::: {.callout-note title = "Learning goals"}

- Apply the procedure for a formal hypothesis test.
- Translate research questions into statistical hypothesis, and vice versa.

:::

::: {.callout-note title = "Additional resources"}

**Required video**

- [Hypothesis testing concepts](https://www.youtube.com/watch?v=c5BUsblTdIE)

**Optional**

- Read: Section 7.3 (stop when you get to Section 7.3.4)  in the [STAT 155 Notes](https://mac-stat.github.io/Stat155Notes/)
- Watch:    
    - [Introduction to Statistical Inference](https://voicethread.com/share/15687646/)
    - [Hypothesis Testing Framework](https://macalester.voicethread.com/myvoice/thread/15687700)
    - [Hypothesis Testing Procedure](https://voicethread.com/share/15713570/)

:::



\
\
\
\

# References {-}

**Hypothesis testing framework**

1. State the hypotheses *using formal notation*.

2. Compare *our* sample results to what we would *expect* if the null hypothesis were
true, i.e. "under the null hypothesis".   
    - Calculate a test statistic: $$\frac{\hat{\beta}_1 - \text{null value}}{s.e.(\hat{\beta}_1)} = \text{ number of s.e. that $\hat{\beta}_1$ falls from the null value}$$
    - Calculate a p-value.
    - Determine an appropriate "significance level" *before* looking at our data.

3. Make a conclusion.
    - Compare the p-value to some significance threshold (or, equivalently, compare the test statistic to some critical value).

4. Consider limitations!




\
\
\
\



# Warm-up {-}


**GOAL**

Understand the *population* model of wages by marital status among 18-25 year-old workers:
    
$$
E(wage | marital) = \beta_0 + \beta_1 maritalsingle
$$



We don't have data on all workers, but we can make *inferences* about this relationship using data on a sample of 1254 workers:

```{r}
# Load packages & data
# Remove outliers that had really high wages
library(tidyverse)
cps <- read.csv("https://mac-stat.github.io/data/cps_2018.csv") %>% 
  filter(age <= 25, wage < 200000)

# Plot the relationship
cps %>% 
  ggplot(aes(x = marital, y = wage)) + 
  geom_boxplot()

# Model the relationship
wage_model_1 <- lm(wage ~ marital, data = cps)
coef(summary(wage_model_1))

# Get CIs
confint(wage_model_1)
```


The *sample model* is

$$
E(wage | marital) = \hat{\beta}_0 + \hat{\beta}_1 maritalsingle = 27399.47 - 6596.13 maritalsingle
$$


Thus:

- We **estimate** that the expected / average wage for single workers is $6596.13 below that of married workers.    

- We expect that this estimate might be off by / have an **error** of $1817.41.

- We're **95% confident** that, in the *broader labor force*, the average wage among single workers is somewhere between $3030.63 and $10161.64 lower than that of married workers. Thus we have *statistically significant* evidence that wage is associated with marital status (a difference of $0 isn't in this interval, hence isn't a plausible value).



\
\
\
\


**EXAMPLE 1: Starting a formal hypothesis test**

$H_0$: Wage is *not* associated with marital status.    
$H_a$: Wage *is* associated with marital status.

Translate this to notation:

$H_0$: 

$H_a$: 

NOTE:

We will evaluate these hypotheses using a 0.05 **significance level**.
It's important to choose our significance level *before* starting our test to avoid post-hoc "tweaks" that suit our narrative.






\
\



**EXAMPLE 2: What would we expect under $H_0$?**

Suppose $H_0$ were true, i.e. $\beta_1$ were 0.
By the Central Limit Theorem, we'd expect the sampling distribution of possible $\hat{\beta}_1$ sample estimates to be Normally distributed around 0 with the standard error obtained from our model summary table:

$$
\hat{\beta}_1 \sim N(\beta_1, s.e.(\hat{\beta}_1)^2)
\;\;\; \Rightarrow \;\;\;
\hat{\beta}_1 \sim N(0, 1817.41^2)
$$

Based on the corresponding plot below, is *our* sample estimate of $\hat{\beta}_1 = -6596.13$ consistent with $H_0$?

```{r fig.width = 8}
# IGNORE THE SYNTAX!!!
# Put OUR sample estimate here
est <- -6596.13

# Put the corresponding standard error here
se <- 1817.41

data.frame(x = 0 + c(-4:4)*se) %>% 
  mutate(y = dnorm(x, sd = se)) %>% 
  ggplot(aes(x = x)) +
  stat_function(fun = dnorm, args = list(mean = 0, sd = se)) +
  geom_segment(aes(x = x, xend = x, y = 0, yend = y), linetype = "dashed") + 
  scale_x_continuous(breaks = c(-4:4)*se) + 
  geom_vline(xintercept = est, color = "blue") +
  labs(y = "density", x = "possible estimates IF H0 were true") 
``` 






\
\






**EXAMPLE 3: Test statistic**

```{r}
coef(summary(wage_model_1))
```


a. Calculate the test statistic for our hypothesis test.



b. Interpret the test statistic: Our sample estimate of $\beta_1$ is ___ standard errors ___ the null value of ___. 


c. So is our sample data consistent with $H_0$?





\
\



**EXAMPLE 4: p-value**

a. Use the 68-95-99.7 Rule with the sketch from Example 2 to *approximate* the p-value:    
    - less than 0.003
    - between 0.003 and 0.05
    - between 0.05 and 0.32
    - bigger than 0.32


b. Report the *exact* p-value from the model `summary()`.    

```{r}
coef(summary(wage_model_1))
```



c. How can we interpret this p-value?   
    - Given our sample data, it's very unlikely that wages are associated with marital status (i.e. that $H_a$ is true). $$P(H_a | data) = \text{prob of $H_a$ being true given our data}$$
    - Given our sample data, it's very unlikely that wages and marital status are unrelated (i.e. that $H_0$ is true). $$P(H_0 | data) = \text{prob of $H_0$ being true given our data}$$
    - If in fact there were no relationship between wages and marital status in the broader labor force (i.e. $H_0$ were true), it's unlikely that we'd have observed such a big wage gap between single and married workers "by chance". $$P(data | H_0) = \text{prob of observing our data given $H_0$ is true}$$
    
d. So is our sample data consistent with $H_0$?




\
\
\
\




**EXAMPLE 5: Conclusion**    

Throughout the test, we've been comparing *our* data to the null hypothesis $H_0$.
Thus at the conclusion of a hypothesis test, either...

- Our results are "statistically significant".
    - We have enough evidence to *reject* $H_0$ in favor of $H_a$
    - NOTE: this does *not* mean that we "accept" $H_a$, we just "have evidence" for it.
- Our results are *not* "statistically significant".
    - We do *not* have enough evidence to reject $H_0$, i.e. we *fail to reject* $H_0$.
    - NOTE: This does *not* mean that we "accept" $H_0$, we just "don't have enough evidence" to to reject it.

a. So, at the 0.05 significance level, what's our conclusion? Is there a **statistically significant** association between wages and marital status? Does this conclusion agree with the one we made using the confidence interval for $\beta_1$, (-10161.64, -3030.627)?


b. *Statistical* significance merely indicates evidence that an association between wages and marital status *exists* (it's non-0). Consider an important follow-up: is the *magnitude* or *size* of this association (a $3,030.63--$10,161.64 wage gap) **practically significant** / **contextually meaningful**?







\
\
\
\

::: {.callout-note title = "CIs, test statistics, and p-values"}

We've observed 3 equivalent approaches to determining whether or not to reject the null hypothesis that $\beta$ equals some null value.
Assuming a 95% confidence level:

- If the **95%** CI for $\beta$ does not include the null value, reject $H_0$!
- If the p-value < **0.05**, reject $H_0$!
- If the |test statistic| > **2**, reject $H_0$! 


:::

\
\
\
\


**EXAMPLE 6: tests for the intercept**

$H_0: \beta_0 = 0$    
$H_a: \beta_0 \ne 0$    

```{r}
coef(summary(wage_model_1))
```

a. The p-value for this test is miniscule (1.85 * 10^(-52)). What do you conclude?    

    - We have statistically significant evidence of a relationship between wages and marriage status.
    - We have statistically significant evidence that there's no relationship between wages and marriage status.
    - We have statistically significant evidence that the average wage of single workers is non-0.
    - We have statistically significant evidence that the average wage of married workers is non-0.    

b. This test isn't very useful. How could we make it more useful? For example, what if we wanted to compare the average wage to $25,000, not $0?







\
\
\
\


::: {.callout-note title = "Hypothesis t-Tests for Model coefficients"}

Consider some population model

$$
E[Y | X_1, X_2, X_3] = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_3
$$

In our model `summary()` table, the reported test statistic (`t value`) and p-value (`Pr(>|t|)`) correspond to the following "**t-test**": 

$H_0$: $\beta = 0$    
$H_a$: $\beta \ne 0$

The *meaning* or *intepretation* of these hypothesis depends upon the meaning of $\beta$, which itself depends upon:

- whether the corresponding $X$ predictor is quantitative or categorical
- what other predictors are in the model

:::

::: {.callout-note title = "Everything is regression! Or, regression is everything!"}

In practice, people often *name* specific types of hypothesis tests.
Though not taught this way in AP Stat or more traditional intro courses, many of these just boil down to specific regression models.
Meaning, *you already know how to do these tests*!
For example:

- **two-sample t-test**   
    - A test for the difference between 2 population group means.
    - Example: Testing hypotheses about $\beta_1$ in the model of `wage` by `marital` status
    
- **one-sample t-test**   
    - A test for the comparison of 1 population mean to some null value.
    - Example: Testing hypotheses about $\beta_0$ in the model of `wage` by `marital` status


:::




\
\
\
\




# Exercises {-}


**GOAL**

Practice and consider some nuances to be aware of in hypothesis testing.


\
\
\
\


## Exercise 1: Controlling for age {-}

Since we haven't controlled for important *confounders*, we can't use the above results to argue that there's a *causal* association between wage and marital status.
For example, age is a potential confounder: it impacts both marital status and wage.
To this end, consider the following model:

$$
E[wage | marital, age] = \beta_0 + \beta_1 maritalsingle + \beta_2 age
$$
    
Some sample insights are below:

```{r}
wage_model_2 <- lm(wage ~ marital + age, data = cps)
coef(summary(wage_model_2))
confint(wage_model_2)
```


### Part a {-}

Consider the following hypotheses about $\beta_1$:    
    
$H_0: \beta_1 = 0$    

$H_a: \beta_1 \ne 0$    
    
How can we interpret these hypotheses?
*Adjust* the interpretations below (which are currently incorrect):
    
$H_0$: Wage is *not* associated with marital status

$H_a$: Wage *is* associated with marital status



### Part b {-}

Report the following (obtained via the above chunk):

- 95% CI for $\beta_1$ = (___, ___)

- test statistic = ___

- p-value = ___


### Part c {-}


At the 0.05 significance level, and using *any* of the 3 pieces of information in Part b, what conclusion do you make?

- We have statistically significant evidence of an association between wage and marital status.
- We do *not* have statistically significant evidence of an association between wage and marital status.
- When controlling for age, we have statistically significant evidence of an association between wage and marital status.
- When controlling for age, we do *not* have statistically significant evidence of an association between wage and marital status.


### Part d {-}

Bringing it all together, summarize what you conclude about the relationship between wages and marital status from the *combined* results from `wage_model_1` (from the warm-up) and `wage_model_2` (from this exercise).


### Part e (OUTSIDE OF CLASS) {-}

*After* class, come back and practice how to interpret the CI, test statistic, and p-value from Part b.




\
\
\
\



## Exercise 2: Interaction? {-}

In some cases, our *research questions* define what model we should construct.
For example, we built `wage_model_2` to explore the relationship between wages and marital status while controlling for age.
In other cases, *hypothesis tests* can aid the model building process.
For example, should we have included an *interaction term* in `wage_model_2`?
That is, does the relationship of wages with marital status significantly differ by age?

### Part a {-}

Address this question using the plot below:

```{r}
cps %>% 
  ggplot(aes(y = wage, x = age, color = marital)) + 
  geom_smooth(method = "lm")
```

### Part b {-}

Now, address this question by building a *new* model.
Support your conclusion with a p-value.

```{r}

```


\
\
\
\


## Exercise 3: Spotify Part 1 {-}    

Let's consider a new context.
Two researchers make a claim: the average duration (in seconds) of a song is different in the `latin` genre than in non-latin genres (pop or edm specifically).
Thus they're interested in the following population model:

E(`duration` | `latin_genre`) = $\beta_0$ + $\beta_1$ `latin_genreTRUE`

And the following hypothesis:

$H_0$: $\beta_1 = 0$

$H_a$: $\beta_1 \ne 0$    


### Part a {-}

To test these hypotheses, Researcher 1 collects data on 20 songs and uses this to model song `duration` (in seconds) by `latin_genre`.
Based on their plot below, what do you anticipate will be Researcher 1's conclusion about the relationship?


```{r}
# Load the data
spotify_small <- read.csv("https://ajohns24.github.io/data/spotify_example_small.csv") %>% 
  select(track_artist, track_name, duration_ms, latin_genre) %>% 
  mutate(duration = duration_ms / 1000)
nrow(spotify_small)

spotify_small %>% 
  ggplot(aes(y = duration, x = latin_genre)) + 
  geom_boxplot()
```


### Part b {-}

Researcher 1's sample model is below.
Report and interpret their sample estimate $\hat{\beta}_1$.
*In the context of song listening, is this a meaningful effect size?*

```{r}
spotify_model_1 <- lm(duration ~ latin_genre, spotify_small)
coef(summary(spotify_model_1))
```   



### Part c {-}

At the 0.05 significance level, does Researcher 1 have statistically significant evidence that the average song duration is different (longer) in the latin genre than in the pop / edm genres?


\
\
\
\



## Exercise 4: Spotify Part 2 {-}    

Researcher 2 has the same research question, but gathers data on *16,216* songs:    
    
```{r}
spotify_big <- read.csv("https://ajohns24.github.io/data/spotify_example_big.csv") %>% 
  select(track_artist, track_name, duration_ms, latin_genre) %>% 
  mutate(duration = duration_ms / 1000)
nrow(spotify_big)
```


### Part a {-}

Based on their below plot of `duration` by `latin_genre`, what do you anticipate will be Researcher 2's conclusion about the relationship?

```{r}
spotify_big %>% 
  ggplot(aes(y = duration, x = latin_genre)) + 
  geom_boxplot()
```



### Part b {-}

Researcher 2's sample model is below.
Report and interpret their sample estimate $\hat{\beta}_1$.
*In the context of song listening, is this a meaningful effect size?*

```{r}
spotify_model_2 <- lm(duration ~ latin_genre, spotify_big)
coef(summary(spotify_model_2))
```   


### Part c {-}

At the 0.05 significance level, does Researcher 2 have statistically significant evidence that the average song duration is different (longer) in the latin genre than in the pop / edm genres?



\
\
\
\



## Exercise 5: Caution!! {-}    

Hypothesis testing is a super, super, super powerful tool!
But only if we use it responsibly.
Consider some cautionary tales.

### Part a {-}

HUH?!?
In the Spotify exercises, you just witnessed the stark difference between **statistical significance** and **practical significance**.
Explain *why* this happened.
That is, when might we observe statistically significant results that aren't practically significant?
(Hence when should we be especially cautious about interpreting the results of a hypothesis test?)


### Part b {-}

Check out this "Still not significant" article: [https://mchankins.wordpress.com/2013/04/21/still-not-significant-2/](https://mchankins.wordpress.com/2013/04/21/still-not-significant-2/)
What cautionary tale is it trying to tell about hypothesis testing?


### Part c {-}

Check out this XKCD comic: [https://imgs.xkcd.com/comics/significant.png](https://imgs.xkcd.com/comics/significant.png)
What cautionary tale is it trying to tell about hypothesis testing?


### Part d {-}

Consider an example that's closer to life.
The article
[https://fivethirtyeight.com/features/science-isnt-broken/](https://fivethirtyeight.com/features/science-isnt-broken/) 
highlights the dangers of fishing for significance and the importance of always questioning: 

- How many things did people test?
- Could I replicate these results if I got a different sample than them?
- What was their motivation?
- How did they measure their variables? Might their results have changed if they used different measurements?

Check out an interactive tool that helps you explore these ideas: [https://stats.andrewheiss.com/hack-your-way/](https://stats.andrewheiss.com/hack-your-way/)
Use this tool to "prove" 2 different claims.
For both, be sure to indicate what response variable (eg: Unemployment rate) and data (eg: senators) you used:

- the economy is significantly better when Republicans are in power
- the economy is significantly better when Democrats are in power


\
\
\
\





# Extra practice {-}

You might not get to these exercises during class.
Rather, they provide IMPORTANT extra practice *outside* of class with both *hypothesis testing* and *logistic regression*.


## Exercise 6 {-}

Can we predict whether or not a mushroom is poisonous based on the shape of its cap?
To answer this and related questions, we'll use data from various species of *gilled* mushrooms in the Agaricus and Lepiota families.
We have information on whether a mushroom is `poisonous` (TRUE if it is, FALSE if it's edible) and its `cap_shape` (a categorical variable with 6 categories):

```{r}
mushrooms <- read_csv("https://Mac-STAT.github.io/data/mushrooms.csv") %>% 
  mutate(cap_shape = relevel(as.factor(cap_shape), ref = "flat")) %>%
  select(poisonous, cap_shape)
```

Note that in the code above, we have forced the reference category for the `cap_shape` predictor to be `flat`, though its not first alphabetically:

```{r}
head(mushrooms)
mushrooms %>% 
  count(cap_shape)
```


### Part a {-}

Fit a logistic regression model to investigate whether a mushroom being `poisonous` can be predicted by its `cap_shape`: 

```{r}
mushroom_mod <- ___(poisonous ~ cap_shape, data = mushrooms, ___)
coef(summary(mushroom_mod))
```

### Part b {-}

Let's start by testing hypotheses about the intercept coefficient.
On the log(odds) scale:

$H_0$: $\beta_0 = 0$

$H_a$: $\beta_0 \ne 0$

Equivalently, on the odds scale:

$H_0$: $e^{\beta_0} = 1$

$H_a$: $e^{\beta_0} \ne 1$

Interpret these hypotheses *in context*.
HINT: It might help to do Part c, then come back!

### Part c {-}

Interpret our sample estimate of the intercept **on the odds scale**.


### Part d {-}

Report and interpret the test statistic for the intercept term (our coefficient of interest).


### Part e {-}

Report and interpret the p-value for the intercept term.


### Part f {-}

At a significance level of 0.05, do we have evidence that mushrooms with flat caps are more likely to be poisonous than edible?



\
\
\
\


## Exercise 7 {-}

Let's continue to explore different elements of our mushroom model:

```{r}
coef(summary(mushroom_mod1))
```

### Part a {-}

Interpret the *exponentiated* estimates of the `cap_shapeconical` and `cap_shapeknobbed` coefficients.
Remember: these are *odds ratios*, comparing to `flat` cap mushrooms.

```{r}

```



### Part b {-}

Based on these odds ratios, which of the 3 mushroom cap shapes would you be the *most* worried about eating: `flat`, `conical`, or `knobbed`?
The *least* worried about eating?


### Part c {-}

Report the p-values for the `cap_shapeconical` and `cap_shapeknobbed` coefficients.
What do you conclude from these?
Specifically, which of the 3 mushroom cap shapes would you be the *most* worried about eating: `flat`, `conical`, or `knobbed`?
The *least* worried about eating?


### Part d {-}

Putting this all together, if you were given a plate of mushrooms with different cap shapes and had to pick one to eat, which one would you choose: `flat`, `conical`, or `knobbed`?
Which would you avoid?
Are your decisions guided by the coefficient estimates, the p-values, or both?

### Part e {-}

Check out the table of counts below:

```{r}
mushrooms %>% 
  count(cap_shape, poisonous)
```

Now, if you were given a plate of mushrooms with different cap shapes and had to pick one shape to eat and one to avoid, would you choose the same shapes? Why or why not?



\
\
\
\



## Exercise 8 {-}

Let's return to the fish dataset:

```{r}
fish <- read_csv("https://Mac-STAT.github.io/data/Mercury.csv")
head(fish)
```

**Research question:** We believe the length of a fish (cm) is *causally* associated with its mercury concentration (ppm).
We suspect that the river a fish is sampled from may be a *confounder*: differences in the river environment may causally influence both the average length of fish (e.g. due to differences in water temperature or food availability) as well as mercury concentration (e.g. due to differences between the two rivers in mercury pollution levels). 

### Part a {-}

Fit a linear regression model that can be used to answer our research question.

```{r}
mod_fish1 <- ___
summary(mod_fish1)
```

### Part b {-}

Interpret the coefficient estimate, test statistic, and p-value for the `RiverWacamaw` coefficient.
Use a significance level of 0.05. 


### Part c {-}

Suppose we now want to determine if the causal effect of fish length on mercury concentration *differs* according to the river a fish was sampled from. 
First, modify the code chunk below to visualize the 3-way relationship between the `Concen`, `Length`, and `River` variables:

```{r}
fish %>% 
  ggplot(aes(x = ___, y = ___, color = ___)) + 
  # [ADDITIONAL GGPLOT LAYER(S)]
```

Next, fit an appropriate linear regression model with an interaction term to investigate this question.

```{r}
mod_fish2 <- ___
summary(mod_fish2)
```

### Part d {-}

Interpret the coefficient estimate, test statistic, and p-value for the `RiverWacamaw:Length` interaction term in this revised model (`mod_fish2`).
Use a significance level of 0.05. 


### Part e {-}

Interpret the coefficient estimate, test statistic, and p-value for the `RiverWacamaw` coefficient in this revised model (`mod_fish2`).


### Part f (CHALLENGE) {-}

Suppose another researcher runs the same model we fit in part c above (`mod_fish2`), but they claim that a more appropriate alternative hypothesis should be $\beta_1 < 0$, (not $\beta_1 \ne 0$ as is assumed by default when running a regression model).
Because of this, they reported a smaller p-value for the coefficient, and claim that the Wacamaw River has a lower baseline mercury concentration (i.e., when `Length = 0cm`).

- What is the p-value they would have reported for the `RiverWacamaw` coefficient in `mod_fish2`?

- What is a potential ethical problem with the other researcher's claim that the alternative hypothesis should be $\beta_1 < 0$? 




## Done!

- Finalize your notes: (1) Render your notes to an HTML file; (2) Inspect this HTML in your Viewer -- check that your work translated correctly; and (3) Outside RStudio, navigate to your 'Activities' subfolder within your 'STAT155' folder and locate the HTML file. Open and inspect this in a web browser.
- Clean up your RStudio session: End the rendering process by clicking the 'Stop' button in the 'Background Jobs' tab in the lower left pane.
- Check the wrap-up announcements and solutions in the online course manual, at the bottom of the corresponding chapter.
- Work on homework!


