Activity: Correlation vs Linear Models

Relationships

Hands-on activity: correlation, simple linear regression, diagnostic plots, and ANOVA for regression on the Nazca booby, lion nose-pigmentation, and prairie-biodiversity datasets.

Author

Bill Perry

Worksheet: Correlation and Linear Regression

How to use this worksheet

Work through each part in order, at your own pace. Type every line of code yourself into a plain R script — do not copy-paste. Blocks marked ▶ Run this are code you should type and execute. Blocks marked ✏️ Your turn ask you to write, modify, or answer something. Boxes marked 🚀 If you finish early are optional bonus material that goes a bit further than what we covered in lecture.

This worksheet demonstrates correlation and regression using three ecological examples from the textbook: Nazca booby aggression (correlation), lion nose pigmentation and age (regression), and prairie plant diversity and stability (a second regression to compare against).

Part 1 · Setup and Data

▶ Run this in your Script:

library(tidyverse)  # For data manipulation and visualization
library(patchwork)  # For combining plots
library(car)         # For regression diagnostics
library(broom)       # For tidy model output

set.seed(123)

▶ Run this — build the three example datasets:

# Lion data from Example 17.1
l_df <- tibble(
  prop_black = c(0.21, 0.14, 0.11, 0.13, 0.12, 0.13, 0.12, 0.18, 0.23, 0.22,
                0.20, 0.17, 0.15, 0.27, 0.26, 0.21, 0.30, 0.42, 0.43, 0.59,
                0.60, 0.72, 0.29, 0.10, 0.48, 0.44, 0.34, 0.37, 0.34, 0.74, 0.79, 0.51),
  age_yr = c(1.1, 1.5, 1.9, 2.2, 2.6, 3.2, 3.2, 2.9, 2.4, 2.1,
          1.9, 1.9, 1.9, 1.9, 2.8, 3.6, 4.3, 3.8, 4.2, 5.4,
          5.8, 6.0, 3.4, 4.0, 7.3, 7.3, 7.8, 7.1, 7.1, 13.1, 8.8, 5.4)
)

# Booby data from Example 16.1
b_df <- tibble(
  visits = c(1, 7, 15, 4, 11, 14, 23, 14, 9, 5, 4, 10,
            13, 13, 14, 12, 13, 9, 8, 18, 22, 22, 23, 31),
  aggression = c(-0.80, -0.92, -0.80, -0.46, -0.47, -0.46, -0.23, -0.16,
                 -0.23, -0.23, -0.16, -0.10, -0.10, 0.04, 0.13, 0.19,
                 0.25, 0.23, 0.15, 0.23, 0.31, 0.18, 0.17, 0.39)
)

# Prairie stability data from Example 17.3
p_df <- tibble(
  spp_n = rep(c(1, 2, 4, 8, 16), times = c(32, 32, 32, 32, 33)),
  log_stability = 1.20 + 0.033 * spp_n + rnorm(161, 0, 0.35)
)

Tip

Package overview

tidyverse — data science toolkit
patchwork — combine multiple ggplot2 plots easily
car — companion to Applied Regression (diagnostic tools)
broom — convert statistical objects into tidy data frames

Tip

🚀 If you finish early: Run glimpse() on all three data frames. Which columns are the predictor (X) and which are the response (Y) in each?

# Write your code here:

Part 2 · Correlation Analysis

Note

Correlation analysis: data types and assumptions

Data types required: both X and Y are continuous numerical, and both should be measured (not manipulated).

Assumptions for Pearson correlation: random sampling, bivariate normality, a linear relationship, and no extreme outliers.

Let’s start with the Nazca booby data.

▶ Run this:

cor(b_df$visits, b_df$aggression)

cor.test(b_df$visits, b_df$aggression)

▶ Run this — save the coefficient and get the variance explained:

booby_corr <- cor(b_df$visits, b_df$aggression)
booby_corr

booby_corr^2  # R-squared (variance explained)

▶ Run this — visualize it:

b_df %>%
  ggplot(aes(x = visits, y = aggression)) +
  geom_point(size = 3, alpha = 0.7)

Important

Activity: interpret the correlation

Based on the output above, answer:

Direction: Is the correlation positive or negative? What does this mean biologically?
Strength: How would you classify this correlation (weak, moderate, strong)?
Significance: Is the correlation statistically significant? What is the p-value?
Variance explained: What percentage of variance in adult aggression is explained by nestling visits?

1. ________________________________________________________
2. ________________________________________________________
3. ________________________________________________________
4. ________________________________________________________

▶ Run this — test the assumptions:

shapiro.test(b_df$visits)
shapiro.test(b_df$aggression)

p1 <- ggplot(b_df, aes(x = visits)) +
  geom_histogram(bins = 10, fill = "lightblue", color = "black") +
  labs(title = "Distribution of Visits", x = "Visits as Nestling", y = "Count") +
  theme_minimal()

p2 <- ggplot(b_df, aes(x = aggression)) +
  geom_histogram(bins = 10, fill = "lightgreen", color = "black") +
  labs(title = "Distribution of Aggression", x = "Future Aggression", y = "Count") +
  theme_minimal()

p3 <- ggplot(b_df, aes(sample = visits)) +
  stat_qq() + stat_qq_line() +
  labs(title = "Q-Q Plot: Visits", x = "Theoretical", y = "Sample") +
  theme_minimal()

p4 <- ggplot(b_df, aes(sample = aggression)) +
  stat_qq() + stat_qq_line() +
  labs(title = "Q-Q Plot: Aggression", x = "Theoretical", y = "Sample") +
  theme_minimal()

(p1 + p2) / (p3 + p4)

⚠️ If normality assumptions are violated (p < 0.05 in the Shapiro-Wilk test), consider Spearman’s rank correlation (non-parametric), a data transformation, or removing outliers (if justified).

▶ Run this — Spearman’s correlation:

cor.test(b_df$visits, b_df$aggression, method = "spearman")

✏️ Your turn: Do the Pearson and Spearman results agree on the direction and significance of the relationship? ________________________

Tip

🚀 If you finish early: Compute Kendall’s tau instead (method = "kendall"). How does its magnitude compare to Pearson’s r and Spearman’s rho for the same data?

# Write your code here:

Part 3 · Simple Linear Regression

Now let’s move from correlation to regression using the lion nose data.

Note

Linear regression: data types and assumptions

Data types required: X (predictor) and Y (response) are both continuous numerical; X can be fixed/controlled, Y is the outcome of interest.

Assumptions: linearity, independence, homoscedasticity, normality of residuals, and no influential outliers.

▶ Run this — fit the model:

lion_model <- lm(age_yr ~ prop_black, data = l_df)
summary(lion_model)

Important

Activity: interpret the regression output

Regression equation: age = ______________ + ______________ × prop_black
Slope interpretation: what does the slope value mean in biological terms?
R-squared: what percentage of variation in age is explained by nose blackness?
Significance: is the relationship statistically significant? How do you know?

1. ________________________________________________________
2. ________________________________________________________
3. ________________________________________________________
4. ________________________________________________________

▶ Run this — visualize the fitted line with its confidence band:

l_df %>%
  ggplot(aes(x = prop_black, y = age_yr)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, color = "red", fill = "pink", alpha = 0.3)

💡 Confidence interval = range for the mean age of all lions with that nose blackness. Prediction interval = range for an individual lion with that nose blackness. Prediction intervals are always wider.

Tip

🚀 If you finish early: Use predict(lion_model, newdata = tibble(prop_black = 0.5), interval = "confidence") and then interval = "prediction" for the same prop_black = 0.5. Confirm the prediction interval is wider.

# Write your code here:

Part 4 · Testing Regression Assumptions

▶ Run this — the four-panel diagnostic plot:

par(mfrow = c(2, 2))
plot(lion_model)
par(mfrow = c(1, 1))

Note

Understanding regression diagnostic plots

Residuals vs Fitted — look for random scatter around the horizontal line at 0; patterns indicate non-linearity or heteroscedasticity.
Q-Q Plot — look for points following the diagonal line; deviations indicate non-normal residuals.
Scale-Location — look for random scatter with a horizontal trend line; increasing spread indicates heteroscedasticity.
Residuals vs Leverage — look for points within Cook’s distance lines; points outside indicate influential observations.

▶ Run this — the formal tests:

shapiro_residuals <- shapiro.test(residuals(lion_model))
shapiro_residuals

library(lmtest)
bp_test <- bptest(lion_model)  # Breusch-Pagan test for homoscedasticity
bp_test

Important

Activity: assess assumption violations

Based on the diagnostic plots and tests:

Linearity: does the relationship appear linear? (Check Residuals vs Fitted)
Normality: are the residuals normally distributed? (Check Q-Q plot and Shapiro test)
Homoscedasticity: is the variance constant? (Check Scale-Location plot and BP test)
Influential points: are there any concerning influential observations?

1. ________________________________________________________
2. ________________________________________________________
3. ________________________________________________________
4. ________________________________________________________

Tip

🚀 If you finish early: bptest()’s null hypothesis is that variance is constant (homoscedastic). Write out, in your own words, what a significant Breusch-Pagan result would mean for this model — and whether that’s what you found.

________________________________________________________

Part 5 · ANOVA for Regression

▶ Run this:

anova_table <- anova(lion_model)
anova_table

▶ Run this — confirm the partition by hand:

ss_total <- sum((l_df$age_yr - mean(l_df$age_yr))^2)
ss_residual <- sum(residuals(lion_model)^2)
ss_regression <- ss_total - ss_residual

print("Manual calculation of sums of squares:")
print(paste("SS Total:", round(ss_total, 2)))
print(paste("SS Regression:", round(ss_regression, 2)))
print(paste("SS Residual:", round(ss_residual, 2)))
print(paste("SS Regression + SS Residual:", round(ss_regression + ss_residual, 2)))

✏️ Your turn: Does ss_regression + ss_residual match ss_total? Does ss_regression match the “prop_black” row’s Sum Sq in anova_table? ________________________

▶ Run this — visualize the variance components for one lion:

l_df$predicted <- predict(lion_model)
mean_age <- mean(l_df$age_yr)
example_point <- 10

variance_plot <- ggplot(l_df, aes(x = prop_black, y = age_yr)) +
  geom_point(size = 3, alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE, color = "blue", linewidth = 1) +
  geom_hline(yintercept = mean_age, linetype = "dashed", color = "darkgreen") +
  geom_segment(aes(x = prop_black[example_point],
                   y = age_yr[example_point],
                   xend = prop_black[example_point],
                   yend = predicted[example_point]),
               color = "red", linewidth = 1) +
  geom_segment(aes(x = prop_black[example_point],
                   y = predicted[example_point],
                   xend = prop_black[example_point],
                   yend = mean_age),
               color = "darkgreen", linewidth = 1) +
  annotate("text", x = 0.15, y = mean_age + 0.5, label = "Mean", color = "darkgreen") +
  annotate("text", x = l_df$prop_black[example_point] + 0.05,
           y = (l_df$age_yr[example_point] + l_df$predicted[example_point])/2,
           label = "Residual", color = "red") +
  annotate("text", x = l_df$prop_black[example_point] + 0.05,
           y = (l_df$predicted[example_point] + mean_age)/2,
           label = "Regression", color = "darkgreen") +
  labs(title = "Variance Components in Regression",
       subtitle = "Total variation = Regression + Residual",
       x = "Proportion Black", y = "Age (years)") +
  theme_minimal()
variance_plot

Tip

🚀 If you finish early: Change example_point to a different row number (try the lion with the largest residual — hint: which.max(abs(residuals(lion_model)))). How does the plot change?

# Write your code here:

Part 6 · Comparing Multiple Datasets

Let’s practice regression with the prairie biodiversity data.

▶ Run this:

prairie_model <- lm(log_stability ~ spp_n, data = p_df)
summary(prairie_model)

prairie_plot <- ggplot(p_df, aes(x = spp_n, y = log_stability)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm")
prairie_plot

Important

Activity: compare the two regressions

Compare the lion and prairie regression models:

Which model explains more variance? (Compare R² values)
Which has a stronger relationship? (Compare standardized slopes or correlation)
Which has more precise estimates? (Compare standard errors relative to estimates)

1. ________________________________________________________
2. ________________________________________________________
3. ________________________________________________________

Tip

🚀 If you finish early: Run plot(prairie_model) (four diagnostic plots, like Part 4). With n = 161, does non-normality of residuals worry you as much as it did for the lion model (n = 32)? Why might sample size change how much assumption violations matter?

# Write your code here:

Summary and Key Takeaways

Tip

What we learned today

Correlation vs. regression — correlation measures association; regression predicts one variable from another.
Assumptions matter — always check assumptions before interpreting results, and use appropriate alternatives when they’re violated.
Interpretation — R² tells us the proportion of variance explained; slopes tell us the rate of change; p-values tell us if relationships are statistically significant.
Practical considerations — correlation ≠ causation, outliers can have major impacts, and sample size affects the power to detect relationships.

Warning

Common mistakes to avoid

Using correlation when you mean regression (or vice versa)
Ignoring assumption violations
Extrapolating beyond the range of the data
Confusing confidence and prediction intervals
Over-interpreting R² values
Forgetting about biological significance vs. statistical significance

Review and checkpoint

At this point you can:

Calculate and interpret a Pearson correlation coefficient and its significance test
Test correlation assumptions and switch to Spearman’s when they fail
Fit a simple linear regression and interpret slope, intercept, and R²
Read the four standard regression diagnostic plots
Explain how ANOVA partitions variance into regression and residual components
Compare R² and precision across two different regression models

Note

📤 What to turn in before next class

Upload both of these to the course management system:

Your code — the scripts/ folder (or just 10_correlation_linear_models.R)
This worksheet, with your written answers

Additional Resources

Whitlock & Schluter, Ch. 16 (Correlation)
Whitlock & Schluter, Ch. 17 (Regression)
Gotelli & Ellison, A Primer of Ecological Statistics, Ch. 9 (Regression)
R for Data Science
Quick-R: Regression

Getting unstuck

When code breaks — and it will, that is normal:

Read the error message out loud. R usually names the line and the problem.
Check the usual suspects: did you run library(tidyverse), library(car), library(broom), and library(lmtest)? Spelling? A missing ) or %>% at the start of a line?
?function_name opens the built-in help page.
Bring the exact error (copy-paste it) to class or office hours.

💡 Key idea: Every working scientist googles error messages daily. Getting stuck is not failing — it is the job.

Other Formats

Worksheet: Correlation and Linear Regression

How to use this worksheet

Part 1 · Setup and Data

Part 2 · Correlation Analysis

Part 3 · Simple Linear Regression

Part 4 · Testing Regression Assumptions

Part 5 · ANOVA for Regression

Part 6 · Comparing Multiple Datasets

Summary and Key Takeaways

Review and checkpoint

Additional Resources

Getting unstuck