Analyze Antarctic Penguin Biodiversity

Use R and tidyverse to analyze real Antarctic penguin research data.

Introduction

30 Second Summary

Data can tell compelling stories about the natural world. But it can also mislead you if you don't know how to look at it from the right angle.

In this project, you will analyze real Antarctic penguin research data using R. You'll explore measurements from three penguin species collected at Palmer Station, Antarctica, clean missing values, compute summary statistics, and build visualizations that reveal hidden patterns.

What You'll Build

Run your finished script in RStudio and see a series of plots and statistics appear, telling the story of how three Antarctic penguin species differ in body size, bill shape, and island habitat.

By the end of this project, you'll have:

  • Scatter plots, boxplots, and histograms built with ggplot2 that show three distinct penguin species clusters by body size and bill shape.
  • Summary statistics from cleaned data revealing which species is the largest, which has the longest flippers, and which islands each species calls home.
  • A misleading trend line you built on purpose and corrected, revealing Simpson's Paradox where a pooled trend completely reverses inside each species.
  • Secret Mission: Investigate sexual dimorphism patterns and quantify how much larger males are than females in each species.

Are there any prerequisites?

This project assumes you already have R and RStudio Desktop installed on your computer. Step 1 handles installing all the packages you need.

Before We Start

Before we dive in, let's lock in what you're about to build and why it matters. This step is your chance to commit to the goal before any hands-on work begins.

Set Up Your R Packages

To analyze Antarctic penguin data, you need the right tools loaded into R. R's real power for data analysis comes from packages, which are collections of pre-built functions that anyone can download for free.

In this step, you'll install two packages that power this entire project. You'll also create the script file where all your analysis code will live.

In this step, get ready to:
  • Install the tidyverse and palmerpenguins packages from CRAN.
  • Create your penguin_analysis.R script and verify both packages load successfully.
Install the tidyverse and palmerpenguins packages

Your two packages are tidyverse (a collection of tools for data wrangling and visualization) and palmerpenguins (real Antarctic penguin research data from Palmer Station). Both are hosted on CRAN, R's central package repository.

  • Install both packages by running these commands one at a time in the RStudio Console pane (the panel at the bottom-left of RStudio with the > prompt):
install.packages("tidyverse")
install.packages("palmerpenguins")
  • Wait for the first command to finish before running the second. Each command prints download and installation messages as R fetches files from CRAN.
  • After each command completes, the > prompt returns in the Console.

What is CRAN?

CRAN is the Comprehensive R Archive Network. It hosts thousands of free R packages that anyone can install.

The install.packages() function downloads a package from CRAN to your computer. You only need to run it once per package.

Seeing an error during installation?

  • Make sure you're connected to the internet. install.packages() downloads files from CRAN's servers.
  • If R asks you to select a CRAN mirror, pick any option from the list or accept the default.
  • Tidyverse is a large collection of packages, so the installation may take a few minutes. Wait for the > prompt to return before running the next command.

Why is my R package installation failing?

Create your analysis script and load the packages

Now that the packages are installed on your computer, you need to load them into your current R session. You also need a script file to save your code so you can rerun it later.

  • Click File in the RStudio menu bar.
  • Click New File, then click R Script.

A blank editor tab opens above the Console. This is where you'll write all your analysis code.

  • Copy the following code into the empty script:
# penguin_analysis.R
# Analyze Antarctic Penguin Biodiversity with R
# Data: Palmer Station LTER, Antarctica

# Load libraries
library(tidyverse)
library(palmerpenguins)

What is the difference between install.packages() and library()?

install.packages() downloads a package from CRAN to your computer. You only run this once per package.

library() loads an already-installed package into your current R session. You need to call library() every time you start a new R session. Think of it like downloading an app once from the app store, then opening it each time you want to use it.

  • Press Cmd+S (macOS) or Ctrl+S (Windows) to save the file.
  • In the save dialog, navigate to your Desktop.
  • Name the file penguin_analysis.R and click Save.
  • Select all the code by pressing Cmd+A (macOS) or Ctrl+A (Windows).

Before you run this, what do you think will happen? Will both packages load silently, or will one of them print something to the Console?

  • Press Cmd+Enter (macOS) or Ctrl+Enter (Windows) to run the selected code.
  • Look at the Console pane below. Tidyverse prints a message listing its attached core packages, including ggplot2, dplyr, tidyr, readr, purrr, tibble, stringr, and forcats. The palmerpenguins package loads silently with no output.

Seeing a 'there is no package' error?

  • This means the install.packages() command from the previous substep did not complete successfully. Run the install command for the missing package again in the Console, then rerun your script.
  • If you see lines mentioning Conflicts in the tidyverse output, that is completely normal. It means some tidyverse functions share names with built-in R functions. Tidyverse's versions take priority automatically.

Why won't my R packages load with library()?

✔️ Awesome, I've got everything!

Make sure you've saved penguin_analysis.R and that both library() calls ran without errors in the Console.

ⓧ I'd like to double check the full code

Your penguin_analysis.R file should contain exactly this:

# penguin_analysis.R
# Analyze Antarctic Penguin Biodiversity with R
# Data: Palmer Station LTER, Antarctica

# Load libraries
library(tidyverse)
library(palmerpenguins)

Your R environment is ready and your script file is set up. Next up, you'll load the penguin dataset, inspect its structure, and build your first scatter plot.

Explore the Penguin Data and Build Your First Plot

Your R environment is ready with tidyverse and palmerpenguins loaded. Now it is time to meet the data you will be analyzing for the rest of this project.

You cannot analyze what you do not understand. Every data analysis starts the same way: look at the structure first. In this step, you will inspect the penguins dataset, understand what each column represents, and build your first visualization.

In this step, get ready to:
  • Inspect the penguins dataset using three functions that reveal its structure, summary statistics, and group sizes.
  • Build a scatter plot of flipper length vs body mass colored by species using ggplot2.
  • Spot a warning about missing values that you will need to clean up in the next step.
Inspect the penguins dataset

Three functions will tell you everything you need to know about a new dataset. You will use these three at the start of every future analysis.

  • In penguin_analysis.R, add the following code below the library() calls:
# ---- Explore the Data ----

# First look at the data structure
glimpse(penguins)

# Summary statistics for every variable
summary(penguins)

# Count penguins by species
penguins %>%
  count(species)

What does this code do?

Each function answers a different question about the dataset.

  • glimpse() shows column names, data types, and the first few values. This tells you what variables exist and whether they are numbers, text, or categories.
  • summary() gives the min, max, median, and mean for numeric columns. It also flags how many NA (missing) values each column has.
  • count() tallies how many rows belong to each group. The %>% pipe operator sends the data on the left into the function on the right, so you can chain operations together.

Running code in RStudio

To run a line of code from your script, place your cursor on that line and press Ctrl+Enter (Windows) or Cmd+Enter (macOS). To run multiple lines, select them first, then press the same shortcut. The output appears in the Console pane below.

  • Place your cursor on the glimpse(penguins) line and press Ctrl+Enter (Windows) or Cmd+Enter (macOS) to run it.

You should see 344 rows and 8 columns listed in the Console. The columns include species, island, bill measurements, flipper length, body mass, sex, and year.

  • Run summary(penguins) the same way.

The summary output shows ranges and averages for every numeric variable. Notice the NA's counts at the bottom of several columns. These are missing values you will deal with later.

  • Select the last two lines (penguins %>% and count(species)) and run them together.

The Console prints a tibble showing three species: Adelie (152), Chinstrap (68), and Gentoo (124).

Not seeing any output?

  • Make sure you are running the code from the script editor (top-left pane), not typing it directly into the Console.
  • Check that library(palmerpenguins) ran successfully in Step 1. If you see object 'penguins' not found, rerun the library() lines at the top of your script first.
  • For the count() pipeline, make sure you select both lines before pressing Ctrl+Enter. Running just the first line will produce an incomplete pipe error.

Still stuck? Help me debug why my penguins data is not loading

Build your first scatter plot

Every ggplot2 plot follows the same pattern. You start with ggplot() to set up a canvas, add a geom_ layer to draw marks on it, then customize with scales and labels. This pattern is the same for every plot you will ever build.

  • Add this scatter plot code at the bottom of penguin_analysis.R, below the count() pipeline:
# First scatter plot: flipper length vs body mass
ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g, color = species, shape = species)) +
  geom_point(size = 2) +
  scale_color_manual(values = c("darkorange", "darkorchid", "cyan4")) +
  labs(
    title = "Penguin Flipper Length vs Body Mass",
    x = "Flipper Length (mm)",
    y = "Body Mass (g)"
  ) +
  theme_minimal()

What does this code do?

  • ggplot(penguins, aes(...)) creates the canvas. The aes() function maps data columns to visual properties: x position, y position, color, and shape.
  • geom_point(size = 2) draws one dot per row in the dataset.
  • scale_color_manual() sets custom colors for each species instead of the defaults.
  • labs() adds a title and axis labels.
  • theme_minimal() applies a clean, minimal background style.

Before you run this, what do you think the plot will look like? Will the three species form distinct clusters, or will the points all overlap?

  • Select all the scatter plot lines (from the comment through theme_minimal()) and press Ctrl+Enter (Windows) or Cmd+Enter (macOS) to run them.

A scatter plot appears in the Plots pane (bottom-right) showing three distinct colored clusters. Gentoo penguins (cyan) are clearly the largest, clustered in the top-right. Adelie (orange) and Chinstrap (purple) are smaller but separated from each other.

Plot not appearing?

  • Make sure you selected ALL the lines of the plot code, from the ggplot() line through theme_minimal(). Each line ends with a + that connects it to the next line. Running only partial lines produces an error.
  • If you see the Plots pane but it is blank, try clicking the Plots tab in the bottom-right pane to make sure it is active.
  • If you see an error about could not find function "ggplot", rerun library(tidyverse) at the top of your script.

Still stuck? Help me figure out why my ggplot scatter plot is not rendering

Spot the missing values warning

Look at the Console after running the plot. You should see a warning about rows being removed due to missing values. This is ggplot2 telling you that 2 rows in the dataset had NA values in the plotted columns, so those points could not be drawn.

Why does this matter?

Those 2 rows are just the tip of the iceberg. The dataset has 19 missing values scattered across several columns. You already saw the NA's counts in your summary() output earlier.

If you ignore missing data, it will silently corrupt your calculations. In the next step, you will clean these missing values out before running any analysis.

✔️ Awesome, I've got everything!

Great. Double check you have saved penguin_analysis.R by pressing Ctrl+S (Windows) or Cmd+S (macOS).

ⓧ I'd like to double check the full code

Here is the complete penguin_analysis.R file so far. Make sure your script matches this exactly.

# penguin_analysis.R
# Analyze Antarctic Penguin Biodiversity with R
# Data: Palmer Station LTER, Antarctica

# Load libraries
library(tidyverse)
library(palmerpenguins)

# ---- Explore the Data ----

# First look at the data structure
glimpse(penguins)

# Summary statistics for every variable
summary(penguins)

# Count penguins by species
penguins %>%
  count(species)

# First scatter plot: flipper length vs body mass
ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g, color = species, shape = species)) +
  geom_point(size = 2) +
  scale_color_manual(values = c("darkorange", "darkorchid", "cyan4")) +
  labs(
    title = "Penguin Flipper Length vs Body Mass",
    x = "Flipper Length (mm)",
    y = "Body Mass (g)"
  ) +
  theme_minimal()
Practice what you learned

Before moving on, try these exercises on your own. They use the same functions you just learned, applied to different variables and questions from the same dataset.

✔️ I'll give it a go!

  • Use count() to tally penguins by island instead of species. Which island has the most penguins?
  • Use count() with two variables: count(species, island). Read the output and identify which species appears on only one island.
  • Modify your scatter plot to show bill_length_mm on the x-axis and bill_depth_mm on the y-axis (still colored by species). Notice how the species separate differently than in the flipper/mass plot.

Try these in the Console rather than adding them to your script. This keeps your script clean while you experiment.

ⓧ Maybe a little help...

Here are hints and solutions for each exercise.

Exercise 1: Replace species with island in the count() call. Biscoe has the most penguins (168).

penguins %>%
  count(island)

Exercise 2: Pass both grouping variables to count(). You will see that Chinstrap only appears on Dream island, and Gentoo only appears on Biscoe.

penguins %>%
  count(species, island)

Exercise 3: Swap the x and y aesthetics in your aes() call. The species form tighter, more overlapping clusters in bill space compared to the flipper/mass plot.

ggplot(penguins, aes(x = bill_length_mm, y = bill_depth_mm, color = species, shape = species)) +
  geom_point(size = 2) +
  scale_color_manual(values = c("darkorange", "darkorchid", "cyan4")) +
  labs(
    title = "Penguin Bill Length vs Bill Depth",
    x = "Bill Length (mm)",
    y = "Bill Depth (mm)"
  ) +
  theme_minimal()

You now know how to meet any new dataset: glimpse() for structure, summary() for statistics, and count() for group sizes. Next up, you will clean those missing values and use dplyr to wrangle the data into grouped summaries that reveal what makes each species unique.

Clean and Wrangle the Data with dplyr

Your scatter plot showed three distinct species clusters. But remember the warning about removed rows? Those 19 missing values are still sitting in your dataset, and if you leave them in, every calculation you run next will either skip them silently or return NA.

Cleaning the data first means your summary statistics will be reliable. In this step, you'll learn the dplyr verbs that form the backbone of every R analysis: drop_na() for cleaning, group_by() and summarize() for grouped statistics, and filter() and select() for drilling into specific rows and columns.

In this step, get ready to:
  • Clean missing values from the dataset using drop_na().
  • Compute species-level summary statistics using group_by() and summarize().
  • Discover which species inhabit which islands and find the heaviest penguin in each species.
Remove missing values with drop_na()

The penguins dataset has 344 rows, but 19 of those rows have at least one missing measurement. drop_na() removes any row that contains a missing value. Since the missing data is a small proportion of the total (about 5.5%), this approach is safe here.

  • In penguin_analysis.R, scroll to the bottom of your file (below the scatter plot code) and add this cleaning code:
# ---- Clean and Wrangle ----

# Remove rows with missing values
penguins_clean <- penguins %>%
  drop_na()

# Verify: 333 rows remain
nrow(penguins_clean)

What does this code do?

  • penguins %>% drop_na() pipes the raw dataset into drop_na(), which removes every row that has at least one NA value.
  • penguins_clean stores the result so you can use the cleaned data in every step going forward.
  • nrow(penguins_clean) prints the row count so you can confirm the cleaning worked.
  • Select these lines and press Ctrl+Enter to run them in the Console.
  • You should see 333 printed in the Console. That confirms 19 incomplete rows were removed (344 minus 19 equals 333).

Seeing NA or a different number?

  • Make sure you ran library(palmerpenguins) earlier in the script. Without it, the penguins object does not exist.
  • Check that you ran all three lines, not just the nrow() call. The penguins_clean object needs to be created first.

Why isn't my drop_na() returning 333 rows?

Compute species-level summary statistics

Now that you have clean data, you can ask the question: how do the three species differ? The group_by() function splits your data into groups. summarize() then collapses each group into one row of computed statistics.

  • Add this pipeline below your nrow() line in penguin_analysis.R:
# Compute summary statistics by species
species_summary <- penguins_clean %>%
  group_by(species) %>%
  summarize(
    count = n(),
    avg_bill_length = mean(bill_length_mm),
    avg_bill_depth = mean(bill_depth_mm),
    avg_flipper_length = mean(flipper_length_mm),
    avg_body_mass = mean(body_mass_g)
  )

print(species_summary)

What does this pipeline do?

  • group_by(species) splits the 333 rows into three groups: one per species.
  • summarize() collapses each group into a single row. Inside it, n() counts the rows and mean() computes the average for each measurement.
  • The result is a tidy 3-row tibble where each row is one species and each column is a summary statistic. This group_by() + summarize() pattern is the single most useful pattern in dplyr.

Before you run this, which species do you think will have the largest average body mass? Think back to the scatter plot you built in the last step.

  • Select the entire pipeline (from species_summary through print(species_summary)) and press Ctrl+Enter to run it.

You should see a 3-row tibble in the Console. Gentoo penguins are the largest species with an average body mass around 5076g. Adelie and Chinstrap are similar in mass, but look closely at the bill measurements. They differ in bill shape.

Seeing an error about penguins_clean?

  • If you see object 'penguins_clean' not found, go back and run the drop_na() chunk from the previous substep first.
  • Make sure you selected all lines of the pipeline, from species_summary <- down through the closing ). Running only part of a pipeline produces an error.

Why is my species summary pipeline returning an error?

Discover species distribution and find the heaviest penguins

You know the average measurements for each species. But which islands do they live on? Using count() with two variables creates a cross-tabulation that answers this question.

  • Add this pipeline below print(species_summary) in your script:
# Which species live on which islands?
penguins_clean %>%
  count(species, island) %>%
  arrange(species)

What does this code do?

  • count(species, island) tallies how many penguins belong to each species-island combination.
  • arrange(species) sorts the output alphabetically by species so it is easy to read.
  • Select these lines and press Ctrl+Enter.

You should see that Adelie penguins appear on all three islands (Biscoe, Dream, Torgersen). Chinstrap are found only on Dream. Gentoo are found only on Biscoe.

Now let's find the heaviest individual penguin in each species. This uses two new dplyr verbs: filter() to keep only rows matching a condition, and select() to pick the columns you want to see.

  • Add this pipeline below the species-island count code:
# Heaviest penguin in each species
penguins_clean %>%
  group_by(species) %>%
  filter(body_mass_g == max(body_mass_g)) %>%
  select(species, island, body_mass_g, sex)

What does this pipeline do?

  • group_by(species) splits the data into species groups, just like before.
  • filter(body_mass_g == max(body_mass_g)) keeps only the row with the maximum body mass within each group.
  • select(species, island, body_mass_g, sex) narrows the output to just the four columns you care about, instead of showing all eight.
  • Select these lines and press Ctrl+Enter.
  • You should see one row per species showing the island, body mass, and sex of each species' heaviest individual.
  • Save your file by pressing Ctrl+S.

Getting unexpected results from filter()?

  • Make sure you used == (double equals, comparison) and not = (single equals, assignment). This is a common R mistake.
  • Check that you have group_by(species) before the filter() line. Without it, you will get only one row total instead of one per species.

Why is my filter + select pipeline not showing the right results?

Practice what you've learned

You've now used six dplyr verbs. Before moving on, try these exercises in the Console to build muscle memory. You do not need to save these into your script.

  • Rewrite the species_summary pipeline but use group_by(island) instead of group_by(species). Which island has the heaviest penguins on average? Think about why that might be.
  • Use filter(species == "Adelie") to create a dataset of only Adelie penguins. Pipe it into summarize() to compute the min(), max(), and sd() of body_mass_g.
  • Try group_by(species) %>% summarize(count = n(), median_flipper = median(flipper_length_mm)) to practice using median() instead of mean(). Compare the median and mean values. Are they close, or is the data skewed?
  • Use mutate(mass_kg = body_mass_g / 1000) on penguins_clean to add a new column that converts grams to kilograms. Pipe the result into head() to see the first few rows.

What does mutate() do?

mutate() adds a new column to your tibble based on a calculation you define. The original columns stay untouched. You will use this verb in the Secret Mission to compute derived metrics like mass ratios.

✔️ Awesome, I've got everything!

Make sure you've saved penguin_analysis.R with all three new sections: the cleaning code, the species summary pipeline, and the species-island and heaviest-penguin queries.

ⓧ I'd like to double check the full code

Your file is now over 30 lines, so here is a checklist of everything that should be in penguin_analysis.R, in order. Scroll through your script and confirm each section is present.

  • File header comments and library(tidyverse) / library(palmerpenguins) at the top.
  • The # ---- Explore the Data ---- section with glimpse(), summary(), count(species), and the scatter plot.
  • The # ---- Clean and Wrangle ---- section header.
  • penguins_clean <- penguins %>% drop_na() and nrow(penguins_clean).
  • The species_summary pipeline with group_by(species) %>% summarize(...) and print(species_summary).
  • count(species, island) %>% arrange(species) for species-island distribution.
  • The heaviest-penguin pipeline with group_by(species) %>% filter(...) %>% select(...).

You now have a clean dataset and clear species-level summaries. You know Gentoo penguins are the largest, Adelie are the most widespread, and Chinstrap live on just one island. Next up, you'll build a plot that tells a convincing lie, and then catch it in the act.

Build a Plot That Lies to You

Your species_summary table just confirmed what the scatter plot hinted at: these three penguin species differ in bill shape, flipper length, and body mass. Every number tells a clean, species-level story.

But what happens if you ignore species and look at all 333 penguins as a single group? In this step, you'll build a scatter plot that tells a compelling story about bill shape. The only problem is that the story is completely wrong.

In this step, get ready to:
  • Build a scatter plot of bill length vs bill depth without species coloring.
  • Add a trend line and discover a misleading relationship.
Create a scatter plot without species grouping

In Step 2, your scatter plot mapped each species to a different color. That color separation revealed three tidy clusters. This time, you'll strip away the species information and see what the raw data looks like without that context.

  • In penguin_analysis.R, scroll to the bottom of the file, below the "Heaviest penguin in each species" section.
  • Add this code to create a new scatter plot of bill length vs bill depth:
# ---- The Misleading Relationship ----

# Plot bill length vs bill depth WITHOUT species grouping
ggplot(penguins_clean, aes(x = bill_length_mm, y = bill_depth_mm)) +
  geom_point(size = 2) +
  labs(
    title = "Bill Length vs Bill Depth",
    subtitle = "All penguins combined",
    x = "Bill Length (mm)",
    y = "Bill Depth (mm)"
  ) +
  theme_minimal()

What does this code do?

  • # ---- The Misleading Relationship ---- is a section comment that labels this part of your script. It matches the same style you used for "Explore the Data" and "Clean and Wrangle."
  • aes(x = bill_length_mm, y = bill_depth_mm) maps bill length to the x-axis and bill depth to the y-axis. Notice there is no color = species this time. Every point will be the same default black.
  • The labs() subtitle says "All penguins combined" to remind you that species are deliberately lumped together.
  • Select the entire new code block (from the section comment through theme_minimal()) and run it by pressing Ctrl+Enter.

A scatter plot appears in the Plots pane. All 333 points are the same color, forming a wide, messy cloud.

Don't see the plot?

  • Make sure you selected all lines from ggplot( through theme_minimal() before pressing Ctrl+Enter. Running only part of the call produces an error.
  • If you see "object 'penguins_clean' not found", you need to re-run the cleaning code from Step 3 first. Select everything from penguins_clean <- penguins %>% through drop_na() and run it.
  • Check that the Plots tab is selected in the bottom-right panel of RStudio.
  • My scatter plot still isn't showing up in RStudio. What else could be wrong?
Add a trend line and question the result

A scatter plot shows raw data points, but sometimes you need a trend line to see the overall direction. The geom_smooth() function in ggplot2 fits a model to your data and draws the result as a line through the points.

  • In penguin_analysis.R, find this line in your new plot code:
  geom_point(size = 2) +
  • Add this line directly below it:
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE, color = "red") +

What does geom_smooth do?

  • geom_smooth() fits a model to your data and draws a trend line through the scatter of points.
  • method = "lm" tells it to fit a linear model. That means a straight line of best fit, the simplest way to summarize the direction of the data.
  • formula = y ~ x makes the formula explicit. Without it, ggplot2 prints an informational message about using a default formula.
  • se = FALSE hides the shaded confidence band around the line, and color = "red" makes the trend line red so it stands out clearly against the black points.

Before you run this, look at the scatter of points in your current plot. Do you think the trend line will slope upward (longer bills go with deeper bills), slope downward (longer bills go with shallower bills), or stay flat?

  • Select the entire plot code block (from ggplot( through theme_minimal()) and run it with Ctrl+Enter.
  • Save penguin_analysis.R by pressing Ctrl+S.

A red line slopes clearly downward across all 333 data points. According to this trend, penguins with longer bills have shallower bills.

Does this make biological sense?

Think about it. If a penguin's bill grows longer, would you expect it to also grow thinner? That seems unlikely. A bill is a 3D structure. Larger bills tend to be larger in every dimension.

Yet the data, plotted this way, says the opposite. If you put this plot in a research paper, your conclusion would be wrong. Something is hiding in the data that this plot does not show.

Seeing a flat or upward trend line?

  • Make sure you are using penguins_clean (not penguins) as the data argument in ggplot().
  • Check that the x-axis is bill_length_mm and the y-axis is bill_depth_mm. Swapping them changes the visual but should still show a negative slope.
  • Confirm there is no color = species inside your aes(). The whole point is that species is NOT mapped here.
  • My trend line doesn't look right. Can you help me debug my ggplot code?
Practice on your own

These exercises are for exploration. Run them in the Console to experiment without changing your saved script.

  • Build a scatter plot of flipper_length_mm vs body_mass_g WITHOUT coloring by species, and add a geom_smooth(method = "lm") trend line. Is the pooled trend misleading here too, or does it tell the truth this time? Compare it to the species-colored version you built in Step 2.
  • In your bill length vs depth plot, try changing method = "lm" to method = "loess". What changes visually? loess fits a curved line instead of a straight one. When might a curved fit be more appropriate?

✔️ I'll give it a go!

Give both exercises a try. You have all the building blocks from the code you just wrote. Swap variables in aes() and change the method argument to see what happens.

ⓧ Maybe a little help...

For the first exercise, paste this into the Console:

# Pooled flipper vs body mass (no species color)
ggplot(penguins_clean, aes(x = flipper_length_mm, y = body_mass_g)) +
  geom_point(size = 2) +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE, color = "red") +
  labs(
    title = "Flipper Length vs Body Mass",
    subtitle = "All penguins combined",
    x = "Flipper Length (mm)",
    y = "Body Mass (g)"
  ) +
  theme_minimal()

The pooled trend here slopes upward, and it matches the within-species trends from Step 2. This relationship is not misleading because all three species follow the same direction.

For the second exercise, change method = "lm" to method = "loess" in the bill length vs depth code and re-run it. The line curves instead of staying straight. A curved fit is useful when the relationship between two variables is not constant across the range of values.

✔️ Awesome, I've got everything!

Great. Double check you've saved penguin_analysis.R.

ⓧ I'd like to double check the full code

Your penguin_analysis.R file is getting long. Scroll through it and confirm these sections appear in order:

  • library(tidyverse) and library(palmerpenguins) at the top.
  • # ---- Explore the Data ---- section with glimpse(), summary(), count(), and the flipper vs body mass scatter plot.
  • # ---- Clean and Wrangle ---- section with drop_na(), species_summary, island counts, and heaviest penguin queries.
  • # ---- The Misleading Relationship ---- section with the bill length vs bill depth scatter plot and the red geom_smooth() trend line.

You just built a plot that would fool anyone. The negative trend looks real, looks clean, and would pass a quick glance in a report. But your biological instinct says something is off. In the next step, you'll add one piece of information back to the plot and watch that entire trend reverse direction.

Reveal the True Pattern and Build Publication-Quality Visualizations

You built a plot that tells a convincing lie. That downward trend line suggests penguins with longer bills have shallower bills. If that finding ended up in a research paper, the conclusion would be wrong.

The fix turns out to be one extra argument in your code. In this step, you will reveal the true pattern hiding inside that misleading trend and learn to recognize Simpson's Paradox. You will also build two more plot types that are polished enough for presentations and reports.

In this step, get ready to:
  • Reveal Simpson's Paradox by adding species coloring to the misleading scatter plot.
  • Build a faceted boxplot comparing body mass across species and islands.
  • Build a layered histogram showing how flipper lengths separate the three species.
Color by species to reveal the true trend

The misleading plot from Step 4 treated all 333 penguins as one group. What happens when you tell ggplot2 to treat each species separately?

Before you add the code below, make a prediction: do you think the trend lines will still slope downward within each species, or will they change direction?

  • Add the following code to the bottom of penguin_analysis.R, below your misleading plot code:
# ---- Reveal the True Pattern ----

# Same plot, but colored by species
ggplot(penguins_clean, aes(x = bill_length_mm, y = bill_depth_mm, color = species)) +
  geom_point(size = 2) +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE) +
  scale_color_manual(values = c("darkorange", "darkorchid", "cyan4")) +
  labs(
    title = "Bill Length vs Bill Depth by Species",
    subtitle = "Simpson's Paradox: the trend reverses within each species",
    x = "Bill Length (mm)",
    y = "Bill Depth (mm)"
  ) +
  theme_minimal()

What does this code do?

  • color = species inside aes() maps each species to a different color. This one change also tells geom_smooth() to fit a separate trend line per group instead of one line across all data.
  • geom_smooth(method = "lm", formula = y ~ x, se = FALSE) draws a straight-line fit for each species. se = FALSE hides the confidence band so the lines are cleaner.
  • scale_color_manual() uses the same three colors from the official palmerpenguins examples so your plots stay visually consistent.
  • Select all the new code and run it by pressing Ctrl+Enter.

Three separate trend lines appear in the Plots pane, and every single one slopes upward. The pooled negative trend was hiding three positive within-species trends.

What is Simpson's Paradox?

Simpson's Paradox is when a trend that appears in combined data reverses or disappears when you split the data into subgroups. Here, the pooled data suggested longer bills go with shallower bills. But within each species, longer bills actually go with deeper bills.

The paradox happens because the three species occupy different regions of the plot. Gentoo penguins have long bills but shallow depth, while Adelie have short bills with deep depth. Mixing them creates a false negative slope. The fix is always the same: break the data down by the relevant grouping variable before drawing conclusions.

Only seeing one trend line?

  • Make sure color = species is inside aes(), not outside it. Placing it outside the aes() call sets a fixed color instead of mapping species to colors.
  • Check that you removed color = "red" from geom_smooth(). If a fixed color is set there, all three lines draw in the same color and look like one.

Why is my species-colored scatter plot only showing one trend line?

Build a faceted boxplot of body mass

Scatter plots show relationships between two continuous variables. Boxplots are the right choice when you want to compare distributions across groups, showing the median, spread, and outliers at a glance.

Faceting takes this further by splitting the plot into separate panels for each level of a categorical variable. You will facet by island to see how species body mass compares across the three islands in the Palmer Archipelago.

  • Add the following code to the bottom of penguin_analysis.R:
# Body mass by species and island
ggplot(penguins_clean, aes(x = species, y = body_mass_g, fill = species)) +
  geom_boxplot(alpha = 0.7, show.legend = FALSE) +
  facet_wrap(~island) +
  scale_fill_manual(values = c("darkorange", "darkorchid", "cyan4")) +
  labs(
    title = "Penguin Body Mass by Species and Island",
    subtitle = "Palmer Archipelago, Antarctica (2007-2009)",
    x = "Species",
    y = "Body Mass (g)",
    caption = "Data: palmerpenguins | Palmer Station LTER"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold"),
    strip.text = element_text(size = 12, face = "bold")
  )

What does this code do?

  • fill = species maps species to the box fill color. show.legend = FALSE hides the legend since the x-axis already labels each species.
  • facet_wrap(~island) creates one panel per island. The ~ (tilde) means "split by."
  • labs() adds a title, subtitle, axis labels, and a data source caption. theme() adjusts font sizes for the title and panel labels to make the plot presentation-ready.
  • Select the new boxplot code and run it by pressing Ctrl+Enter.

Three panels appear in the Plots pane, one per island. Notice that not all species appear on every island. Chinstrap penguins only show up on Dream, and Gentoo only on Biscoe.

Seeing an error about fill?

  • Check that you used fill = species (not color = species) inside aes(). Boxplots use fill for the box interior. color only controls the border.
  • Make sure scale_fill_manual() matches the aesthetic you used. If you used fill in aes(), the scale function must be scale_fill_manual(), not scale_color_manual().

My faceted boxplot isn't displaying correctly. Can you help?

Build a layered histogram of flipper lengths

Histograms show the shape of a single variable's distribution. When you layer multiple species on the same histogram, you can see exactly where their ranges overlap and where they separate.

  • Add the following code to the bottom of penguin_analysis.R:
# Flipper length distribution by species
ggplot(penguins_clean, aes(x = flipper_length_mm, fill = species)) +
  geom_histogram(alpha = 0.6, position = "identity", binwidth = 2) +
  scale_fill_manual(values = c("darkorange", "darkorchid", "cyan4")) +
  labs(
    title = "Flipper Length Distribution by Species",
    x = "Flipper Length (mm)",
    y = "Count",
    fill = "Species"
  ) +
  theme_minimal()

What does this code do?

  • position = "identity" tells ggplot2 to place each species' bars at their true position. Without it, the bars stack on top of each other. Stacking hides the shape of each distribution.
  • alpha = 0.6 makes the bars 60% opaque so overlapping regions are visible. A value of 1 is fully solid and 0 is invisible.
  • binwidth = 2 sets each bar to span 2mm of flipper length. Smaller bins show more detail. Larger bins smooth out the shape.
  • Save penguin_analysis.R by pressing Ctrl+S.
  • Select the new histogram code and run it by pressing Ctrl+Enter.

The histogram appears in the Plots pane. Gentoo flipper lengths sit clearly to the right of Adelie and Chinstrap, which overlap in the middle. The three species separate into distinct clusters.

Bars look stacked instead of overlapping?

  • Make sure you included position = "identity" inside geom_histogram(). The default position is "stack", which piles bars on top of each other.
  • If you can see the bars overlapping but one species is hidden, check that alpha = 0.6 is set. Without transparency, the front bars completely cover the ones behind.

Help me fix my layered histogram so the species overlap correctly.

✔️ Awesome, I've got everything!

Great. Double check you have saved penguin_analysis.R.

ⓧ I'd like to double check the full code

Your file is long enough that a single code block would be hard to compare. Scroll through your penguin_analysis.R and confirm these sections appear in order:

  • Library loading: library(tidyverse) and library(palmerpenguins).
  • Explore the Data: glimpse(), summary(), count(species), and the first scatter plot of flipper length vs body mass.
  • Clean and Wrangle: drop_na(), nrow() verification, species_summary, island counts, and the heaviest-penguin query.
  • The Misleading Relationship: bill length vs bill depth scatter plot with a single red trend line.
  • Reveal the True Pattern: the same scatter plot now colored by species with three upward-sloping trend lines.
  • Faceted boxplot: body mass by species, split into three island panels with a title, subtitle, and caption.
  • Layered histogram: flipper length distribution filled by species with position = "identity" and transparency.
Practice on your own

Try these exercises before moving on. Each one uses a tool you just learned, applied to a different variable or question from the same dataset.

  • Build a boxplot of bill_length_mm by species with no faceting. Add a title and custom fill colors. Write it from scratch without looking at the previous boxplot code.
  • Build a histogram of body_mass_g filled by species. Try three different binwidth values: 50, 200, and 500. Notice how the bin size changes the story the plot tells.
  • Take your faceted boxplot and change facet_wrap(~island) to facet_wrap(~year). Now you see body mass by species across the three study years (2007, 2008, 2009). Does body mass look stable across years?

✔️ I'll give it a go!

You have all the building blocks. The ggplot2 formula is always the same: ggplot(data, aes()) + geom_*() + scale_*() + labs() + theme_*(). Swap in the variable names and you are set.

ⓧ Maybe a little help...

For the boxplot, start with the same structure as your body mass boxplot but replace y = body_mass_g with y = bill_length_mm. Remove the facet_wrap() line entirely.

For the histogram, copy your flipper length histogram and change x = flipper_length_mm to x = body_mass_g. Then change binwidth = 2 to each of the three values and run it each time to compare.

For the faceted boxplot change, just swap ~island for ~year inside facet_wrap() and run the code again.

Secret mission

Investigate Sexual Dimorphism

You've compared species, but what about differences within each species? In this mission, you'll investigate whether male and female penguins differ in body mass. You'll visualize the gap with a faceted boxplot and quantify exactly how much larger males are using a new data reshaping technique called pivot_wider().

Clean Up Your Resources

Clean Up Your Resources

Decide whether to keep your resources running, pause them to come back later, or delete them entirely. This project runs entirely locally with no cloud services, so there are no ongoing costs.

Resources you used:

  • penguin_analysis.R script file.
  • tidyverse R package (installed from CRAN).
  • palmerpenguins R package (installed from CRAN).

Keep everything running

No action needed. Choose this if you plan to keep exploring the penguin data or use your script as a template for other datasets.

  • Your R packages stay installed and available in any future R session.
  • Your penguin_analysis.R script stays saved on your machine.
  • The tidyverse package is a general-purpose R toolkit. You will likely reuse it in future R projects.

Pause - I'll come back to this later

Close RStudio to free up memory. Your script file and installed packages remain on your machine for when you return.

  • Close RStudio. Your penguin_analysis.R file stays saved wherever you stored it.
  • The tidyverse and palmerpenguins packages remain installed. You just need to call library() again when you open a new R session.

Delete - I don't want to use this again

Remove all project resources and start fresh.

  • Remove the palmerpenguins data package by running this in the RStudio Console:
remove.packages("palmerpenguins")
  • Remove the tidyverse metapackage by running:
remove.packages("tidyverse")

What about ggplot2, dplyr, and the other tidyverse packages?

The remove.packages() function only removes the specific package you name. It does not remove dependencies.

Running remove.packages("tidyverse") removes the tidyverse wrapper, but the individual packages it installed (like ggplot2, dplyr, and tidyr) stay on your machine. If you want a complete removal, run remove.packages() on each individual package as well.

  • Delete your penguin_analysis.R file. Open File Explorer, navigate to the folder where you saved the file, right-click it, and select Delete.

Nice Work!

Nice Work!

You just built a complete data analysis of real Antarctic penguin research data, from raw observations all the way to publication-quality visualizations. You caught a trend that looked real but was actually backwards, and you know the fix that makes sure it never fools you again.

You've learned how to:

  • Run the full R data analysis workflow (load, inspect, clean, summarize, visualize) using dplyr verbs like filter(), group_by(), summarize(), and mutate() to wrangle real research data from 333 penguin observations.
  • Build three types of ggplot2 visualizations from scratch: scatter plots for relationships between continuous variables, boxplots for comparing distributions across groups, and histograms for revealing the shape of a single variable. You know how to layer on color, facets, titles, captions, and custom themes to make them presentation-ready.
  • Spot Simpson's Paradox in action. You built a pooled trend line that pointed the wrong direction, then revealed the true pattern by grouping by species. You will always check for hidden subgroup reversals before trusting an aggregated trend.
  • Secret Mission: Investigate sexual dimorphism by reshaping data with pivot_wider() and computing male-to-female mass ratios to quantify how much larger males are in each species.

Ready to quiz yourself?