---
title: "Lab 8: Confidence Intervals and Hypothesis Testing with Infer"
author: "Insert Name"
output:
  pdf_document
 
---

# Instructions


# Due: 12am night before lab

# Instructions
Work through the problems below (it may be helpful to have both the .rmd and .pdf files open side-by-side), and submit this document as a knitted .pdf file with answers to the Lab 8 assignment on Gradescope.




```{r, include=FALSE}
# Do not edit this code block/chunk
knitr::opts_chunk$set(echo = TRUE, message=FALSE, warning = FALSE, fig.width = 5, fig.height = 3)
```


# Confidence Intervals

## Background

If you have access to data on an entire population, say the opinion of every adult in the United States on whether or not they think climate change is affecting their local community, it's straightforward to answer questions like, "What percent of US adults think climate change is affecting their local community?". 

Similarly, if you had demographic information on the population you could examine how, if at all, this opinion varies among young and old adults and adults with different leanings.

But if you only have access to a sample of the population, as is often the case, the task becomes more complicated. What is your best guess for this proportion if you only have data from a small sample of adults? This type of situation requires that you use your sample to make inference on what your population looks like.

## An Important Note on Reproducibility

You will take random samples and build sampling distributions in this lab, which means you should set a seed using `set.seed()` in each code chunk where you generate random values to ensure that each time you knit your document, you obtain an identical random sample.


### Load packages

In this lab, we will explore and visualize the data using the **tidyverse** suite of packages, and perform statistical inference using **infer**.

Let's load the packages.

```{r load-packages, message=FALSE}
library(tidyverse)
library(openintro)
library(moderndive)
library(infer)
``` 
 
### The data


To keep our computation simple, we will assume a total population size of 100,000 (even though that's smaller than the population size of all US adults).

A 2019 Pew Research report states the following:

> Roughly six-in-ten U.S. adults (62%) say climate change is currently affecting their local community either a great deal or some, according to a new Pew Research Center survey.
>
>**Source:** [Most Americans say climate change impacts their community, but effects vary by region](https://www.pewresearch.org/fact-tank/2019/12/02/most-americans-say-climate-change-impacts-their-community-but-effects-vary-by-region/)
In this lab, you will assume this 62% is a true population proportion and learn about how sample proportions can vary from sample to sample by taking smaller samples from the population. 
We will first create our population assuming a population size of 100,000. 
This means 62,000 (62%) of the adult population think climate change impacts their community, and the remaining 38,000 does not think so.

```{r}
us_adults <- tibble(
  climate_change_affects = c(rep("Yes", 62000), rep("No", 38000))
)
```

The name of the data frame is `us_adults` and the name of the variable that contains responses to the question *"Do you think climate change is affecting your local community?"* is `climate_change_affects`.

We can quickly visualize the distribution of these responses using a bar plot.

```{r bar-plot-pop, fig.height=2.5, fig.width=10}
ggplot(us_adults, aes(x = climate_change_affects)) +
  geom_bar() +
  labs(
    x = "", y = "",
    title = "Do you think climate change is affecting your local community?"
  ) +
  coord_flip() 
```

We can also obtain summary statistics to confirm we constructed the data frame correctly.

```{r}
us_adults %>%
  count(climate_change_affects) %>%
  mutate(p = n /sum(n))
```

In this lab, you'll start with a simple random sample of size 50 from the population.  Make sure you change the seed before running the code.

```{r sample}
set.seed(1)

samp <- us_adults %>%
  rep_sample_n(size = 50, reps = 1)
```


### Exercise 1
What percent of the adults in your sample think climate change affects their local community? **Hint:** Just like we did with the population, we can calculate the proportion of those **in this sample** who think climate change affects their local community. Be sure to change the value of the seed in the previous chunk to a number of your choice.

****************************************************

****************************************************

### Exercise 2
Would you expect another student's sample proportion to be identical to yours? Would you expect it to be similar? Why or why not? Check with another student in the class to verify your hypotheses.

****************************************************

****************************************************

## Confidence intervals

Return for a moment to the question that first motivated this lab: based on this sample, what can you infer about the population? 

With just one sample, the best estimate of the proportion of US adults who think climate change affects their local community would be the sample proportion, usually denoted as $\hat{p}$ (here we are calling it `p_hat`). 

That serves as a good **point estimate**, but it would be useful to also communicate how uncertain you are of that estimate. This uncertainty can be quantified using a **confidence interval**.

One way of calculating a 95% confidence interval for a population proportion is to use the rule of thumb that in a bell-shaped distribution, 95% of all values are within 2 standard error units of the mean. For the sampling distribution of the sample proportion, this produces the following confidence interval:
$$\hat{p} \pm 2 \sqrt{ \frac{\hat{p} (1-\hat{p})}{n} }$$

Another way is using simulation, or to be more specific, using **bootstrapping**. The general bootstrapping procedure is the following:

- **Step 1.** Take a bootstrap sample - a random sample taken **with replacement** from the original sample, of the same size as the original sample.
- **Step 2.** Calculate the bootstrap statistic - a statistic such as mean, median, proportion, slope, etc. computed on the bootstrap samples.
- **Step 3.** Repeat steps (1) and (2) many times to create a bootstrap distribution - a distribution of bootstrap statistics.
- **Step 4.** Calculate the bounds of the C% confidence interval as the middle C% of the bootstrap distribution.

Instead of coding up each of these steps, we will construct confidence intervals using the **infer** package.

Below is an overview of the functions we will use to construct this confidence interval:

Function     | Purpose
------------:|:------------------------------------------
`specify`    | Identify your variable of interest
`generate`   | The number of samples you want to generate
`calculate`  | The sample statistic you want to do inference with, or you can also think of this as the population parameter you want to do inference for
`get_ci`     | Find the confidence interval
`visualize`  | Create a histogram of the bootstrap statistics

This code will create a bootstrap distribution based on your sample. 

```{r}
set.seed(1)

bootstrap_dist<- samp %>%
  specify(response = climate_change_affects, success = "Yes") %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "prop")  
```

The code in the next chunk will compute the 95% confidence interval using the percentile method for the proportion of US adults who think climate change affects their local community.

```{r}
percentile_ci<-bootstrap_dist %>% 
  get_confidence_interval(level = 0.95, type = "percentile")
percentile_ci
```


We can also compute a confidence interval using the standard error method, **provided the bootstrap distribution appears symmetric, unimodal and bell-shaped.** To do so, we first need to store the sample statistic as a variable, as shown in the following code chunk. 

```{r}
sample_statistic<- samp %>%
  specify(response = climate_change_affects, success = "Yes") %>%
  calculate(stat = "prop")  
sample_statistic
```
Note that we use the same code as we did for the bootstrap distribution, except we **removed** the line with the `generate` verb.

The next code chunk will compute the 95% confidence interval using the standard error method, which follows a similar pattern to the 95% interval using the percentile method, except we need to additionally specify the statistic we are using as our `point estimate`.

```{r}
se_ci<-bootstrap_dist %>% 
  get_confidence_interval(level = 0.95, type = "se", point_estimate = sample_statistic)
se_ci
```


The code in the following chunk will create a histogram of the bootstrap statistics.

```{r}
bootstrap_dist %>% 
  visualise( )
```

We can add an additional layer to `visualize` in order to shade the region corresponding to the confidence interval, as demonstrated in the following code:

```{r}
bootstrap_dist %>% 
  visualise( )+ 
  shade_confidence_interval(endpoints = percentile_ci)
```



To summarize:
- In `specify` we specify the `response` variable and the level of that variable we are calling a `success`.
- In `generate` we provide the number of resamples we want from the population in the `reps` argument (this should be a reasonably large number) as well as the type of resampling we want to do, which is `"bootstrap"` in the case of constructing a confidence interval.
- Then, we `calculate` the sample statistic of interest for each of these resamples (`prop`, which is short for proportion)
- Finally, we may either `visualize` the bootstrap distribution by creating a histogram of the statistic's values, or `get_confidence_interval` and specify the type (either `"percentile"` or `"se"`) and confidence `level`.




## Confidence levels


### Exercise 3a

Follow the procedures outlined above to create a new sample called `samp60` of size 60. Create a bootstrap distribution called `bootstrap_dist60`, plot it, and calculate the 95% confidence interval using the percentile method.

**************************************************

**************************************************


### Exercise 3b

Use the sample you created in exercise 3a (`samp60`) to find a 95% confidence interval using the **standard error method** for the proportion of US Adults who think climate change is affecting their local community. Compare the intervals produced by the methods in parts 3a and 3b.


**************************************************

**************************************************

### Exercise 4a

In this case, we have the rare luxury of knowing the true population proportion (62%) since we have data on the entire population. Does your confidence interval capture the true population proportion of US adults who think climate change affects their local community? If you are working on this lab in a classroom, does your neighbor's interval capture this value? 

**************************************************

**************************************************

### Exercise 4b

Each student should have gotten a slightly different confidence interval. What proportion of those intervals would you expect to capture the true population proportion? Why?

**************************************************

**************************************************


### Exercise 5

Since we know the population parameter, we can compute the sampling distribution. Plot a sampling distribution of 1000 samples of size 60 taken from the population, using `rep_sample_n`. Compare it to the plot of the bootstrap distribution you created in Exercise 3a. Do the distributions have the same center? Spread? Shape? 

**************************************************

**************************************************

In the next part of the lab, you will collect many samples to learn more about how sample proportions and confidence intervals constructed based on those samples vary from one sample to another.

Find the .Rmd file titled "Confidence_Intervals_App.Rmd" on the [Labs page](https://reed-statistics.github.io/math141s22-wells-website/lab.html) of the course website, upload this file to your RStudio server, open the .Rmd file, and then click the **Run Document** button at the top (where the **Knit Document** button is usually located). Use the interactive app that opens to answer some of the following question.


### Exercise 6

Given a sample size of 60, as well as 1000 bootstrap samples for each interval, and 50 confidence intervals constructed (the default values for the above app), what proportion of your confidence intervals include the true population proportion? Is this proportion exactly equal to the confidence level? If not, explain why. 

**************************************************

**************************************************




### Exercise 7a

Choose a different confidence level than 95%. Would you expect a confidence interval at this level to be wider or narrower than the confidence interval you calculated at the 95% confidence level? Explain your reasoning.

**************************************************

**************************************************



### Exercise 7b

Using the app, experiment with different sample sizes and comment on how the widths of intervals change as sample size changes (increases and decreases).

**************************************************

**************************************************




### Exercise 8
    
Finally, given a fixed sample size (say, 60) and fixed confidence level (say, 95%), how does the width of the interval change as you increase the number of bootstrap samples? Support your answer with graphics and/or calculations. **Hint:** Does changing the number of bootstrap samples affect the standard error?

**************************************************

**************************************************

# Hypothesis Testing
 
## Background
  
As the novel coronavirus spread throughout 2020, the benefit of facemasks and other protective gear for preventing the transmission of airborne viruses received renewed attention. Historically, two types of disposable masks see frequent use: cloth/cotton surgical masks and heavier N95 respirators. In this lab, we will use hypothesis testing to answer the following question: is there a difference in effectiveness between the two types of masks?

## Data

C. MacIntyre et al compare the efficacy of cloth masks to medical masks in hospital healthcare workers in <https://bmjopen.bmj.com/content/5/4/e006577.long>

A total of 1607 healthcare workers were recruited in the study from 14 hospitals in Hanoi, Vietnam. Wards in participating hospitals were randomly sorted into 3 groups: a group instructed to wear `medical` masks at all times on shift, a group instructed to wear `cloth` masks at all times on shift, and a `control` group instructed to follow their standard routine for mask usage. Participants were followed for 4 weeks of facemask use and for one week after to determine the whether each exhibited symptoms of a clinical respiratory illness (`CRI`).


A subset of data from the study is presented in the following `masks` data frame.

```{r}
masks<-tibble(
  infected = c(
    replicate(n = 28,"yes" ),
    replicate(n = 552,"no" ), 
    replicate(n = 43,"yes" ), 
    replicate(n = 526,"no" ),
    replicate(n = 32,"yes" ), 
    replicate(n = 426,"no" )),
  mask_type = c(
    replicate(580, "medical"), 
    replicate(569, "cloth"), 
    replicate(458, "control"))
)
```

# Hypotheses and Data Exploration

One questions researchers are interested in resolving is whether the two masks were equally effective in preventing influenza-like illnesses. Without any additional information, the null hypothesis is that there is no difference between medical masks and cloth masks, and the alternative hypothesis is that there is a difference.

### Exercise 9

Using p_m and p_c to denote the theoretical infection rate among all health-care workers wearing **m**edical masks and **c**loth masks (respectively), write out the null and alternative hypotheses in mathematical notation.


********************************************

********************************************


### Exercise 10

 Use `ggplot` to create an appropriate data visualization that displays the proportion of health care workers that were infected in the 3 groups.


********************************************

********************************************


# Testing the hypothesis

We'll use `infer` to run this hypothesis test. 

## Step 1: Calculate the observed difference

While we could use `dplyr` verbs to calculate the difference in infection proportions between the `medical` and `cloth` masks groups, the `infer` package makes this calculation straightforward.


```{r}
masks_med_cloth<-masks %>% 
  filter(mask_type %in% c("medical", "cloth"))

obs_diff_masks <- masks_med_cloth %>%
  specify(infected ~ mask_type, success = "yes") %>% 
  calculate( stat = "diff in props", order = c("medical", "cloth"))
```

Notice in the above code that we needed to first restrict our attention only to those health care workers who wore cloth or medical masks (what group did we remove?)

We also needed to chose an order for the levels of the `mask_type` variable. The particular order we used doesn't matter, as long as we are consistent in all subsequent calculations.

### Exercise 11

What is the observed difference in proportion calculated above? Before you conduct any further analysis, make a guess: is this difference statistically significant?

********************************************

********************************************


## Step 2: Generate the permuted samples

This step involves generating simulated values of differences in proportions *as if* we lived in a world where there truly is no difference in effectiveness between the two mask types. Returning to the idea of permutation, this is akin to shuffling the `mask_type` label among all healthy care workers. 

```{r}
null_masks_samples <- masks_med_cloth %>% 
  specify(infected ~ mask_type, success = "yes") %>% 
  hypothesize(null = "independence") %>% 
  generate( reps = 2000, type = "permute")
```

### Exercise 12

What was the size of each shuffled sample? How many times did we shuffle? How many rows are there in the `null_masks_samples` data frame?

********************************************

********************************************

## Step 3: Create a null distribution of differences in proportions

The following code calculates the differences in infection rate between `medical` and `cloth` mask health-care workers.

```{r}
null_masks_proportions <- null_masks_samples %>% 
  calculate( stat = "diff in props", order = c("medical", "cloth"))
```

### Exercise 13

How many rows are there in the `null_masks_proportions` data frame? Why does this differ from the number of rows in the `null_masks_samples` data frame?

********************************************

********************************************

## Step 4: Visualize how the observed difference compares to the null distribution.

The following code plots the differences in proportions for each of the shuffled replicates.

```{r}
null_masks_proportions %>% 
  visualize()
```

**Note that zero is the center of this null distribution.** The null hypothesis is that there is no difference in infection rates between mask types. In the null distribution, zero was the most common difference, because we broke any existing relationship in the sample---infected healthcare workers were sorted into `medical` and `cloth` groups at random.

### Exercise 14

Based on this graphic, what do you think is a plausible range of values for the difference in proportions under the null hypothesis?

********************************************

********************************************


We can also visualize where the observed difference lies in the null distribution using the following code:

```{r warning = F}
null_masks_proportions %>% 
  visualize( ) +
  shade_p_value(obs_stat = obs_diff_masks, direction = "both")
```

**Note that since we used a two-sided alternate hypothesis, we set `direction = "both"` to indicate we should look for extreme values in both tails of the distribution. **

### Exercise 15

Based on this graphic, does the observed difference seem plausible if the null hypothesis were true? Explain.

********************************************

********************************************


## Step 5: Calculate a p-value

Use the following code to calculate the p-value for the observed difference:

```{r}
null_masks_proportions %>% 
  get_p_value(obs_stat = obs_diff_masks, direction = "both")
```

The p-value indicates the probability that we would see a difference in proportion equal to, or more extreme, as the one we did if there were truly no difference in effectiveness between mask types.


### Exercise 16

Compare the p-value above to the standard significance level of alpha = 0.05. Does this test give evidence to reject the null hypothesis. Does this imply that the null hypothesis is false?

********************************************

********************************************

## Exercise 17

The researchers in <https://bmjopen.bmj.com/content/5/4/e006577.long> actually drew the opposite conclusion, arguing that the data indeed showed a statistically significant difference in effectiveness between cloth and medical masks. This is because the researchers performed further statistical analysis to correct for the fact that the sample was generated using **cluster sampling** rather than **simple random sampling** (i.e. each health worker in the same ward was required to wear the same type of mask).

Theorize two reasons why the researchers may have decided to use cluster sampling rather than simple random sampling.

********************************************

********************************************


