Problem Set 3: General(ized) Linear Models

PLS 206 — Applied Multivariate Modeling

Author

Grey Monroe

Published

January 1, 2026

Note

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

Parts 1–2 use the first half of the week (groups as predictors, interactions); Parts 3–4 use the second half (logistic regression). The week’s code page has working code for every step; the hints point to its sections.

Submission Instructions

Upload two files to Canvas:

Written answers — a .pdf document with your written responses 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: HW3emailID.R (code) and HW3emailID.pdf (written answers), where emailID is the part of your UC Davis email before the @.

Example: Grey Monroe would submit HW3gmonroe.R and HW3gmonroe.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: disease and yield in pistachio

Pistachios are an important crop for California agriculture. As in most crops, genotypes can differ in traits like yield and susceptibility to disease. This simulated dataset has 300 trees of three pistachio genotypes (“kerman”, “chaparillo”, and “t41”; 100 trees each) planted in a common garden. At the end of the growing season (September), each tree was scored for disease (disease: 1 = diseased, 0 = healthy; a disease that causes leaf die-off) and its yield was measured (yield: total harvested nuts). Drones also measured each tree’s greenness (NDVI) in March, June, and August (NDVI_march, NDVI_june, NDVI_august); larger NDVI means a greener tree.

We want to know (a) whether yield differs between genotypes, (b) whether the disease affects yield, and differently in different genotypes, (c) whether rates of disease differ between genotypes, and (d) whether the drone data can flag diseased trees before the end of the season.

Read it straight from the course website:

pistachio <- read.csv("https://greymonroe.github.io/PLS_206/data/pistachio.csv")
str(pistachio)

Part 1: Genotype and Yield

Q1. Plot yield by genotype (a box or violin plot with the points on top). Fit lm(yield ~ genotype). In one sentence each: what does the (Intercept) estimate mean, and what does the genotypekerman estimate mean? Give both in yield units.

Hint

ggplot(pistachio, aes(genotype, yield)) + geom_boxplot() + geom_jitter(width = 0.1); summary(lm(yield ~ genotype, data = pistachio)); tapply(pistachio$yield, pistachio$genotype, mean) to check. Code page sections 2 and 6; slides Two groups, one line and Three groups, two dummies.

Q2. Do the genotypes differ in yield? Answer with one test of the whole factor (anova()), then say which pairs (if any) differ with TukeyHSD().

Hint

anova(model); TukeyHSD(aov(model)). Code page sections 7 and 8; slides Do the varieties differ at all? anova() and Which pairs? TukeyHSD().

Part 2: Disease, Yield, and Genotype

Q3. Plot yield by genotype with a separate box for diseased and healthy trees (fill = factor(disease)). Fit lm(yield ~ genotype * factor(disease)). Use emmeans to get the yield lost to disease within each genotype. Which genotype loses the least, and how much does it lose?

Hint

pistachio$disease_f <- factor(pistachio$disease) first; then m2 <- lm(yield ~ genotype * disease_f, data = pistachio) and pairs(emmeans(m2, ~ disease_f | genotype)) (package emmeans). Code page sections 11 and 12; slide The drought effect in each variety: emmeans.

Q4. In summary() of your Q3 model, genotypekerman has p < 0.001. In anova() of the same model, the genotype row has p = 0.11. Both are correct. Explain in two or three sentences what each one is testing.

Hint

Which trees does the genotypekerman coefficient compare, given that disease_f = 0 is the reference level? Code page section 11b; slides Two tables, two questions, What each coefficient compares, and Change the reference level.

Q5. Is there evidence that the disease affects the genotypes differently? Which genotype would you recommend if the disease is expected to stay in the orchard, and which if it could be prevented (for example, with a fungicide)? Two or three sentences, using your Q3 numbers.

Hint

The genotype:disease_f row of anova(m2) tests the interaction. For the recommendation, compare the emmeans(m2, ~ disease_f | genotype) means for healthy and for diseased trees. Slide Is the interaction real?

Part 3: Disease Rates

Q6. What proportion of each genotype is diseased? Fit a logistic regression, glm(disease ~ genotype, family = binomial). Using plogis() on the coefficients, calculate the model’s probability of disease for a chaparillo tree and for a t41 tree, and check them against your proportions.

Hint

prop.table(table(pistachio$genotype, pistachio$disease), 1); chaparillo is the reference: plogis(coef(g)[1]); t41: plogis(coef(g)[1] + coef(g)[3]). Code page section 19; slides A group predictor: glm(rust ~ water) and Back to probabilities: plogis().

Q7. Calculate the odds ratio of disease for t41 relative to chaparillo, with its 95% confidence interval. Write it as one sentence for a results section.

Hint

exp(coef(g)) and exp(confint(g)). Code page section 19; slides Odds ratio for a group and Reporting a GLM.

Part 4: An Early Warning from Drones

The disease was scored at the end of the season, in September. Could the drone images tell a farmer which trees are diseased earlier, while there is still time to treat them? Treatment works better the earlier it starts, but it is expensive, so it can’t be applied to every tree. The farmer needs to know when the NDVI data become reliable enough to act on.

Q8. Fit glm(disease ~ NDVI_march, family = binomial). Classify each tree as diseased if its fitted probability is above 0.5 and calculate the proportion classified correctly. Compare it with the accuracy of simply guessing “diseased” for every tree. Is March NDVI useful?

Hint

pred <- as.numeric(fitted(model) > 0.5); mean(pred == pistachio$disease); the baseline is mean(pistachio$disease == 1). Code page section 22; slide How often is it right?

Q9. Add June NDVI (disease ~ NDVI_march + NDVI_june), then August NDVI as well, and calculate the accuracy of each model. In the March + June model, the two slopes have opposite signs. What kind of tree does the model flag as likely diseased?

Hint

Same accuracy code as Q8 for each model. For the signs: a positive March slope and a negative June slope mean a high March NDVI raises the log-odds of disease and a high June NDVI lowers them. What does a tree that is green in March but brown in June look like? Optional check: pistachio$change <- pistachio$NDVI_june - pistachio$NDVI_march, then fit disease ~ change.

Q10. After which month would you use the NDVI predictions to choose trees to treat? Two or three sentences weighing accuracy against how early you can act.

Hint

No new code: compare your Q8 and Q9 accuracies with each other and with the baseline. There is no single right answer; justify yours.

--- title: "Problem Set 3: General(ized) Linear Models" 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 19, 1:10 PM (start of class) — submitted through Canvas. Parts 1–2 use the first half of the week (groups as predictors, interactions); Parts 3–4 use the second half (logistic regression). The [week's code page](Week3_GLM.html) has working code for every step; the hints point to its sections. ::: --- ## Submission Instructions Upload **two files** to Canvas: 1. **Written answers** — a `.pdf` document with your written responses 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:** `HW3emailID.R` (code) and `HW3emailID.pdf` (written answers), where `emailID` is the part of your UC Davis email before the `@`. > **Example:** Grey Monroe would submit `HW3gmonroe.R` and `HW3gmonroe.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: disease and yield in pistachio Pistachios are an important crop for California agriculture. As in most crops, genotypes can differ in traits like yield and susceptibility to disease. This simulated dataset has 300 trees of three pistachio genotypes ("kerman", "chaparillo", and "t41"; 100 trees each) planted in a common garden. At the end of the growing season (September), each tree was scored for disease (`disease`: 1 = diseased, 0 = healthy; a disease that causes leaf die-off) and its yield was measured (`yield`: total harvested nuts). Drones also measured each tree's greenness (NDVI) in March, June, and August (`NDVI_march`, `NDVI_june`, `NDVI_august`); larger NDVI means a greener tree. We want to know (a) whether yield differs between genotypes, (b) whether the disease affects yield, and differently in different genotypes, (c) whether rates of disease differ between genotypes, and (d) whether the drone data can flag diseased trees before the end of the season. Read it straight from the course website: ```{r} pistachio <- read.csv("https://greymonroe.github.io/PLS_206/data/pistachio.csv") str(pistachio) ``` --- ## Part 1: Genotype and Yield **Q1.** Plot yield by genotype (a box or violin plot with the points on top). Fit `lm(yield ~ genotype)`. In one sentence each: what does the `(Intercept)` estimate mean, and what does the `genotypekerman` estimate mean? Give both in yield units. ::: {.callout-tip collapse="true"} ## Hint `ggplot(pistachio, aes(genotype, yield)) + geom_boxplot() + geom_jitter(width = 0.1)`; `summary(lm(yield ~ genotype, data = pistachio))`; `tapply(pistachio$yield, pistachio$genotype, mean)` to check. Code page sections 2 and 6; slides *Two groups, one line* and *Three groups, two dummies*. ::: **Q2.** Do the genotypes differ in yield? Answer with one test of the whole factor (`anova()`), then say which pairs (if any) differ with `TukeyHSD()`. ::: {.callout-tip collapse="true"} ## Hint `anova(model)`; `TukeyHSD(aov(model))`. Code page sections 7 and 8; slides *Do the varieties differ at all? anova()* and *Which pairs? TukeyHSD()*. ::: --- ## Part 2: Disease, Yield, and Genotype **Q3.** Plot yield by genotype with a separate box for diseased and healthy trees (`fill = factor(disease)`). Fit `lm(yield ~ genotype * factor(disease))`. Use `emmeans` to get the yield lost to disease **within each genotype**. Which genotype loses the least, and how much does it lose? ::: {.callout-tip collapse="true"} ## Hint `pistachio$disease_f <- factor(pistachio$disease)` first; then `m2 <- lm(yield ~ genotype * disease_f, data = pistachio)` and `pairs(emmeans(m2, ~ disease_f | genotype))` (package `emmeans`). Code page sections 11 and 12; slide *The drought effect in each variety: emmeans*. ::: **Q4.** In `summary()` of your Q3 model, `genotypekerman` has p < 0.001. In `anova()` of the same model, the `genotype` row has p = 0.11. Both are correct. Explain in two or three sentences what each one is testing. ::: {.callout-tip collapse="true"} ## Hint Which trees does the `genotypekerman` coefficient compare, given that `disease_f` = 0 is the reference level? Code page section 11b; slides *Two tables, two questions*, *What each coefficient compares*, and *Change the reference level*. ::: **Q5.** Is there evidence that the disease affects the genotypes differently? Which genotype would you recommend if the disease is expected to stay in the orchard, and which if it could be prevented (for example, with a fungicide)? Two or three sentences, using your Q3 numbers. ::: {.callout-tip collapse="true"} ## Hint The `genotype:disease_f` row of `anova(m2)` tests the interaction. For the recommendation, compare the `emmeans(m2, ~ disease_f | genotype)` means for healthy and for diseased trees. Slide *Is the interaction real?* ::: --- ## Part 3: Disease Rates **Q6.** What proportion of each genotype is diseased? Fit a logistic regression, `glm(disease ~ genotype, family = binomial)`. Using `plogis()` on the coefficients, calculate the model's probability of disease for a chaparillo tree and for a t41 tree, and check them against your proportions. ::: {.callout-tip collapse="true"} ## Hint `prop.table(table(pistachio$genotype, pistachio$disease), 1)`; chaparillo is the reference: `plogis(coef(g)[1])`; t41: `plogis(coef(g)[1] + coef(g)[3])`. Code page section 19; slides *A group predictor: glm(rust ~ water)* and *Back to probabilities: plogis()*. ::: **Q7.** Calculate the odds ratio of disease for t41 relative to chaparillo, with its 95% confidence interval. Write it as one sentence for a results section. ::: {.callout-tip collapse="true"} ## Hint `exp(coef(g))` and `exp(confint(g))`. Code page section 19; slides *Odds ratio for a group* and *Reporting a GLM*. ::: --- ## Part 4: An Early Warning from Drones The disease was scored at the end of the season, in September. Could the drone images tell a farmer which trees are diseased **earlier**, while there is still time to treat them? Treatment works better the earlier it starts, but it is expensive, so it can't be applied to every tree. The farmer needs to know when the NDVI data become reliable enough to act on. **Q8.** Fit `glm(disease ~ NDVI_march, family = binomial)`. Classify each tree as diseased if its fitted probability is above 0.5 and calculate the proportion classified correctly. Compare it with the accuracy of simply guessing "diseased" for every tree. Is March NDVI useful? ::: {.callout-tip collapse="true"} ## Hint `pred <- as.numeric(fitted(model) > 0.5)`; `mean(pred == pistachio$disease)`; the baseline is `mean(pistachio$disease == 1)`. Code page section 22; slide *How often is it right?* ::: **Q9.** Add June NDVI (`disease ~ NDVI_march + NDVI_june`), then August NDVI as well, and calculate the accuracy of each model. In the March + June model, the two slopes have opposite signs. What kind of tree does the model flag as likely diseased? ::: {.callout-tip collapse="true"} ## Hint Same accuracy code as Q8 for each model. For the signs: a positive March slope and a negative June slope mean a high March NDVI raises the log-odds of disease and a high June NDVI lowers them. What does a tree that is green in March but brown in June look like? Optional check: `pistachio$change <- pistachio$NDVI_june - pistachio$NDVI_march`, then fit `disease ~ change`. ::: **Q10.** After which month would you use the NDVI predictions to choose trees to treat? Two or three sentences weighing accuracy against how early you can act. ::: {.callout-tip collapse="true"} ## Hint No new code: compare your Q8 and Q9 accuracies with each other and with the baseline. There is no single right answer; justify yours. :::