Activity: Correlation vs Linear Models

Relationships

In-class activity in two halves: code the correlation on the Nazca booby and Indian rope trick data, then fit the lion regression, build its sums of squares by hand, and test its assumptions.

Author

Bill Perry

Worksheet: Correlation and Linear Regression

📥 Booby data 📥 Rope trick data 📥 Lion data 📥 Tools script (.R) 📥 Script skeleton (.R)

This worksheet runs in two halves, matching the lecture. Parts 1–4 are correlation: two measured variables, one number, no direction. Parts 5–8 are regression: a fitted line, the sums of squares it is made of, and the diagnostics that decide whether to believe it.

Three datasets, all from Whitlock & Schluter:

booby_aggression.csv — Example 16.1. Nestling visits and adult aggression in Nazca boobies
rope_trick.csv — Example 16.5. Years since seeing the Indian rope trick, and how impressive it is remembered as
lion_noses.csv — Example 17.1. Proportion of black on a male lion’s nose, and its age

Part 1 · Setup and the data

▶ Run this:

# ── Part 1.1 · Setup and the data ───────────────────────────────────────
library(tidyverse)
library(patchwork)

source("themes_functions/r_themes_and_functions.R")

booby_df <- read_csv("data/booby_aggression.csv")
rope_df  <- read_csv("data/rope_trick.csv")
lion_df  <- read_csv("data/lion_noses.csv")

booby_df

Part 2 · Correlation

▶ Run this — look at the data before you test it:

# ── Part 2.1 · Correlation ──────────────────────────────────────────────
ggplot(booby_df, aes(visits_as_nestling, future_aggression)) +
  geom_point()

✏️ Your turn: Is the trend straight, or is there a curve in it? Any point far off on its own?



____________________________________________________________

▶ Run this — the coefficient, then the test:

# ── Part 2.2 · Correlation ──────────────────────────────────────────────
cor(booby_df$visits_as_nestling, booby_df$future_aggression)

cor.test(booby_df$visits_as_nestling, booby_df$future_aggression)

▶ Run this — square it for the variance explained:

# ── Part 2.3 · Correlation ──────────────────────────────────────────────
r_booby <- cor(booby_df$visits_as_nestling, booby_df$future_aggression)

r_booby^2

✏️ Your turn: Fill these in from the output above.



r = ____________   p = ____________   r squared = ____________

What share of the variation do the two variables share? ____________ %

✏️ Your turn: Write one sentence stating what the correlation means biologically — without using the word “cause.”


____________________________________________________________


____________________________________________________________

Part 3 · Correlation assumptions

Pearson’s r assumes both variables are normal, the relationship is linear, and no outlier is driving it.

▶ Run this — the formal test on each variable:

# ── Part 3.1 · Correlation assumptions ──────────────────────────────────
shapiro.test(booby_df$visits_as_nestling)

shapiro.test(booby_df$future_aggression)

▶ Run this — and the picture:

# ── Part 3.2 · Correlation assumptions ──────────────────────────────────
h1 <- ggplot(booby_df, aes(visits_as_nestling)) + geom_histogram(bins = 10)
h2 <- ggplot(booby_df, aes(future_aggression)) + geom_histogram(bins = 10)
q1 <- ggplot(booby_df, aes(sample = visits_as_nestling)) + geom_qq() + geom_qq_line()
q2 <- ggplot(booby_df, aes(sample = future_aggression)) + geom_qq() + geom_qq_line()

(h1 / h2) | (q1 / q2)

✏️ Your turn: Remember H₀ here is “these data ARE normal.” Do both variables pass?



visits: p = __________ pass? ______    aggression: p = __________ pass? ______

✏️ Your turn: One of them lands very close to 0.05 while its QQ plot looks fine. At n = 24, which do you believe, and does it change the coefficient you report?



____________________________________________________________

Part 4 · When normality fails

Same workflow, a dataset that does not cooperate: an impressiveness score, which is bounded and lumpy.

▶ Run this:

# ── Part 4.1 · When normality fails ─────────────────────────────────────
ggplot(rope_df, aes(years_since_seen, impressiveness)) +
  geom_point()

shapiro.test(rope_df$impressiveness)

▶ Run this — both coefficients on the same data:

# ── Part 4.2 · When normality fails ─────────────────────────────────────
cor.test(rope_df$years_since_seen, rope_df$impressiveness)

cor.test(rope_df$years_since_seen, rope_df$impressiveness,
         method = "spearman")

✏️ Your turn: Which coefficient do you report for these data, and why?


Pearson r = __________   Spearman rho = __________

I report ______________ because ______________________________

✏️ Your turn: Spearman works on ranks. What did you give up by switching to it?

____________________________________________________________

That is the correlation half. The rest of the worksheet fits a line.

Part 5 · Fitting a regression

▶ Run this:

# ── Part 5.1 · Fitting a regression ─────────────────────────────────────
ggplot(lion_df, aes(proportion_black, age_years)) +
  geom_point()

lion_model <- lm(age_years ~ proportion_black, data = lion_df)

summary(lion_model)

✏️ Your turn: Write the fitted equation, then the slope as a sentence in the units of the data.


age = ____________ + ____________ x proportion_black


For every 0.1 more black on the nose, predicted age goes up by ______ years.

✏️ Your turn: The intercept is the predicted age at proportion_black = 0. The youngest lion measured was 1.1 years old. Is the intercept a statement about newborn lions? Why not?


____________________________________________________________

▶ Run this — the line with its confidence band:

# ── Part 5.2 · Fitting a regression ─────────────────────────────────────
ggplot(lion_df, aes(proportion_black, age_years)) +
  geom_point() +
  geom_smooth(method = "lm")

▶ Run this — predict a lion with a half-black nose, two ways:

# ── Part 5.3 · Fitting a regression ─────────────────────────────────────
new_lion <- tibble(proportion_black = 0.50)

predict(lion_model, new_lion, interval = "confidence")

predict(lion_model, new_lion, interval = "prediction")

✏️ Your turn: Which interval is wider, and what question does each one answer?

____________________________________________________________

Part 6 · Sums of squares

The lecture’s centrepiece: build the partition yourself, then let R hand you the same numbers.

▶ Run this:

# ── Part 6.1 · Sums of squares ──────────────────────────────────────────
fitted_age <- predict(lion_model)
mean_age   <- mean(lion_df$age_years)

ss_total <- sum((lion_df$age_years - mean_age)^2)
ss_reg   <- sum((fitted_age - mean_age)^2)
ss_resid <- sum((lion_df$age_years - fitted_age)^2)

c(total = ss_total, regression = ss_reg, residual = ss_resid)

ss_reg + ss_resid

✏️ Your turn: Does the partition add up, and what fraction of the total does the line account for?

ss_reg + ss_resid = __________   ss_total = __________

ss_reg / ss_total = __________   Where have you seen that number before? ______

▶ Run this — the same three numbers from R:

# ── Part 6.2 · Sums of squares ──────────────────────────────────────────
anova(lion_model)

✏️ Your turn: Match the table to what you computed.

Sum Sq, proportion_black row = __________   Sum Sq, Residuals row = __________

Df = ______ and ______      Mean Sq = Sum Sq / ______

▶ Run this — the t on the slope, squared:

# ── Part 6.3 · Sums of squares ──────────────────────────────────────────
summary(lion_model)$coefficients

summary(lion_model)$coefficients[2, 3]^2

✏️ Your turn: Compare that to the F value in the ANOVA table. What does the match tell you about the two tests?

t squared = __________   F = __________

____________________________________________________________

Part 7 · Regression assumptions

Every check from here is on the residuals, not the raw data.

▶ Run this — the four standard diagnostic plots:

# ── Part 7.1 · Regression assumptions ───────────────────────────────────
par(mfrow = c(2, 2))
plot(lion_model)
par(mfrow = c(1, 1))

▶ Run this — the formal test on the residuals:

# ── Part 7.2 · Regression assumptions ───────────────────────────────────
shapiro.test(residuals(lion_model))

✏️ Your turn: One line per assumption. Say what you looked at, not just pass or fail.

Linearity (residuals vs fitted): ____________________________

Equal variance (scale-location): ____________________________

Normality (QQ + Shapiro p = ______): ________________________

Influential points (Cook's distance): _______________________

✏️ Your turn: Would you report this model as it stands? One sentence.

____________________________________________________________

Part 8 · The figure you hand in

Everything above was for looking. This one is for showing, so it gets the full treatment.

▶ Run this:

# ── Part 8.1 · The figure you hand in ───────────────────────────────────
final_plot <- ggplot(lion_df, aes(proportion_black, age_years)) +
  geom_point(size = 2.5) +
  geom_smooth(method = "lm", colour = "blue") +
  AND MAKE IT YOUR OWN

final_plot

ggsave(
  "figures/yourname_lion_regression.pdf",
  plot   = final_plot,
  width  = 6,
  height = 4,
  units  = "in"
)

✏️ Your turn: Open the PDF and check the axis labels are readable before you submit it.

Review and checkpoint

At this point you can:

Run and interpret cor.test(), and square r for the variance explained
Check correlation assumptions on the raw variables and switch to Spearman when they fail
Fit a line with lm() and state the slope in the units of the data
Build \(SS_{regression}\), \(SS_{residual}\) and \(SS_{total}\) by hand and match them to anova()
Show that t² = F for a one-predictor model
Read the four regression diagnostics off the residuals
Export a labelled figure to PDF

Note

📤 What to turn in

10_correlation_linear_models.R
figures/lion_regression.pdf
This worksheet, with your written answers

Homework

Homework 5 — Correlation & Regression: The Transformation You Have to Find

One hundred stream invertebrate species, body mass against population density. You run Parts 5–8 of this worksheet again — fit, interpret, partition into sums of squares, check the residuals — in the same order and with the same functions.

The difference is that here the diagnostics do not confirm the model. On the raw data the regression says there is nothing there: a slope near zero, R² of 0.001, p = 0.73. The residuals say otherwise, and finding the transformation that fixes it is the assignment. Once you have it, there is a second question waiting: the slope has a predicted value from theory, so “is the slope different from zero” is not the interesting test.

Data: stream_invertebrates.csv — Whitlock & Schluter, The Analysis of Biological Data, Chapter 16, Figure 16.4-1. Note that the book ships these data already log-transformed; ours are the back-transformed raw values, so look the figure up after you have run your own diagnostics, not before.

Hand in: your 10_correlation_regression_homework.R script, your final figure as a PDF, and the worksheet with your written answers.

Other Formats

Worksheet: Correlation and Linear Regression

Part 1 · Setup and the data

Part 2 · Correlation

Part 3 · Correlation assumptions

Part 4 · When normality fails

Part 5 · Fitting a regression

Part 6 · Sums of squares

Part 7 · Regression assumptions

Part 8 · The figure you hand in

Review and checkpoint

Homework