Chapter 6
Charm Meskin
2026-10-28
rm(list = ls()) # clean startlibrary(dplyr)
library(ggplot2)# Diagnostic panel for any lm object: Residuals vs Fitted, Normal Q-Q,
# Scale-Location, Residuals vs Leverage. Base graphics, so no extra packages.
diag_plots <- function(model) {
op <- par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
on.exit(par(op)) # restore settings even if plot() errors
plot(model)
}PART 1: The two-way ANOVA
1.1 Import and prepare
growth.moo <- read.csv("growth.csv")
# read.csv leaves these as <chr> under R >= 4.0; they must be factors
growth.moo <- growth.moo %>%
mutate(diet = as.factor(diet),
supplement = as.factor(supplement))
glimpse(growth.moo)## Rows: 48
## Columns: 3
## $ supplement <fct> supergain, supergain, supergain, supergain, control, contro…
## $ diet <fct> wheat, wheat, wheat, wheat, wheat, wheat, wheat, wheat, whe…
## $ gain <dbl> 17.37125, 16.81489, 18.08184, 15.78175, 17.70656, 18.22717,…1.2 Levels and the reference level
levels() shows the order. Alphabetical order puts
agrimore first, so control must be set
explicitly or every contrast in summary() is read against
the wrong baseline.
levels(growth.moo$diet) # reference = barley (first alphabetically)## [1] "barley" "oats" "wheat"levels(growth.moo$supplement)## [1] "agrimore" "control" "supergain" "supersupp"growth.moo <- growth.moo %>%
mutate(supplement = relevel(supplement, ref = "control"))
stopifnot(levels(growth.moo$supplement)[1] == "control")
levels(growth.moo$supplement)## [1] "control" "agrimore" "supergain" "supersupp"1.3 Summary data: all 12 combinations
sumMoo <- growth.moo %>%
group_by(diet, supplement) %>%
summarise(meanGrow = mean(gain))
sumMoo## # A tibble: 12 × 3
## # Groups: diet [3]
## diet supplement meanGrow
## <fct> <fct> <dbl>
## 1 barley control 23.3
## 2 barley agrimore 26.3
## 3 barley supergain 22.5
## 4 barley supersupp 25.6
## 5 oats control 20.5
## 6 oats agrimore 23.3
## 7 oats supergain 19.7
## 8 oats supersupp 21.9
## 9 wheat control 17.4
## 10 wheat agrimore 19.6
## 11 wheat supergain 17.0
## 12 wheat supersupp 19.71.4 Interaction plot
Points only (Figure 6.1), then join by diet (Figure
6.2). group = diet draws the lines;
colour = diet colours them.
ggplot(sumMoo, aes(x = supplement, y = meanGrow)) +
geom_point() +
theme_bw()
Figure 6.1 Cow weight gain, unadorned.
ggplot(sumMoo, aes(x = supplement, y = meanGrow,
colour = diet, group = diet)) +
geom_point() +
geom_line() +
theme_bw()
Figure 6.2 Does the effect of supplement type on cow weight gain depend on the diet?
Parallel lines indicate no interaction. Here they are not parallel, so test the interaction.
1.5 The model
* expands to
diet + supplement + diet:supplement.
model_cow <- lm(gain ~ diet * supplement, data = growth.moo)1.6 Assumptions
diag_plots(model_cow)
Figure 6.3 Diagnostic plots for the cow weight gain two-way ANOVA.
Look for: no pattern in residuals vs fitted; near-flat scale–location; no point both high-leverage and large-residual; an approximately straight Q–Q line. Mild departures from normality are tolerable.
1.7 anova() table
anova(model_cow)## Analysis of Variance Table
##
## Response: gain
## Df Sum Sq Mean Sq F value Pr(>F)
## diet 2 287.171 143.586 83.5201 2.999e-14 ***
## supplement 3 91.881 30.627 17.8150 2.952e-07 ***
## diet:supplement 6 3.406 0.568 0.3302 0.9166
## Residuals 36 61.890 1.719
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Variation explained by diet, then by
supplement after diet, then whether letting
the supplement effect vary by diet explains additional variation (the
interaction row). Each F is a Mean Sq divided by the residual Mean
Sq.
The row order is irrelevant only because this is a balanced,
orthogonal, fully factorial designed experiment. For unbalanced or
observational data, order of testing matters — anova()
alone won’t do.
1.8 summary() table
summary(model_cow)##
## Call:
## lm(formula = gain ~ diet * supplement, data = growth.moo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.48756 -1.00368 -0.07452 1.03496 2.68069
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 23.2966499 0.6555863 35.536 < 2e-16 ***
## dietoats -2.8029851 0.9271390 -3.023 0.00459 **
## dietwheat -5.8911317 0.9271390 -6.354 2.34e-07 ***
## supplementagrimore 3.0518277 0.9271390 3.292 0.00224 **
## supplementsupergain -0.8305263 0.9271390 -0.896 0.37631
## supplementsupersupp 2.2786527 0.9271390 2.458 0.01893 *
## dietoats:supplementagrimore -0.2471088 1.3111726 -0.188 0.85157
## dietwheat:supplementagrimore -0.8182729 1.3111726 -0.624 0.53651
## dietoats:supplementsupergain -0.0001351 1.3111726 0.000 0.99992
## dietwheat:supplementsupergain 0.4374395 1.3111726 0.334 0.74060
## dietoats:supplementsupersupp -0.9120830 1.3111726 -0.696 0.49113
## dietwheat:supplementsupersupp -0.0158299 1.3111726 -0.012 0.99043
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.311 on 36 degrees of freedom
## Multiple R-squared: 0.8607, Adjusted R-squared: 0.8182
## F-statistic: 20.22 on 11 and 36 DF, p-value: 3.295e-12(Intercept) is the reference combination
(barley–control); every other row is a difference from it (treatment
contrasts). Interaction terms are non-significant, agreeing with the
ANOVA table. A priori contrasts would be better practice than reading
this table; multcomp, contrast and
rms support that (not used here).
1.9 Mean ± SE figure
sumMoo <- growth.moo %>%
group_by(diet, supplement) %>%
summarise(meanGrow = mean(gain),
# n() counts rows including NAs; use sum(!is.na(gain)) if any are present
seGrow = sd(gain) / sqrt(n()))
sumMoo## # A tibble: 12 × 4
## # Groups: diet [3]
## diet supplement meanGrow seGrow
## <fct> <fct> <dbl> <dbl>
## 1 barley control 23.3 0.703
## 2 barley agrimore 26.3 0.919
## 3 barley supergain 22.5 0.771
## 4 barley supersupp 25.6 1.06
## 5 oats control 20.5 0.506
## 6 oats agrimore 23.3 0.613
## 7 oats supergain 19.7 0.349
## 8 oats supersupp 21.9 0.413
## 9 wheat control 17.4 0.460
## 10 wheat agrimore 19.6 0.710
## 11 wheat supergain 17.0 0.485
## 12 wheat supersupp 19.7 0.475geom_errorbar() needs its own
ymin/ymax; width sets the cap
size.
ggplot(sumMoo, aes(x = supplement, y = meanGrow,
colour = diet, group = diet)) +
geom_point() +
geom_line() +
geom_errorbar(aes(ymin = meanGrow - seGrow,
ymax = meanGrow + seGrow), width = 0.1) +
theme_bw()
Figure 6.4 A quite nice graph of the cow growth data, finishing off our two-way ANOVA example.
PART 2: Analysis of covariance (ANCOVA)
One categorical explanatory variable (SEASON) plus one
continuous (DENSITY); response EGGS. Data from
Quinn & Keough (2002).
2.1 Import and prepare
SEASON must be a factor — character data makes
levels() return NULL, which
expand.grid() silently drops and predict()
then fails on.
limp <- read.csv("limpet.csv")
limp <- limp %>% mutate(SEASON = as.factor(SEASON))
stopifnot(is.factor(limp$SEASON), !is.null(levels(limp$SEASON)))
glimpse(limp)## Rows: 24
## Columns: 3
## $ DENSITY <int> 8, 8, 8, 8, 8, 8, 15, 15, 15, 15, 15, 15, 30, 30, 30, 30, 30, …
## $ SEASON <fct> spring, spring, spring, summer, summer, summer, spring, spring…
## $ EGGS <dbl> 2.875, 2.625, 1.750, 2.125, 1.500, 1.875, 2.600, 1.866, 2.066,…levels(limp$SEASON) # spring, summer — spring is the reference level## [1] "spring" "summer"2.2 Always start with a picture
ggplot(limp, aes(x = DENSITY, y = EGGS, colour = SEASON)) +
geom_point() +
scale_colour_manual(values = c(spring = "green", summer = "red")) +
theme_bw()
Figure 6.5 Always begin analyses with a picture. Egg production by limpets in two seasons and at four limpet densities.
A straight line is \(y = b + m x\): intercept \(b\), slope \(m\). Ask whether the intercept (eggs at zero density) differs between seasons, and whether the slope (density dependence) differs. Candidate models run from no effect → common slope with different intercepts (parallel, additive) → different intercepts and slopes (interaction).
Hypothesis: the effect of density on egg production depends on the season.
# Figure 6.6 competing models, for reference only
lm(EGGS ~ 1, data = limp) # no effect
lm(EGGS ~ SEASON, data = limp) # different intercepts, no slope
lm(EGGS ~ DENSITY, data = limp) # one common line
lm(EGGS ~ DENSITY + SEASON, data = limp) # parallel lines (additive)
lm(EGGS ~ DENSITY * SEASON, data = limp) # different intercepts and slopes2.3 The model
limp.mod <- lm(EGGS ~ DENSITY * SEASON, data = limp)* gives the main effect of density, the main effect of
season, and the possibility that the density effect depends on season.
lm() stores a lot on the fitted object — see
names(limp.mod) and ?lm (Value section).
names(limp.mod)## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "contrasts" "xlevels" "call" "terms"
## [13] "model"2.4 Assumptions
diag_plots(limp.mod)
Figure 6.7 Diagnostic graphs for the linear model of the limpet data.
2.5 anova() table
anova(limp.mod)## Analysis of Variance Table
##
## Response: EGGS
## Df Sum Sq Mean Sq F value Pr(>F)
## DENSITY 1 5.0241 5.0241 30.1971 2.226e-05 ***
## SEASON 1 3.2502 3.2502 19.5350 0.0002637 ***
## DENSITY:SEASON 1 0.0118 0.0118 0.0711 0.7925333
## Residuals 20 3.3275 0.1664
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Sequential sums of squares: common slope (DENSITY, 1
df), different intercepts (SEASON, 1 df), then varying
slopes (DENSITY:SEASON, 1 df). The interaction row is the
formal test of whether the density effect depends on season.
2.6 summary() table
summary(limp.mod)##
## Call:
## lm(formula = EGGS ~ DENSITY * SEASON, data = limp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.65468 -0.25021 -0.03318 0.28335 0.57532
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.664166 0.234118 11.380 3.45e-10 ***
## DENSITY -0.033650 0.008259 -4.074 0.000591 ***
## SEASONsummer -0.812282 0.331092 -2.453 0.023450 *
## DENSITY:SEASONsummer 0.003114 0.011680 0.267 0.792533
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4079 on 20 degrees of freedom
## Multiple R-squared: 0.7135, Adjusted R-squared: 0.6705
## F-statistic: 16.6 on 3 and 20 DF, p-value: 1.186e-05(Intercept) and DENSITY are the spring line
(spring sorts first). SEASONsummer is the difference in
intercepts and DENSITY:SEASONsummer the difference in
slopes, so:
coef(limp.mod)## (Intercept) DENSITY SEASONsummer
## 2.664166462 -0.033649651 -0.812282493
## DENSITY:SEASONsummer
## 0.003113571# EGGS_spring = 2.66 - 0.033 * DENSITY
# EGGS_summer = (2.66 - 0.812) + (-0.033 + 0.003) * DENSITY2.7 Prediction grid
Build the x values with expand.grid(), then get y and
95% CIs with predict().
expand.grid(FIRST = c("A", "B"), SECOND = c(1, 2)) # first factor varies fastest## FIRST SECOND
## 1 A 1
## 2 B 1
## 3 A 2
## 4 B 2predict(limp.mod) # fitted values at the observed x## 1 2 3 4 5 6 7 8
## 2.3949692 2.3949692 2.3949692 1.6075953 1.6075953 1.6075953 2.1594217 2.1594217
## 9 10 11 12 13 14 15 16
## 2.1594217 1.3938428 1.3938428 1.3938428 1.6546769 1.6546769 1.6546769 0.9358016
## 17 18 19 20 21 22 23 24
## 0.9358016 0.9358016 1.1499321 1.1499321 1.1499321 0.4777604 0.4777604 0.4777604new.x <- expand.grid(
DENSITY = seq(from = 8, to = 45, length.out = 10),
SEASON = levels(limp$SEASON))
# guard: expand.grid() drops NULL arguments without warning
stopifnot(all(c("DENSITY", "SEASON") %in% names(new.x)),
is.factor(new.x$SEASON),
nrow(new.x) == 20)
str(new.x) # 20 obs. of 2 variables## 'data.frame': 20 obs. of 2 variables:
## $ DENSITY: num 8 12.1 16.2 20.3 24.4 ...
## $ SEASON : Factor w/ 2 levels "spring","summer": 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "out.attrs")=List of 2
## ..$ dim : Named int [1:2] 10 2
## .. ..- attr(*, "names")= chr [1:2] "DENSITY" "SEASON"
## ..$ dimnames:List of 2
## .. ..$ DENSITY: chr [1:10] "DENSITY= 8.00000" "DENSITY=12.11111" "DENSITY=16.22222" "DENSITY=20.33333" ...
## .. ..$ SEASON : chr [1:2] "SEASON=spring" "SEASON=summer"new.y <- predict(limp.mod, newdata = new.x, interval = "confidence")
head(new.y) # fit, lwr, upr## fit lwr upr
## 1 2.394969 2.019285 2.770654
## 2 2.256632 1.931230 2.582034
## 3 2.118294 1.834274 2.402315
## 4 1.979957 1.724062 2.235852
## 5 1.841619 1.595998 2.087241
## 6 1.703282 1.447918 1.958646# bind them and rename fit -> EGGS so the names match the raw data;
# that name match is what lets geom_smooth() inherit the x/y aesthetics
addThese <- data.frame(new.x, new.y) %>% rename(EGGS = fit)
head(addThese)## DENSITY SEASON EGGS lwr upr
## 1 8.00000 spring 2.394969 2.019285 2.770654
## 2 12.11111 spring 2.256632 1.931230 2.582034
## 3 16.22222 spring 2.118294 1.834274 2.402315
## 4 20.33333 spring 1.979957 1.724062 2.235852
## 5 24.44444 spring 1.841619 1.595998 2.087241
## 6 28.55556 spring 1.703282 1.447918 1.958646Optional: hold DENSITY at its mean to get predicted eggs
per season at average density — the same recipe generalises.
new.x2 <- expand.grid(DENSITY = mean(limp$DENSITY),
SEASON = levels(limp$SEASON))
data.frame(new.x2, predEgg = predict(limp.mod, newdata = new.x2))## DENSITY SEASON predEgg
## 1 24.5 spring 1.83975
## 2 24.5 summer 1.103752.8 Final figure
DENSITY and EGGS are inherited from the
first ggplot() call because addThese uses the
same names. stat = "identity" tells
geom_smooth() to draw the supplied values instead of
fitting its own smoother.
ggplot(limp, aes(x = DENSITY, y = EGGS, colour = SEASON)) +
geom_point(size = 5) +
geom_smooth(data = addThese,
aes(ymin = lwr, ymax = upr, fill = SEASON),
stat = "identity") +
scale_colour_manual(values = c(spring = "green", summer = "red")) +
scale_fill_manual(values = c(spring = "green", summer = "red")) +
theme_bw()
Figure 6.8 Fitted lines for both seasons, added to the raw data for egg production by limpets.