Worksheet 11 — Linear Regression

Building a calibration curve to predict leaf area from paper tracing mass

Doing linear regressions and plotting

Regression — Predicting Leaf Area

Recap from Worksheets 07, 09–10

• Quarto — one .qmd file holds your writing and your code and renders to Word, HTML, or PowerPoint
• Three ingredients — YAML front matter, Markdown text, and uniquely named code chunks
• Reproducibility — every number and figure is computed, never copied by hand
• Inline code — a live result drops straight into a sentence with `r `
• Carried forward (Worksheets 09–10) — shady leaves were significantly heavier than sunny leaves (Welch’s t-test, p < 0.001)
• But: weight in grams tells us about mass — not directly about leaf surface area

Today’s Objectives

1. Load and inspect the paper calibration data
2. Make a scatter plot and visually assess linearity
3. Fit a linear regression with lm() and read every line of summary()
4. Extract slope, intercept, and R² from the model
5. Check all four regression assumptions (linearity, equal variance, normality of residuals)
6. Use the calibration equation to predict leaf area from a tracing mass
7. Use predict() with confidence and prediction intervals
8. Write a complete results paragraph following W&S Chapter 17 conventions

How to use this worksheet

• Work through each part in order. Type the code into a new R script in Positron and run it line by line.

  • Code blocks marked ▶ Run this should be executed as written.

  • Blocks marked ✏️ Your turn ask you to write, modify, or interpret.

  • The Going further section is optional.

Note🔮 Predict before you run — and type, don’t paste

Before each ▶ Run this block, cover the output and predict what R will print. Then type the code in yourself rather than pasting it.

A regression handout has a lot of numbers in it — slope, intercept, R², a t-value, a p-value — and it’s tempting to just read them off the screen. Committing to a guess about one of them first (is R² closer to 0.5 or 1? is the slope positive?) gives you a reason to actually look at the output instead of skimming past it. Typing the formula by hand also has a side benefit: it’s how you catch that you flipped area_cm2 ~ mass_g backwards, or typed newdata with the wrong column name — mistakes that copy-paste would have hidden.

Tip📂 The same analysis, two ways — download both and compare

Today you get to feel the difference the Quarto intro talked about. These two files do the exact same regression — one as a plain script, one as a report:

• regression_analysis_script.R — the analysis as a classic R script. Great for running the numbers. To share results you would copy them into a document by hand.
• regression_analysis_report.qmd — the analysis as a Quarto report. Same code, but wrapped in writing, with figures and numbers that render straight into a Word document.

Open them side by side in Positron. As you work through this worksheet, notice which one you would rather hand in.

---

🧩 Chunk 1 — Set up the data (after lecture Chunk 1)

Parts 1–3: load the paper data, inspect it, and make the scatter plot.

Part 1 · Load libraries and data

Libraries

▶ Run this at the top of your script:

# Load packages at the top — always ----------------------
library(readxl)      # reading Excel files
library(tidyverse)   # data manipulation + ggplot2
library(skimr)       # fast dataset overview

The calibration dataset

The file paper_area_weights.xlsx contains 118 paper squares of known area (1–567 cm²) with their measured masses in grams. We will fit a regression to convert mass → area, then apply it to leaf tracings.

▶ Run this:

# Load the paper calibration data -----------------------
paper_df <- read_excel("data/paper_area_weights.xlsx")

✏️ Your turn: Run dim(paper_df). How many rows and columns does the dataset have?

# Write your code here:

Rows:
Columns:
Column names and their roles:
  area_cm2 →  (response / explanatory / circle one)
  mass_g   →  (response / explanatory / circle one)

---

Part 2 · Inspect the data

▶ Run this:

# Quick overview with skim --------------------------------
skim(paper_df)

# How many squares were measured at each area level? -----
paper_df %>%
  group_by(area_cm2) %>%
  summarize(
    n          = n(),
    mean_mass  = round(mean(mass_g), 5),
    min_mass   = round(min(mass_g),  5),
    max_mass   = round(max(mass_g),  5)
  )

✏️ Your turn: Fill in the table from the output above:

area_cm2   mean mass (g)
--------   -------------
1
4
16
64
100
400

✏️ Your turn: Is there a clear pattern — as area increases, what happens to mass? Is this what you expected from your knowledge of paper?

Your answer:

---

Part 3 · First scatter plot

▶ Run this:

# Scatter plot — check the relationship visually ---------
scatter_09_plot <- paper_df %>%
  ggplot(aes(x = mass_g, y = area_cm2)) +
  geom_point(alpha = 0.5, size = 2) +
  labs(
    title = "Paper Mass vs. Known Area",
    x     = "Mass (g)",
    y     = "Area (cm²)"
  ) +
  theme_minimal()

scatter_09_plot

✏️ Your turn: Describe the relationship. Is it linear? Is the variance (spread of points) roughly constant across the range of mass values?

Shape of relationship (linear / curved):
Variance (constant / increases / decreases with X):
Any obvious outliers?  Y / N

✏️ Your turn: In this study, which variable is the explanatory (X) and which is the response (Y)? Why does that direction make sense for our goal of predicting leaf area?

Explanatory (X):
Response (Y):
Reason for this direction:

💡 Key idea: You could measure leaf area directly with a planimeter or a scanner, but a balance is what’s sitting on the lab bench. That’s the whole point of a calibration curve — trade a hard measurement (area) for an easy one (mass) once you’ve established that the two move together.

---

🧩 Chunk 2 — Fit & read the model (after lecture Chunk 2)

Parts 4–7: fit lm(), decode summary(), pull the equation and R².

Part 4 · Fit the regression with lm()

🔮 Predict first: area rises with mass. Before you read summary() — will the slope be positive or negative, and roughly how many cm² per gram?

▶ Run this:

# Fit the least-squares linear model --------------------
# lm(response ~ predictor, data = data_frame)
paper_lm_model <- lm(area_cm2 ~ mass_g, data = paper_df)

# Full summary — every number matters -------------------
summary(paper_lm_model)

✏️ Your turn: Copy the Coefficients table output here:

Paste the Coefficients section of the summary() output:

---

Part 5 · Decode the summary() output line by line

✏️ Your turn: Match each piece of output to its meaning. Fill in the blanks:

(Intercept) estimate = _____
  → This is "a" in Ŷ = a + bX.
    It is the predicted area when mass = _____.
    Does this make physical sense? (a piece of paper with zero mass has zero area)  Y / N

mass_g estimate = _____
  → This is "b" (the slope).
    It means: for every additional 1 gram of paper, area increases by _____ cm².
    Biologically, 1 gram of this paper stock has _____ cm² of surface area.

Std. Error for mass_g = _____
  → The uncertainty in the slope estimate.

t value for mass_g = _____
  → slope / SE.  This tests H₀: β = _____.

Pr(>|t|) for mass_g = _____
  → The p-value for the slope test. Is it significant at α = 0.05?  Y / N

Multiple R-squared = _____
  → _____ % of the variation in area is explained by mass.

F-statistic = _____  on _____ and _____ df
  → Overall test of the model.  p = _____

📖 Whitlock & Schluter §17.3 p. 551–553: The slope is tested with a t-statistic: t = b / SE_b, with df = n − 2. If H₀: β = 0 is rejected, mass significantly predicts area.

---

Part 6 · Extract the calibration equation

▶ Run this:

# Pull intercept and slope from the model ---------------
a_intercept <- coef(paper_lm_model)[1]
b_slope     <- coef(paper_lm_model)[2]

cat("Intercept (a) =", round(a_intercept, 4), "cm²\n")
cat("Slope (b)     =", round(b_slope, 2), "cm²/g\n")

✏️ Your turn: Write the calibration equation for predicting area from mass:

area_cm2 = _____ × mass_g + _____

✏️ Your turn: Use the equation by hand to predict the area of a paper tracing that weighs 0.120 g. Show your work:

area = _____ × 0.120 + _____
     = _____  cm²

---

Part 7 · Extract and interpret R²

🔮 Predict first: For uniform paper, predict R² — near 0, 0.5, or 1? Write your guess before you run it.

▶ Run this:

# Pull R-squared from the model summary -----------------
r2_val <- summary(paper_lm_model)$r.squared
cat("R² =", round(r2_val, 6), "\n")

✏️ Your turn: What does R² = 0.9999 mean in plain language for our calibration?

Your answer:

✏️ Your turn: In ecology, you might see regression results with R² = 0.45. Is that a useful result? What does it mean?

Your answer:

📖 Whitlock & Schluter §17.4 p. 555: “R² measures the fraction of variation in Y that is explained by X.” A high R² indicates a tight fit; a low R² indicates high scatter around the line — but a low R² does not mean the relationship is not real or useful.

---

🧩 Chunk 3 — Visualize & check assumptions (after lecture Chunk 3)

Parts 8–10: plot the curve, then check the residual and QQ/Shapiro assumptions.

Part 8 · Plot the calibration curve

▶ Run this:

# Regression line with 95% confidence band --------------
calibration_plot <- paper_df %>%
  ggplot(aes(x = mass_g, y = area_cm2)) +
  geom_point(alpha = 0.5, size = 1.8) +
  geom_smooth(method = "lm", se = TRUE,
              color = "darkblue", fill = "lightblue") +
  labs(
    title    = "Paper Calibration Curve",
    subtitle = paste0("area = ", round(b_slope, 2),
                      " × mass + ", round(a_intercept, 3),
                      "  |  R² = ", round(r2_val, 4)),
    x        = "Paper Mass (g)",
    y        = "Area (cm²)"
  ) +
  theme_minimal()

calibration_plot

✏️ Your turn: Where is the shaded confidence band narrowest? Where is it widest? Why does the width vary?

Narrowest at:
Widest at:
Why it varies:

Save the plot

▶ Run this:

ggsave("figures/paper_calibration_curve.png",
       plot   = calibration_plot,
       width  = 5, height = 5, units = "in", dpi = 300)

---

Part 9 · Check assumptions — the base R diagnostic plots

🔮 Predict first: plot() on a fitted model gives four panels at once. One is a QQ plot of the residuals — what do you predict the other three will show?

▶ Run this:

# plot() on an lm object gives all four diagnostics -----
par(mfrow = c(2, 2))
plot(paper_lm_model)

✏️ Your turn: Describe panel 1 (Residuals vs Fitted). Do the points form a random cloud around zero, or is there a pattern?

Pattern or random cloud?:
Any evidence of non-linearity (curved pattern)?  Y / N
Any evidence of unequal variance (funnel shape in panel 1 or 3)?  Y / N

⚠️ Watch out! A funnel shape (residuals getting larger at higher fitted values) means variance is not constant — a violation of the equal variance assumption. A curved pattern means the relationship is not linear.

📖 Whitlock & Schluter §17.5 p. 559 (Figure 17.5-4): “A residual plot should show a roughly symmetric cloud of points with no pattern.”

---

Part 10 · Confirming normality — Shapiro-Wilk

▶ Run this:

# Shapiro-Wilk test on residuals — NOT on raw data -----
shapiro.test(residuals(paper_lm_model))

✏️ Your turn: Record the Shapiro-Wilk result:

W = _____   p = _____
Decision (normal / not normal):

✏️ Your turn: Does the Shapiro-Wilk decision agree with what panel 2 (Normal Q-Q) of plot(paper_lm_model) showed you?

Your answer:

✏️ Your turn: An important distinction — in regression, normality is checked on the residuals, not on Y or X directly. Why does that matter?

Your answer:

---

🧩 Chunk 4 — Predict & report (after lecture Chunk 4)

Parts 11–12: predict leaf area from tracing mass, then write the results paragraph.

Part 11 · Predict leaf area from tracing mass

🔮 Predict first: sunny tracing 0.092 g, shady 0.138 g. Which predicts the larger area? Does that match “shady leaves are bigger”? Predict before running.

Method 1 — Direct calculation

▶ Run this:

# Direct calculation using the calibration equation ----
mass_sunny <- 0.092   # g — example sunny leaf tracing
mass_shady <- 0.138   # g — example shady leaf tracing

area_sunny <- b_slope * mass_sunny + a_intercept
area_shady <- b_slope * mass_shady + a_intercept

cat("Sunny tracing:", mass_sunny, "g →",
    round(area_sunny, 2), "cm²\n")
cat("Shady tracing:", mass_shady, "g →",
    round(area_shady, 2), "cm²\n")

✏️ Your turn: Does the shady leaf have a larger predicted area? In Worksheets 09–10 we found shady leaves were significantly heavier. Does a larger predicted area match that finding? What does this tell you about the relationship between leaf weight and leaf area?

Sunny predicted area:
Shady predicted area:
Larger side:
Does larger area match the heavier weight from Worksheets 09–10?  Y / N
What does this tell you about the relationship between weight and area in leaves?

Method 2 — Using predict() with intervals

▶ Run this:

# predict() gives uncertainty around predictions ------
new_masses <- tibble(mass_g = c(0.092, 0.138))

# 95% confidence interval for the MEAN area
predict(paper_lm_model,
        newdata  = new_masses,
        interval = "confidence") %>%
  as_tibble() %>%
  mutate(mass_g = c(0.092, 0.138), .before = 1) %>%
  round(3)

# 95% prediction interval for a SINGLE observation
predict(paper_lm_model,
        newdata  = new_masses,
        interval = "prediction") %>%
  as_tibble() %>%
  mutate(mass_g = c(0.092, 0.138), .before = 1) %>%
  round(3)

✏️ Your turn: The prediction interval is wider than the confidence interval. In your own words, explain why — what additional source of uncertainty does a prediction interval capture?

Your answer:

📖 Whitlock & Schluter §17.2 p. 549: “Confidence bands measure the precision of the predicted mean Y for each value of X. Prediction intervals measure the precision of the predicted single Y-values for each X.”

---

Part 12 · Write a scientific results paragraph

✏️ Your turn: Write a complete results paragraph using the information below. Follow the format from Lecture 11.

Information to include:

• What the regression tested
• F-statistic, df₁, df₂, and p-value (from the F-statistic line in summary())
• R²
• The calibration equation (slope and intercept)
• 95% confidence interval for the slope (use confint(paper_lm_model))
• Brief mention of assumption checks

# Get the 95% CI for the slope -------------------------
confint(paper_lm_model)

Write your results paragraph here:

Tip✍️ Now put it in a report

Open regression_analysis_report.qmd. Its Discussion section already pulls the slope, R², and p-value into the sentences with inline code. Paste your results paragraph in, press Render, and you have a Word document — no numbers copied by hand.

---

Part 13 · Review and checkpoint

At this point you should be able to:

• Identify the response (Y) and explanatory (X) variables for a regression problem
• Make a scatter plot and assess linearity visually
• Fit a linear regression with lm(Y ~ X, data) and call summary()
• Extract slope, intercept, and R² from the model
• Interpret slope, intercept, t-value, p-value, and R² in biological terms
• Run plot(model) and read all four diagnostic panels
• Check the residuals vs fitted panel for linearity and equal variance
• Check the Normal Q-Q panel and Shapiro-Wilk on residuals for normality
• Use coef(model) to write the prediction equation
• Use predict(model, newdata, interval = "prediction") for predictions with uncertainty
• Write a complete scientific results paragraph with F, df, p, R², and equation

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

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

---

Part 14 · Going further

This section is optional — work through it if you finish early or want to push deeper.

Explore the fit at each area group

▶ Try this:

# How well does the model fit within each area group? --
paper_df %>%
  mutate(
    predicted = fitted(paper_lm_model),
    residual  = residuals(paper_lm_model)
  ) %>%
  group_by(area_cm2) %>%
  summarize(
    n             = n(),
    mean_resid    = round(mean(residual), 4),
    max_abs_resid = round(max(abs(residual)), 4)
  )

✏️ Your turn: Are the residuals larger at bigger area values? What does that pattern (or lack of it) tell you about the equal-variance assumption?

Your observation:

Prediction for a real leaf tracing

Suppose you traced a leaf onto the same type of paper used in the calibration, cut it out, and weighed it. The tracing weighs 0.175 g.

▶ Try this:

# Predict area for your leaf tracing -------------------
my_tracing <- tibble(mass_g = 0.175)

predict(paper_lm_model,
        newdata  = my_tracing,
        interval = "prediction") %>%
  round(2)

✏️ Your turn: What is the predicted leaf area? What is the 95% prediction interval? How confident are you in this prediction?

Predicted area:
95% prediction interval:  _____ to _____ cm²
Confidence in prediction:

Extrapolation warning

The calibration data ranges from 1 to 567 cm². Suppose a very large leaf tracing weighs 5.2 g.

✏️ Your turn: Would you trust a prediction for a tracing mass of 5.2 g? Why or why not? (Hint: what is the largest mass in the calibration data?)

Max mass in calibration data:
Is 5.2 g within the calibration range?  Y / N
Should you predict at 5.2 g?  Y / N
Why:

📖 Whitlock & Schluter §17.2 p. 550: “Extrapolation is the prediction of Y at values of X beyond the range of X-values in the data. Extrapolation is problematic because there is no way to ensure the relationship between X and Y continues to be linear.”

Connect to Worksheets 09–10

In Worksheets 09–10 you ran a Welch’s t-test comparing leaf weight (grams) between sunny and shady sides and found a significant difference (p < 0.001). Now you have a calibration curve that converts tracing weight to area (cm²) — a more biologically meaningful measurement. You could use the regression equation as a mutate() step to add predicted area to a leaf data frame, then redo the t-test on area rather than weight.

▶ Sketch the code (do not necessarily run it — you would need actual tracing masses):

# How you WOULD apply the calibration to the leaf data
# (assuming trace_df has columns: side, tracing_mass_g)

# trace_df <- trace_df %>%
#   mutate(predicted_area_cm2 = b_slope * tracing_mass_g + a_intercept)
#
# t.test(predicted_area_cm2 ~ side, data = trace_df,
#        var.equal = FALSE, alternative = "two.sided")

✏️ Your turn: What is the advantage of working in units of area (cm²) rather than raw weight (g) when comparing sunny vs. shady leaves?

Your answer:

---

What your figures/ folder should contain after this worksheet

figures/
├── paper_calibration_curve.png     ← from Part 8

---

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

Turn this in with your worksheet. paper_df is the same shared dataset for everyone — this task can’t be copied, because the object is one you make.

Important

E2 and E3 must be handwritten on paper, photographed, and embedded (![caption](my_photo.jpg)). Typed answers to E2/E3 get at most half credit, even if correct — I want your own eyeballed prediction and your own reasoning about your own numbers.

E1 · Calibrate an object only you have (4 pts)

1. Cut any shape you like from ordinary printer/notebook paper — trace your hand, your initials in bubble letters, a rough state outline.
2. Weigh it in grams (kitchen or postal scale). No scale? Estimate by comparing to a US penny (2.5 g) or a sheet of printer paper (~4.5 g) and say so.
3. In R, use paper_lm_model and predict(..., interval = "prediction") to estimate its area from its mass, with a 95% prediction interval.
4. Measure the real area low-tech: trace onto graph paper and count squares, or approximate with rectangles/triangles you compute by hand.

Report: the shape, its mass, the predicted area + interval, your hand-measured area, and whether the real area fell inside the predicted interval.

E2 · Predict, then check — ✍️ by hand (3 pts)

Before running E1’s predict(), look only at your calibration plot from the in-class part. By hand: mark roughly where your shape’s mass sits on the x-axis, sketch where the line puts the area, and write your guessed area range and why (steep/shallow line? near the middle of the data or an edge?). Photograph it. Then run predict() and add a line: were you close? What did you misjudge?

E3 · Explain it — ✍️ by hand, with YOUR numbers (3 pts)

1. Using your R² from the in-class fit, explain in plain language what fraction of the variation in area is explained by mass — and what R² does not tell you about whether the relationship is real or causal.
2. Your prediction interval in E1 is wider than a confidence interval for the mean area at that mass. Explain why predicting one new object is riskier than estimating an average.
3. Using max(paper_df$mass_g), explain why extrapolating to a much heavier object (a 5 g cardboard cutout) is a different kind of risk than the interpolation you did in E1.

---

Getting unstuck

1. lm() error: the formula must be lm(Y ~ X, data) — response on the left, predictor on the right. Check column names with names(paper_df).
2. coef() gives two values: coef(model)[1] is the intercept, coef(model)[2] is the slope for the first predictor.
3. residuals() vs raw data: always apply shapiro.test() and plot(model) to the fitted model — the residuals they check come from residuals(model), not the raw Y column.
4. predict() needs a tibble: newdata must be a data frame or tibble with the exact same column name as the predictor (mass_g). A typo here is the most common error.
5. geom_smooth(se = TRUE) shows confidence band, not prediction interval. For prediction intervals, use predict() manually.
6. Cheat sheets — https://posit.co/resources/cheatsheets/

💡 Key idea: Plot it, fit it, check the residuals, interpret the coefficients, then predict — you followed those same five steps for this one calibration curve, and you’ll follow them again, mostly unchanged, the first time you run an ANOVA next module.

---

End of the Regression worksheet. Next: downloading and summarizing real climate data.