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: any variable in male vs female Adelie or 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% 
## 699.8286 913.1289

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:

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?

To answer this, we begin by filtering the penguins data to Adelie penguins only, keeping only rows with complete data for body mass and sex.

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

nrow(adelie)

## [1] 146

table(adelie$sex)

## 
## female   male 
##     73     73

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

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")

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.

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_a <- mean((adelie %>% filter(sex == "male"))$body_mass_g) -
                 mean((adelie %>% filter(sex == "female"))$body_mass_g)

observed_diff_a

## [1] 674.6575

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

The Randomization Test

one_shuffle_diff_a <- 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_a <- replicate(1000, one_shuffle_diff_a())

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

null_df_a <- data.frame(sim_diff = null_distribution_a)

ggplot(null_df_a, aes(x = sim_diff)) +
  geom_histogram(bins = 40, fill = "steelblue", color = "white", alpha = 0.8) +
  geom_vline(xintercept = observed_diff_a, color = "firebrick", linewidth = 1.2,
             linetype = "dashed") +
  annotate("text", x = observed_diff_a - 200, y = 60,
           label = paste0("Observed\nDifference\n", round(observed_diff_a, 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_a <- sum(null_distribution_a >= observed_diff_a) / length(null_distribution_a)
p_value_a

## [1] 0

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

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

The Bootstrap Confidence Interval

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

one_boot_diff_a <- 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_a <- replicate(1000, one_boot_diff_a())

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_a) %>%
  ggplot(aes(x = boot_diff)) +
  geom_histogram(bins = 40, fill = "steelblue", color = "white") +
  geom_vline(xintercept = observed_diff_a,    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_a <- quantile(boot_diffs_a, probs = c(0.025, 0.975))
ci_diff_a

##     2.5%    97.5% 
## 577.3887 774.3151

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

Unit 7 Project

Samson Buras

03/06/2026

Scenario: Inference for Penguin Data using Randomization and Bootstrapping

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

Setting up the data

The Observed Difference

The Randomization Test

The Bootstrap Confidence Interval

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

The Observed Difference

The Randomization Test

The Bootstrap Confidence Interval