Scenario: Inference for Penguin Data using Randomization and Bootstrapping

We are interested in the Palmer Penguins data set, and we will use two inferential tools to study group differences. In this project, we will work with a randomization test to ask whether an observed difference between two groups is real or plausibly due to chance, and a bootstrap confidence interval to estimate how large that difference plausibly is.

We will do this twice. First, we will work through a complete example together: studying the difference in body mass between male and female Gentoo penguins. Then, you will carry out the same analysis yourself for a second comparison: bill depth in Adelie versus Chinstrap penguins.

Part 1: Body Mass in Male vs. Female Gentoo Penguins

Setting up the data

We begin by filtering the penguins data to Gentoo penguins only, keeping only rows with complete data for body mass and sex.

gentoo <- penguins %>%
  filter(species == "Gentoo", !is.na(body_mass_g), !is.na(sex))

nrow(gentoo)
## [1] 119
table(gentoo$sex)
## 
## female   male 
##     58     61

We can visualize the raw data to get a sense of whether a difference exists before doing any inference.

ggplot(gentoo, aes(x = sex, y = body_mass_g, fill = sex)) +
  geom_boxplot(alpha = 0.7) +
  labs(title = "Body mass of Gentoo penguins by sex",
       x = "Sex", y = "Body mass (g)") +
  theme_minimal() +
  theme(legend.position = "none")

There appears to be a visible difference. But could this have arisen by chance? We use a randomization test to find out.

The Observed Difference

Our test statistic is the difference in mean body mass: male mean minus female mean.

gentoo %>%
  group_by(sex) %>%
  summarise(mean_body_mass = mean(body_mass_g))
## # A tibble: 2 × 2
##   sex    mean_body_mass
##   <fct>           <dbl>
## 1 female          4680.
## 2 male            5485.
observed_diff <- mean((gentoo %>% filter(sex == "male"))$body_mass_g) -
                 mean((gentoo %>% filter(sex == "female"))$body_mass_g)

observed_diff
## [1] 805.0947

The observed difference in mean body mass is approximately 805 grams, with males heavier on average.

The Randomization Test

To ask whether this difference is real, we simulate a world in which sex has no relationship to body mass. We do this by shuffling the sex labels and recomputing the difference in means, repeating the process 1,000 times to build a null distribution.

one_shuffle_diff <- function() {
  gentoo %>%
    mutate(sex_shuffle = sample(sex)) %>%
    group_by(sex_shuffle) %>%
    summarise(mean_mass = mean(body_mass_g)) %>%
    summarise(diff = mean_mass[sex_shuffle == "male"] -
                mean_mass[sex_shuffle == "female"]) %>%
    pull(diff)
}

null_distribution <- replicate(1000, one_shuffle_diff())

We can visualize the null distribution and mark where our observed difference falls.

null_df <- data.frame(sim_diff = null_distribution)

ggplot(null_df, aes(x = sim_diff)) +
  geom_histogram(bins = 40, fill = "steelblue", color = "white", alpha = 0.8) +
  geom_vline(xintercept = observed_diff, color = "firebrick", linewidth = 1.2,
             linetype = "dashed") +
  annotate("text", x = observed_diff - 200, y = 60,
           label = paste0("Observed\nDifference\n", round(observed_diff, 0), "g"),
           color = "firebrick", hjust = 0, size = 3.5) +
  labs(title = "Null distribution: 1,000 shuffles of sex labels",
       x = "Simulated difference in mean body mass (g)",
       y = "Count") +
  theme_minimal()

The observed difference falls far out in the tail of the null distribution — none of the 1,000 shuffles came close to producing a difference as large as what we observed.

We can quantify this with an informal p-value: the proportion of simulated differences that were as large or larger than what we observed.

p_value <- sum(null_distribution >= observed_diff) / length(null_distribution)
p_value
## [1] 0

The p-value is 0. Out of 1,000 random shuffles, zero produced a difference as large as the observed difference of 805 grams. This is strong evidence that the difference in body mass between male and female Gentoo penguins is not due to chance.

The Bootstrap Confidence Interval

The randomization test tells us the difference is real. Now we ask: how large is it? We use bootstrapping to estimate a plausible range for the true difference. Unlike the randomization test, bootstrapping resamples from the real data with replacement.

gentoo_males   <- gentoo %>% filter(sex == "male")
gentoo_females <- gentoo %>% filter(sex == "female")

one_boot_diff <- function() {
  boot_males   <- sample(gentoo_males$body_mass_g,   size = nrow(gentoo_males),   replace = TRUE)
  boot_females <- sample(gentoo_females$body_mass_g, size = nrow(gentoo_females), replace = TRUE)
  mean(boot_males) - mean(boot_females)
}

boot_diffs <- replicate(1000, one_boot_diff())

Once we have our bootstrapped distribution, we can look at it, and calculate the 95% confidence interval by chopping off the most extreme values (2.5% on each side).

data.frame(boot_diff = boot_diffs) %>%
  ggplot(aes(x = boot_diff)) +
  geom_histogram(bins = 40, fill = "steelblue", color = "white") +
  geom_vline(xintercept = observed_diff,    color = "firebrick",   linewidth = 1.2) + 
  labs(title = "Bootstrap distribution: difference in mean body mass",
       subtitle = "Red = observed difference",
       x = "Bootstrap difference in mean body mass (male - female)",
       y = "Count") +
  theme_minimal()

ci_diff <- quantile(boot_diffs, probs = c(0.025, 0.975))
ci_diff
##     2.5%    97.5% 
## 707.6578 908.9335

The 95% confidence interval runs from approximately 710g to 910g (although exact values will change depending on the randomness in our specific simulation). We are 95% confident that male Gentoo penguins are, on average, somewhere between 710 and 910 grams heavier than female Gentoo penguins. Notice that this interval does not include zero, which is consistent with our randomization test conclusion that the difference is real.

Part 2: Do a Similar Analysis for another Penguin Species!

Now it is your turn. You will carry out the same two-part analysis for a different comparison. Pick a different species (Adelie or Chinstrap) and a quantitative variable (body mass is fine to do, but you may find more interesting results with bill depth, bill length, or flipper length) and address the following questions:

adelie <- penguins %>%
  filter(species == "Adelie", !is.na(body_mass_g), !is.na(sex))

nrow(adelie)
## [1] 146
table(adelie$sex)
## 
## female   male 
##     73     73
ggplot(adelie, aes(x = sex, y = body_mass_g, fill = sex)) +
  geom_boxplot(alpha = 0.7) +
  labs(title = "Body mass of Adelie penguins by sex",
       x = "Sex", y = "Body mass (g)") +
  theme_minimal() +
  theme(legend.position = "none")

Question A: Consider the difference in mean between the male and female penguins of your chosen species. Is that difference real, or is it plausibly due to chance?

adelie %>%
  group_by(sex) %>%
  summarise(mean_body_mass = mean(body_mass_g))
## # A tibble: 2 × 2
##   sex    mean_body_mass
##   <fct>           <dbl>
## 1 female          3369.
## 2 male            4043.
observed_diff <- mean((adelie %>% filter(sex == "male"))$body_mass_g) -
                 mean((adelie %>% filter(sex == "female"))$body_mass_g)

observed_diff
## [1] 674.6575

The observed difference in mean body mass is approximately 675 grams, with males heavier on average.

The Randomization Test

To ask whether this difference is real, we simulate a world in which sex has no relationship to body mass. We do this by shuffling the sex labels and recomputing the difference in means, repeating the process 1,000 times to build a null distribution.

one_shuffle_diff <- function() {
  adelie %>%
    mutate(sex_shuffle = sample(sex)) %>%
    group_by(sex_shuffle) %>%
    summarise(mean_mass = mean(body_mass_g)) %>%
    summarise(diff = mean_mass[sex_shuffle == "male"] -
                mean_mass[sex_shuffle == "female"]) %>%
    pull(diff)
}

null_distribution <- replicate(1000, one_shuffle_diff())

We can visualize the null distribution and mark where our observed difference falls.

null_df <- data.frame(sim_diff = null_distribution)

ggplot(null_df, aes(x = sim_diff)) +
  geom_histogram(bins = 40, fill = "steelblue", color = "white", alpha = 0.8) +
  geom_vline(xintercept = observed_diff, color = "firebrick", linewidth = 1.2,
             linetype = "dashed") +
  annotate("text", x = observed_diff - 200, y = 60,
           label = paste0("Observed\nDifference\n", round(observed_diff, 0), "g"),
           color = "firebrick", hjust = 0, size = 3.5) +
  labs(title = "Null distribution: 1,000 shuffles of sex labels",
       x = "Simulated difference in mean body mass (g)",
       y = "Count") +
  theme_minimal()

Question B: What is a plausible range of values for the size of that difference?

To do this, you will want to set up and filter your data similar to what we’ve done already; shuffle the sexes of the penguins to create a randomization test; and create a bootstrapped confidence interval. Most (all?) of your code can be drawn from Script 2 and Script 3, with key values changed so that we’re looking at a different penguin species and, possibly, a different quantitative variable.

Include the following:

(1) An explicit calculation of the observed statistic. (This is a difference in mean between males and females of your species. For the Gentoo body mass example, this was roughly 805.)

adelie %>%
  group_by(sex) %>%
  summarise(mean_body_mass = mean(body_mass_g))
## # A tibble: 2 × 2
##   sex    mean_body_mass
##   <fct>           <dbl>
## 1 female          3369.
## 2 male            4043.
observed_diff <- mean((adelie %>% filter(sex == "male"))$body_mass_g) -
                 mean((adelie %>% filter(sex == "female"))$body_mass_g)

observed_diff
## [1] 674.6575

The observed difference in mean body mass is approximately 675 grams, with males heavier on average.

(2) A histogram showing the difference in mean that is created using the simulated randomization test. (This comes from Script 2) ## The Randomization Test

To ask whether this difference is real, we simulate a world in which sex has no relationship to body mass. We do this by shuffling the sex labels and recomputing the difference in means, repeating the process 1,000 times to build a null distribution.

one_shuffle_diff <- function() {
  adelie %>%
    mutate(sex_shuffle = sample(sex)) %>%
    group_by(sex_shuffle) %>%
    summarise(mean_mass = mean(body_mass_g)) %>%
    summarise(diff = mean_mass[sex_shuffle == "male"] -
                mean_mass[sex_shuffle == "female"]) %>%
    pull(diff)
}

null_distribution <- replicate(1000, one_shuffle_diff())

We can visualize the null distribution and mark where our observed difference falls.

null_df <- data.frame(sim_diff = null_distribution)

ggplot(null_df, aes(x = sim_diff)) +
  geom_histogram(bins = 40, fill = "steelblue", color = "white", alpha = 0.8) +
  geom_vline(xintercept = observed_diff, color = "firebrick", linewidth = 1.2,
             linetype = "dashed") +
  annotate("text", x = observed_diff - 200, y = 60,
           label = paste0("Observed\nDifference\n", round(observed_diff, 0), "g"),
           color = "firebrick", hjust = 0, size = 3.5) +
  labs(title = "Null distribution: 1,000 shuffles of sex labels",
       x = "Simulated difference in mean body mass (g)",
       y = "Count") +
  theme_minimal()

(3) A histogram showing the bootstrapped distribution via resampling. (This comes from Script 3)

adelie_males   <- adelie %>% filter(sex == "male")
adelie_females <- adelie %>% filter(sex == "female")

one_boot_diff <- function() {
  boot_males   <- sample(adelie_males$body_mass_g,   size = nrow(adelie_males),   replace = TRUE)
  boot_females <- sample(adelie_females$body_mass_g, size = nrow(adelie_females), replace = TRUE)
  mean(boot_males) - mean(boot_females)
}

boot_diffs <- replicate(1000, one_boot_diff())
data.frame(boot_diff = boot_diffs) %>%
  ggplot(aes(x = boot_diff)) +
  geom_histogram(bins = 40, fill = "steelblue", color = "white") +
  geom_vline(xintercept = observed_diff,    color = "firebrick",   linewidth = 1.2) + 
  labs(title = "Bootstrap distribution: difference in mean body mass",
       subtitle = "Red = observed difference",
       x = "Bootstrap difference in mean body mass (male - female)",
       y = "Count") +
  theme_minimal()

ci_diff <- quantile(boot_diffs, probs = c(0.025, 0.975))
ci_diff
##     2.5%    97.5% 
## 578.0822 776.0274

(4) Appropriate analysis and conclusion. (How do these two graphs compare to each other? What can you conclude about the difference between male and female penguins of your species – is the observed difference likely due to an underlying difference, or is it possible that it’s due to random chance?)

The difference between the Null distribution and the Bootstrap distribution is that the Bootstrap histogram has a much wider range while the Null is much more condensed.