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

# Instructions

**Due: Before lab Thursday or Friday March 18 or 19**

Read through the following document (it may be helpful to have both the .rmd and .pdf files open side-by-side).  Work through the problems below, and submit this document as a knitted .pdf file with answers to the Lab 6 on [Gradescope](https://www.gradescope.com/courses/227442/). 


```{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 imitate the sampling activity we did in class on Monday. 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 `types of conifer trees`Condition` of Portland parks Western Redcedars:

```{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 people. 

```{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 1a

Ask two of your classmates what their estimate of $\hat{p}$ was. How do the $\hat{p}$ estimates from different samples compare? 

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

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

# Exercise 1b

**Why** did everyone get a different estimate?

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

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

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

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

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


# 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
-------------- | ---------------------------------- | --------------
15             | x                                  | a
50             | y                                  | b 

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

 

# 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)

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

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

## Exercise 12

The quintessential activity that Reed alumni do upon returning to campus is head up to the thesis tower to check to see if anyone has checked out their senior thesis. They’re wondering: what is the chance that someone was interested in my work? To answer that question and learn more generally about the way in which senior theses are checked out, you will explore a sample of theses from the thesis tower.

Before exploring the dataset, let's consider sampling theses from the thesis tower.  The the tower has two floors and the thesis are organized in alphabetical order by author's last name, across both floors.

Imagine that we would like to collect a sample with the following variables:

* barcode number
* year of publication
* number of checkouts
* number of pages 
* division (all departments are grouped into divisions, found at <http://www.reed.edu/academics.html>)


A. First, estimate the number of theses stored in the tower.  It is ok if you have never been to the thesis tower, just make your best guess. Enter your estimate in the space below.

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

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

B. Provide a detailed plan of how you would collect a sample of data that is representative and efficient.  Imagine that you have the resources to collect data on around 40 theses.  You might consider simple random sampling, stratified sampling, cluster sampling, a hybrid, or another random technique all together.  Your plan should be explicit enough that someone else could follow in your footsteps next semester and know how to collect a sample in an identical manner. 


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

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

C. Justify why you think your sampling method is representative and efficient.


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

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

Unfortunately, due to the pandemic, we will not visit the thesis tower and implement our plan.  We will however, look at previously obtained data from the thesis tower.

Run the following code chunk to load the data set into R, make sure to change the seed value to obtain a unique sample:

```{r echo = FALSE, eval = TRUE, message = FALSE}

theses <- read_csv("https://reed-statistics.github.io/math141s21/data/theses.csv")

set.seed(42)

sample_theses <- theses %>% 
  rep_sample_n(size = 40, reps = 1)
```


The first time you run this code chunk, you will need to cache OAuth access credentions. Navigate to the console and type "1". Your webbrowser will then provide an authorization code. Copy and paste the code into the console. After this, you will not need to do anything special to run the code in the future.

In the space below, run the `glimps(theses)` chunk to verify that you've loaded your data into R.

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

```{r}
glimpse(sample_theses)
```

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


D.  Create a plot to visualize the relationship between the number of checkouts and the year of publication. Describe the structure that you find. How many checkouts can you expect in your thesis if you come back for your 10th reunion? Your 20th?

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

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

E.  Create a plot that addresses the question: “Is there a trend in length of theses over time?”. Answer that question also in words.

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

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

F.  Create a plot that addresses the question: “Are theses from the different divisions checked out at different rates?”. Answer that question also in words.

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

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

G.  Combine the variables from the plots in parts D. and F. into a single plot that shows the trends in checkouts over time, separated out for each division. Describe any structure that you see.
    
*************************************************

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




## Exercise 13

Suppose Jonathan offers an amazing divination service, where you ask him a yes or no question and he responds with an appropriate answer. Amazingly, it turns out that Jonathan offers incorrect answers only 5% of the time on average! Of course, this a random process, so if you ask him 100 questions, he may not answer exactly 95 of them correctly. 

To study the variability in the number of answers he gets correct, you first ask him a sequence of 20 questions, repeated 10000 times, and then a sequence of 100 questions, repeated 10000 times.

Run the following code to generate data frames recording your findings.


```{r}
## This is the correct code
set.seed(42)

divination<-tibble(Answer = c("Incorrect Answer", replicate(19, "Correct Answer")))  

questions20<-rep_sample_n(divination, size = 20, reps = 10000, replace = TRUE)
questions100<-rep_sample_n(divination, size = 100, reps = 10000, replace = TRUE)
```

A. Use the data above to create sampling distributions for the proportion correct answers Jonathan gives when asked (a) twenty questions, and (b) one-hundred questions.


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

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




B. Based on the histograms, what appears to be the center of both distributions? Verify by computing the exact means of each data sets.

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

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

C. Approximate the standard error of each distribution. Which standard error is smaller?


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

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

D. Suppose you take one more sample in each case. Would a sample proportion of 100% questions correct be reasonably likely from a sample of size 20? Would it be reasonably likely from a sample of size 100? Explain.

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

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


## Exercise 14

In each of the following scenarios, identify: the implied population and population size (if listed), the sample and sample size, the parameter of interest, the statistic used as a point estimate, and whether the sample represents an SRS of the population.

A. A recent article in the Quest stated that Reed students get an average of 5.5 hours of sleep each night. A student was skeptical of this claim and decided to survey the 25 members of a Math 141 class. The student found that these Math 141 students sleep 6.25 hours on average per night.

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

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

B. Suppose you want to determine the percentage of videos on YouTube that are cat videos. Since 82.2 years worth of content is uploaded to YouTube every day, there is no way to perform an exhuastive search. Instead, you use a random video picker (like <https://randomtube.net/>) to select 200 videos for you, and you find that 9 are cat videos.

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

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

C. The Gallup Poll uses a procedure called random digit dialing, which creates phone numbers based on a list of all area codes in America in conjunction with the associated number of residential households in each area code. Based on this method, the Gallup Poll calls 5440 phone numbers and finds that only 1221 are actually assigned to residential households. Of these 1221 numbers, the Gallup poll successfully gathers yearly household income information from 912 individuals, finding that the yearly household income for these individuals was \$32,000. 

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

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



 