knitr::opts_chunk$set(eval = TRUE, message = FALSE, warning = FALSE)
set.seed(1234)Foundations for statistical inference - Sampling distributions
In this lab, you will investigate the ways in which the statistics from a random sample of data can serve as point estimates for population parameters. We’re interested in formulating a sampling distribution of our estimate in order to learn about the properties of the estimate, such as its distribution.
Setting a seed: We will take some random samples and build sampling distributions in this lab, which means you should set a seed at the start of your lab. If this concept is new to you, review the lab on probability.
Getting Started
Load packages
In this lab, we will explore and visualize the data using the tidyverse suite of packages. We will also use the infer package for resampling.
Let’s load the packages.
library(tidyverse)
library(openintro)
library(infer)
library(shiny)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.
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))
)The name of the data frame is global_monitor and the name of the variable that contains responses to the question “Do you believe that the work scientists do benefit people like you?” is scientist_work.
We can quickly visualize the distribution of these responses using a bar plot.
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() 
We can also obtain summary statistics to confirm we constructed the data frame correctly.
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.2The unknown sampling distribution
In this lab, you have access to the entire population, but this is rarely the case in real life. Gathering information on an entire population is often extremely costly or impossible. Because of this, we often take a sample of the population and use that to understand the properties of the population.
If you are interested in estimating the proportion of people who don’t think the work scientists do benefits them, you can use the sample_n command to survey the population.
samp1 <- global_monitor %>%
sample_n(50)This command collects a simple random sample of size 50 from the global_monitor dataset, and assigns the result to samp1. This is similar to randomly drawing names from a hat that contains the names of all in the population. Working with these 50 names is considerably simpler than working with all 100,000 people in the population.
- Describe the distribution of responses in this sample. How does it compare to the distribution of responses in the population. Hint: Although the
sample_nfunction 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 summarizing the population data will still be useful for the sample, however be careful to not label your proportionpsince 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.
The distribution of responses in this sample is similar but not exactly the same as the distribution of responses in the population. This is expected based on our understanding of the distribution from a sample of large enough size (30 or above) from the population works. The general ‘shape’ of the distribution of responses is the same, which makes sense since the sample is literally from the population: There are more people who say science ‘Benefits’ than people who say it ‘Doesn’t benefit’ in both the population and the sample. However, the exact numbers aren’t exactly the same in terms of the sample proportion and the population proportion. The sample resembles approximately the population: in the sample, 74% of respondents represent scientist_work == 'Benefits' while in the population, 80% of respondents represent scientist_work == 'Benefits'. In the sample, 26% of respondents represent scientist_work == 'Doesn't benefit' compared to 20% in the population.
samp1 |>
count(scientist_work) |>
mutate(proportion_samp1 = n / sum(n) )# A tibble: 2 × 3
scientist_work n proportion_samp1
<chr> <int> <dbl>
1 Benefits 37 0.74
2 Doesn't benefit 13 0.26If you’re interested in estimating the proportion of all people who do not believe that the work scientists do benefits them, but you do not have access to the population data, your best single guess is the sample mean.
samp1 %>%
count(scientist_work) %>%
mutate(p_hat = n /sum(n))# A tibble: 2 × 3
scientist_work n p_hat
<chr> <int> <dbl>
1 Benefits 37 0.74
2 Doesn't benefit 13 0.26Depending on which 50 people you selected, your estimate could be a bit above or a bit below the true population proportion of 0.26. In general, though, the sample proportion turns out to be a pretty good estimate of the true population proportion, and you were able to get it by sampling less than 1% of the population.
- 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 - it’d be unlikely to get exactly the same sample proportion in any two samples, but it would certainly be possible. After all, we’re sampling from the same population – so it’s possible, just not likely – especially with a population so large and a sample size representing less than 1% of the population. However, we can expect any sample proportion will approximate the population proportion – as noted above, the sample proportion is a pretty good estimate of the true population proportion, so we could expect a similar approximation of the population in any sample of the sample population.
- Take a second sample, also of size 50, and call it
samp2. How does the sample proportion ofsamp2compare with that ofsamp1? 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)samp2 |>
count(scientist_work) |>
mutate(proportion_samp2 = n / sum(n) )# A tibble: 2 × 3
scientist_work n proportion_samp2
<chr> <int> <dbl>
1 Benefits 42 0.84
2 Doesn't benefit 8 0.16The samp2 is different from samp1 but this samp1 proportion still represents an approximation of the population proportion, as expected. While the sample proportions in samp1 and samp2 are different, they are both accurate approximations of the population proportion. samp2 represents a different approximation, but both are still ‘approximations’ of the population: samp2 represents 84% ‘Benefits’ compared to samp1 74% ‘Benefits and 80% ’Benefits’ for the population; samp2 represents 16% ‘Doesn’t benefit’ compared to samp1 26% ‘Doesn’t benefit’ and 20% ‘Doesn’t benefit’ for the population. It is not at all surprising, and in fact, expected to see this approximation of the population proportion.
If we took two more samples, one of size 100 and one of size 1000, we could also expect to see approximations of the population proportion. We can expect that as the sample size n increases, so does the level of approximating the population proportion. So, the larger the sample size, the more likely that the sample proportion is closer / more accurately representative of the population proportion.
Not surprisingly, every time you take another random sample, you might get a different sample proportion. It’s useful to get a sense of just how much variability you should expect when estimating the population mean this way. The distribution of sample proportions, called the sampling distribution (of the proportion), can help you understand this variability. In this lab, because you have access to the population, you can build up the sampling distribution for the sample proportion by repeating the above steps many times. Here, we use R to take 15,000 different samples of size 50 from the population, calculate the proportion of responses in each sample, filter for only the Doesn’t benefit responses, and store each result in a vector called sample_props50. Note that we specify that replace = TRUE since sampling distributions are constructed by sampling with replacement.
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")And we can visualize the distribution of these proportions 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"
)Next, you will review how this set of code works.
- 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.
sample_props50 # A tibble: 15,000 × 4
# Groups: replicate [15,000]
replicate scientist_work n p_hat
<int> <chr> <int> <dbl>
1 1 Doesn't benefit 11 0.22
2 2 Doesn't benefit 8 0.16
3 3 Doesn't benefit 8 0.16
4 4 Doesn't benefit 6 0.12
5 5 Doesn't benefit 10 0.2
6 6 Doesn't benefit 11 0.22
7 7 Doesn't benefit 10 0.2
8 8 Doesn't benefit 14 0.28
9 9 Doesn't benefit 11 0.22
10 10 Doesn't benefit 11 0.22
# ℹ 14,990 more rowsThere are 15000 total ‘elements’ or observations that represent the proportion of “Doesn’t benefit” for 15000 samples of size n = 50. In sample_props50 each row represents the proportion of ‘Doesn’t benefit’ in each sample of n = 50 taken (so, each sample had 50 elements). So we have 15000 sample proportions of “Doesn’t benefit” from the 15000 samples taken from the population. There were 15000 samples of size n = 50. What we see is a sampling distribution – a distribution of sample proportions (p-hat of ‘Doesn’t benefit’).
The sample distribution, as expected based on the Central Limit Theorem, is approximately normal, meaning it’s approximately symmetrical about the mean, and the mean of the sampling distribution is an accurate approximation of the population mean. We can see that, like in the population proportion of 0.2 (20%) for “Doesn’t benefit”, so too the sampling distribution represents a center of 0.2 – which is expected based on the Central Limit Theorem, which includes that the mean of the sample means is equal to the population mean.
ggplot(data = sample_props50, aes(x = p_hat)) +
geom_histogram(binwidth = 0.02, fill = "lightgrey", color = "darkgrey") +
labs(
x = "p_hat (Doesn't benefit)",
title = "Sampling distribution of p_hat 'Doesn't benefit'",
subtitle = "Sample size = 50, Number of samples = 15000"
)
We can see similar sampling distribution to represent the population proportion if we look at the “Benefits” results. As expected, the center of 0.8 – the mean of sample means in this sampling distribution – is equal to the population mean of 0.8. Nice!
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 == "Benefits")ggplot(data = sample_props50, aes(x = p_hat)) +
geom_histogram(binwidth = 0.02, fill = "lightgrey", color = "darkgrey") +
labs(
x = "p_hat (Benefits)",
title = "Sampling distribution of p_hat 'Benefits'",
subtitle = "Sample size = 50, Number of samples = 15000"
)
Interlude: Sampling distributions
The idea behind the rep_sample_n function is repetition. Earlier, you took a single sample of size n (50) from the population of all people in the population. With this new function, you can repeat this sampling procedure rep times in order to build a distribution of a series of sample statistics, which is called the sampling distribution.
Note that in practice one rarely gets to build true sampling distributions, because one rarely has access to data from the entire population.
Without the rep_sample_n function, this would be painful. We would have to manually run the following code 15,000 times
global_monitor %>%
sample_n(size = 50, replace = TRUE) %>%
count(scientist_work) %>%
mutate(p_hat = n /sum(n)) %>%
filter(scientist_work == "Doesn't benefit")# A tibble: 1 × 3
scientist_work n p_hat
<chr> <int> <dbl>
1 Doesn't benefit 8 0.16as well as store the resulting sample proportions each time in a separate vector.
Note that for each of the 15,000 times we computed a proportion, we did so from a different sample!
- To make sure you understand how sampling distributions are built, and exactly what the
rep_sample_nfunction 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 namedsample_props_small. Print the output. How many observations are there in this object calledsample_props_small? What does each observation represent?
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")print(sample_props_small)# A tibble: 22 × 4
# Groups: replicate [22]
replicate scientist_work n p_hat
<int> <chr> <int> <dbl>
1 1 Doesn't benefit 1 0.1
2 2 Doesn't benefit 2 0.2
3 3 Doesn't benefit 6 0.6
4 4 Doesn't benefit 1 0.1
5 5 Doesn't benefit 1 0.1
6 6 Doesn't benefit 3 0.3
7 8 Doesn't benefit 1 0.1
8 9 Doesn't benefit 3 0.3
9 10 Doesn't benefit 2 0.2
10 12 Doesn't benefit 1 0.1
# ℹ 12 more rowsInterestingly, and likely related to the small sample size of n = 10, it appears that 3 of the 25 samples included 0 as a proportion for “Doesn’t benefit”. There are 22 rows in sample_props_small which represent the sample proportions where scientist_work == "Doesn't benefit". This means that for 3 of the 25 samples, 100% (proportion of 1) represented “Benefit”. Because we filtered for “Doesn’t benefit”, we see a result of 22 rows rather than 25 – the proportion for 3 of the samples for “Doesn’t benefit” was 0 so they’re not included. The sample size of n = 10 is not sufficiently large (the recommendation is n = 30 or above).
Sample size and the sampling distribution
Mechanics aside, let’s return to the reason we used the rep_sample_n function: to compute a sampling distribution, specifically, the sampling distribution of the proportions from samples of 50 people.
ggplot(data = sample_props50, aes(x = p_hat)) +
geom_histogram(binwidth = 0.02)The sampling distribution that you computed tells you much about estimating the true proportion of people who think that the work scientists do doesn’t benefit them. Because the sample proportion is an unbiased estimator, the sampling distribution is centered at the true population proportion, and the spread of the distribution indicates how much variability is incurred by sampling only 50 people at a time from the population.
In the remainder of this section, you will work on getting a sense of the effect that sample size has on your sampling distribution.
- 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.)
For a sample size of 10 for 5000 simulations: the mean is 0.22, the standard error is 0.11, and the shape is right skewed.
For a sample size of 50 for 5000 simulations: the mean is 0.2, the standard error is 0.06, and the shape is approximately normal.
For a sample size of 100 for 5000 simulations: the mean is 0.2, the standard error is 0.04, and the shape is approximately normal.
Like with the sampling distributions we were looking at previously, each observation represents the proportion of “Doesn’t benefit” in each sample size of n=10, n=50, and n=100, each with 5000 simulations (so 5000 samples of n=10, 5000 samples of n=50, 5000 samples of n=100). And in each case we’re looking at the sampling distribution of sample means (sample proportions) of “Doesn’t benefit”. We can see that as the sample size increases, the shape becomes more normal and the standard error decreases.
More Practice
So far, you have only focused on estimating the proportion of those you think the work scientists doesn’t benefit them. Now, you’ll try to estimate the proportion of those who think it does.
Note that while you might be able to answer some of these questions using the app, you are expected to write the required code and produce the necessary plots and summary statistics. You are welcome to use the app for exploration.
- 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 enhances their lives?
samp1_benefits <- global_monitor |>
sample_n(15)samp1_benefits |>
count(scientist_work) |>
mutate(p_hat = n /sum(n))# A tibble: 2 × 3
scientist_work n p_hat
<chr> <int> <dbl>
1 Benefits 10 0.667
2 Doesn't benefit 5 0.333samp1_benefits_p_hat <- samp1_benefits |>
count(scientist_work) |>
mutate(p_hat = n /sum(n)) |>
filter(scientist_work == "Benefits") |>
pull(p_hat) |>
round(2)
samp1_benefits_p_hat[1] 0.67Interestingly, we see the exact same sample proportion of 0.8 as the population proportion of 0.8 which represents the proportion of people in this sample (and in the population) who think that the work scientists do benefits or enhances their lives. The proportion in the sample size of n=15 is an estimation that exactly represents the population proportion of scientist_work == "Benefits".
- Since you have access to the population, simulate the sampling distribution of proportion of those who think the work scientists do encances 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 enhances their lives to be? Finally, calculate and report the population proportion.
sample_props15 <- global_monitor |>
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.08, fill = "lightgrey", color = "darkgrey") +
labs(
x = "p_hat (Benefits)",
title = "Sampling distribution of p_hat 'Benefits'",
subtitle = "Sample size = 15, Number of samples = 2000"
)
The shape of this sampling distribution is left skewed. Based on this sampling distribution, I would guess the true proportion of those who think the work scientists do enhances their lives to be close to .8 or 80%. This is also what we see in the population proportion.
population <- global_monitor |>
group_by(scientist_work) |>
summarize(popcount = n()) |>
mutate(population_proportion = popcount / sum(popcount))
population# A tibble: 2 × 3
scientist_work popcount population_proportion
<chr> <int> <dbl>
1 Benefits 80000 0.8
2 Doesn't benefit 20000 0.2- 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 enhances their lives?
sample_props150 <- global_monitor |>
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_props150, aes(x = p_hat)) +
geom_histogram(binwidth = 0.02, fill = "lightgrey", color = "darkgrey") +
labs(
x = "p_hat (Benefits)",
title = "Sampling distribution of p_hat 'Benefits'",
subtitle = "Sample size = 150, Number of samples = 2000"
)
Once again, based on this sampling distribution, I would guess around .8 to be the true proportion of those who think the work scientists do enhances their lives. This is an approximation of the population proportion as reported above. As expected, the shape of the distribution appeared normal about the mean, and there is closer approximation to the population proportion with a greater sample size.
- 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?
mean_15 <- round(mean(sample_props15$p_hat), 4)
se_15 <- round(sd(sample_props15$p_hat) / sqrt(nrow(sample_props15)), 4)
mean_150 <- round(mean(sample_props150$p_hat), 4)
se_150 <- round(sd(sample_props150$p_hat) / sqrt(nrow(sample_props150)), 4)
summary_15_150 <- data.frame(
sampling_distribution = c("sample_props15", "sample_props150"),
mean = c(mean_15, mean_150),
standard_error = c(se_15, se_150)
)
summary_15_150 sampling_distribution mean standard_error
1 sample_props15 0.8019 0.0023
2 sample_props150 0.8005 0.0007