Lab 5: Inference Lab

library(tidyverse)
library(openintro)
library(infer)

global_monitor <- tibble(
  scientist_work = c(rep("Benefits", 80000), rep("Doesn't benefit", 20000))
)

In this lab we will look at a dataset containing responses to the question “Do you believe that the work scientists do benefit people like you?”

Visualize the distribution of responses with a barplot:

ggplot(global_monitor, aes(x = scientist_work)) +
  geom_bar() +
  labs(
    x = "", y = "",
    title = "Do you believe that the work scientists do benefit people like you?"
  ) +
  coord_flip() 

Compare with summary statistics:

global_monitor %>%
  count(scientist_work) %>%
  mutate(p = n /sum(n))

## # A tibble: 2 x 3
##   scientist_work      n     p
## * <chr>           <int> <dbl>
## 1 Benefits        80000   0.8
## 2 Doesn't benefit 20000   0.2

Going forward, we will assume these responses represent the entire population, but also that we do not have access to the full set of responses. The idea is to examine what we can infer using random sampling from a larger population which we would generally not have access to. We’ll use sample_n to generate a random sample from the dataset:

#Default is sampling without replacement
samp1 <- global_monitor %>%
  sample_n(50)
samp1

## # A tibble: 50 x 1
##    scientist_work 
##    <chr>          
##  1 Benefits       
##  2 Benefits       
##  3 Benefits       
##  4 Benefits       
##  5 Benefits       
##  6 Benefits       
##  7 Doesn't benefit
##  8 Benefits       
##  9 Benefits       
## 10 Benefits       
## # … with 40 more rows

table(samp1$scientist_work)

## 
##        Benefits Doesn't benefit 
##              41               9

Exercise 1

Describe the distribution of responses in this sample. How does it compare to the distribution of responses in the population. Hint: Although the sample_n function takes a random sample of observations (i.e. rows) from the dataset, you can still refer to the variables in the dataset with the same names. Code you presented earlier for visualizing and summarising the population data will still be useful for the sample, however be careful to not label your proportion p since you’re now calculating a sample statistic, not a population parameters. You can customize the label of the statistics to indicate that it comes from the sample.

samp_monitor <- tibble(samp1 = c(rep("Benefits", 37), rep("Doesn't benefit", 13)))

Visualize the sample distribution

ggplot(samp_monitor, aes(x = samp1)) +
  geom_bar() +
  labs(
    x = "", y = "",
    title = "Do you believe that the work scientists do benefit people like you?"
  ) +
  coord_flip() 

How does this plot correspond with the summary data?

samp_monitor %>%
  count(samp1) %>%
  mutate(p_hat = n /sum(n))

## # A tibble: 2 x 3
##   samp1               n p_hat
## * <chr>           <int> <dbl>
## 1 Benefits           37  0.74
## 2 Doesn't benefit    13  0.26

This is a pretty good estimate of the population mean, considering the relatively small sample size

Exercise 2

Would you expect the sample proportion to match the sample proportion of another student’s sample? Why, or why not? If the answer is no, would you expect the proportions to be somewhat different or very different? Ask a student team to confirm your answer.

Not necessarily, because this is a relatively small random sample of the population. I would expect such samples to form a distribution of their own.

Exercise 3

Take a second sample, also of size 50, and call it samp2. How does the sample proportion of samp2 compare with that of samp1? Suppose we took two more samples, one of size 100 and one of size 1000. Which would you think would provide a more accurate estimate of the population proportion?

samp2 <- global_monitor %>%
  sample_n(50)
table(samp2$scientist_work)

## 
##        Benefits Doesn't benefit 
##              38              12

I would expect the larger sample to more accurately reflect the greater population, since the sample itself more closely resembles the population as the size increases.

Now we’ll perform repeated random sampling with the rep_sample_n function:

sample_props50 <- global_monitor %>%
                    infer::rep_sample_n(size = 50, reps = 15000, replace = TRUE) %>%
                    count(scientist_work) %>%
                    mutate(p_hat = n /sum(n)) %>%
                    filter(scientist_work == "Doesn't benefit")

and visualize it with a histogram:

ggplot(data = sample_props50, aes(x = p_hat)) +
  geom_histogram(binwidth = 0.02) +
  labs(
    x = "p_hat (Doesn't benefit)",
    title = "Sampling distribution of p_hat",
    subtitle = "Sample size = 50, Number of samples = 15000"
  )

Exercise 4

How many elements are there in sample_props50? Describe the sampling distribution, and be sure to specifically note its center. Make sure to include a plot of the distribution in your answer.

There are 15,000 elements in sample_props50, each of which is a measure of the proportion of “Doesn’t benefit” replies in a random sample of 50. From the plot, we can quickly see that the mean sample proportion is about 0.2, which matches our “population”.

sample_props50

## # A tibble: 15,000 x 4
## # Groups:   replicate [15,000]
##    replicate scientist_work      n p_hat
##        <int> <chr>           <int> <dbl>
##  1         1 Doesn't benefit    10  0.2 
##  2         2 Doesn't benefit     9  0.18
##  3         3 Doesn't benefit    14  0.28
##  4         4 Doesn't benefit     9  0.18
##  5         5 Doesn't benefit    16  0.32
##  6         6 Doesn't benefit     6  0.12
##  7         7 Doesn't benefit    13  0.26
##  8         8 Doesn't benefit     9  0.18
##  9         9 Doesn't benefit    15  0.3 
## 10        10 Doesn't benefit     9  0.18
## # … with 14,990 more rows

Exercise 5

To make sure you understand how sampling distributions are built, and exactly what the rep_sample_n function does, try modifying the code to create a sampling distribution of 25 sample proportions from samples of size 10, and put them in a data frame named sample_props_small. Print the output. How many observations are there in this object called sample_props_small? What does each observation represent?

sample_props_small <- global_monitor %>%
                    infer::rep_sample_n(size = 10, reps = 25, replace = TRUE) %>%
                    count(scientist_work) %>%
                    mutate(p_hat = n /sum(n)) %>%
                    filter(scientist_work == "Doesn't benefit")
sample_props_small

## # A tibble: 23 x 4
## # Groups:   replicate [23]
##    replicate scientist_work      n p_hat
##        <int> <chr>           <int> <dbl>
##  1         1 Doesn't benefit     2   0.2
##  2         2 Doesn't benefit     2   0.2
##  3         3 Doesn't benefit     2   0.2
##  4         5 Doesn't benefit     2   0.2
##  5         6 Doesn't benefit     2   0.2
##  6         7 Doesn't benefit     2   0.2
##  7         8 Doesn't benefit     2   0.2
##  8         9 Doesn't benefit     3   0.3
##  9        10 Doesn't benefit     1   0.1
## 10        11 Doesn't benefit     2   0.2
## # … with 13 more rows

There are 25 observations in the sample_props_small object, each one
representing the proportion of “Doesn’t benefit” replies in a sample of size 10.

Exercise 6

Use the app below to create sampling distributions of proportions of Doesn’t benefit from samples of size 10, 50, and 100. Use 5,000 simulations. What does each observation in the sampling distribution represent? How does the mean, standar error, and shape of the sampling distribution change as the sample size increases? How (if at all) do these values change if you increase the number of simulations? (You do not need to include plots in your answer.)

The standard deviation of the sample mean decreases with sample size. The mean stays approximately consistent.

Exercise 7

Take a sample of size 15 from the population and calculate the proportion of people in this sample who think the work scientists do enchances their lives. Using this sample, what is your best point estimate of the population proportion of people who think the work scientists do enchances their lives?

samp3 <- global_monitor %>%
  sample_n(15)
table(samp3$scientist_work)

## 
##        Benefits Doesn't benefit 
##              12               3

This sample suggests that 12/15 = 4/5 of people believe that scientists’ work does benefit them.

Exercise 8

Since you have access to the population, simulate the sampling distribution of proportion of those who think the work scientists do enchances their lives for samples of size 15 by taking 2000 samples from the population of size 15 and computing 2000 sample proportions. Store these proportions in as sample_props15. Plot the data, then describe the shape of this sampling distribution. Based on this sampling distribution, what would you guess the true proportion of those who think the work scientists do enchances their lives to be? Finally, calculate and report the population proportion.

sample_props15 <- global_monitor %>%
                    infer::rep_sample_n(size = 15, reps = 2000, replace = TRUE) %>%
                    count(scientist_work) %>%
                    mutate(p_hat = n /sum(n)) %>%
                    filter(scientist_work == "Benefits")

ggplot(data = sample_props15, aes(x = p_hat)) +
  geom_histogram(binwidth = 0.02) +
  labs(
    x = "p_hat (Benefits)",
    title = "Sampling distribution of p_hat",
    subtitle = "Sample size = 15, Number of samples = 2000"
  )

mean(sample_props15$p_hat)

## [1] 0.7961333

Exercise 9

Change your sample size from 15 to 150, then compute the sampling distribution using the same method as above, and store these proportions in a new object called sample_props150. Describe the shape of this sampling distribution and compare it to the sampling distribution for a sample size of 15. Based on this sampling distribution, what would you guess to be the true proportion of those who think the work scientists do enchances their lives?

sample_props150 <- global_monitor %>%
                    infer::rep_sample_n(size = 150, reps = 2000, replace = TRUE) %>%
                    count(scientist_work) %>%
                    mutate(p_hat = n /sum(n)) %>%
                    filter(scientist_work == "Benefits")

ggplot(data = sample_props15, aes(x = p_hat)) +
  geom_histogram(binwidth = 0.02) +
  labs(
    x = "p_hat (Benefits)",
    title = "Sampling distribution of p_hat",
    subtitle = "Sample size = 150, Number of samples = 2000"
  )

mean(sample_props150$p_hat)

## [1] 0.79992

Exercise 10

Of the sampling distributions from 2 and 3, which has a smaller spread? If you’re concerned with making estimates that are more often close to the true value, would you prefer a sampling distribution with a large or small spread?

We would expect for the larger sample size of 150 to lead to a smaller spread. Theoretically, s = sigma/(sqrt(n)), implying that the sample standard deviation decreases with increasing n. Let’s check this anyways:

sd(sample_props15$p_hat) > sd(sample_props150$p_hat)

## [1] TRUE

All else being equal, we would prefer the larger sampler size.