Foundations for statistical inference - Sampling distributions


Getting Started


Load packages

We are operating in the tidyverse.

# Load packages ----------------------------------------------------------------
library(tidyverse)
library(openintro)
library(infer)  # Needed for resampling
library(shiny)  # Needed for the shiny app


The data

A 2019 Gallup report states the following:

The premise that scientific progress benefits people has been embodied in discoveries throughout the ages – from the development of vaccinations to the explosion of technology in the past few decades, resulting in billions of supercomputers now resting in the hands and pockets of people worldwide. Still, not everyone around the world feels science benefits them personally.

Source: World Science Day: Is Knowledge Power?

The Wellcome Global Monitor finds that 20% of people globally do not believe that the work scientists do benefits people like them. In this lab, you will assume this 20% 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 20,000 (20%) of the population think the work scientists do does not benefit them personally and the remaining 80,000 think it does.

global_monitor <- tibble(
  scientist_work = c(rep("Benefits", 80000), rep("Doesn't benefit", 20000))
)
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()


Exercises


Exercise 1

Describe the distribution of responses in this sample. How does it compare to the distribution of responses in the population.

The proportion of people in the sample who don’t believe Science benefits them is 18% which is close the total dataset’s 20%. I hope that 20% have cell phones. Thank you, Science!

# Generate statistics for the data
global_monitor %>%
  count(scientist_work) %>%
  mutate(p = n /sum(n))
## # A tibble: 2 × 3
##   scientist_work      n     p
##   <chr>           <int> <dbl>
## 1 Benefits        80000   0.8
## 2 Doesn't benefit 20000   0.2
# Create the sample
set.seed(7585)
samp1 <- global_monitor %>%
  sample_n(50)

# Generate statistics for the sample
samp1_p_hat <- samp1 %>% 
  count(scientist_work) %>% 
  mutate(p_hat = n /sum(n)) %>% 
  filter(scientist_work == "Doesn't benefit") %>% 
  pull(p_hat) %>% 
  round(2)
samp1_p_hat
## [1] 0.18


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.

Based on the formula for the standard error below, the error is inversely proportional to the square root of the sample size.

\[ σ /\sqrt{n} \]

With a sample size of 50, that’s equivalent to dividing the standard error of the sample mean from the true mean by roughly 7. I would expect my sample proportion to be close another student’s sample proportion but not exactly the same. I didn’t ask another student but ran a second sample in the next exercise where I got a p_hat of 80% instead of the 82% from sample 1; close but not necessarily exactly the same.


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?

With the second sample I get a p_hat of 20% instead of the 18% from sample 1; close but not necessarily exactly the same. As we increase the sample size the standard error of the sample mean will go down so size 1000 will tend to be more accurate than 100, than 50.

# Create a second sample
set.seed(7686)
samp2 <- global_monitor %>%
  sample_n(50)

# Generate statistics for sample 2
samp2_p_hat <- samp2 %>% 
  count(scientist_work) %>% 
  mutate(p_hat = n /sum(n)) %>% 
  filter(scientist_work == "Doesn't benefit") %>% 
  pull(p_hat) %>% 
  round(2)
samp2_p_hat
## [1] 0.2


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.

Tibble-wise there are 15,000 rows x 4 columns or 60,000 elements in sample_props50. Experiment-wise there are 15000 repetitions x 50 samples or 750,000 individuals sampled, 7.5 times the total population! The sampling distribution looks normally distributed with a center at the true population mean of 20%.

# Take a sample distribution of the data
set.seed(4)
sample_props50 <- global_monitor %>%
                    rep_sample_n(size = 50, reps = 15000, replace = TRUE) %>%
                    count(scientist_work) %>%
                    mutate(p_hat = n /sum(n)) %>%
                    filter(scientist_work == "Doesn't benefit")

# Plot results ------------------------------------------------
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 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?

There are 25 observations in this sampling distribution where each observation is a sample of size 10 where we count how many of the ten don’t believe science benefits them and represent that as a proportion.

# Take a sample distribution of the data
set.seed(5)
sample_props_small <- global_monitor %>%
                    rep_sample_n(size = 10, reps = 25, replace = TRUE) %>%
                    count(scientist_work) %>%
                    mutate(p_hat = n /sum(n)) %>%
                    filter(scientist_work == "Doesn't benefit")

# Plot results ------------------------------------------------
ggplot(data = sample_props_small, aes(x = p_hat)) +
  geom_histogram(binwidth = 0.05) +
  labs(
    x = "p_hat (Doesn't benefit)",
    title = "Sampling distribution of p_hat",
    subtitle = "Sample size = 50, Number of samples = 15000"
  )

sample_props_small
## # A tibble: 18 × 4
## # Groups:   replicate [18]
##    replicate scientist_work      n p_hat
##        <int> <chr>           <int> <dbl>
##  1         2 Doesn't benefit     3   0.3
##  2         4 Doesn't benefit     2   0.2
##  3         5 Doesn't benefit     1   0.1
##  4         6 Doesn't benefit     1   0.1
##  5         8 Doesn't benefit     1   0.1
##  6        10 Doesn't benefit     1   0.1
##  7        11 Doesn't benefit     1   0.1
##  8        13 Doesn't benefit     2   0.2
##  9        14 Doesn't benefit     3   0.3
## 10        15 Doesn't benefit     4   0.4
## 11        17 Doesn't benefit     4   0.4
## 12        18 Doesn't benefit     3   0.3
## 13        20 Doesn't benefit     3   0.3
## 14        21 Doesn't benefit     1   0.1
## 15        22 Doesn't benefit     3   0.3
## 16        23 Doesn't benefit     5   0.5
## 17        24 Doesn't benefit     1   0.1
## 18        25 Doesn't benefit     1   0.1


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, standard 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.)

Each observation represents a sampling from the total population of the specified size. With 5,000 simulations the mean with sample size 50 and 100 were consistently 20%. The mean with sample size 10 was consistently 22%, even with 50,000 simulations. This is probably due to by using sample size 10 the results were forced into discrete steps of 10%. The standard error changed from .11 to .06 to .04 with the change in sample size from 10 to 50 to 100, and was independent of number of simulations. Going from 5,000 to 50,000 did not change any of these values.


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 enhances 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?

My best point estimate of the population proportion of people who think the work scientists do enhances their lives is 60%.

# Create a second sample
set.seed(7)
samp7 <- global_monitor %>%
  sample_n(15)

# Generate statistics for sample 2
samp7_p_hat <- samp7 %>% 
  count(scientist_work) %>% 
  mutate(p_hat = n /sum(n)) %>% 
  filter(scientist_work == "Benefits") %>% 
  pull(p_hat) %>% 
  round(2)
samp7_p_hat
## [1] 0.6


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.

The sampling distribution shape looks almost normal, but discretized where it can only take results in increments of 1/15th or 6.67%, and it’s skewed to the left with relatively large spread. Since it skews left, while the mode is at 80% (12/15ths), the mean is probably lower. Based on this sampling distribution I would guess the true proportion of those who think the science benefits them would be 80% but I would want to rerun and increase the sample size so the results are less discretized and the mean and mode get closer. If I want my proportion to be accurate to 1%, I can’t be measuring in increments of 6.67%! The population proportion among all of the samples is 79.36%.

# Take a sample distribution of the data
set.seed(8378)
samplesize = 15
numberofreps = 2000
sample_props15 <- global_monitor %>%
                    rep_sample_n(samplesize, numberofreps, replace = TRUE) %>%
                    count(scientist_work) %>%
                    mutate(p_hat = n /sum(n)) %>%
                    filter(scientist_work == "Benefits")

# Plot results ------------------------------------------------
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 = 50, Number of samples = 15000"
  )

# Take the population proportion among all of the samples
sum(sample_props15$n)/(samplesize*numberofreps)
## [1] 0.7936


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?

Now the sampling distribution looks like a normal distribution and doesn’t appear discretized or skewed, the spread is smaller and the mean and the mode are virtually one and the same. Now I would guess the true proportion of those who think the work scientists do enhances their lives to be 80%.

# Take a sample distribution of the data
set.seed(9753)
samplesize = 150
numberofreps = 2000
sample_props150 <- global_monitor %>%
                    rep_sample_n(samplesize, numberofreps, replace = TRUE) %>%
                    count(scientist_work) %>%
                    mutate(p_hat = n /sum(n)) %>%
                    filter(scientist_work == "Benefits")

# Plot results ------------------------------------------------
ggplot(data = sample_props150, aes(x = p_hat)) +
  geom_histogram(binwidth = 0.02) +
  labs(
    x = "p_hat (Benefits)",
    title = "Sampling distribution of p_hat",
    subtitle = "Sample size = 50, Number of samples = 15000"
  )


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?

The distribution from exercise 9 with a sample size of 150 has a smaller spread than the distribution from exercise 8 with a sample size of 15. When we’re concerned with making estimates that are more often close to the true value, we prefer a sampling distribution with a small spread.