Activity: Multiple Regression

Several predictors

Hands-on activity: multicollinearity, variable transformations, model diagnostics, model comparison, and predictor importance on the Michaletz et al. (2014) global forest NPP dataset.

Author

Bill Perry

Worksheet: Multiple Regression — Forest Net Primary Production

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 examines the relationship between Net Primary Production (npp) and various climate and forest characteristics across global forest sites, based on Michaletz et al. (2014). We’ll explore multicollinearity, model selection, and variable transformations — the real-data, many-predictor counterpart to the lecture’s two-predictor ant example.

Part 1 · Setup and Data

▶ Run this in your Script:

library(tidyverse)    # For data manipulation and visualization
library(car)           # For regression diagnostics (VIF, etc.)
library(corrplot)      # For correlation plots
library(GGally)        # For pairs plots
library(broom)         # For tidy model outputs

forest_df <- read_csv("data/michaletz_etal_2014_clean.csv")
head(forest_df)

✏️ Your turn: Run glimpse(forest_df). Which columns look continuous, and which look categorical? ________________________

Part 2 · Exploring Multicollinearity

▶ Run this — a correlation matrix, and a visual version of it:

num_vars <- forest_df %>%
  select_if(is.numeric)

cor_matrix <- cor(num_vars, use = "complete.obs")
cor_matrix

corrplot(cor_matrix, method = "color", type = "upper",
         addCoef.col = "grey45", tl.cex = 0.8, number.cex = 0.7)

▶ Run this — a pairs plot for a fuller picture:

forest_df %>%
  select(-leaf) %>%  # Exclude categorical variable for pairs plot
  ggpairs(
    upper = list(continuous = wrap("cor", size = 5)),
    lower = list(continuous = wrap("points", alpha = 0.6, size = 0.8))
  ) +
  theme_minimal()

✏️ Your turn: From the correlation matrix, which two predictors look most strongly correlated with each other? ________________________

▶ Run this — fit the full model with every predictor:

model_init <- lm(npp ~ age + biomass + season + temp +
             precip + teb + leaf, data = forest_df)
summary(model_init)

▶ Run this — check Variance Inflation Factors:

vif_values <- vif(model_init)
vif_values

vif_df <- data.frame(
  Variable = names(vif_values),
  VIF = as.numeric(vif_values)
) %>%
  arrange(desc(VIF))

ggplot(vif_df, aes(x = reorder(Variable, VIF), y = VIF)) +
  geom_col(fill = "darkblue", alpha = 0.7) +
  geom_hline(yintercept = 10, color = "red", linetype = "dashed", linewidth = 1) +
  geom_hline(yintercept = 5, color = "orange", linetype = "dashed", linewidth = 1) +
  coord_flip() +
  labs(
    title = "Variance Inflation Factors",
    subtitle = "Red line: VIF = 10 (serious concern), Orange line: VIF = 5 (moderate concern)",
    x = "Variables", y = "VIF"
  ) +
  theme_minimal()

✏️ Your turn: Which predictor(s) have VIF > 10? Does that match what you flagged from the correlation matrix? ________________________

Tip

🚀 If you finish early: Recompute VIF using only 3 predictors of your choice (drop the rest from the model formula). Do the remaining VIFs change? Why would dropping a correlated predictor lower the VIFs of the ones that are left?

# Write your code here:

Part 3 · Addressing Multicollinearity

Season and temperature are highly correlated. Let’s remove season.

▶ Run this:

model_2 <- lm(npp ~ age + biomass + temp + precip + teb + leaf,
             data = forest_df)
summary(model_2)

vif_values2 <- vif(model_2)
vif_values2

✏️ Your turn: Did removing season fix the VIF problem? ________________________

Part 4 · Checking the Shape of Relationships

▶ Run this — Added Variable (partial regression) plots:

par(mfrow = c(2, 3))
avPlots(model_2, main = "Partial Regression Plots")
par(mfrow = c(1, 1))

The original analysis found age showed a curvy relationship. Let’s try a log transformation. (Really, you should try transforming the response variable first — we’re doing age here to see the AV-plot difference clearly.)

▶ Run this:

forest_df <- forest_df %>%
  mutate(log_age = log10(age))

model_3 <- lm(npp ~ log_age + biomass + temp + precip + teb + leaf,
             data = forest_df)
summary(model_3)

par(mfrow = c(2, 3))
avPlots(model_3, main = "Partial Regression Plots (Log age)")
par(mfrow = c(1, 1))

✏️ Your turn: Compare the AV plot for age before and after the log transform. Does the relationship look straighter now? ________________________

Part 5 · Model Diagnostics

▶ Run this — the four standard diagnostic plots:

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

▶ Run this — a closer look at residuals vs. fitted:

residuals_data <- data.frame(
  Fitted = fitted(model_3),
  Residuals = residuals(model_3),
  Standardized_Residuals = rstandard(model_3)
)

ggplot(residuals_data, aes(x = Fitted, y = Residuals)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", color = "red") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "Residuals vs Fitted Values", x = "Fitted Values", y = "Residuals") +
  theme_minimal()

Let’s also try a cube-root transformation of the response, npp:

▶ Run this:

forest_df <- forest_df %>%
  mutate(npp_cuberoot = npp^(1/3))

model_4 <- lm(npp_cuberoot ~ log_age + biomass + temp + precip +
             teb + leaf, data = forest_df)
summary(model_4)

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

✏️ Your turn: Do the residual plots for model_4 look better than for model_3? What specifically improved (or didn’t)? ________________________

Tip

🚀 If you finish early: Try a square-root transform of npp instead of cube-root (sqrt(npp)). Refit the model and compare its residual plots to model_4’s.

# Write your code here:

Part 6 · Model Simplification and Comparison

▶ Run this — drop the non-significant precip term:

model_5 <- lm(npp_cuberoot ~ log_age + biomass + temp + teb + leaf,
             data = forest_df)
summary(model_5)

▶ Run this — compare model_4 and model_5 by AIC:

model_comparison <- data.frame(
  Model = c("Model 4 (Full)", "Model 5 (No precip)"),
  AIC = c(AIC(model_4), AIC(model_5)),
  R_squared = c(summary(model_4)$r.squared, summary(model_5)$r.squared),
  Adj_R_squared = c(summary(model_4)$adj.r.squared, summary(model_5)$adj.r.squared)
)
model_comparison

summary(model_5, conf.int = TRUE)

✏️ Your turn: Does model_5 have a lower AIC than model_4? Given the “penalize complexity” logic from lecture, what does that tell you about dropping precip? ________________________

Part 7 · Predictor Importance

▶ Run this — partial R² for every predictor at once, using sensemakr:

library(sensemakr)
partial_r2_sensemakr <- partial_r2(model_5)
partial_r2_sensemakr

▶ Run this — the same idea, computed by hand from a Type III ANOVA table:

anova_type3 <- Anova(model_5, type = "III")
anova_type3

f_stats <- anova_type3$`F value`[!is.na(anova_type3$`F value`)]
df_num <- anova_type3$Df[!is.na(anova_type3$`F value`)]
df_den <- anova_type3$Df[nrow(anova_type3)]  # Residual df

partial_r2_from_f <- f_stats * df_num / (f_stats * df_num + df_den)

results_table <- data.frame(
  Variable = rownames(anova_type3)[!is.na(anova_type3$`F value`)],
  F_statistic = f_stats,
  p_value = anova_type3$`Pr(>F)`[!is.na(anova_type3$`F value`)],
  Partial_R_squared = partial_r2_from_f
)
results_table

✏️ Your turn: Do the sensemakr partial R² values roughly match your by-hand Partial_R_squared column? ________________________

Tip

🚀 If you finish early: Which single predictor has the highest partial R² in model_5? Does that match which predictor had the smallest p-value in summary(model_5)?

# Write your code here:

Part 8 · An Alternative: Standardized Variables

▶ Run this — refit with every numeric predictor standardized (mean 0, SD 1), including a temp × precip interaction:

forest_standardized <- forest_df %>%
  mutate(
    npp_sqrt_scaled = scale(sqrt(npp))[,1],
    log_age_scaled = scale(log10(age))[,1],
    biomass_scaled = scale(biomass)[,1],
    temp_scaled = scale(temp)[,1],
    precip_scaled = scale(precip)[,1],
    teb_scaled = scale(teb)[,1]
  )

model_std <- lm(npp_sqrt_scaled ~ log_age_scaled + biomass_scaled +
                temp_scaled * precip_scaled + teb_scaled,
                data = forest_standardized)
summary(model_std)

✏️ Your turn: In model_std, which predictor has the largest absolute standardized coefficient? Is that the same predictor that had the highest partial R² in Part 7? ________________________

Key Findings

Multicollinearity — growing season length and temperature were highly correlated; removed season to address it.
Variable transformations — log transformation of age improved model fit; cube-root transformation of npp addressed assumption violations.
Final model results — significant predictors: age (negative), biomass (positive), temp (positive); teb had a negative effect; leaf type differences were significant.
Biological interpretation — younger stands had higher npp (for a given biomass); higher biomass was associated with higher npp; temperature was positively related to npp; coniferous forests had lower npp than broadleaf forests.

Note

Reference

Michaletz, S.T., Cheng, D., Kerkhoff, A.J. & Enquist, B.J. (2014). Convergence of terrestrial plant production across global climate gradients. Nature, 512, 39–43.

Review and checkpoint

At this point you can:

Detect multicollinearity with a correlation matrix and VIFs, and address it by dropping a redundant predictor
Use Added Variable plots to check whether a predictor’s relationship with Y is linear
Apply and evaluate predictor and response transformations
Read the standard four-panel regression diagnostic plots
Compare nested models with AIC and adjusted R²
Quantify predictor importance with partial R² and standardized coefficients

Note

📤 What to turn in before next class

Upload both of these to the course management system:

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

Part 9 · Take-Home Extension — Nutrients and Chlorophyll in New Zealand Lakes

Due Monday, October 19 — before the Analysis of Variance class.

Same workflow as Parts 1–6 — correlation/multicollinearity check, fit the model, diagnose, transform if needed, interpret — new dataset, new question. This time you decide what (if anything) needs fixing before you trust the model; nobody walks you through it step by step.

Background

Abell et al. (2010) studied nutrient limitation in New Zealand lakes. Chlorophyll-a is a standard proxy for algal biomass, and both phosphorus and nitrogen can drive it up. Below are measurements from 111 New Zealand lakes.

▶ Run this:

nz_df <- read_csv("data/abel_et_al_lakes_data.csv")
head(nz_df)

Your Task

Question: How do total phosphorus and total nitrogen jointly predict chlorophyll-a concentration across these lakes?

✏️ Your turn: State your hypotheses. Check for multicollinearity between the two predictors and decide whether it’s a problem here. Look at the shape of each relationship and decide whether a transformation is needed before you trust a linear model (hint: nutrient–chlorophyll relationships are rarely linear on the raw scale). Justify whatever choices you make, then fit the model, interpret it, and report the result properly.

H0 =
________________________________________________________
Ha =
________________________________________________________
Model/transformation I will use and why:
________________________________________________________
________________________________________________________

# Write your code here:

Interpretation:
________________________________________________________
________________________________________________________

Final Figure

Produce one publication-quality figure showing the relationship between chlorophyll-a and the predictor(s) that mattered most — proper axis labels, no default ggplot grey background — and export it with ggsave().

# Write your code here:

Note

📤 What to turn in — due Oct 19

Your write-up should include:

The question
Your hypotheses
Which model/transformation you used, and your justification for choosing it
The code and its output
Your interpretation of the result
A final figure showing the result

Note

Reference

Abell, J.M., Özkundakci, D. & Hamilton, D.P. (2010). Nitrogen and phosphorus limitation of phytoplankton growth in New Zealand lakes. Ecosystems, 13, 966–977.

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(corrplot), library(GGally), and library(sensemakr)? 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: Multiple Regression — Forest Net Primary Production

How to use this worksheet

Part 1 · Setup and Data

Part 2 · Exploring Multicollinearity

Part 3 · Addressing Multicollinearity

Part 4 · Checking the Shape of Relationships

Part 5 · Model Diagnostics

Part 6 · Model Simplification and Comparison

Part 7 · Predictor Importance

Part 8 · An Alternative: Standardized Variables

Key Findings

Review and checkpoint

Part 9 · Take-Home Extension — Nutrients and Chlorophyll in New Zealand Lakes

Background

Your Task

Final Figure

Getting unstuck