---
title: "Hypothesis testing details and practice"
subtitle: "Notes and in-class exercises"
format: 
  html:
    embed-resources: true
    toc: true
eval: true       
---

# Notes

- **File organization:** Save this file in the "Activities" subfolder of your "STAT155" folder.

## Learning goals

By the end of this lesson, you should be able to:

- Apply the procedure for a formal hypothesis test

- Articulate how we can formalize a research question as a testable, statistical hypothesis

## Readings and videos

Please complete the following reading or videos **before** class:

- Reading: Section 7.3 (stop when you get to Section 7.3.4)  in the [STAT 155 Notes](https://mac-stat.github.io/Stat155Notes/)

- Video 1: [Introduction to Statistical Inference](https://voicethread.com/share/15687646/)
- Video 2: [Hypothesis Testing Framework](https://macalester.voicethread.com/myvoice/thread/15687700)
- Video 3: [Hypothesis Testing Procedure](https://voicethread.com/share/15713570/)


# Exercises

## Exercise 1

**Research Question**: Can we predict whether or not a mushroom is poisonous based on the shape of its cap?

For this exercise, we will look at data from various species of gilled mushrooms in the Agaricus and Lepiota Family. We have information on whether a mushroom is `poisonous` (TRUE if it is, FALSE if it's edible), the shape of its cap (`cap_shape`, a categorical variable with 6 categories), the texture of its cap surface (`cap_surface`, a categorical variable with 4 categories), and the size of its gills (`gill_size`, a categorical variable with two categories)

```{r warning=FALSE, message=FALSE}
# Load the data & packages
library(tidyverse)
library(readr)
library(broom)

mushrooms <- read_csv("https://Mac-STAT.github.io/data/mushrooms.csv")

mushrooms <- mushrooms %>%
  mutate(cap_shape = relevel(as.factor(cap_shape), ref="flat")) %>%
  dplyr::select(poisonous, cap_shape)

head(mushrooms)
```

### Part a

One of the most poisonous species of mushrooms is the *Amanita phalloides* or ["Death Cap" mushroom](http://www.bccdc.ca/health-info/prevention-public-health/death-cap-mushrooms), which typically has a flat cap shape when mature. Based on this anecdote, we hypothesize that species of mushrooms with flat caps in general may be *more* likely to be poisonous than edible.

First, let's translate this question to an appropriate null and alternative hypothesis that we can compare with a formal hypothesis test. Remember that `poisonous` is a binary outcome, so we need to frame our null and alternative hypotheses in terms of odds (i.e., Odds(poisionous | flat cap) = P(poisonous|flat cap)/P(edible | flat cap)).

> **Your answer**

### Part b

- Fit a logistic regression model to investigate whether `cap_shape` is associated with a mushroom being `poisonous`. (Note that in the setup code above, we have forced the reference category for the `cap_shape` predictor to be `flat`; without this, the reference category by default would be set as `bell`, which is the first category when sorted alphabetically).

```{r}
mushroom_mod1 <- ()

coef(mushroom_mod1)
```



### Part c

Provide an appropriate interpretation of the intercept coefficient **on the odds scale**. Based on this interpretation, do you believe mushrooms with flat caps are more likely to be poisonous, or more likely to be edible?

> Your answer here


### Part d

Let's look at the full model summary:

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

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

> **Your answer**

### Part e

- Report and interpret the p-value for the intercept term.
- Based on this p-value and a significance level of 0.05, do we have evidence that mushrooms with flat caps are more likely to be poisonous than edible?

> **Your answer**

### Part f

Now suppose we are interested in whether the odds of being poisonous are different for mushrooms with other cap shapes. 

By hand, calculate the odds of being poisonous for mushrooms with knobbed caps, conical caps, and sunken caps (remember that the non-exponentiated coefficients represent a difference in log-odds compared to the reference category):

> odds(poisonous | knobbed cap) =

> odds(poisonous | conical cap) =

> odds(poisonous | sunken cap) =

### Part g

Based on these odds, which of the 4 mushroom cap shapes we've investigated (flat, knobbed, conical, and sunken) do you believe is the best indicator that it's edible? Which cap shape do you expect is most likely to be poisonous?

> **Your answer**

### Part h

Let's get the full model summary again:

```{r}
tidy(mushroom_mod1) %>% 
    mutate(exp_estimate = exp(estimate)) %>% 
    select(term, estimate, exp_estimate, everything())
```

Now report and interpret the p-values for the coefficients corresponding to `cap_shapeknobbed`, `cap_shapeconical`, and `cap_shapesunken`:

> **Your answer**

### Part i

Based on the model summary output in **part h** above, if you were given a plate of mushrooms with different cap shapes and had to pick one to eat, which one would you choose? Which cap shape would you absolutely avoid at all costs? Are your decisions guided by the coefficient estimates, the p-values, or both?

> **Your answer**

### Part j

Let's look at the data a slightly different way, using a 6x2 table of counts:

```{r}
mushrooms %>% 
  mutate(cap_shape=as.factor(cap_shape),
         poisonous=as.factor(poisonous)) %>%
  dplyr::count(cap_shape, poisonous, .drop=FALSE) %>% 
  pivot_wider(names_from=poisonous, values_from=n, names_prefix="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 absolutely avoid, would you choose the same shapes? Why or why not?

> **Your answer**



## Exercise 2

For this exercise, let's return to the fish dataset from a previous activity.

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

head(fish)
```

**Research question:** We believe the length of a fish (measured in centimeters) is causally associated with its mercury concentration (measured in parts per million [ppm]). We suspect that the river a fish is sampled from may be a confounder, since 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. Assume we have specified a significance level of 0.05. 

> **Response**

### 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 = ___, colour = ___)) + 
  # [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`). Assume we've set a significance level of 0.05. 

> **Response**

### Part e

Interpret the coefficient estimate, test statistic, and p-value for the `RiverWacamaw` coefficient in this revised model (`mod_fish2`). (again, you can assume we've set a significance level of 0.05).

### 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, (and not Beta_1 ≠ 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`?

> **Response**

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

> **Response**

### Part g (CHALLENGE)

You point out to the other researcher that the intercept and `RiverWacamaw` coefficients are both negative, so whatever difference in mercury concentration between the two rivers your model predicts "at baseline" is not useful or meaningful--you cannot have a fish that is 0cm long, nor a mercury concentration <0ppm.

You propose that a more appropriate model should transform the `Length` variable in some way to make the intercept more interpretable. Create a new variable named `Length_adj` with this transformation and use it to re-fit the model:

```{r}
mod_fish3 <- lm(Concen ~ Length_adj*River, data=fish)
summary(mod_fish3)
```

Compare the output of this model to that of `mod_fish2`. What happened to the estimate, test statistic, and p-value for the `RiverWacamaw` coefficient? How does this affect your conclusion? How about the other researcher's conclusion? 

> **Response**


