---
title: "Lab 7: Sampling"
author: "Insert Name"
output:
  pdf_document
 
---

# Due: 12am night before lab

# Instructions
Work through the problems below, and submit this document as a knitted .pdf file with answers to the Lab 6 assignment on Gradescope. 

**In some problems, an incomplete R chunk has been provided. To evaluate the chunk, follow the given instructions, and then change `eval=FALSE` to `eval=TRUE`.

**For each problem, put your solution between the bars of red stars.**


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


 


## Background

In this lab, we will gain experience sampling from a population. We will use some data from the `pdx_trees` dataset that we have been using.

For this lab, the **population** of interest will be **ALL** 964 Western Redcedar trees in Portland parks. As an analogy to the tactile sampling you did in class,  Portland parks are the "deck" and the 964 trees are the 52 cards.

If you get stuck as you are working through this Problem Set, you may find it helpful to re-read Chapter 9 in ModernDive, in particular section 7.2 on "Virtual Sampling" and section 7.3 on "Sampling Framework". 


### Symbols and Notation 

Symbol           | POPULATION PARAMETER               | SAMPLE PARAMETER
---------------- | ---------------------------------- | --------------
Number of cases  | $N$                                | $n$
Proportion       | $p$                                | $\hat{p}$
Standard error   | $SE$                               | $\widehat{SE}$

***

## Setup

First load the necessary packages:

```{r}
library(ggplot2)
library(dplyr)
library(forcats)
library(moderndive)
```

You can load and take a `glimpse()` of the  `pdxTrees_parks` data set in the `forcats` package like so: 

```{r message = F, results = "hide"}
# Grab trees near Reed
library(pdxTrees)
pdxTrees_parks <- get_pdxTrees_parks()
glimpse(pdxTrees_parks)
```


Be sure to take a look at the data in the **viewer**. And type `?get_pdxTrees_parks()` into the **console** to see a description of the variables in this data set. 

 

## Exploratory data wrangling

This data set includes many variables. To start, we will restrict our analysis to only Western Redcedar trees, and to only the variable indicating the `Condition` of each tree 


```{r}
trees <- pdxTrees_parks %>% 
  filter(Common_Name == "Western Redcedar") %>% 
  select(Condition)
```


The following shows the different levels of `Condition` of the Western Redcedars in Portland parks:
```{r}
trees  %>% 
  distinct(Condition) 
```

 

### Setting a seed

**Setting a seed:** We will take some random samples in this Problem Set. In order to make sure R takes the same random sample every time you run your code, you can do what is called "setting a seed". Do this in any code chunk that you take a random sample! 

You can set a seed like so. Any number will do. (You do not need to run this right now. This chunk show shows you how.)
```{r}
set.seed(42)
```

 

## The true population proportion $p$ of Western Redcedars

Again, for this exercise, the **population** of interest will be **ALL** 964 Western Redcedar trees in Portland parks.  Since we have data on **ALL** 964 Western Redcedar trees in Portland parks, we can compute the **exact population proportion $p$ of Western Redcedars in good condition directly** like so, using **ALL** the data: 

```{r}
trees %>% 
  summarize(good = sum(Condition == "Good"), 
            N = n()) %>% 
  mutate(p = good / N)
```

 

> Note: that we use $N$ for the size of the full population of 964 trees, and $p$ because we are calculating the TRUE population proportion $p$

> Note that no inference to the population is needed. We do not need to use a **sample** to try to infer something about the **true population proportion $p$** of Western Redcedars in good condition in Portland parks. We know that $p$ = 0.13. 
Of course, this lab is not a realistic reflection of a real life problem. However, we are *simulating* the act of sampling from this population to understand and study how factors like sample size influence **sampling variation**.

 

## Demo: sampling 50 trees

### a) Estimating $\hat{p}$ from a single sample

We are first going to use random sampling to **ESTIMATE** the true **population** proportion $p$ of the Western Redcedars in good condition with only a **sample** of 50 trees. 

> This will represent a situation of only having the resources to record 50 Western Redcedar trees in Portland parks.
Be sure to look at the results in the viewer. Remember, you can set the seed to whatever value you like. 

```{r}
set.seed(42)
n50_1rep <- trees %>% 
  rep_sample_n(size = 50, reps = 1)
```

 

Next, let's calculate the **sample proportion** $\hat{p}$ of good condition Western Redcedars in our sample of 50 trees. 

```{r}
n50_1rep %>% 
  summarize(Good = sum(Condition == "Good"), 
            n = n()) %>% 
  mutate(p_hat = Good / n)
```


This sample proportion $\hat{p}$ is an **ESTIMATE**, and our **best guess** of what the **true population** proportion $p$ of good condition Western Redcedars in Portland parks, based on a sample of only 50 trees.  It is reasonably close to the true population proportion $p = 0.12$ we calculated from the full population. 

# Exercise 1

Change the seed value in the `set.seed()` in the code chunk above that takes a sample so it is something different from "42", and rerun the code to estimate $\hat{p}$.  Ask two of your classmates what their estimate of $\hat{p}$ was. How do the $\hat{p}$ estimates from different samples compare? 

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

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


### b) Estimating $\widehat{SE}$ from a single sample

Typically we only have the opportunity to collect **one sample** for our study, and so we have to use the amount of variability in our **single sample** as an estimate of the amount of variability we might expect in our results if we had taken a random sample of 50 different trees. The $\widehat{SE}_{\hat{p}}$ serves as an **ESTIMATE**  of **sampling variability** if you only have a **single sample**. The formula for estimating the standard error of $\hat{p}$ is the following:

$$\widehat{SE}_{\hat{p}}  \approx  \sqrt{\frac{\hat{p} \times (1-\hat{p})}{n}}$$

> Note that we use $n$ for the size of the sample, that p "wears a hat", like so: $\hat{p}$ because we are ESTIMATING a proportion based on only a sample, and that the SE "wears a hat" as well because we ESTIMATING $\widehat{SE}$ based on only a sample. 
 

The standard error of $\hat{p}$ can be estimated in R like so:

```{r}
n50_1rep %>% 
  summarize(Good = sum(Condition == "Good"), 
            n = n()) %>% 
  mutate(p_hat = Good / n,
         se_hat = sqrt(p_hat * (1 - p_hat) / n))
```

 

## Demo: Generating a sampling distribution of $\hat{p}$


If you ran the code chunk that takes a random sample of 50 cases a thousand more times....and wrote down every $\hat{p}$ you got, you would have what is called a "sampling distribution". 

<br> 

> A sampling distribution shows every [or nearly every!] possible result a sampling statistic can have under every [or nearly every!] possible sample **of a given sample size** from a population.  
<br> 

### a) Sampling distribution of $\hat{p}$ for n = 50

Instead of running the sampling code chunk for n = 50 over and over, we can "collect" 1000 samples of n = 50 really easily in R. The following code chunk takes 1000 **different** samples of n = 50 and stores them in the data frame `n50_1000rep`:

```{r}
set.seed(42)
n50_1000rep <- trees %>% 
  rep_sample_n(size = 50, reps = 1000)
```

 
Be sure to look at `n50_rep1000` in the data viewer to get a sense of these 1000 samples look like.

 

# Exercise 2a

What is the name of the column that identifies which of the 1000 samples each row is from?

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

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

# Exercise 2b

What is the sample size n for each of the 1000 samples we took? (i.e. how many trees are sampled in each replicate)?

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

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

The following code chunk calculates the sample proportion $\hat{p}$ of good condition Western Redcedars for each of the **1000 samples** 

```{r}
p_hat_n50_1000rep <- n50_1000rep %>% 
  group_by(replicate) %>% 
  summarize(Good = sum(Condition == "Good"), 
            n = n()) %>% 
  mutate(p_hat = Good / n)
```

Take a look at the first five rows of the results:

```{r}
p_hat_n50_1000rep %>%
  slice(1:5)
```

### b) Visualize the sampling distribution of $\hat{p}$ for n = 50

We can plot the **sampling distribution** of these 1000 $\hat{p}$ estimates of good Western Redcedars with a histogram, like so:

```{r}
ggplot(p_hat_n50_1000rep, aes(x = p_hat)) +
  geom_histogram(binwidth = 0.02, color = "black", fill = "aquamarine3") +
  labs(x = "Sample proportion of Western Redcedars in Good Condition", 
       title = "Sampling distribution of p-hat based on n = 50") 
```

 

# Exercise 3 

Based on your histogram, what appeared to be a very common value of $\hat{p}$? What was a very uncommon value?

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

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

### c) Mean and standard error of the sampling distribution of $\hat{p}$ for n = 50

Finally we can estimate  the mean of the sampling distribution by calculating the mean of all 1000 $\hat{p}$ estimates, and the standard error of the sampling distribution by calculating the standard deviation of all 1000 $\hat{p}$ values like so: 

```{r}
p_hat_n50_1000rep %>% 
  summarize(M_p_hat = mean(p_hat), 
            SE_p_hat = sd(p_hat))
```

 

# Exercise 4

How do these values compare to the estimates we got above of $\hat{p}$ and $\widehat{SE}$ of $\hat{p}$ for good condition Western Redcedars based on your **single** sample of 50 trees earlier in this Problem Set?  


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

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


# Exercise 5a

Use the `rep_sample_n` function to collect 1000 virtual samples of size *n* = 15. **BE SURE TO NAME YOUR SAMPLE SOMETHING NEW, TO ENSURE YOU CAN DISTINGUISH IT FROM THE n = 50 SAMPLE ABOVE!** Use seed 910 to generate your sample. Include your code, but you don't need to print the output of the code.

****************************************
 
```{r}
set.seed(910)
n15_1000rep <- trees %>% 
  rep_sample_n(size = 15, reps = 1000)
```

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


# Exercise 5b

Calculate sample proportion $\hat{p}$ of good condition Western Redcedars for each replicate of your n = 15 sampling. Include your code, but you don't need to print the output of the code.

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

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

# Exercise 5c

Visualize the sampling distribution of $\hat{p}$ from your n = 15 sampling with a histogram.

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

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

# Exercise 5d

Calculate the mean of the n = 15 sampling distribution, and the standard error of the n = 15 sampling distribution.

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

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

# Exercise 6a

How does the standard error of the n = 15 sampling distribution compare to the standard error of the n = 50 sampling distribution?

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

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

# Exercise 6b

Explain any observed differences from Exercise 6a.

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

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

# Exercise 7a

Use the `rep_sample_n` function to collect 1000 virtual samples of size *n* = 600. **Note: BE SURE TO NAME YOUR SAMPLE SOMETHING NEW, TO ENSURE YOU CAN DISTINGUISH IT FROM THE n = 50, and n = 15 SAMPLES ABOVE!** Use seed 84 to generate your sample. Include your code, but you don't need to print the output of the code.

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

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

# Exercise 7b

Calculate the proportion $\hat{p}$ of good condition Western Redcedar trees for each replicate of your n = 600 sampling. Include your code, but you don't need to print the output of the code.

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

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

# Exercise 7c

Calculate the mean of the n = 600 sampling distribution, and the standard error of the n = 600 sampling distribution.

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

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



# Exercise 7d

Was there more **variability** from sample to sample when we took a sample size of 600 or a sample size of 50? **Explain what evidence you have for assessing this**

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

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


# Exercise 8

Which sampling distribution looked more normally distributed (bell shaped and symmetric); the one built on n = 15, 50 or 600? **Why?**

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

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


 
 
## Estimating $\hat{p}$ and the standard error of $\hat{p}$ from a single sample (revisited)

In most instances, we do not have access to the full population as we did in this Portland trees data; instead we have to take a **sample** to try to say something about the **larger population**. Furthermore, in the real world, we typically only take a **single** sample from the population, due to time or money constraints. 

So how do we **ESTIMATE** a $\hat{p}$ and a standard error of $\hat{p}$ when we only have a single sample, and not 1000 repeated samples? As demonstrated at the very beginning of the Problem Set we:

* estimate $\hat{p}$ from the sample 
* use the formula for the standard error of $\hat{p}$ below, to estimate SE based on a single sample 


$$\widehat{SE}_{\hat{p}}  \approx  \sqrt{\frac{\hat{p} \times (1-\hat{p})}{n}}$$

 
# Exercise 9

Imagined we collected only a single small sample of 15 respondents like so: 

```{r}
set.seed(53)
n15_1rep <- trees %>% 
  rep_sample_n(size = 15, reps = 1)
```


Following the example from the beginning of the Problem Set (roughly line 140), estimate the **sample proportion** $\hat{p}$ of good condition Western Redcedar trees based on `n15_1rep`, as well as the **standard error of $\hat{p}$**

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

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

> You should have gotten a value reasonably close to the estimate we made earlier from our sampling distribution for n = 15! Note that when you must estimate a standard error from **only a single sample**, the formula **contains the sample size, n**. The larger the sample size n, the larger the number in the denominator of the SE formula. 


# Exercise 10 

Fill in the R Markdown table below with all the standard errors you computed for this problem set. In other words:

1. Replace `x` with the standard error you obtained by taking the standard deviation of the n = 15 sampling distribution 
2. Replace `a` with the standard error you obtained for a single sample of n = 15 using the mathematical formula
3. etc. 

When you are done, make sure all the `|` in the table still line up so your results print out in a table! 

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

Sample size n  | SE via sd of sampling distribution | SE via one sample and formula
-------------- | ---------------------------------- | --------------
x              | x.                                 | x
x              | x                                  | x

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

 

# Exercise 11

Based on what you observed for Exercise 10, **IF** you collected a single sample from 600 trees, do you think the standard error will be smaller or larger than the one you calculated for n = 15. **Explain your reasoning** (Note: if you are not sure you can collect a sample and calculate the standard error)

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

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


As we've seen so far, the size of the sample determines how well we can estimate a parameter of a population.  However, what role does the number of sample repetitions play in dertiming how well we understand the sampling distribution?  Consider the sample of 50 trees obtained in the demo at the beginning of this lab.  

## Exercise 12a

Using a seed value of 42, take 100 **different** samples of n = 50, calculate the mean fo the sampling distribution, the standard deviation of the sampling distribution, and visualize the sampling distribution.

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

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


## Exercise 12b

Using a seed value of 42, take 10,000 **different** samples of n = 50, calculate the mean fo the sampling distribution, the standard deviation of the sampling distribution, and visualize the sampling distribution.

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

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

## Exercise 12c

Do the mean and standard deviation of the sampling distributions change as the number of repetitions change? Why or why not?

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

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

## Exercise 12d

Compare and contrast the visualizations for 100 repetitions, 1000 repetitions and 10,000 repetitions.  What do you see.

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

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



