library(car) # vif()
wheat <- read.csv("wild_wheat.csv")
str(wheat) # 54 obs. of 23 variables: Accession, Flowering, Latitude, Longitude, bio1-bio19Problem Set 2 Key: Linear Regression
PLS 206 — Applied Multivariate Modeling
Tip
Try the problems yourself before consulting this key. Every number below was checked in R against data/wild_wheat.csv (54 genotypes, 23 columns) on 2026-10-02.
Part 1: Look at the Data
Q1. Mean flowering time 105.2 days, SD 9.87 days.
mean(wheat$Flowering) # 105.17
sd(wheat$Flowering) # 9.87
hist(wheat$Flowering, xlab = "Days to flowering", main = "")Q2. 19 climate variables. bio1 = annual mean temperature (°C here); bio13 = precipitation of the wettest month.
bios <- grep("^bio", names(wheat), value = TRUE)
length(bios) # 19Part 2: One Predictor
Q3. Any reasoned hypothesis is fine. Example: genotypes from hotter, drier sites flower earlier to finish their life cycle before summer drought (drought escape), so flowering should decline with bio1.
Q4. With bio1 (the key uses this; other variables are fine if the numbers are reported correctly):
m1 <- lm(Flowering ~ bio1, data = wheat)
summary(m1)
confint(m1)- slope −2.77 days per °C, 95% CI −3.37 to −2.17
- p = 1.3 × 10⁻¹² (t = −9.28 on 52 df)
- R² = 0.62
Genotypes from warmer sites flower earlier: each extra degree of annual mean temperature at home goes with flowering about 2.8 days earlier, and temperature explains 62% of the variation among genotypes.
Q5. cor.test() gives r = −0.790, t = −9.28, p = 1.3 × 10⁻¹². The t and p are identical to the slope’s t-test; r² = 0.790² = 0.62 = R². With one predictor, the correlation test and the regression are the same test.
cor.test(wheat$bio1, wheat$Flowering)
(-0.7895)^2Q6. Example: “Flowering time declined with annual mean temperature at the site of origin (−2.77 days per °C, 95% CI −3.37 to −2.17; t₅₂ = −9.28, p < 0.001; R² = 0.62, n = 54).”
Part 3: Screen All the Climate Variables
Q7. Strongest: bio1, r = −0.79 (negative). Next: bio10 (−0.73), bio11 (−0.68), bio5 (−0.67).
cors <- sapply(wheat[, bios], function(v) cor(v, wheat$Flowering))
round(sort(cors), 2)
names(which.max(abs(cors))) # "bio1"
# the same with a for() loop
cors <- c()
for (b in bios) cors[b] <- cor(wheat[[b]], wheat$Flowering)Q8. Same model as Q4. Plain English: wild wheat from warmer places flowers earlier; across the collection, a genotype from a site 10 °C warmer flowers about four weeks earlier on average, and temperature accounts for most, but not all, of the differences among genotypes.
plot(wheat$bio1, wheat$Flowering, xlab = "Annual mean temperature (°C)", ylab = "Days to flowering")
abline(m1, col = "red", lwd = 2)
library(ggplot2)
ggplot(wheat, aes(bio1, Flowering)) + geom_point() + geom_smooth(method = "lm") + theme_classic()Part 4: Check the Fit
Q9.
par(mfrow = c(2, 2)); plot(m1); par(mfrow = c(1, 1))
shapiro.test(resid(m1)) # W = 0.958, p = 0.058- Residuals vs fitted: no curve; the spread is similar across fitted values (a mild hint of less spread at the warm end, nothing alarming).
- Q–Q: roughly on the line with slight tail departures; Shapiro p = 0.058, so no strong evidence against normality.
- No genotype falls near the Cook’s distance contours; no single point drives the line.
Part 5: Predictions
Q10. −2.77 days per °C: the bio1 coefficient in the Estimate column of summary(m1) (or coef(m1)[2]).
Q11. New mean 99.6 days, i.e. 5.55 days earlier than the observed mean of 105.2. Shifting every genotype by +2 °C shifts every prediction by 2 × (−2.77) = −5.55 days, so the mean shifts by exactly that.
warmer <- predict(m1, newdata = data.frame(bio1 = wheat$bio1 + 2))
mean(warmer) # 99.62
mean(warmer) - mean(wheat$Flowering) # -5.55
2 * coef(m1)[2] # -5.55Q12. The prediction is 144.9 − 2.77 × 30 = 61.6 days, but bio1 in the data runs from 8.0 to 18.9 °C. 30 °C is far outside the range: the line is an assumption there, not an estimate (no wild einkorn grows at an annual mean of 30 °C, and the relationship need not stay linear).
predict(m1, newdata = data.frame(bio1 = 30)) # 61.6
range(wheat$bio1) # 8.0 to 18.9Part 6: A Reviewer’s Objection
Q13. r(bio1, Latitude) = −0.60: sites farther north are colder. The reviewer’s diagram: Latitude → Flowering (via day length) and Latitude → Temperature, with the observed Temperature–Flowering link possibly borrowed from latitude. Latitude is the proposed confounder.
cor(wheat$bio1, wheat$Latitude) # -0.596Q14. Flowering increases with latitude: +1.14 days per degree of latitude, p = 2 × 10⁻⁴, R² = 0.23. On its own, latitude looks like a real predictor, consistent with the reviewer’s worry.
summary(lm(Flowering ~ Latitude, data = wheat))Q15. With both in the model, temperature keeps its slope (−2.74 days per °C, p < 0.001) while latitude collapses to 0.04 days per degree, p = 0.88. Latitude’s apparent effect in Q14 was borrowed from temperature, not the other way round: once temperature is held fixed, latitude adds nothing. VIF = 1.55 for both: moderate correlation, no collinearity problem, so the comparison is trustworthy.
m3 <- lm(Flowering ~ Latitude + bio1, data = wheat)
summary(m3)
vif(m3) # 1.55, 1.55Part 7: The Kitchen-Sink Model
Q16. bio7 (temperature annual range) is defined as bio5 − bio6 (max temperature of the warmest month minus min temperature of the coldest month). It is an exact linear combination of two other predictors, so the model has no unique solution for its slope and R drops it (NA, “not defined because of singularities”): perfect multicollinearity.
full <- lm(Flowering ~ . - Accession, data = wheat)
summary(full)
all.equal(wheat$bio7, wheat$bio5 - wheat$bio6) # TRUEQ17. R² rises from 0.62 to 0.74, but adjusted R² is 0.59, below the one-predictor model’s 0.62 (adjusted 0.62). Twenty-one predictors for 54 genotypes: R² can only go up as predictors are added, and the adjustment charges for each one; the extra 20 predictors are not earning their keep.
full2 <- lm(Flowering ~ . - Accession - bio7, data = wheat)
summary(full2)$r.squared # 0.744
summary(full2)$adj.r.squared # 0.588
summary(m1)$adj.r.squared # 0.616Q18. The 19 bioclim variables are different summaries of the same two things (temperature and precipitation), so they are highly correlated with one another: bio1 is nearly a function of bio5, bio6, bio10, bio11. The model cannot tell whose slope the temperature effect belongs to, so every standard error balloons (bio1’s SE goes from 0.30 to 15.7) and nothing is significant even though the model as a whole fits. VIFs are enormous (bio10 ≈ 19,000; most in the hundreds to thousands; the rule of thumb is > 10 is a problem). Remedies: keep one representative variable per climate axis (e.g. bio1 and bio12), drop the highest-VIF variables one at a time, or combine them into composite axes (PCA, Week 7).
round(vif(full2), 1)Part 8: Reporting
Q19. The one-predictor model (or the two-predictor model, which gives the same temperature slope and answers the reviewer): a precise, interpretable effect (−2.77 days per °C, CI −3.37 to −2.17), the best adjusted R², and no collinearity (VIF 1.55). The kitchen-sink model has a higher raw R² but a lower adjusted R², uninterpretable slopes with huge SEs, and VIFs in the thousands; nothing in it can be reported as an effect.
Q20. In SD units, temperature matters far more: about −0.78 SD of flowering per SD of bio1 versus 0.02 SD per SD of latitude. Raw slopes can’t be compared because they are in different units (days per °C vs days per degree of latitude); standardizing puts both on the same scale.
m_std <- lm(scale(Flowering) ~ scale(Latitude) + scale(bio1), data = wheat)
round(coef(m_std), 2)