Activity 06 - Summary statistics

Mean, median, spread, and describing groups the tidy way

tidyverse

descriptive-stats

Hands-on companion to the summary statistics lecture. Handle NAs correctly with sum(!is.na()), compute mean, median, SD, and SE by group, then summarize many columns at once with across() and skimr.

Author

Bill Perry

Published

September 10, 2026

Describing Our Leaf Data

Recap from Activity 05

Used filter(), select(), mutate(), and arrange() to wrangle the leaf data
Chained all four verbs into a single pipeline
Used the pipe %>% to chain steps

Today’s Objectives

Understand why NA values require special handling with sum(!is.na())
Calculate mean, median, SD, and SE using filter() and pull()
Use group_by() + summarize() to compute stats the tidy way
Summarize many columns at once with across()
Explore data quickly with skimr

How this activity works

You are building one R script this whole class, and you turn it in. Create it now: scripts/06_summary_stats.R.

Every line you run goes in that script — not in the Console, not typed into this page. Type it, don’t paste it.

Start every code chunk in your script with a short # comment that says what it does. The comments are part of the grade.

The top of your script, in this order: a title comment, then your library() calls, then the line that loads the data into leaf_df.

Code marked ▶ Run this is typed into your script exactly as shown. Code marked ✏️ Your turn is a change you make in that same script and run.

🔮 Predict before you run

Before you run any ▶ Run this block, predict the answer — a row count, a mean, which group comes out larger. Predicting is what turns “I saw it on a slide” into “I can write it.”

Part 1 · Load libraries and data

📂 Get the data

Same leaf file you have used since Activity 02 — it should already be in your project’s data/ folder. If not: 2026_09_03_data_sci_leaf_area.xlsx → put it in data/.

▶ Type this at the very top of scripts/06_summary_stats.R:

# ---- Activity 06: Summary statistics ---------------------
# your name, today's date

# ---- Libraries -------------------------------------------
library(readxl)      # read Excel files
library(tidyverse)   # dplyr + ggplot2
library(janitor)     # clean_names()
library(skimr)       # fast descriptive summaries

# ---- Load data ------------------------------------------
leaf_df <- read_excel("data/2026_09_03_data_sci_leaf_area.xlsx") %>%
  clean_names()

glimpse(leaf_df)     # look at the data right after loading

Part 2 · The NA problem — counting observations correctly

What is an NA?

NA stands for “Not Available” — it marks a missing value. Missing data are common in real ecological studies. Our leaf data has some: not every leaf got a paper_mass_g measurement.

The problem with `length()`

▶ Run this:

# a small vector with two missing values, for demonstration
x <- c(0.4, 0.3, NA, 0.7, NA)

# length() counts ALL positions — including the NAs
length(x)

length(x) returned:
Number of real (non-NA) values:

Breaking down `is.na()` step by step

▶ Run each line one at a time:

is.na(x)          # TRUE where a value is missing
!is.na(x)         # ! flips it: TRUE where a value is REAL
sum(!is.na(x))    # sum() counts the TRUEs = the real values

✏️ Your turn — in your script: Now do the same on the real data. Compare nrow(leaf_df) to sum(!is.na(leaf_df$paper_mass_g)).

nrow(leaf_df):
sum(!is.na(leaf_df$paper_mass_g)):
How many paper_mass_g values are missing?

⚠️ Watch out! Any stats function on a vector containing NA returns NA — not an error. R will not tell you something went wrong. Always include na.rm = TRUE.

Part 3 · Descriptive statistics in base R

Extract one group with `filter()` and `pull()`

▶ Run this:

# pull the sunny leaf masses out as a plain vector
sunny_mass <- leaf_df %>%
  filter(shade == "sunny") %>%
  pull(mass_g)

sunny_mass    # print the raw values

✏️ Your turn — in your script: Write the code to extract the shady leaf masses as a vector called shady_mass.

Mean

▶ Run this:

# mean mass for the sunny side
mean_sunny <- mean(sunny_mass, na.rm = TRUE)
cat("Mean sunny mass:", round(mean_sunny, 3), "g\n")

✏️ Your turn — in your script: Calculate the mean for the shady side. Store it as mean_shady and print it.

Mean shady leaf mass (g):
Is the shady mean larger than sunny?  Y / N

Median

▶ Run this:

# median mass for each side
med_sunny <- median(sunny_mass, na.rm = TRUE)
med_shady <- median(shady_mass, na.rm = TRUE)

cat("Median sunny:", round(med_sunny, 3), "\n")
cat("Median shady:", round(med_shady, 3), "\n")

✏️ Your turn: Compare the mean and median for the sunny side. If they are about equal, what does that tell you about the shape of the distribution?

Mean sunny:
Median sunny:
What it suggests about distribution shape:

Standard Deviation

▶ Run this:

# SD for each side
sd_sunny <- sd(sunny_mass, na.rm = TRUE)
sd_shady <- sd(shady_mass, na.rm = TRUE)

cat("SD sunny:", round(sd_sunny, 3), "g\n")
cat("SD shady:", round(sd_shady, 3), "g\n")

Which side has more variable leaf mass?
What does a larger SD tell you biologically?

Standard Error

🔮 Predict first: SE = SD / √n. Since n is larger than 1, will the SE be larger or smaller than the SD you just calculated?

▶ Run this:

# SE = SD / sqrt(n) — n must be the count of REAL values
n_sun <- sum(!is.na(sunny_mass))
n_sha <- sum(!is.na(shady_mass))

se_sunny <- sd_sunny / sqrt(n_sun)
se_shady <- sd_shady / sqrt(n_sha)

cat("n sunny  =", n_sun, "   SE sunny =", round(se_sunny, 3), "g\n")
cat("n shady  =", n_sha, "   SE shady =", round(se_shady, 3), "g\n")

✏️ Your turn: In your own words, what does the SE tell you that the SD does not? And if you collected 40 leaves per side instead of ~25, would the SE get larger or smaller?

SE vs SD:
40 leaves per side → SE gets:

Part 4 · Tidy stats with `group_by()` + `summarize()`

🔮 Predict first: How many rows will stats_df have? (How many values does shade take?)

▶ Run this:

# group_by() + summarize() — the core tidyverse pattern
stats_df <- leaf_df %>%
  group_by(shade) %>%
  summarize(
    n         = sum(!is.na(mass_g)),
    mean_mass = round(mean(mass_g,   na.rm = TRUE), 3),
    med_mass  = round(median(mass_g, na.rm = TRUE), 3),
    sd_mass   = round(sd(mass_g,     na.rm = TRUE), 3),
    se_mass   = round(sd_mass / sqrt(n), 3)
  )

stats_df

💡 Key idea: group_by(shade) splits the data by the shade column, and every summarize() calculation runs separately within each group — that’s why stats_df comes back with one row per shade instead of one row overall.

✏️ Your turn: Do the numbers in stats_df match what you calculated by hand in Part 3?

Y / N — if not, explain any differences:

Summary across multiple variables

▶ Run this:

# mean of all three leaf measurements, by shade
size_stats_df <- leaf_df %>%
  group_by(shade) %>%
  summarize(
    n         = n(),
    mean_mass = round(mean(mass_g,       na.rm = TRUE), 3),
    mean_pet  = round(mean(petiole_mm,   na.rm = TRUE), 1),
    mean_thick = round(mean(thickness_mm, na.rm = TRUE), 3)
  )

size_stats_df

✏️ Your turn: Looking at size_stats_df, fill in the table:

Measurement    | Sunny mean | Shady mean | Shady larger? (Y/N)
---------------|------------|------------|---------------------
mass_g         |            |            |
petiole_mm     |            |            |
thickness_mm   |            |            |

Part 5 · The fast way — `across()`

▶ Run this:

# apply the SAME function to every measurement column at once
mean_all_df <- leaf_df %>%
  group_by(shade) %>%
  summarize(
    across(
      c(mass_g, petiole_mm, thickness_mm),
      ~ round(mean(.x, na.rm = TRUE), 3)
    )
  )

mean_all_df

✏️ Your turn — in your script: Change mean inside across() to sd and run it.

Same numbers as size_stats_df?  Y / N
What changed:

Part 6 · Fast overview with `skimr`

▶ Run this:

# skim() grouped by shade — n_missing, mean, SD, and percentiles at once
leaf_df %>%
  group_by(shade) %>%
  skim()

✏️ Your turn: Look at the n_missing row for paper_mass_g, and the mean vs p50 (median) for mass_g.

n_missing paper_mass_g (sunny / shady):
Are mean and median of mass_g close for each group? What does that suggest about skew?

💡 Key idea: one line of skim() gets you n_missing, mean, SD, and percentiles for every column, grouped by shade — the same numbers you computed by hand in Part 3, without writing filter() or pull() at all.

Part 7 · Review and checkpoint

At this point you should be able to:

Explain why length() is unreliable with NAs and use sum(!is.na()) instead
Use filter() and pull() to extract a group’s values as a vector
Calculate mean, median, SD, and SE with na.rm = TRUE
Use group_by() + summarize() to compute stats for all groups at once
Use across() to apply one function to many columns at once
Use skim() for a fast full-dataset overview

✏️ Your turn — before you move on: Run your entire script top to bottom with Ctrl/Cmd + Shift + Enter (Source). Does it complete without errors?

Ran cleanly?  Y / N
If not, what error appeared:

Extension — out of class (~30–40 min)

Add this to the bottom of scripts/06_summary_stats.R and turn it in with the rest. Put your written answers in # comments right under the code they go with. In class you computed n, mean, SD, and SE. Now you add three spread statistics — one of them new — and reason about what each one tells you.

E1 · Three more spread statistics (4 pts)

Extend your grouped-by-shade summary of mass_g to also include:

iqr_mass — the interquartile range (IQR(), mind na.rm)
range_mass — max(mass_g, na.rm = TRUE) - min(mass_g, na.rm = TRUE)
cv_mass — the coefficient of variation, computed by hand as sd(mass_g) / mean(mass_g) (a unitless measure of relative spread — you have not been shown this one; work it out from the definition)

Show the full table (n, mean, sd, iqr_mass, range_mass, cv_mass) by shade.

E2 · Predict, then check (3 pts)

Before running E1, write in comments:

Which side do you predict has the larger sd? The larger cv_mass? (They need not be the same side — think about why.)
One reason, based on the raw Part-3 numbers you already saw.

Then run E1 and say whether both predictions held. If the bigger-SD side did not have the bigger CV, explain why that can happen.

E3 · Explain it, with YOUR values (3 pts)

In plain language, what does your cv_mass tell you about the relative variability of each side that the raw sd alone does not?
Why is the CV useful for comparing spread between two groups with different means, when the raw SD can mislead?

Getting unstuck

Read the error message out loud. R usually names the line and the problem.
Check the usual suspects: Did you load library(readxl), library(tidyverse), and library(skimr)?
Spelling? R is case-sensitive — "Sunny" ≠ "sunny". Check with names(leaf_df).
na.rm = TRUE missing? Any stat function on a column with NAs returns NA without it — no error, just a silent wrong answer.
Cheat sheets — https://posit.co/resources/cheatsheets/
Bring the exact error (copy-paste it) to class, Canvas, or office hours.

💡 Key idea: a missing na.rm = TRUE, a typo’d column name — these are the errors this activity is built around, and reading the message before you panic will solve most of them.

End of the Summary Statistics activity. Next: Quarto — writing reproducible reports that mix prose, code, and results in one document.

Other Formats

Describing Our Leaf Data

Recap from Activity 05

Today’s Objectives

How this activity works

🔮 Predict before you run

Part 1 · Load libraries and data

Part 2 · The NA problem — counting observations correctly

What is an NA?

The problem with length()

Breaking down is.na() step by step

Part 3 · Descriptive statistics in base R

Extract one group with filter() and pull()

Mean

Median

Standard Deviation

Standard Error

Part 4 · Tidy stats with group_by() + summarize()

Summary across multiple variables

Part 5 · The fast way — across()

Part 6 · Fast overview with skimr

Part 7 · Review and checkpoint

Extension — out of class (~30–40 min)

E1 · Three more spread statistics (4 pts)

E2 · Predict, then check (3 pts)

E3 · Explain it, with YOUR values (3 pts)

Getting unstuck

The problem with `length()`

Breaking down `is.na()` step by step

Extract one group with `filter()` and `pull()`

Part 4 · Tidy stats with `group_by()` + `summarize()`

Part 5 · The fast way — `across()`

Part 6 · Fast overview with `skimr`