---
title: "Hypothesis testing: additional considerations (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')

# Use a color blind friendly color palette throughout doc
library(tidyverse)
cb_palette <- c("#0072B2", "#D55E00", "black", "#E69F00", "#56B4E9", "#009E73", "#F0E442",  "#CC79A7")
scale_colour_discrete <- function(...) scale_colour_manual(values = cb_palette, ...)
scale_fill_discrete   <- function(...) scale_fill_manual(values = cb_palette, ...)
theme_set(theme_bw())
```


## 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 eval=TRUE}
# 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))
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: Hypothesis test review

Consider the following hypotheses:

$H_0$: $\beta_1 = 0$ (wage is *not* associated with marital status)

$H_a$: $\beta_1 \ne 0$ (wage *is* associated with marital status)

Evaluate these hypotheses using a 0.05 **significance level**.

Be sure to state a clear conclusion, and support that conclusion with appropriate evidence.



### Example 2: p-value

The p-value for the above test was 0.0003 (which is very small!).
How can we interpret this p-value?   

a. 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}$$

b. 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}$$

c. 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}$$



### Example 3: Just some considerations

Hypothesis tests are a powerful tool, but need to be considered carefully!

a. *Causal misinterpretations*. The above results do *not* establish a causal association between wage and marital status, i.e. that being single leads to lower wages. Why?

b. *Statistical vs practical significance*. We observed a *statistically significant* association between wages and marital status. BUT *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**?

```{r}
confint(wage_model_1)
```

c. *Type I & Type II errors*. Our test might have led us to an incorrect conclusion, not because we did anything wrong, but because we got unlucky data! Check out the errors below. What type of error might we have made in our conclusion about wages and marital status?


Conclusion            $H_0$ actually true            $H_0$ actually false
--------------------- ------------------------------ --------------------------------
fail to reject $H_0$  correct!                       Type II error (false negative)
reject $H_0$          Type I error (false positive)  correct!




### Example 4: Type I & Type II error rate

Type I and Type II error *rates* measure the *chance* that we end up making one of these errors. These can be expressed as conditional probabilities:

- Type I Error (false positive) rate    
    - P(reject $H_0$ | $H_0$ actually true)
    - P(conclude an effect *does* exist when it actually *doesn't*)
    
- Type II Error (false negative) rate    
    - P(fail to reject $H_0$ | $H_0$ actually false)
    - P(conclude an effect *doesn't* exist when it actually *does*)

Type II error rates depend upon various factors, that we only have some control over (e.g. actual effect size and sample size).
However, we *do* have control Type I error rate.
How?





## Exercises {-}

**Goal**

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



### Exercise 1: Conceptual understanding

a. Suppose that you and a friend are in two different sections (each with the same number of students) of the same class. On your respective midterm exams, you each obtained 85%, and the class average in both of your classes was 80%. Could one of you or your friend still be considered further above the class average than the other? Briefly explain.

b. Now suppose that your section's test scores were more tightly packed around 80%: maybe the standard deviation of your section's scores was 2.5, whereas the standard deviation of your friend's section's scores was 5. Which of you or your friend was further above the class average? Explain/justify your answer.

c. Broadly speaking, does a p-value measure the chance of a hypothesis being true, or, the chance of the data having occurred?

d. Why can't a p-value measure the other quantity that you didn't choose in Part c?

e. Explain in words why, in calculating a p-value, we need to assume that the null hypothesis is true.

f. Suppose that a hypothesis test yields a p-value of 1e-6 ($1\times 10^{-6}$). What can you tell about the magnitude of the effect or the uncertainty of the effect from this p-value? (i.e., What can you tell about the coefficient estimate or the standard error?)



### Exercise 2: Statistical vs. practical significance

Music researchers compiled information on 16,216 Spotify songs. They examined the relationship between a song's genre (latin vs. not latin) and song duration in seconds. Their modeling code and output is below:

```{r}
spotify <- 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)
```

a. Based on their below plot of `duration` by `latin_genre`, does song duration appear to be associated with genre?

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

b. Check out a model of this relationship and interpret the `latin_genreTRUE` coefficient. *In the context of song listening, is this a meaningful effect size?*

```{r}
spotify_model <- lm(duration ~ latin_genre, spotify)
coef(summary(spotify_model))
```   

c. Report and interpret the p-value for the `latin_genreTRUE` coefficient.

d. Use the p-value to make a yes/no decision about the evidence for a relationship between genre and song duration.

e. This exercise highlights the 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?)



### Exercise 3: Full model context

When interpreting hypothesis tests, we must be consider the full model context -- what else is / is not in our model? Consider our question about a wage gap and marital status. Compare and contrast our conclusions from the `maritalsingle` coefficient in the 2 models below, one of which controls for `age`, a potential *confounder*.

IMPORTANT: These results are *not* contradictory! Rather, together they provide a more complete picture about the relationship of interest.

```{r}
wage_model_1 <- lm(wage ~ marital, data = cps)
coef(summary(wage_model_1))

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







### Exercise 4: Power

Recall: Via our chosen significance level (typically 0.05), we control the Type I error (aka false positive) rate:

P(Type I error) = P(reject $H_0$ | $H_0$ is actually true) = 0.05

We have less control over the Type II error (aka false negative) rate:

P(Type II error) = P(fail to reject $H_0$ | $H_0$ actually false) = ???

An important related concept is that of statistical **power**, the *true positive rate*. Specifically, power is the probability of (correctly!) rejecting $H_0$ when $H_0$ is actually false:

power = P(reject $H_0$ | $H_0$ actually false) = 1 - P(Type II error)

Less formally, statistical power is the probability of detecting a relationship when there truly is a relationship. Thus power is a good thing!! We want it to be as high as possible!! Mainly, if an effect exists, we want there to be a high chance that we detect it.

Unfortunately, we can't control power. It's influenced by various factors. Navigate to [this page](https://rpsychologist.com/d3/nhst/) to check out an interactive visualization of the factors that influence statistical power.

Under "Settings", next to the "Solve for?" text, click "Power". You will vary the 3 different parameters (significance level, sample size, and effect size) one at a time to understand how these factors affect power.

Some context behind this interactive visualization:

- Visualization is based on a one sample Z-test:    
    This is a test for whether the true population mean equals a particular value. (e.g., true mean = 30)
- The effect size slider is measured with a metric called Cohen's d:
    - Cohen's d = magnitude of effect/standard deviation of response variable
    - Here: how far is the true mean from the null value in units of SD?
    - e.g., If the null value is 30, true mean is 40, and the true population SD of the quantity is 5, the Cohen’s d effect size is (40-30)/5 = 2.

a. What is your intuition about how changing the significance level will change power? Check your intuition with the visualization and explain why this happens.

b. Repeat Part a for the sample size.

c. Repeat Part a for the effect size.



### Exercise 5: Ethical considerations

a. Check out this ["Still not significant" article](https://mchankins.wordpress.com/2013/04/21/still-not-significant-2/). What ethical consideration is arising here? (Just for fun: a related [xkcd comic](https://xkcd.com/1478/))

b. Check out the xkcd comic [here](https://xkcd.com/882/). What ethical consideration is arising here?

c. 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


### Exercise 6: Multiple testing

In the previous exercise, parts b and c illustrate the dangers of **multiple testing**: the more things we test, the more likely something will appear to be "significant" just by chance (even if it's not). Thus fishing around for significance can lead to misleading conclusions.

a. Let's consider some math. Using a 0.05 significance level, if we perform 1 test the probability of getting a Type I error is 0.05. If we perform *2* tests (with $H_0$ being true in each case), the probability of getting *at least 1* Type I error is:   

- Pr(at least one Type I error) = 1 - P(no Type I errors) = 1 - (1 - 0.05)^2 = 0.0975
    
```{r}
1 - (1 - 0.05)^2
```

This is called the *overall or family-wise Type I error rate*. If we perform *3* tests (with $H_0$ being true in each case), the probability of getting *at least 1* Type I error is:   

- Pr(at least one Type I error) = 1 - P(no Type I errors) = 1 - (1 - 0.05)^3 = 0.1426
    
```{r}
1 - (1 - 0.05)^3
```

Calculate the probability of getting *at least 1* Type I error if we perform 20 tests (with $H_0$ being true in each case). This is how many tests were completed in the jelly bean example!


b. What's the takeaway from part a?

c. If you are a Bio major, read on: To counter the phenomenon above, one solution when doing multiple testing is to inflate the observed p-value from each test. For example, the *Bonferroni method* penalizes for testing 'too many' hypotheses by artificially inflating the p-value:

Bonferroni-adjusted p-value = observed p-value * total # of tests

Equivalently, we artificially deflate the significance level:

Bonferroni-adjusted significance level = 0.05 / total # of tests

For example, if we do 20 tests at the 0.05 level, we reject each $H_0$ if its p-value is less than 0.0025 (i.e. 0.05/20). This makes it tougher to reject $H_0$, thus to make a Type I error. BUT what are the trade-offs? NOTE: There are less conservative approaches out there! Bonferroni is just easy to use.




## 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 `inclass_activities` subfolder within your `stat155` folder and locate the HTML file -- you can open it again in your browser.
- Clean up your RStudio session: End the rendering process by clicking the 'Stop' button in the 'Background Jobs' pane.
- Check the solutions in the course website, at the bottom of the corresponding chapter.
- Work on homework and/or any extra practice exercises!


