population <- tibble(id = 1:1e5) |> # initialize with 100K rows
mutate(parent = if_else(id <= 70000, "Yes", "No")) # first 70K "yes"2 Surveys and samples
2.1 Populations and samples
Studying populations is nice:
- all countries in the world
- all states in the US
- all cities in a state
However, we cannot study (say) all adults in a country. So we usually work with samples. This raises the issue of using samples to make inferences about populations.
2.1.1 Simple random sampling
Sampling where every eligible case has an equal probability of selection.
Note
In real-life surveys, simple random sampling is pretty uncommon. But it’s important as a baseline! The GSS, for example, uses a more complicated multi-stage design because it would be far too expensive to send interviewers to randomly scattered addresses all over the country. We’ll treat samples as simple random samples throughout this book anyway. The basic ideas are the same, but the formulas need some adjusting.
2.2 Simulations
We can use simulations to build intuition about sampling.
A simulation is when we make up “true data”, hide it from ourselves, and see how well we can figure out the truth using some procedure.
2.2.1 A simple survey
Imagine a city of 100,000 adults. Of these, 70,000 (i.e., 70%) have at least one child.
How close could we get to this number by drawing different random samples?
Let’s set up the “true” population:
This simple code makes the first 70,000 rows “yes” and the next 30,000 “no.” We now know the “truth”, which we can use for comparison.
Show code
plt(~ parent, data = population, type = "bar", xlab = "", ylab = "Count")
Show code
ggplot(population, aes(x = parent)) +
geom_bar() +
labs(x = "", y = "Count")
2.2.2 Data types
Because of the way we created it, parent will be a character <chr> variable. We often use 1 to mean “yes” and 0 to mean “no” in statistics. We could add a numeric version of parent as follows.
population <- population |>
mutate(parent_num = if_else(parent == "Yes", 1L, 0L))
Tip
Assigning an object to its “old” name allows you to add things to the original object. In this case, we are adding a new column, parent_num.
We can use glimpse() to easily see what type of variable things are.
glimpse(population)Rows: 100,000
Columns: 3
$ id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, ~
$ parent <chr> "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Ye~
$ parent_num <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,~
Tip
glimpse() is very useful. It shows you your data “sideways” so you can see information about all the columns (which are shown as rows).
Data type affects what we can do to a variable (or column). For example, we can take the mean (average) of a set of numbers but we can’t take the mean of a set of characters.
population |>
summarize(mean1 = mean(parent),
mean2 = mean(parent_num))Warning: There was 1 warning in `summarize()`.
i In argument: `mean1 = mean(parent)`.
Caused by warning in `mean.default()`:
! argument is not numeric or logical: returning NA# A tibble: 1 x 2
mean1 mean2
<dbl> <dbl>
1 NA 0.7
Note
In R, NA stands for “not available” and means that the data are missing. This cell is missing because what we asked for couldn’t be calculated.
NoteAside: doing stuff to a variable
We can pick out a column to operate on in two ways.
Tidyverse:
population |>
pull(parent_num) |>
mean()[1] 0.7Base R:
mean(population$parent_num)[1] 0.7Now back to our story…
2.3 Parameters and statistics
In our city, 70% is the population parameter because exactly 70,000 out of 100,000 people actually have at least one child.
We can’t afford to ask everyone, though. So what if we asked, say, 1000 randomly selected adults. Then we could compute the proportion of the sample that has a child. This would be a sample statistic. We use sample statistics to make inferences about population parameters.
This is the Roman vs. Greek distinction from Chapter 1. The population mean \(\mu\) is a parameter; the sample mean \(\bar{x}\) is a statistic. For proportions, the population proportion is \(p\) and the sample proportion is \(\hat{p}\), or “p-hat,” where the hat means “estimated from data.”
2.3.1 Drawing a sample
set.seed(722)
my_sample <- population |>
slice_sample(n = 1000,
replace = FALSE)Sampling is a random process. We will get a different result every time. By using set.seed(), we ensure we get the same “random” result every time the code is run.
Note
Sampling theory is based on sampling with replacement. However, to make it more straightforward, we will use sampling without replacement here.
2.3.2 The sample statistic
To estimate the population proportion, we will use the sample proportion.
my_sample |>
group_by(parent) |>
summarize(n = n())# A tibble: 2 x 2
parent n
<chr> <int>
1 No 318
2 Yes 682In our sample, 682 people are parents. This is 68.2%, which isn’t exactly the 70% in the population. This is because we randomly sampled from our population. It could be higher or lower.
2.4 Repeating the experiment
In real life, we only get to sample once. Sampling is expensive! But since this is just a simulation, we can ask what would happen if we sampled 1000 people many, many times.
2.4.1 A custom function
We can first make a function that does what we want once. This is hard at first but usually pays off.
get_count <- function(n = 1000) { # default n = 1000
slice_sample(population, n = n) |> # take a sample
summarize(sum = sum(parent_num)) |> # count the parents
as.integer() # save the number
}
Tip
If you run all the code up to here, you can call get_count() interactively in the console many times to get a feel for it.
2.4.2 Iterating
This would seem to make sense, but it doesn’t work.
set.seed(722)
my_bad_samples <- tibble(
sample_id = 1:100,
samp_count = get_count(n = 1000))
head(my_bad_samples)# A tibble: 6 x 2
sample_id samp_count
<int> <int>
1 1 682
2 2 682
3 3 682
4 4 682
5 5 682
6 6 682Every row is the same! That’s because tibble() (like mutate()) calls get_count() just once and then copies the answer all the way down the column. We want the function to run once for every row.
2.4.3 Iterating with rowwise()
set.seed(722)
my_samples <- tibble(
sample_id = 1:100) |>
rowwise() |>
mutate(samp_count = get_count(n = 1000))
head(my_samples, n = 3)# A tibble: 3 x 2
# Rowwise:
sample_id samp_count
<int> <int>
1 1 682
2 2 716
3 3 694
Tip
I don’t need n = 1000 because I set it as the default when I made my function.
2.4.4 Plotting the results
set.seed(722)
my_many_samples <- tibble(
sample_id = 1:2500) |>
rowwise() |>
mutate(samp_count = get_count())Show code
plt(~ samp_count, data = my_many_samples, type = "hist",
freq = FALSE, breaks = seq(600, 800, 5),
xlab = "Count of parents in the sample", ylab = "Density")
plt_add(~ samp_count, data = my_many_samples, type = "density", lwd = 2)
Show code
ggplot(my_many_samples, aes(x = samp_count)) +
geom_histogram(aes(y = after_stat(density)), # for overlay
boundary = 697.5, # why would I choose this?
binwidth = 5, # somewhat arbitrary
color = "white",
fill = "gray") +
scale_x_continuous(breaks = seq(600, 800, 10)) +
geom_density(linewidth = 1) # overlay density
WarningSample size and number of samples
When we are doing simulations like this, it can be easy to confuse the sample size (1000) with the number of samples in our simulation (2500). They are not the same thing!
The sample size is the number of people we would survey “in the real world”.
The number of samples is how many times we want to run our simulated experiment.
2.4.5 How accurate are we?
We will do this formally later. But now we can quantify how accurately a sample proportion of 1000 people might estimate this population proportion by using the interquartile range. This is how wide the middle half of the data is.
my_many_samples |> pull(samp_count) |> quantile(c(.25, .75))25% 75%
691 710 my_many_samples |> pull(samp_count) |> IQR()[1] 19
TipReminder
Remember: in real life we only get one of these samples.
Show code
plt(~ samp_count, data = my_many_samples, type = "density", lwd = 2,
xlab = "Count of parents in the sample", ylab = "Density")
abline(v = quantile(my_many_samples$samp_count, c(.25, .75)), lty = 3, lwd = 2)
Show code
ggplot(my_many_samples, aes(x = samp_count)) +
geom_density(fill = "gray", color = NA) +
scale_x_continuous(breaks = seq(600, 800, 10)) +
geom_vline(xintercept = quantile(my_many_samples$samp_count, c(.25, .75)),
linetype = "dotted", linewidth = 1)
2.5 Sample size
Remember that we drew a sample of 1000 people to estimate our sample proportions. What if we had different sample sizes? Let’s compare the following:
- \(n\) = 60
- \(n\) = 250
- \(n\) = 1000
Data prep
set.seed(722)
my_n60_samples <- tibble(
sample_id = 1:2500,
sample_size = "n = 60") |>
rowwise() |>
mutate(samp_count = get_count(n = 60),
samp_prop = samp_count / 60) # proportion
my_n250_samples <- tibble(
sample_id = 1:2500,
sample_size = "n = 250") |>
rowwise() |>
mutate(samp_count = get_count(n = 250),
samp_prop = samp_count / 250)
my_many_samples <- my_many_samples |> # adding the group var and prop
mutate(sample_size = "n = 1000",
samp_prop = samp_count / 1000)
samp_size_compare <-
bind_rows(my_n60_samples,
my_n250_samples,
my_many_samples)Show code
plt(~ samp_prop | sample_size, data = samp_size_compare, type = "density",
lwd = 2, xlab = "Sample proportion", ylab = "Density")
Show code
ggplot(samp_size_compare,
aes(x = samp_prop, group = sample_size, color = sample_size)) +
geom_density() +
labs(x = "Sample proportion", y = "Density", color = "")
TipProportions
The x-axis is now proportion because we can no longer compare raw counts.
The interquartile ranges (i.e., widths of the middle half of the data) decrease a lot with sample size.
samp_size_compare |>
group_by(sample_size) |>
summarize(IQR = IQR(samp_prop))# A tibble: 3 x 2
sample_size IQR
<chr> <dbl>
1 n = 1000 0.0190
2 n = 250 0.0360
3 n = 60 0.0667
Warning
We will explore these issues more formally very soon using the concepts sampling distribution and standard error. For now, the goal is to understand how to use simulations to build qualitative intuition about sample size.
You might notice that each sample size is about four times the one before it, and each IQR is about half the one before it. Precision improves with the square root of the sample size, so you need four times as many people to cut the spread in half. (This is one reason most national surveys have a couple thousand respondents rather than tens of thousands.)
NoteBias is not the same as noise
All three distributions above are centered on 0.7. A bigger sample doesn’t make an estimate less biased. Random samples of any size are unbiased; bigger ones are just less noisy. But if your sample is biased (say, because some kinds of people are less likely to answer), making it bigger won’t help.
2.6 Thinking with real data: GSS
The General Social Survey is a repeated cross-sectional survey that has been fielded every year or other year since 1972. It is the “Hubble Telescope” of sociology!
We’ll use the 2024 file you saved in the setup chapter.
gss2024 <- readRDS(here::here("data", "gss2024.rds")) |>
haven::zap_labels()2.6.1 Introducing abany
The GSS abany item asks “Please tell me whether or not you think it should be possible for a pregnant woman to obtain a legal abortion if the woman wants it for any reason?” The answers are “yes” (1) and “no” (2).
Warning
We do not want variables coded 1 and 2. As we will see later, it’s better if (almost) all variables have a meaningful 0 value.
d <- gss2024 |>
select(abany)
d |> group_by(abany) |>
summarize(n = n())# A tibble: 3 x 2
abany n
<dbl> <int>
1 1 1281
2 2 859
3 NA 1169d <- d |>
drop_na() |> # drop NA values
mutate(abany = if_else(abany == 1, 1, 0))
d |> pull(abany) |> table()
0 1
859 1281 2.6.2 abany sample proportion
We can use mean() to calculate the sample proportion.
mean(d$abany)[1] 0.5985981We find that 59.9% of our sample supports abortion rights for any reason.
2.6.3 Inference
How close is this sample statistic to the population parameter? We’ll never know.
We can use a simulation to give us a sense of what kind of accuracy is possible with a sample of 2140 respondents.
We could build an imaginary US adult population where 59.9% of adults support abortion rights. But we can instead use random number functions to draw a sample from an infinite population instead.
set.seed(722)
rbinom(n = nrow(d), size = 1, prob = mean(d$abany)) |> # sample from inf. pop.
mean() # take the mean[1] 0.5682243
TipWhy that value?
Why assume that the population parameter is the sample proportion? Because we are interested in how widely spread the simulations are and there isn’t a more reasonable value to choose.
Let’s do this 5000 times and collect the results.
set.seed(722)
gss_sims <- tibble(
sim_id = 1:5000) |>
rowwise() |>
mutate(samp_prop = mean(rbinom(nrow(d), 1, mean(d$abany))))The IQR tells us how spread out the middle half of the estimates are.
gss_sims |> pull(samp_prop) |> IQR()[1] 0.01495327Thus, with 2140 cases, half of the sample proportions will be within approximately 0.7 points of the true value.
Show code
plt(~ samp_prop, data = gss_sims, type = "density", lwd = 2,
xlab = "Simulated sample proportion", ylab = "Density")
abline(v = quantile(gss_sims$samp_prop, c(.025, .25, .75, .975)),
lty = 3, lwd = 2)
Show code
ggplot(gss_sims, aes(x = samp_prop)) +
geom_density(fill = "gray", color = NA) +
geom_vline(xintercept = quantile(gss_sims$samp_prop, c(.025, .25, .75, .975)),
linetype = "dotted") +
labs(x = "Simulated sample proportion", y = "Density",
caption = "lines at 2.5th, 25th, 75th, and 97.5th percentiles")
We still don’t know the true value, of course. We are probably within a couple of points of the true value. But we could be 3 (or possibly more) points away.
WarningWhat this doesn’t tell you
This simulation assumes simple random sampling, that everyone we selected actually responded, and that they all answered honestly. None of these is exactly true in a real survey. Sampling variability is just the part of the problem we can calculate.
2.7 Recap
- we want to know about populations
- we end up having to use samples
- samples are random subsets of the population that are expensive to collect
- the larger the sample, the more accurately we can infer the population proportion
- we can use simulations to understand how this works