Problem Set 2: Linear Regression

PLS 206 — Applied Multivariate Modeling

Author

Grey Monroe

Published

January 1, 2026

Note

Due: Monday, October 12, 1:10 PM (start of class) — submitted through Canvas.

Parts 1–5 use the first half of the week’s lectures (one predictor); Parts 6–8 use the second half (more than one predictor, confounding, multicollinearity). The Week 2 script, Week2_Linear_Regression.R, has working code for every step: Parts A–C of the script cover Parts 1–5 here, Parts D–F cover Parts 6–8.

Submission Instructions

Upload two files to Canvas:

Written answers — a .pdf document with your written responses, screenshots, and interpretations.
R script — a .R file with all the code you used to produce your answers.

List the first and last names of any collaborators in both files.

Naming convention: HW2emailID.R (code) and HW2emailID.pdf (written answers), where emailID is the part of your UC Davis email before the @.

Example: Grey Monroe would submit HW2gmonroe.R and HW2gmonroe.pdf.

Feel free to use Google, StackOverflow, and other resources if you are stuck with code, and post questions on Slack if anything is unclear.

The data: flowering time in wild wheat

In plants, the timing of flowering is one of the most important ecophysiological traits determining which environments a genotype will thrive in. This week’s dataset is flowering time (days from germination to flowering) for 54 genotypes of wild einkorn wheat (Triticum monococcum) collected across the Middle East and Europe. Flowering time was measured in a common garden: seeds from every collection site were grown together in a greenhouse, so differences among genotypes are genetic, not environmental.

Using the latitude and longitude of each genotype’s origin, we extracted the 19 bioclimatic variables of its home environment from WorldClim (bio1 to bio19). We will use linear regression to test hypotheses about climate drivers of natural variation in flowering time.

Download wild_wheat.csv from the course data folder and read it into R:

wheat <- read.csv("wild_wheat.csv")
str(wheat)

Part 1: Look at the Data

Q1. What is the mean flowering time, and its standard deviation? Plot a histogram of flowering time and paste it in your written answers.

Hint

mean(), sd(), hist(wheat$Flowering). Script Part A, section 1 (reading and looking at a data frame); Week 1 script for the summaries.

Q2. How many climate variables (bio1 … bio19) are there? Look up the WorldClim definitions: what do bio1 and bio13 measure? (bio1 is in °C.)

Hint

grep("^bio", names(wheat), value = TRUE) lists the climate columns; length() counts them. The definitions are on the WorldClim page linked above.

Part 2: One Predictor

Q3. Pick a climate variable you hypothesize could affect flowering time (for example: genotypes from drier places flower earlier). State the hypothesis and your biological reasoning in two or three sentences, before running anything.

Hint

No code: reasoning only. The data description above and the WorldClim variable list are all you need.

Q4. Test your hypothesis with a linear model: lm(Flowering ~ yourvariable, data = wheat). Report, from summary() and confint():

1. the slope, in its units (days of flowering per unit of your variable), with its 95% confidence interval
1. the p-value for the slope
1. R²: how much of the variation in flowering time your variable explains

Hint

model <- lm(Flowering ~ bio1, data = wheat); summary(model) for the slope, p, and R²; confint(model) for the CI. Script Part A section 4 and Part B section 7; slides lm() finds the same line, summary(model), Confidence interval for the slope.

What do you conclude?

Q5. Run cor.test() on the same two variables. How do its r, t statistic, and p-value relate to the numbers from your linear model in Q4? (Hint: square r.)

Hint

cor.test(wheat$bio1, wheat$Flowering); compare its t and p-value lines with the slope row of summary(model), and square its cor estimate. Script Part B section 6 (cor()^2); slides R²: the share explained and The F-test (F = t²).

Q6. Write the result from Q4 as one sentence for the results section of a paper: effect size in units with its CI, the test statistic, the p-value, R², and n.

Hint

Slide Reporting a regression has the template sentence; swap in your numbers from Q4.

Part 3: Screen All the Climate Variables

Q7. Calculate the correlation between flowering time and every climate variable (bio1 … bio19). Kudos if you do it with a for() loop or an apply()-family function instead of 19 copies of cor(). Which variable has the strongest correlation with flowering time (largest absolute value)? Is it positive or negative?

Hint

Put the column names in a vector (bios from Q2), then either for (b in bios) cors[b] <- cor(wheat[[b]], wheat$Flowering) or sapply(wheat[, bios], function(v) cor(v, wheat$Flowering)). which.max(abs(cors)) finds the strongest. The birds section at the end of Linear Regression 2.R shows the same loop-vs-apply pattern; the Friday crash-course script covers for() and apply().

Q8. Build a linear model between flowering time and that variable. Make a scatterplot with the fitted line on it (abline(model) in base R, or geom_smooth(method = "lm") in ggplot) and paste it in your written answers. In one or two plain-English sentences, how would you describe this result to someone who has not taken a statistics course?

Hint

plot(x, y) then abline(model, col = "red"), or ggplot(...) + geom_point() + geom_smooth(method = "lm"). Script Part A section 4 and Part B section 7; slide The band: geom_smooth(method = “lm”).

Part 4: Check the Fit

Q9. Run plot(model) on your Q8 model (use par(mfrow = c(2, 2)) to see all four panels). Paste the panels and answer, one sentence each:

1. Residuals vs fitted: does the relationship look linear, and is the variance roughly equal across the fitted values?
1. Q–Q plot: do the residuals look normal? Confirm with shapiro.test(resid(model)).
1. Residuals vs leverage: is any genotype pulling the line unusually hard?

Hint

par(mfrow = c(2, 2)); plot(model); shapiro.test(resid(model)). Script Part C section 9; slides plot(model): four diagnostic plots, Four patterns, four fixes, Finding it: residuals vs leverage.

Part 5: Predictions

Q10. Based on your Q8 model, what is the predicted change in flowering time for a one degree increase in annual mean temperature (bio1)? Where in the model output does that number live?

Hint

coef(model); the slope is the second number. Slide lm() finds the same line (reading a slope in units).

Q11. Use predict() to calculate the flowering time of every genotype if its home climate warmed by 2 °C (newdata = data.frame(bio1 = wheat$bio1 + 2)). What is the new mean flowering time, and how much does it differ from the observed mean in Q1? Explain why the difference is exactly what the slope implies.

Hint

predict(model, newdata = data.frame(bio1 = wheat$bio1 + 2)) gives 54 predictions; mean() them and subtract the Q1 mean. Script Part A section 5; slide Predictions: predict().

Q12. Predict the flowering time of a genotype from a site with bio1 = 30. Why should you distrust this number? (range(wheat$bio1) helps.)

Hint

predict(model, newdata = data.frame(bio1 = 30)); range(wheat$bio1). Script section 5 (the “nonsense plants” lines); slide Interpolate, don’t extrapolate.

Part 6: A Reviewer’s Objection (more than one predictor)

Congratulations: you have submitted your findings. A reviewer writes back:

How do you know the relationship between flowering time and temperature isn’t driven by a variable that is correlated with temperature? Latitude affects day length, which could select on flowering time directly, and latitude is also correlated with temperature.

Q13. Calculate the correlation between bio1 and Latitude. Draw (by hand or in words) the causal diagram the reviewer has in mind: which variable is the proposed confounder?

Hint

cor(wheat$bio1, wheat$Latitude). For the diagram: three boxes and two or three arrows; slides Reveal: a confounder and Confounder: adjust for it show the shape.

Q14. Build a linear model of flowering time as a function of latitude alone. What does it suggest?

Hint

summary(lm(Flowering ~ Latitude, data = wheat)).

Q15. Build a linear model with flowering time as a function of both latitude and bio1. What do you conclude, and why? Justify with the slopes and p-values from the model. Check vif() (package car) on this two-predictor model: is collinearity a problem here?

Hint

m3 <- lm(Flowering ~ Latitude + bio1, data = wheat); summary(m3); library(car); vif(m3). Script Part D section 12 (two predictors, why the height slope changed) and Part E section 18 (VIF); slides Why the height slope changed, Detect it: the variance inflation factor.

Part 7: The Kitchen-Sink Model

Q16. Build a linear model of flowering time as a function of all climate variables plus latitude and longitude (Flowering ~ . - Accession is a shortcut). Look at summary(): bio7 has NA for every value. Why? Look at the WorldClim definitions of bio5, bio6, and bio7 again, and check your explanation in R.

Hint

full <- lm(Flowering ~ . - Accession, data = wheat) (the dot means every other column); summary(full). To test your explanation: all.equal(wheat$bio7, wheat$bio5 - wheat$bio6). Script Part E section 18 (the microbiome2 example); slide Perfect collinearity: R gives up.

Q17. Remove bio7 and rebuild the model. Compare its R² and adjusted R² with the one-predictor model from Q8. Why are they so far apart here?

Hint

full2 <- lm(Flowering ~ . - Accession - bio7, data = wheat); summary(full2)$r.squared and summary(full2)$adj.r.squared, same for the Q8 model. Script Part D section 13; slide R² only goes up: use adjusted R².

Q18. In this model, temperature (bio1) is no longer significant (nor is anything else). Why not? Check with vif(). Do the results suggest excessive collinearity, and what could you do about it? Two or three sentences.

Hint

round(vif(full2), 1); the rule of thumb is on slide Detect it: the variance inflation factor. Script Part E sections 18–19 (dropping the worst offender, update(model, . ~ . - x)).

Part 8: Reporting

Q19. Of the one-predictor model (Q8), the two-predictor model (Q15), and the kitchen-sink model (Q17), which would you report in a paper, and why? Two or three sentences; mention at least two of: effect size, confidence interval, adjusted R², VIF, interpretability.

Hint

No new code. Slides Reporting a regression and Rough summary and guide; think effect size with CI, adjusted R², VIF, and whether you could explain the slopes to someone.

Q20. Standardize both predictors in your Q15 model (scale()) and refit. Which predictor matters more in standard-deviation units? Why is this comparison not possible from the raw slopes?

Hint

lm(scale(Flowering) ~ scale(Latitude) + scale(bio1), data = wheat). Script Part F section 21; slide Reporting a regression (slopes in SD units).

--- title: "Problem Set 2: Linear Regression" subtitle: "PLS 206 — Applied Multivariate Modeling" author: "Grey Monroe" date: "Fall 2026" format: html: theme: flatly toc: true toc-depth: 2 toc-title: "Parts" highlight-style: github code-fold: false execute: echo: true eval: false --- ::: {.callout-note} **Due:** Monday, October 12, 1:10 PM (start of class) — submitted through Canvas. Parts 1–5 use the first half of the week's lectures (one predictor); Parts 6–8 use the second half (more than one predictor, confounding, multicollinearity). The Week 2 script, [Week2_Linear_Regression.R](Week2_Linear_Regression.R), has working code for every step: Parts A–C of the script cover Parts 1–5 here, Parts D–F cover Parts 6–8. ::: --- ## Submission Instructions Upload **two files** to Canvas: 1. **Written answers** — a `.pdf` document with your written responses, screenshots, and interpretations. 2. **R script** — a `.R` file with all the code you used to produce your answers. List the first and last names of any collaborators in both files. **Naming convention:** `HW2emailID.R` (code) and `HW2emailID.pdf` (written answers), where `emailID` is the part of your UC Davis email before the `@`. > **Example:** Grey Monroe would submit `HW2gmonroe.R` and `HW2gmonroe.pdf`. Feel free to use Google, StackOverflow, and other resources if you are stuck with code, and post questions on Slack if anything is unclear. --- ## The data: flowering time in wild wheat In plants, the timing of flowering is one of the most important ecophysiological traits determining which environments a genotype will thrive in. This week's dataset is flowering time (days from germination to flowering) for 54 genotypes of wild einkorn wheat (*Triticum monococcum*) collected across the Middle East and Europe. Flowering time was measured in a common garden: seeds from every collection site were grown together in a greenhouse, so differences among genotypes are genetic, not environmental. Using the latitude and longitude of each genotype's origin, we extracted the 19 bioclimatic variables of its home environment from [WorldClim](https://www.worldclim.org/data/bioclim.html) (`bio1` to `bio19`). We will use linear regression to test hypotheses about climate drivers of natural variation in flowering time. Download `wild_wheat.csv` from the [course data folder](../data/wild_wheat.csv) and read it into R: ```{r} wheat <- read.csv("wild_wheat.csv") str(wheat) ``` --- ## Part 1: Look at the Data **Q1.** What is the **mean** flowering time, and its **standard deviation**? Plot a histogram of flowering time and paste it in your written answers. ::: {.callout-tip collapse="true"} ## Hint `mean()`, `sd()`, `hist(wheat$Flowering)`. Script Part A, section 1 (reading and looking at a data frame); Week 1 script for the summaries. ::: **Q2.** How many climate variables (`bio1` … `bio19`) are there? Look up the WorldClim definitions: what do `bio1` and `bio13` measure? (`bio1` is in °C.) ::: {.callout-tip collapse="true"} ## Hint `grep("^bio", names(wheat), value = TRUE)` lists the climate columns; `length()` counts them. The definitions are on the WorldClim page linked above. ::: --- ## Part 2: One Predictor **Q3.** Pick a climate variable you hypothesize could affect flowering time (for example: genotypes from drier places flower earlier). State the hypothesis and your biological reasoning in two or three sentences, **before** running anything. ::: {.callout-tip collapse="true"} ## Hint No code: reasoning only. The data description above and the WorldClim variable list are all you need. ::: **Q4.** Test your hypothesis with a linear model: `lm(Flowering ~ yourvariable, data = wheat)`. Report, from `summary()` and `confint()`: - a) the **slope**, in its units (days of flowering per unit of your variable), with its **95% confidence interval** - b) the **p-value** for the slope - c) **R²**: how much of the variation in flowering time your variable explains ::: {.callout-tip collapse="true"} ## Hint `model <- lm(Flowering ~ bio1, data = wheat)`; `summary(model)` for the slope, p, and R²; `confint(model)` for the CI. Script Part A section 4 and Part B section 7; slides *lm() finds the same line*, *summary(model)*, *Confidence interval for the slope*. ::: What do you conclude? **Q5.** Run `cor.test()` on the same two variables. How do its r, t statistic, and p-value relate to the numbers from your linear model in Q4? (Hint: square r.) ::: {.callout-tip collapse="true"} ## Hint `cor.test(wheat$bio1, wheat$Flowering)`; compare its `t` and `p-value` lines with the slope row of `summary(model)`, and square its `cor` estimate. Script Part B section 6 (`cor()^2`); slides *R²: the share explained* and *The F-test* (F = t²). ::: **Q6.** Write the result from Q4 as **one sentence for the results section of a paper**: effect size in units with its CI, the test statistic, the p-value, R², and n. ::: {.callout-tip collapse="true"} ## Hint Slide *Reporting a regression* has the template sentence; swap in your numbers from Q4. ::: --- ## Part 3: Screen All the Climate Variables **Q7.** Calculate the correlation between flowering time and **every** climate variable (`bio1` … `bio19`). Kudos if you do it with a `for()` loop or an `apply()`-family function instead of 19 copies of `cor()`. Which variable has the **strongest** correlation with flowering time (largest absolute value)? Is it positive or negative? ::: {.callout-tip collapse="true"} ## Hint Put the column names in a vector (`bios` from Q2), then either `for (b in bios) cors[b] <- cor(wheat[[b]], wheat$Flowering)` or `sapply(wheat[, bios], function(v) cor(v, wheat$Flowering))`. `which.max(abs(cors))` finds the strongest. The `birds` section at the end of *Linear Regression 2.R* shows the same loop-vs-apply pattern; the Friday crash-course script covers `for()` and `apply()`. ::: **Q8.** Build a linear model between flowering time and that variable. Make a scatterplot with the fitted line on it (`abline(model)` in base R, or `geom_smooth(method = "lm")` in ggplot) and paste it in your written answers. In one or two plain-English sentences, how would you describe this result to someone who has not taken a statistics course? ::: {.callout-tip collapse="true"} ## Hint `plot(x, y)` then `abline(model, col = "red")`, or `ggplot(...) + geom_point() + geom_smooth(method = "lm")`. Script Part A section 4 and Part B section 7; slide *The band: geom_smooth(method = "lm")*. ::: --- ## Part 4: Check the Fit **Q9.** Run `plot(model)` on your Q8 model (use `par(mfrow = c(2, 2))` to see all four panels). Paste the panels and answer, one sentence each: - a) Residuals vs fitted: does the relationship look linear, and is the variance roughly equal across the fitted values? - b) Q–Q plot: do the residuals look normal? Confirm with `shapiro.test(resid(model))`. - c) Residuals vs leverage: is any genotype pulling the line unusually hard? ::: {.callout-tip collapse="true"} ## Hint `par(mfrow = c(2, 2)); plot(model)`; `shapiro.test(resid(model))`. Script Part C section 9; slides *plot(model): four diagnostic plots*, *Four patterns, four fixes*, *Finding it: residuals vs leverage*. ::: --- ## Part 5: Predictions **Q10.** Based on your Q8 model, what is the predicted change in flowering time for a **one degree** increase in annual mean temperature (`bio1`)? Where in the model output does that number live? ::: {.callout-tip collapse="true"} ## Hint `coef(model)`; the slope is the second number. Slide *lm() finds the same line* (reading a slope in units). ::: **Q11.** Use `predict()` to calculate the flowering time of every genotype if its home climate warmed by **2 °C** (`newdata = data.frame(bio1 = wheat$bio1 + 2)`). What is the new mean flowering time, and how much does it differ from the observed mean in Q1? Explain why the difference is exactly what the slope implies. ::: {.callout-tip collapse="true"} ## Hint `predict(model, newdata = data.frame(bio1 = wheat$bio1 + 2))` gives 54 predictions; `mean()` them and subtract the Q1 mean. Script Part A section 5; slide *Predictions: predict()*. ::: **Q12.** Predict the flowering time of a genotype from a site with `bio1 = 30`. Why should you distrust this number? (`range(wheat$bio1)` helps.) ::: {.callout-tip collapse="true"} ## Hint `predict(model, newdata = data.frame(bio1 = 30))`; `range(wheat$bio1)`. Script section 5 (the "nonsense plants" lines); slide *Interpolate, don't extrapolate*. ::: --- ## Part 6: A Reviewer's Objection *(more than one predictor)* Congratulations: you have submitted your findings. A reviewer writes back: > *How do you know the relationship between flowering time and temperature isn't driven by a variable that is correlated with temperature? Latitude affects day length, which could select on flowering time directly, and latitude is also correlated with temperature.* **Q13.** Calculate the correlation between `bio1` and `Latitude`. Draw (by hand or in words) the causal diagram the reviewer has in mind: which variable is the proposed **confounder**? ::: {.callout-tip collapse="true"} ## Hint `cor(wheat$bio1, wheat$Latitude)`. For the diagram: three boxes and two or three arrows; slides *Reveal: a confounder* and *Confounder: adjust for it* show the shape. ::: **Q14.** Build a linear model of flowering time as a function of **latitude alone**. What does it suggest? ::: {.callout-tip collapse="true"} ## Hint `summary(lm(Flowering ~ Latitude, data = wheat))`. ::: **Q15.** Build a linear model with flowering time as a function of **both** latitude and `bio1`. What do you conclude, and why? Justify with the slopes and p-values from the model. Check `vif()` (package `car`) on this two-predictor model: is collinearity a problem here? ::: {.callout-tip collapse="true"} ## Hint `m3 <- lm(Flowering ~ Latitude + bio1, data = wheat)`; `summary(m3)`; `library(car); vif(m3)`. Script Part D section 12 (two predictors, *why the height slope changed*) and Part E section 18 (VIF); slides *Why the height slope changed*, *Detect it: the variance inflation factor*. ::: --- ## Part 7: The Kitchen-Sink Model **Q16.** Build a linear model of flowering time as a function of **all** climate variables plus latitude and longitude (`Flowering ~ . - Accession` is a shortcut). Look at `summary()`: `bio7` has `NA` for every value. Why? Look at the WorldClim definitions of `bio5`, `bio6`, and `bio7` again, and check your explanation in R. ::: {.callout-tip collapse="true"} ## Hint `full <- lm(Flowering ~ . - Accession, data = wheat)` (the dot means every other column); `summary(full)`. To test your explanation: `all.equal(wheat$bio7, wheat$bio5 - wheat$bio6)`. Script Part E section 18 (the `microbiome2` example); slide *Perfect collinearity: R gives up*. ::: **Q17.** Remove `bio7` and rebuild the model. Compare its **R²** and **adjusted R²** with the one-predictor model from Q8. Why are they so far apart here? ::: {.callout-tip collapse="true"} ## Hint `full2 <- lm(Flowering ~ . - Accession - bio7, data = wheat)`; `summary(full2)$r.squared` and `summary(full2)$adj.r.squared`, same for the Q8 model. Script Part D section 13; slide *R² only goes up: use adjusted R²*. ::: **Q18.** In this model, temperature (`bio1`) is no longer significant (nor is anything else). Why not? Check with `vif()`. Do the results suggest excessive collinearity, and what could you do about it? Two or three sentences. ::: {.callout-tip collapse="true"} ## Hint `round(vif(full2), 1)`; the rule of thumb is on slide *Detect it: the variance inflation factor*. Script Part E sections 18–19 (dropping the worst offender, `update(model, . ~ . - x)`). ::: --- ## Part 8: Reporting **Q19.** Of the one-predictor model (Q8), the two-predictor model (Q15), and the kitchen-sink model (Q17), which would you report in a paper, and why? Two or three sentences; mention at least two of: effect size, confidence interval, adjusted R², VIF, interpretability. ::: {.callout-tip collapse="true"} ## Hint No new code. Slides *Reporting a regression* and *Rough summary and guide*; think effect size with CI, adjusted R², VIF, and whether you could explain the slopes to someone. ::: **Q20.** Standardize both predictors in your Q15 model (`scale()`) and refit. Which predictor matters more in standard-deviation units? Why is this comparison not possible from the raw slopes? ::: {.callout-tip collapse="true"} ## Hint `lm(scale(Flowering) ~ scale(Latitude) + scale(bio1), data = wheat)`. Script Part F section 21; slide *Reporting a regression* (slopes in SD units). :::