Common Code 11 — Correlation and Simple Linear Regression
Penguin flipper length and body mass — from r to a fitted line
How to do a linear regreession in R and interpret the output and test assumptions
Correlation and linear regression
These two analyses are closely related but answer different questions:
- Correlation — how strongly are two variables associated? (no cause assumed)
- Linear regression — how well does one variable predict the other? (X causes or precedes Y)
We use penguin flipper length and body mass throughout. Both are measured variables with no manipulation — making correlation appropriate first, and regression appropriate when we want a predictive equation.
⬇️ Download the companion R script:
10_correlation_regression.R
Packages and data
library(tidyverse)
library(palmerpenguins)
library(skimr)
library(broom) # tidy(), glance() — clean model output
source("themes/r_themes_for_3_sizes.R")
# Complete cases only
pg <- penguins |> drop_na(flipper_length_mm, body_mass_g, species)Before the improved version, two issues worth flagging:
In the DEET regression file: The bites column in chap17q30DEETMosquiteBites.csv contains values like 4.062, 3.535, 2.549 — these are already square-root transformed bite counts (√16.5, √12.5, √6.5). The original script treats them as raw counts, fits a linear model, then writes a methods section saying “number of mosquito bites” — but the response variable is actually √bites. The regression equation, R², intercept interpretation and results text are all on the wrong scale. Always check skim() or head() and compare to the paper before fitting any model.
In both correlation and regression files: n() is used inside summarise() for sample sizes rather than sum(!is.na()), and the results sections contain unfilled [value] placeholders rather than extracted model values.
In both files: par(mfrow = c(2,2)); plot(model) uses base R diagnostics that cannot be saved with ggsave() and do not match course style. Replaced throughout with ggplot versions built from fitted(), residuals(), and rstandard().
PART 1 · Correlation
1 · What correlation measures
The Pearson correlation coefficient r quantifies the strength and direction of a linear association between two continuous variables. It ranges from −1 to +1.
\[H_0: \rho = 0 \qquad \text{(no linear association in the population)}\] \[H_A: \rho \neq 0 \qquad \text{(there is a linear association)}\]
| r | Interpretation |
|---|---|
| 0.00 – ±0.30 | Weak |
| ±0.30 – ±0.70 | Moderate |
| ±0.70 – ±1.00 | Strong |
Correlation measures association only. Neither variable is designated as cause or effect. Both flipper length and body mass are responses to the same underlying biology (species, age, sex, food availability). Correlation is the right tool when neither variable was experimentally manipulated.
2 · Explore the data first
Plot the relationship before computing anything. The LOESS curve tells you whether a straight line is a reasonable description.
ggplot(pg, aes(x = flipper_length_mm, y = body_mass_g)) +
geom_point(alpha = 0.5, color = "steelblue", size = 2) +
geom_smooth(method = "loess", se = FALSE,
color = "tomato", linetype = "dashed", linewidth = 0.8) +
labs(x = "Flipper length (mm)", y = "Body mass (g)",
title = "Does the relationship look linear?") +
theme_regular()- Does the LOESS curve follow a roughly straight path? → linear relationship is plausible
- Is there a clear positive or negative direction?
- Are there unusual clusters or outliers?
- Does the spread of points look similar across the x-axis? (similar spread = equal variance)
If the curve bends strongly, consider a log transformation before correlating.
3 · Check assumptions
Pearson correlation requires both variables to be approximately normally distributed. Check each one separately with a Q-Q plot and Shapiro-Wilk test.
# Q-Q plot for flipper length
ggplot(pg, aes(sample = flipper_length_mm)) +
stat_qq(color = "steelblue", alpha = 0.6) +
stat_qq_line(color = "tomato") +
labs(title = "Q-Q: Flipper length",
x = "Theoretical quantiles", y = "Sample quantiles") +
theme_regular()
# Q-Q plot for body mass
ggplot(pg, aes(sample = body_mass_g)) +
stat_qq(color = "steelblue", alpha = 0.6) +
stat_qq_line(color = "tomato") +
labs(title = "Q-Q: Body mass",
x = "Theoretical quantiles", y = "Sample quantiles") +
theme_regular()
# Shapiro-Wilk test for each
shapiro.test(pg$flipper_length_mm)
shapiro.test(pg$body_mass_g)Spearman’s rank correlation is the non-parametric alternative to Pearson. It tests whether the ranks of the two variables are associated — suitable for non-normal data or when outliers are present.
cor.test(pg$flipper_length_mm, pg$body_mass_g, method = "spearman")Report as ρ (rho) rather than r.
4 · Run the correlation
cor_result <- cor.test(pg$flipper_length_mm, pg$body_mass_g,
method = "pearson")
cor_result
# Clean single-row output with broom
tidy(cor_result)t = 32.7, df = 340, p-value < 2.2e-16
95% CI: [0.843, 0.886]
cor = 0.871- cor — the Pearson r value
- t — the test statistic (r converted to a t-distribution)
- df — n − 2
- p-value — probability of r this large if ρ = 0 in the population
- 95% CI — plausible range for the true population correlation; if it contains 0 the correlation is not significant
r = 0.871 is a strong positive association — larger flippers go with heavier penguins. r² = 0.759, meaning flipper length accounts for about 76% of the variance in body mass across all species combined.
Extract values for reporting
r_val <- round(cor_result$estimate, 3)
r2_val <- round(cor_result$estimate^2, 3)
p_val <- round(cor_result$p.value, 4)Correlation by species
pg |>
group_by(species) |>
summarise(
n = sum(!is.na(flipper_length_mm) & !is.na(body_mass_g)),
r = round(cor(flipper_length_mm, body_mass_g, use = "complete.obs"), 3),
.groups = "drop"
)The overall r across all species (≈ 0.87) is much higher than the r within any individual species (≈ 0.58–0.73). This is Simpson’s paradox — species differences inflate the overall correlation. Always check correlations within groups as well as overall, especially in ecological data.
5 · Publication correlation plot
eq_label <- paste0("r = ", r_val, "\np = ", p_val)
cor_plot <- ggplot(pg, aes(x = flipper_length_mm, y = body_mass_g,
color = species)) +
geom_point(alpha = 0.6, size = 2) +
stat_ellipse(linewidth = 0.8, alpha = 0.7) +
annotate("text", x = 175, y = 6000,
label = eq_label, hjust = 0, size = 3.5, color = "grey30") +
labs(x = "Flipper length (mm)", y = "Body mass (g)",
color = "Species",
title = "Flipper length vs body mass — Palmer penguins") +
theme_regular()
cor_plot
ggsave("figures/penguin_correlation.pdf",
plot = cor_plot, width = 6, height = 5, units = "in")stat_ellipse() draws a 95% confidence ellipse for the bivariate distribution. A narrow, tilted ellipse indicates strong correlation; a nearly circular ellipse indicates weak correlation.
Figure 1. Relationship between flipper length (mm) and body mass (g) for three species of Palmer Archipelago penguins (n = 342). Ellipses show 95% confidence regions. There was a significant positive correlation between flipper length and body mass across all species (Pearson r = 0.871, 95% CI: [0.843, 0.886], p < 0.001).
PART 2 · Simple Linear Regression
6 · Correlation vs regression
| Correlation | Regression | |
|---|---|---|
| Goal | Measure association | Predict Y from X |
| Variables | Neither is response/predictor | X predicts Y |
| Output | r and p-value | Equation + R² + predictions |
| Direction | Symmetric (r of X,Y = r of Y,X) | Asymmetric (X causes/predicts Y) |
Use regression when one variable is logically the predictor — body size predicts metabolic rate, temperature predicts growth rate, DEET dose predicts bites. Use correlation when both are measured responses with no causal direction.
The regression model:
\[\text{body mass} = \alpha + \beta \times \text{flipper length} + \varepsilon\]
- α (intercept) — predicted body mass when flipper length = 0 (often biologically meaningless)
- β (slope) — change in body mass (g) for each 1 mm increase in flipper length
- ε — residual error
7 · Fit the model
reg_model <- lm(body_mass_g ~ flipper_length_mm, data = pg)
# Full summary
summary(reg_model)
# Clean coefficients table with 95% CI
tidy(reg_model, conf.int = TRUE)
# Model-level fit statistics
glance(reg_model)Coefficients:
Estimate Std.Error t value Pr(>|t|)
(Intercept) -5780.8 288.1 -20.1 < 2e-16 ***
flipper_length_mm 49.7 1.5 33.0 < 2e-16 ***
R² = 0.759 F(1,340) = 1070 p < 2e-16- Intercept (−5780.8) — predicted body mass at flipper length = 0; biologically meaningless here (no penguin has 0 mm flippers) but mathematically required
- Slope (49.7) — each additional 1 mm of flipper length is associated with 49.7 g more body mass
- R² — the model explains 75.9% of the variance in body mass
- F-statistic — overall test that the slope ≠ 0 (same result as the t-test on the slope in a simple regression)
Equation: body_mass_g = −5780.8 + 49.7 × flipper_length_mm
Extract key values cleanly
intercept <- round(coef(reg_model)[1], 1)
slope <- round(coef(reg_model)[2], 1)
r2 <- round(summary(reg_model)$r.squared, 3)broom::tidy() and broom::glance()
tidy() gives you coefficients as a tibble. glance() gives you R², F, df, AIC, and other model-level statistics — all in one row, easy to extract. Both work on lm, aov, t.test, and many other model types.
8 · Check regression assumptions
The four ANOVA assumptions apply here too, but the focus shifts slightly — you check the residuals, not the raw variables.
| Assumption | Check | What you want |
|---|---|---|
| Linearity | Residuals vs Fitted | Random scatter around 0, no curves |
| Equal variance | Residuals vs Fitted | Similar spread left to right |
| Normality of residuals | Q-Q plot + Shapiro-Wilk | Points on the line; p > 0.05 |
| Independence | Study design | Not testable statistically |
diag_df <- tibble(
fitted = fitted(reg_model),
residuals = residuals(reg_model),
std_resid = rstandard(reg_model),
cooks_d = cooks.distance(reg_model)
)Residuals vs Fitted
ggplot(diag_df, aes(x = fitted, y = residuals)) +
geom_point(alpha = 0.5, color = "steelblue") +
geom_hline(yintercept = 0, linetype = "dashed", color = "tomato") +
geom_smooth(method = "loess", se = FALSE,
color = "grey40", linewidth = 0.8) +
labs(title = "Residuals vs Fitted",
x = "Fitted values (g)", y = "Residuals") +
theme_regular()Q-Q plot of residuals
ggplot(diag_df, aes(sample = std_resid)) +
stat_qq(alpha = 0.5, color = "steelblue") +
stat_qq_line(color = "tomato") +
labs(title = "Normal Q-Q of residuals",
x = "Theoretical quantiles", y = "Standardised residuals") +
theme_regular()Cook’s distance — influential observations
ggplot(diag_df, aes(x = seq_along(cooks_d), y = cooks_d)) +
geom_col(fill = "steelblue", alpha = 0.7) +
geom_hline(yintercept = 4 / nrow(pg),
linetype = "dashed", color = "tomato") +
labs(title = "Cook's distance",
x = "Observation", y = "Cook's D",
caption = "Dashed line = 4/n threshold") +
theme_regular()Cook’s D measures how much each observation influences the fitted regression line. Points above the 4/n threshold are potentially influential — worth examining, though not automatically removed. Check whether they are data entry errors or genuinely extreme biology.
# Formal normality test on residuals
shapiro.test(residuals(reg_model))9 · Predictions
new_data <- tibble(flipper_length_mm = c(180, 200, 210, 220))
# Confidence interval — uncertainty in the MEAN response at that X
predict(reg_model, new_data, interval = "confidence") |> round(1)
# Prediction interval — uncertainty for a SINGLE new observation
predict(reg_model, new_data, interval = "prediction") |> round(1)Both widen as X moves away from its mean, but they are different questions:
- Confidence interval — where does the average body mass fall for all penguins with flipper length = 200 mm? (narrower)
- Prediction interval — where would a single new penguin’s body mass fall if its flipper length = 200 mm? (wider — individual variation adds uncertainty)
Always use prediction intervals when making a forecast about one individual. Use confidence intervals when describing where the population mean lies.
The equation is only valid within the range of flipper lengths in your sample (~170–230 mm). Plugging in flipper_length_mm = 100 or 300 produces a number, but that number is meaningless — you have no data to support the model at those values.
10 · Publication regression plot
eq_label <- paste0("y = ", intercept, " + ", slope, "x\n",
"R² = ", r2)
reg_plot <- ggplot(pg, aes(x = flipper_length_mm, y = body_mass_g)) +
geom_point(aes(color = species), alpha = 0.5, size = 2) +
geom_smooth(method = "lm", color = "grey20",
fill = "grey70", alpha = 0.2, linewidth = 1) +
annotate("text", x = 172, y = 6000,
label = eq_label, hjust = 0, size = 3.5, color = "grey20") +
labs(x = "Flipper length (mm)", y = "Body mass (g)",
color = "Species",
title = "Body mass predicted by flipper length") +
theme_regular()
reg_plot
ggsave("figures/penguin_regression.pdf",
plot = reg_plot, width = 6, height = 5, units = "in")11 · How to report the results
“There was a significant positive correlation between flipper length and body mass in Palmer Archipelago penguins (Pearson r = 0.871, 95% CI: [0.843, 0.886], n = 342, p < 0.001). Flipper length accounted for 75.9% of the variance in body mass (r² = 0.759).”
“Flipper length was a significant positive predictor of body mass in Palmer Archipelago penguins (simple linear regression: F₁,₃₄₀ = 1070, p < 0.001, R² = 0.759). For each additional millimetre of flipper length, body mass increased by 49.7 g (95% CI: [46.8, 52.6] g mm⁻¹). The fitted equation was: body mass (g) = −5780.8 + 49.7 × flipper length (mm).”
A high R² (even 0.99) does not guarantee the model assumptions are met. Always check the residual plots — a curved residual pattern with R² = 0.95 still means you have fitted the wrong model. Check assumptions first, R² second.
Quick reference
| Task | Code |
|---|---|
| Pearson correlation | cor.test(x, y, method = "pearson") |
| Spearman correlation | cor.test(x, y, method = "spearman") |
| Tidy correlation output | broom::tidy(cor_result) |
| Grouped correlation | group_by() \|> summarise(r = cor(x, y, use = "complete.obs")) |
| Ellipse on scatter | stat_ellipse() |
| Fit regression | lm(y ~ x, data = df) |
| Coefficients + CI | tidy(model, conf.int = TRUE) |
| R², F, AIC | glance(model) |
| Extract residuals | tibble(fitted = fitted(m), resid = residuals(m), std_resid = rstandard(m)) |
| Cook’s distance | cooks.distance(model) |
| Residuals vs Fitted | ggplot on diag_df with geom_point() + geom_hline(y=0) |
| Normality of residuals | shapiro.test(residuals(model)) |
| Confidence interval | predict(model, new_data, interval = "confidence") |
| Prediction interval | predict(model, new_data, interval = "prediction") |
| Regression line on plot | geom_smooth(method = "lm", se = TRUE) |
End of Common Code 10 — Correlation and Simple Linear Regression. Next: Common Code 11 — One-Way ANOVA.