########################################################################
# PLS 206 — Week 2: Linear regression
# Grey Monroe, Fall 2026
#
# Every block is tagged:
#   [CORE]   the analysis: code you would run on your own data
#   [HANDY]  useful, not essential: wrangling, plotting, small helpers
#   [DEMO]   for the lecture only: simulated data, sums by hand, and the
#            code behind slide figures
########################################################################

# [CORE]
# install.packages(c("ggplot2", "car", "scatterplot3d"))
library(ggplot2)

########################################################################
# PART A. Fitting a line
########################################################################

# ---- 1) The data: 40 wheat plants ------------------------------------

# [CORE]
wheat <- read.csv("https://greymonroe.github.io/PLS_206/data/wheat_plants.csv")
head(wheat)
nrow(wheat)

# Always plot first.

# [CORE]
plot(wheat$height, wheat$biomass, xlab = "Height (cm)", ylab = "Biomass (g)")

ggplot(wheat, aes(height, biomass)) +
  geom_point() +
  labs(x = "Height (cm)", y = "Biomass (g)") +
  theme_classic()

# ---- 2) Scoring a line: residuals and the sum of squares -------------

# Any line through the middle of the data: biomass = b0 + slope *
# height. Score each candidate line by its residual sum of squares
# (RSS). lm() does all of this for you in section 4.

# [DEMO: a function that scores a line by its RSS]
# score a line through the middle of the data by its RSS
model_line <- function(slope) {
  b0 <- mean(wheat$biomass) - slope * mean(wheat$height)
  predicted <- b0 + slope * wheat$height
  sum((wheat$biomass - predicted)^2)
}
model_line(0.12)    # line A
model_line(0.279)   # line B
model_line(0.45)    # line C

# [DEMO: the three candidate lines]
# the three candidate lines on the plot (all through the middle of the data)
plot(wheat$height, wheat$biomass, xlab = "Height (cm)", ylab = "Biomass (g)")
for (slope in c(0.12, 0.279, 0.45)) {
  b0 <- mean(wheat$biomass) - slope * mean(wheat$height)
  abline(b0, slope, col = c("blue", "red", "darkgreen")[match(slope, c(0.12, 0.279, 0.45))], lwd = 2)
}

# [DEMO: residuals from line B, and why we square them]
# residuals: the vertical distance from each point to line B
b0 <- mean(wheat$biomass) - 0.279 * mean(wheat$height)
predicted <- b0 + 0.279 * wheat$height
residual  <- wheat$biomass - predicted
plot(wheat$height, wheat$biomass, xlab = "Height (cm)", ylab = "Biomass (g)")
abline(b0, 0.279, col = "red", lwd = 2)
segments(wheat$height, wheat$biomass, wheat$height, predicted, col = "blue", lty = 2)
sum(residual)       # they cancel out...
sum(residual^2)     # ...so we square them: the residual sum of squares

# [DEMO: the RSS drawn as squares]
# squares whose side is the residual: the RSS is their total area
plot(wheat$height, wheat$biomass, xlab = "Height (cm)", ylab = "Biomass (g)", asp = 1)
abline(b0, 0.279, col = "red", lwd = 2)
rect(wheat$height, pmin(wheat$biomass, predicted), wheat$height + abs(residual),
     pmax(wheat$biomass, predicted), col = adjustcolor("blue", 0.2), border = "blue")

# [DEMO: the RSS for every slope is a bowl]
# the RSS for every slope from 0 to 0.6: a bowl with its bottom at 0.279
slopes <- seq(0, 0.6, by = 0.005)
rss <- sapply(slopes, model_line)
plot(slopes, rss, type = "l", xlab = "Slope (g per cm)", ylab = "Residual sum of squares")
slopes[which.min(rss)]

# [DEMO: both knobs at once (intercept and slope)]
# both knobs at once: the RSS over intercept and slope is a bowl
rss2 <- function(b0, b1) sum((wheat$biomass - (b0 + b1 * wheat$height))^2)
b0s <- seq(-20, 10, length.out = 60); b1s <- seq(0.1, 0.45, length.out = 60)
bowl <- outer(b0s, b1s, Vectorize(rss2))
persp(b0s, b1s, bowl, theta = 35, phi = 25, xlab = "intercept", ylab = "slope", zlab = "RSS", col = "lightblue")
# the bottom of the bowl:
c(b0s[which(bowl == min(bowl), arr.ind = TRUE)[1]], b1s[which(bowl == min(bowl), arr.ind = TRUE)[2]])

# ---- 3) The least-squares slope and intercept by hand ----------------

# [DEMO: slope and intercept from cov(), var(), and cor()]
# slope
b1 <- cov(wheat$height, wheat$biomass) / var(wheat$height)
b1
# intercept
b0 <- mean(wheat$biomass) - b1 * mean(wheat$height)
b0
# the same slope, via r
r <- cor(wheat$height, wheat$biomass)
r * sd(wheat$biomass) / sd(wheat$height)

# ---- 4) lm(): the same line ------------------------------------------

# [CORE]
# lm(response ~ predictor, data = data.frame)
model <- lm(biomass ~ height, data = wheat)
model
coef(model)

plot(wheat$height, wheat$biomass, xlab = "Height (cm)", ylab = "Biomass (g)")
abline(model, col = "red", lwd = 2)

# ---- 5) Predictions --------------------------------------------------

# [CORE]
predict(model, newdata = data.frame(height = 100))
-5.65 + 0.279 * 100     # the same thing from the equation
new <- data.frame(height = c(60, 80, 100))
predict(model, newdata = new)

# [CORE]
# fitted values and residuals for our own 40 plants
wheat$predicted <- predict(model)   # fitted values
wheat$residual  <- resid(model)
head(wheat[, c("height", "biomass", "predicted", "residual")])

# [DEMO: observed = predicted + residual]
# observed = predicted + residual
head(wheat$predicted + wheat$residual)
head(wheat$biomass)
sum(wheat$residual)        # residuals sum to zero (up to rounding)

# [CORE]
# interpolate, don't extrapolate: predict() will answer anything
predict(model, newdata = data.frame(height = c(0, 200)))   # nonsense plants
range(wheat$height)                                        # where the data are

# ---- 5b) The line is a conditional mean ------------------------------

# At each height, biomass is a distribution: mean on the line, spread
# sigma.

# [CORE]
sigma(model)                        # residual standard error: the typical miss

# [DEMO: the normal curve of biomass at three heights]
mu <- predict(model, data.frame(height = c(60, 80, 100)))
cbind(height = c(60, 80, 100), mean = mu, low = mu - 2 * sigma(model), high = mu + 2 * sigma(model))
# the model's assumption, drawn: a normal curve of biomass at each of three heights
plot(wheat$height, wheat$biomass, xlab = "Height (cm)", ylab = "Biomass (g)", xlim = c(45, 115))
abline(model, col = "red", lwd = 2)
for (x0 in c(60, 80, 100)) {
  y <- seq(mu[1] - 10, mu[3] + 10, length.out = 100)
  lines(x0 + 40 * dnorm(y, predict(model, data.frame(height = x0)), sigma(model)), y, col = "blue")
}

########################################################################
# PART B. How good is the line?
########################################################################

# ---- 6) R-squared ----------------------------------------------------

# [DEMO: R2 by hand from the sums of squares]
ss_total <- sum((wheat$biomass - mean(wheat$biomass))^2)
ss_resid <- sum(resid(model)^2)
ss_total
ss_resid
1 - ss_resid / ss_total     # R2 by hand

# [CORE]
summary(model)$r.squared
cor(wheat$height, wheat$biomass)^2   # = r2

# [DEMO: SS total vs SS residual, drawn]
# the two pictures: variation around the mean vs around the line
par(mfrow = c(1, 2))
plot(wheat$height, wheat$biomass, main = "SS total", xlab = "Height (cm)", ylab = "Biomass (g)")
abline(h = mean(wheat$biomass), lwd = 2)
segments(wheat$height, wheat$biomass, wheat$height, mean(wheat$biomass), col = "purple", lty = 2)
plot(wheat$height, wheat$biomass, main = "SS residual", xlab = "Height (cm)", ylab = "Biomass (g)")
abline(model, col = "red", lwd = 2)
segments(wheat$height, wheat$biomass, wheat$height, fitted(model), col = "blue", lty = 2)
par(mfrow = c(1, 1))

# [DEMO: what different R2 values look like (simulated)]
# what different R2 values look like (simulated: same slope, more or less noise)
set.seed(7)
par(mfrow = c(1, 4))
for (r2_target in c(0.9, 0.6, 0.3, 0.05)) {
  x <- runif(60, 0, 10)
  y <- x + rnorm(60, sd = sqrt((1 - r2_target) / r2_target * var(x)))
  plot(x, y, main = paste("R2 =", round(summary(lm(y ~ x))$r.squared, 2)))
  abline(lm(y ~ x), col = "red")
}
par(mfrow = c(1, 1))

# ---- 7) summary(): the full output -----------------------------------

# [CORE]
summary(model)

# the slope's t-test: estimate / SE, 38 df
coef(summary(model))

# [DEMO: t and p by hand]
0.27908 / 0.03737                     # t
2 * pt(7.467, df = 38, lower.tail = FALSE)   # two-sided p-value

# [CORE]
# confidence interval for the slope: estimate +/- t * SE
confint(model)

# [DEMO: the CI by hand, and where the SE comes from (a bootstrap)]
0.27908 + c(-1, 1) * qt(0.975, df = 38) * 0.03737

# where the SE comes from: re-sample the plants and refit (a bootstrap)
set.seed(206)
slopes <- replicate(1000, {
  i <- sample(40, replace = TRUE)   # re-sample the plants
  coef(lm(biomass ~ height, data = wheat[i, ]))[2]
})
sd(slopes)        # spread of the slope
summary(model)$coefficients[2, "Std. Error"]
hist(slopes, xlab = "Slope (g per cm)", main = "1000 re-sampled slopes")

# [CORE]
# the line with its confidence band
ggplot(wheat, aes(x = height, y = biomass)) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_classic()

# [CORE]
# two intervals: for the LINE (mean biomass at that height) and for ONE NEW PLANT
new <- data.frame(height = c(60, 80, 100))
predict(model, newdata = new, interval = "confidence")   # the line
predict(model, newdata = new, interval = "prediction")   # a new plant: much wider

# [HANDY]
# both bands on one plot: predict over a grid of heights, then draw the lines
grid <- data.frame(height = 45:110)
ci <- predict(model, grid, interval = "confidence")
pi <- predict(model, grid, interval = "prediction")
plot(wheat$height, wheat$biomass, xlab = "Height (cm)", ylab = "Biomass (g)", ylim = c(0, 36))
lines(grid$height, ci[, "fit"], col = "red", lwd = 2)
lines(grid$height, ci[, "lwr"], col = "red", lty = 2); lines(grid$height, ci[, "upr"], col = "red", lty = 2)
lines(grid$height, pi[, "lwr"], col = "blue", lty = 3); lines(grid$height, pi[, "upr"], col = "blue", lty = 3)

# ---- 8) The F-test ---------------------------------------------------

# [DEMO: F by hand from R2]
r2 <- summary(model)$r.squared
k  <- 1            # predictors
n  <- nrow(wheat)  # plants
f  <- (r2 / k) / ((1 - r2) / (n - k - 1))
f
pf(f, k, n - k - 1, lower.tail = FALSE)   # p-value

# [CORE]
summary(model)$fstatistic
7.467^2            # one predictor: F = t^2

########################################################################
# PART C. Assumptions: what can go wrong
########################################################################

# ---- 9) Diagnostic plots ---------------------------------------------

# [CORE]
par(mfrow = c(2, 2))   # four panels in one window
plot(model)
par(mfrow = c(1, 1))
# 1 residuals vs fitted: a shapeless cloud = linear
# 2 Q-Q: residuals on the line = normal
# 3 scale-location: flat = equal variance
# 4 residuals vs leverage: points with outsized pull on the line

# [DEMO: the four residual-vs-fitted patterns to recognize (simulated)]
# the four residual-vs-fitted patterns to recognize (simulated)
set.seed(11)
x <- runif(80, 0, 10)
par(mfrow = c(1, 4))
plot(x, rnorm(80), main = "fine", xlab = "fitted", ylab = "residual"); abline(h = 0, lty = 2)
plot(x, -((x - 5)^2) / 6 + 2 + rnorm(80, sd = 0.7), main = "curve", xlab = "fitted", ylab = "residual"); abline(h = 0, lty = 2)
plot(x, rnorm(80, sd = 0.2 + 0.3 * x), main = "fan", xlab = "fitted", ylab = "residual"); abline(h = 0, lty = 2)
plot(x, c(rnorm(79, sd = 0.8), 8), main = "outlier", xlab = "fitted", ylab = "residual"); abline(h = 0, lty = 2)
par(mfrow = c(1, 1))

# ---- 10) Violation 1: a curve ----------------------------------------

# [DEMO: simulate growth that curves with water]
set.seed(1)
water  <- runif(100, 0, 60)   # mm per week
growth <- 20 + 2.2 * water - 0.03 * water^2 +
  rnorm(100, sd = 4)
pots   <- data.frame(water, growth)

# [CORE]
model_w <- lm(growth ~ water, data = pots)
summary(model_w)$r.squared

plot(pots$water, pots$growth, xlab = "Water (mm/week)", ylab = "Growth (g)")
abline(model_w, col = "red", lwd = 2)
plot(model_w, which = 1)   # residuals vs fitted: an arch

# [CORE]
# fix: add a squared term (I() protects the ^ inside a formula)
model_w2 <- lm(growth ~ water + I(water^2), data = pots)
round(coef(model_w2), 3)
summary(model_w2)$r.squared

# [HANDY]
# draw a curved fit: predict over a grid, then lines()
plot(pots$water, pots$growth, xlab = "Water (mm/week)", ylab = "Growth (g)")
lines(0:60, predict(model_w2, data.frame(water = 0:60)), col = "red", lwd = 2)

# [CORE]
plot(model_w2, which = 1)

# ---- 11) Violation 2: a fan ------------------------------------------

# [DEMO: simulate yield whose spread grows with dose]
set.seed(2)
dose  <- runif(150, 0, 60)   # fertilizer, kg/ha
yield <- exp(1.5 + 0.03 * dose + rnorm(150, sd = 0.25))
plots <- data.frame(dose, yield)

# [CORE]
model_f <- lm(yield ~ dose, data = plots)
plot(plots$dose, plots$yield, xlab = "Fertilizer (kg/ha)", ylab = "Yield (t/ha)")
abline(model_f, col = "red", lwd = 2)
plot(model_f, which = 1)   # the spread grows with the fitted value
plot(model_f, which = 3)   # scale-location: rising = unequal variance

# fix: log the response. The slope becomes a percent change.
model_log <- lm(log(yield) ~ dose, data = plots)
round(coef(model_log), 3)
exp(0.030)      # x 1.03 per kg/ha: +3%
plot(model_log, which = 1)

# ---- 11b) Violation 3: residuals that are not normal (a count response) ---

# [DEMO: simulate seed counts]
set.seed(4)
size  <- runif(80, 1, 10)   # plant size
# seeds per plant: a count, many zeros
seeds <- rpois(80, lambda = exp(-1 + 0.35 * size))
pods  <- data.frame(size, seeds)

# [CORE]
model_s <- lm(seeds ~ size, data = pods)
shapiro.test(resid(model_s))
par(mfrow = c(1, 2)); plot(model_s, which = 2); hist(resid(model_s)); par(mfrow = c(1, 1))
# the cure is a model for counts: a Poisson GLM (next week)
# glm(seeds ~ size, data = pods, family = poisson)

# ---- 11c) Violation 4: one bad point (leverage, Cook's distance) -----

# [DEMO: plant a typo in the data]
typo <- wheat
typo$biomass[7] <- 164   # 16.4 typed as 164

# [CORE]
model_t <- lm(biomass ~ height, data = typo)
coef(model)      # the real slope
coef(model_t)    # with the typo
plot(typo$height, typo$biomass); abline(model, col = "darkgreen", lwd = 2); abline(model_t, col = "red", lty = 2, lwd = 2)
which.max(cooks.distance(model_t))
max(cooks.distance(model_t))
plot(model_t, which = 5)   # residuals vs leverage, with Cook's distance contours
# the fix is not statistical: go back to the data sheet

# ---- 11d) Violation 5: not independent (pseudoreplication) -----------

# [DEMO: simulate 5 tillers in each of 8 pots]
set.seed(22)
pot_h <- rnorm(8, sd = 4)     # each pot has its own soil...
pot_b <- rnorm(8, sd = 2.5)   # ...and its own luck
pot   <- rep(1:8, each = 5)   # 5 tillers per pot, 8 pots
height  <- 70 + pot_h[pot] + rnorm(40, sd = 3)
biomass <- 2 + 0.2 * height + pot_b[pot] + rnorm(40, 0, 1)
pots <- data.frame(pot, height, biomass)

# [CORE]
# 40 tillers treated as 40 plants
model_p <- lm(biomass ~ height, data = pots)
round(coef(summary(model_p)), 3)
# the honest version: one value per pot, n = 8
pot_means <- aggregate(cbind(height, biomass) ~ pot, pots, mean)
round(coef(summary(lm(biomass ~ height, pot_means))), 3)
# the giveaway: residuals clump by pot
plot(fitted(model_p), resid(model_p), col = pots$pot, pch = 16); abline(h = 0, lty = 2)
# (a mixed model with pot as a random effect is the Week 10 bonus)

# ---- 11e) Missing data: lm() drops rows silently ---------------------

# [DEMO: knock out three nitrogen values]
wheat_na <- wheat
# three plants were not measured
wheat_na$nitrogen[c(3, 11, 25)] <- NA

# [CORE]
colSums(is.na(wheat_na))
model_na <- lm(biomass ~ height + nitrogen, data = wheat_na)
nobs(model_na)                 # rows actually used
length(resid(model_na))        # shorter than the data frame!
sum(complete.cases(wheat_na))
model_na <- lm(biomass ~ height + nitrogen, wheat_na,
  na.action = na.exclude)
# padded with NA: lines up with the data
length(resid(model_na))
wheat_na$residual <- resid(model_na)

########################################################################
# PART D. Multiple regression
########################################################################

# ---- 12) Two predictors ----------------------------------------------

# [CORE]
model2 <- lm(biomass ~ height + nitrogen, data = wheat)
summary(model2)
# each slope: the effect of that predictor holding the other fixed

# [HANDY]
# the fitted plane (needs scatterplot3d)
library(scatterplot3d)
s3 <- scatterplot3d(wheat$height, wheat$nitrogen, wheat$biomass, pch = 16,
                    xlab = "Height (cm)", ylab = "Nitrogen (kg/ha)", zlab = "Biomass (g)")
s3$plane3d(model2)

# [HANDY]
# "holding the other fixed": the height slope at three nitrogen levels
grid <- expand.grid(height = c(50, 100), nitrogen = c(25, 75, 125))
grid$biomass <- predict(model2, grid)
grid

# [DEMO: parallel lines, one per nitrogen level]
# the same slope for height at each nitrogen level: parallel lines
plot(wheat$height, wheat$biomass, xlab = "Height (cm)", ylab = "Biomass (g)",
     col = c("salmon", "darkgreen", "steelblue")[cut(wheat$nitrogen, c(-1, 50, 100, 151))], pch = 16)
for (nit in c(25, 75, 125)) {
  abline(coef(model2)[1] + coef(model2)[3] * nit, coef(model2)[2],
         col = c("salmon", "darkgreen", "steelblue")[match(nit, c(25, 75, 125))], lwd = 2)
}
legend("topleft", legend = c("low N", "mid N", "high N"), col = c("salmon", "darkgreen", "steelblue"), lwd = 2, bty = "n")

# [CORE]
# why the height slope changed
coef(model)    # biomass ~ height
coef(model2)   # biomass ~ height + nitrogen
cor(wheat$height, wheat$nitrogen)

# ---- 13) R2 only goes up: adjusted R2 --------------------------------

# [DEMO: add a predictor that is pure noise]
set.seed(206)
wheat$noise <- rnorm(40)   # a predictor that is pure noise
r2 <- function(m) summary(m)$r.squared
r2(lm(biomass ~ height, wheat))
r2(lm(biomass ~ height + nitrogen, wheat))
r2(lm(biomass ~ height + nitrogen + noise, wheat))
adj <- function(m) summary(m)$adj.r.squared
adj(lm(biomass ~ height + nitrogen + noise, wheat))

# ---- 13b) Does nitrogen earn its place? Nested models ----------------

# [CORE]
anova(model, model2)   # does nitrogen improve the fit?

# [DEMO: F = t^2 when one term is added]
coef(summary(model2))["nitrogen", "t value"]^2   # = F (one term added)

# [CORE]
# the F-test earns its keep when several terms are added at once:
# a 3-level factor, x + I(x^2), a block of variables
anova(model2)   # sequential: height, then nitrogen
# the other order
anova(lm(biomass ~ nitrogen + height, wheat))
# each term given the others (= the t-tests)
drop1(model2, test = "F")

# ---- 14) The F-test with k = 2 ---------------------------------------

# [DEMO: F by hand with two predictors]
r2 <- summary(model2)$r.squared
k <- 2; n <- 40
f <- (r2 / k) / ((1 - r2) / (n - k - 1))
f
pf(f, k, n - k - 1, lower.tail = FALSE)

# ---- 15) Interactions ------------------------------------------------

# [CORE]
# y ~ x1 * x2  is  x1 + x2 + x1:x2
model_int <- lm(biomass ~ height * nitrogen, data = wheat)
round(coef(summary(model_int)), 4)
# height:nitrogen p = 0.93: no interaction here; keep the simpler model

# [DEMO: what an interaction looks like (simulated)]
# what an interaction looks like (simulated): the slope of x1 depends on x2
set.seed(5)
x1 <- runif(60, 0, 10); x2 <- rep(c("low", "high"), each = 30)
y_no  <- 2 + 1 * x1 + ifelse(x2 == "high", 3, 0) + rnorm(60)             # parallel
y_int <- 2 + ifelse(x2 == "high", 2, 0.4) * x1 + rnorm(60)               # not parallel
par(mfrow = c(1, 2))
plot(x1, y_no, col = ifelse(x2 == "high", "darkcyan", "salmon"), pch = 16, main = "no interaction")
plot(x1, y_int, col = ifelse(x2 == "high", "darkcyan", "salmon"), pch = 16, main = "interaction")
par(mfrow = c(1, 1))
round(coef(summary(lm(y_int ~ x1 * x2))), 3)   # the x1:x2 term is the interaction

# A real one: a fertilizer x irrigation trial.

# [DEMO: simulate the fertilizer x irrigation trial]
set.seed(18)
trial <- data.frame(nitrogen = runif(60, 0, 150),
                    water    = runif(60, 10, 50))
trial$yield <- 5 + 0.01 * trial$nitrogen + 0.05 * trial$water +
  0.002 * trial$nitrogen * trial$water + rnorm(60, sd = 1.2)

# [CORE]
model_nw <- lm(yield ~ nitrogen * water, data = trial)
round(coef(summary(model_nw)), 4)
# the nitrogen slope at three water levels: b1 + b3 * water
coef(model_nw)["nitrogen"] + coef(model_nw)["nitrogen:water"] * c(15, 30, 45)

# [HANDY]
# always plot an interaction: predict over a grid, one line per water level
g <- expand.grid(nitrogen = c(0, 150), water = c(15, 30, 45))
g$yield <- predict(model_nw, g)
plot(trial$nitrogen, trial$yield, xlab = "Nitrogen (kg/ha)", ylab = "Yield (t/ha)")
for (w in c(15, 30, 45)) lines(g$nitrogen[g$water == w], g$yield[g$water == w], lwd = 2, col = c("salmon", "darkgreen", "steelblue")[match(w, c(15, 30, 45))])

# [CORE]
# centring: the product term overlaps its parts
library(car)    # vif()
vif(model_nw)
trial$nitrogen_c <- trial$nitrogen - mean(trial$nitrogen)
trial$water_c    <- trial$water - mean(trial$water)
model_nwc <- lm(yield ~ nitrogen_c * water_c, data = trial)
vif(model_nwc)   # centred: the overlap is gone
# same interaction; main effects now at the mean
round(coef(model_nwc), 4)

########################################################################
# PART E. Hidden relationships and multicollinearity
########################################################################

# ---- 16) Case 1: two paths that cancel (a suppressor) ----------------

# [DEMO: simulate temperature, insects, and yield]
set.seed(1)
temperature <- rnorm(100)
insects     <- rnorm(100) - temperature   # heat: fewer
yield       <- temperature + insects + rnorm(100)

# [CORE]
plot(temperature, yield)
round(coef(summary(lm(yield ~ temperature))), 3)            # nothing
round(coef(summary(lm(yield ~ temperature + insects))), 3)  # both

# [HANDY]
# colour points by a third variable
agri <- data.frame(temperature, insects, yield)
ggplot(agri, aes(temperature, yield, colour = insects)) +
  geom_point(size = 2) +
  scale_colour_gradient2() +
  theme_classic()

# ---- 17) Case 2: a confounder ----------------------------------------

# [DEMO: simulate rainfall, nutrients, and yield]
set.seed(2)
rainfall  <- rnorm(200)
nutrients <- -rainfall + rnorm(200)   # rain leaches them
yield     <- rainfall + rnorm(200)    # nutrients: no effect

# [CORE]
plot(nutrients, yield)
round(coef(summary(lm(yield ~ nutrients))), 3)              # "significant"
round(coef(summary(lm(yield ~ rainfall + nutrients))), 3)   # gone

# ---- 17b) Which variables belong in the model? A collider ------------

# x and y are independent; both cause c. Adjusting for c invents an x-y
# association.

# [DEMO: simulate a collider]
set.seed(6)
x <- rnorm(300)
y <- rnorm(300)
c <- x + y + rnorm(300, sd = 0.5)    # e.g. selecting plants by overall size

# [CORE]
round(coef(summary(lm(y ~ x))), 3)       # nothing, correctly
round(coef(summary(lm(y ~ x + c))), 3)   # a strong negative x effect from nowhere
# confounders: adjust. mediators: depends on the question. colliders: leave out.

# ---- 18) Case 3: multicollinearity -----------------------------------

# [DEMO: simulate a predictor that is almost the sum of two others]
set.seed(4)
temperature   <- rnorm(200)
precipitation <- rnorm(200)
microbiome    <- temperature + precipitation +
  rnorm(200, sd = 0.1)                      # almost their sum
growth        <- temperature + precipitation + rnorm(200)

# [CORE]
model_c <- lm(growth ~ temperature + precipitation +
  microbiome)
summary(model_c)   # F huge, R2 fine, no coefficient significant

# variance inflation factors (car package)
library(car)
vif(model_c)
model_c2 <- update(model_c, . ~ . - microbiome)   # drop it
vif(model_c2)
round(coef(summary(model_c2)), 3)
# rule of thumb: VIF > 5 look closely; VIF > 10 drop or combine

# [DEMO: perfect collinearity, R drops the redundant term (NA)]
# perfect collinearity: R drops the redundant term (NA)
microbiome2 <- temperature + precipitation   # their exact sum
coef(lm(growth ~ temperature + precipitation + microbiome2))

# ---- 19) Real data: bird density in forest patches -------------------

# bd = bird density; grzinv = grazing (inverse); dist = distance to
# town; height = forest height; peri = patch perimeter; mammal = large
# mammals; area = patch area; leg = cover by leguminous trees
# (predictors standardized).

# [CORE]
birds <- read.csv("https://greymonroe.github.io/PLS_206/data/birds.csv")[, -1]   # drop column 1
head(round(birds, 2), 4)
model_birds <- lm(bd ~ ., data = birds)   # "." = all others
summary(model_birds)
round(vif(model_birds), 1)

model_birds2 <- update(model_birds, . ~ . - dist)
round(vif(model_birds2), 1)
round(coef(summary(model_birds2)), 3)
# the same model, two other ways
model_birds3 <- lm(bd ~ . - dist, data = birds)
model_birds4 <- lm(bd ~ grzinv + height + peri + mammal + area + leg, data = birds)

########################################################################
# PART F. Under the hood, and reporting
########################################################################

# ---- 20) The design matrix -------------------------------------------

# [DEMO: least squares in one line of matrix algebra]
head(model.matrix(model2), 4)   # a column of 1s for the intercept
# least squares in one line of matrix algebra: (X'X)^-1 X'y
X <- model.matrix(model2); y <- wheat$biomass
solve(t(X) %*% X) %*% t(X) %*% y
coef(model2)

# ---- 21) Reporting: effect size in units, CI, R2, n ------------------

# [CORE]
coef(model)[2]; confint(model)[2, ]; summary(model)$r.squared; nrow(wheat)
# "Biomass increased with height (0.28 g per cm, 95% CI 0.20 to 0.35;
#  t38 = 7.5, p < 0.001). Height explained 59% of the variation in
#  biomass (R2 = 0.59, n = 40)."

# [CORE]
# which predictor matters more? standardize both: slopes in SD units
model_std <- lm(scale(biomass) ~ scale(height) +
  scale(nitrogen), data = wheat)
round(coef(model_std), 2)   # slopes in SD units