Exam II Review

Weeks 7–11 | AGST 50104 Experimental Design

Dr. Samuel B Fernandes

2026-01-01

What You Should Be Able to Do

  • Factorial: identify when interaction changes interpretation; run Type III SS correctly
  • Split-plot: determine whole-plot vs. sub-plot factor; name the two error strata
  • Strip-plot: name the three error strata; fit the correct Error() in R
  • Mixed models: distinguish fixed from random effects; fit with lme4
  • ASReml-R: use asreml() for split-plot and repeated measures; read wald() output

Complex Factorial Designs

Interaction: When It Changes Everything

When interaction is NOT significant:

  • Interpret main effects directly
  • Effects of A are the same at all levels of B

\[Y_{ijk} = \mu + \alpha_i + \beta_j + \rho_k + \varepsilon_{ijk}\]

When interaction IS significant:

  • Do not report main effects alone
  • Use simple effects: effect of A within each level of B

\[Y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \rho_k + \varepsilon_{ijk}\]

Warning

Most common mistake: Reporting main effects and ignoring a significant interaction.

Complex Factorials: R Analysis Pattern

ExpDes::ex4 — 4 rotation speeds × 2 manure treatments (CRD, 3 reps per cell). Response: leaf nitrogen content (%).

library(dplyr); library(car); library(emmeans)
library(ExpDes)

data(ex4)
ex4 <- ex4 |> mutate(revol = factor(revol))

# Step 1: sum-to-zero contrasts for correct Type III SS
options(contrasts = c("contr.sum", "contr.poly"))
m_fac <- lm(n ~ revol * esterco, data = ex4)
car::Anova(m_fac, type = "III")
#> Anova Table (Type III tests)
#> 
#> Response: n
#>               Sum Sq Df  F value    Pr(>F)    
#> (Intercept)   66.667  1 815.3282 3.729e-15 ***
#> revol          0.956  3   3.8993   0.02883 *  
#> esterco        0.101  1   1.2401   0.28190    
#> revol:esterco  0.314  3   1.2783   0.31555    
#> Residuals      1.308 16                       
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
options(contrasts = c("contr.treatment", "contr.poly"))  # always restore!

Simple Effects Follow-Up

# Interaction NOT significant → report main effects
emmeans(m_fac, ~ revol) |> pairs(adjust = "tukey")
#>  contrast          estimate    SE df t.ratio p.value
#>  revol5 - revol10     0.215 0.165 16   1.302  0.5746
#>  revol5 - revol15     0.185 0.165 16   1.121  0.6826
#>  revol5 - revol20    -0.287 0.165 16  -1.736  0.3383
#>  revol10 - revol15   -0.030 0.165 16  -0.182  0.9978
#>  revol10 - revol20   -0.502 0.165 16  -3.039  0.0354
#>  revol15 - revol20   -0.472 0.165 16  -2.857  0.0504
#> 
#> Results are averaged over the levels of: esterco 
#> P value adjustment: tukey method for comparing a family of 4 estimates
library(ggplot2)
ex4 |>
  group_by(revol, esterco) |>
  summarise(mean_n = mean(n), se = sd(n)/sqrt(n()), .groups = "drop") |>
  ggplot(aes(x = revol, y = mean_n, color = esterco, group = esterco)) +
  geom_line(linewidth = 1) + geom_point(size = 3) +
  geom_errorbar(aes(ymin = mean_n - se, ymax = mean_n + se), width = 0.2) +
  labs(x = "Rotation speed (rpm)", y = "Leaf N (%)",
       color = "Manure", title = "Non-significant interaction → parallel lines") +
  theme_minimal(base_size = 14)

Tip

Rule: Non-significant interaction → report main effects. Significant → use ~ A | B in emmeans(), not ~ A alone.

Split-Plot & Strip-Plot Designs

agridat::yates.oats — Rothamsted 1931. 3 oat varieties (whole-plot) × 4 nitrogen rates (sub-plot), 6 blocks.

Split-Plot: Design Decision

The key question:

Which factor is hard to change in the field?

Scenario Whole-plot Sub-plot
Central pivot irrigation Irrigation freq. N rate
Tillage × herbicide Tillage Herbicide
Soil fumigation × variety Fumigation Variety

The hard-to-change factor gets the larger experimental unit.

Two-stage randomization:

  1. Randomly assign WP factor to whole-plots within each block
  2. Randomly assign SP factor to sub-plots within each whole-plot

\(\Rightarrow\) Two separate randomizations = two separate error terms.

Split-Plot: Two Error Strata

\[Y_{ijk} = \mu + \rho_k + \alpha_i + \underbrace{\delta_{ik}}_{\text{WP error}} + \beta_j + (\alpha\beta)_{ij} + \underbrace{\varepsilon_{ijk}}_{\text{SP error}}\]

Source Tested against df (example: 4 blocks, 2 WP, 3 SP)
Block — \(b-1 = 3\)
Whole-plot (A) Block × A (\(\sigma^2_{WP}\)) \(a-1 = 1\)
Sub-plot (B) Residual (\(\sigma^2_{SP}\)) \(b-1 = 2\)
A × B Residual (\(\sigma^2_{SP}\)) \((a-1)(b-1) = 2\)

Warning

⚠ Never test A against the residual MS in a split-plot. Whole-plot error is always larger.

Split-Plot: R Code

Classical ANOVA (two error strata):

library(agridat)
dat <- yates.oats |>
  mutate(nitro = factor(nitro),
         gen   = factor(gen),
         block = factor(block))

# gen = whole-plot; nitro = sub-plot
model_sp <- aov(yield ~ gen * nitro +
                  Error(block / gen),
                data = dat)
summary(model_sp)
#> 
#> Error: block
#>           Df Sum Sq Mean Sq F value Pr(>F)
#> Residuals  5  15875    3175               
#> 
#> Error: block:gen
#>           Df Sum Sq Mean Sq F value Pr(>F)
#> gen        2   1786   893.2   1.485  0.272
#> Residuals 10   6013   601.3               
#> 
#> Error: Within
#>           Df Sum Sq Mean Sq F value   Pr(>F)    
#> nitro      3  20020    6673  37.686 2.46e-12 ***
#> gen:nitro  6    322      54   0.303    0.932    
#> Residuals 45   7969     177                     
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mixed-model for EMMs:

library(lme4); library(emmeans); library(multcomp)

model_lme <- lmer(yield ~ gen * nitro +
                    (1 | block / gen),
                  data = dat, REML = TRUE)

emm_sp <- emmeans(model_lme, ~ nitro)
multcomp::cld(emm_sp, Letters = letters,
              adjust = "tukey")
#>  nitro emmean   SE   df lower.CL upper.CL .group
#>  0       79.4 7.17 6.79     55.3      103  a    
#>  0.2     98.9 7.17 6.79     74.8      123   b   
#>  0.4    114.2 7.17 6.79     90.2      138    c  
#>  0.6    123.4 7.17 6.79     99.3      147    c  
#> 
#> Results are averaged over the levels of: gen 
#> Degrees-of-freedom method: kenward-roger 
#> Confidence level used: 0.95 
#> Conf-level adjustment: sidak method for 4 estimates 
#> P value adjustment: tukey method for comparing a family of 4 estimates 
#> significance level used: alpha = 0.05 
#> NOTE: If two or more means share the same grouping symbol,
#>       then we cannot show them to be different.
#>       But we also did not show them to be the same.

Strip-Plot: Three Error Strata

\[Y_{ijk} = \mu + \rho_k + \alpha_i + \underbrace{\delta_{ik}}_{\text{Block×H}} + \beta_j + \underbrace{\gamma_{jk}}_{\text{Block×V}} + (\alpha\beta)_{ij} + \underbrace{\varepsilon_{ijk}}_{\text{Residual}}\]

Source Tested against Notes
Horizontal strips (H) Block × H error Strips run in one direction
Vertical strips (V) Block × V error Strips run perpendicularly
H × V interaction Residual Intersection plots

Warning

Strip-plot ≠ split-plot. Both factors have their own block interaction as error. The interaction is tested against the residual.

Strip-Plot: R Code

Classical ANOVA:

# H = horizontal strips, V = vertical strips  
model_strip <- aov(
  y ~ H * V +
    Error(Block + Block:H + Block:V),
  data = d)
summary(model_strip)

Three error strata explicitly named.

Mixed-model (for EMMs):

library(lme4)
# Three random terms = three error strata
model_strip_lme <- lmer(
  y ~ H * V +
    (1 | Block) +       # block main effect
    (1 | Block:H) +     # WP error for H
    (1 | Block:V),      # WP error for V
  data = d, REML = TRUE)

Design at a Glance

Factorial (RCBD) Split-Plot Strip-Plot
EU One plot = one combination Whole-plot for A; sub-plot for A×B Column strip for H; row strip for V; cell for H×V
Error strata 1 (residual) 2 (WP error, SP error) 3 (Block×H, Block×V, residual)
R Error() term Error(block) Error(block/A) Error(block + block:H + block:V)
When to use All treatments equally randomizable One factor hard to change in field Both factors hard to change in different directions
Precision Equal for all effects Subplot factors more precise than WP H and V factors each have own precision
Common example Herbicide × dose RCBD Irrigation × N-rate with pivot Tillage direction × herbicide band

Fixed vs. Random: The Decision

Fixed effects

  • Specific levels you chose and will repeat
  • Your research question is about these levels
  • Estimable: treatment means, contrasts

Examples: always depends on the context, but often: fertilizer type, irrigation frequency, tillage method

Random effects

  • Random sample from a larger population
  • Your question is about variance, not specific levels
  • Estimable: variance component (\(\sigma^2_b\))

Examples: always depends on the context, but often: block, field location, year, individual plant/plot

In split-plot (ANOVA framework using Error()): Block and Block×A are always random. Treatments (A, B) are always fixed.

lme4: Key Syntax Patterns

library(lme4); library(emmeans)

dat <- yates.oats |>
  mutate(nitro = factor(nitro), gen = factor(gen), block = factor(block))

# 1. RCBD: block as random intercept (equivalent to classical RCBD ANOVA)
m_rcbd <- lmer(yield ~ gen * nitro + (1 | block),
               data = dat, REML = TRUE)

# 2. Split-plot: block/gen nesting = WP error + SP error
m_sp <- lmer(yield ~ gen * nitro + (1 | block / gen),
             data = dat, REML = TRUE)

(1 | block / gen) expands to (1 | block) + (1 | block:gen) — one random intercept per block, one per whole-plot. This is the split-plot error structure in mixed-model form.

ASReml-R vs lme4

Use lme4 when:

  • Random intercepts and slopes
  • Nested / crossed random effects
  • Standard split-plot or RCBD

ASReml-R needed when:

  • Residual covariance structures (AR1, CS, diagonal)
  • Spatial row × column adjustment
  • Heterogeneous residual variances by environment

ASReml-R: Split-Plot Example

library(asreml)

# Split-plot: irrigation = whole-plot, cv = sub-plot
split_model <- asreml(
  fixed  = resp ~ irrigation * cv,   # fixed effects
  random = ~ block + block:irrigation, # WP error terms
  data   = alfaSPLIT
)

# If convergence issues:
split_model <- update.asreml(split_model)

# F-tests for fixed effects (equivalent to ANOVA table)
wald(split_model)

# Variance components: WP error, SP error
summary(split_model)$varcomp

random = ~ block + block:irrigation encodes the two random error strata: one for blocks, one for the Block×Irrigation whole-plot error.

ASReml-R: Repeated Measures

#> Online License checked out Wed Apr  1 14:24:37 2026
# Independent (baseline — no correlation)
m_id <- asreml(y ~ Tmt * Time,
               data = grassUV,
               trace = F)

# Compound Symmetry: constant correlation across all time points
m_cs <- asreml(y ~ Tmt * Time,
               residual = ~ id(Plant):cor(Time),
               data = grassUV,
               trace = F)

# AR(1): correlation decays with time lag
m_ar1 <- asreml(y ~ Tmt * Time,
                residual = ~ id(Plant):corh(Time),
                data = grassUV,
               trace = F)

Choose the model with lowest AIC/BIC. CS assumes equal correlation; AR(1) assumes closer time points are more correlated.

# Compare models with AIC/BIC
bind_rows(
  asremlPlus::infoCriteria(m_id),
  asremlPlus::infoCriteria(m_cs),
  asremlPlus::infoCriteria(m_ar1)
)
#>   fixedDF varDF NBound      AIC      BIC    loglik
#> 1       0     1      0 420.8836 422.9779 -209.4418
#> 2       0     2      0 397.7535 401.9422 -196.8768
#> 3       0     6      1 364.7580 377.3241 -176.3790
lrt(m_id, m_cs)  # CS better than independent
#> Likelihood ratio test(s) assuming nested random models.
#> (See Self & Liang, 1987)
#> 
#>           df LR-statistic Pr(Chisq)    
#> m_cs/m_id  1        25.13  2.68e-07 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Exam II: Key Skills Checklist

Step 1 — Identify the design:

☐ How many factors? Are any hard to change?
☐ Factorial → one error term
☐ Split-plot → two error strata
☐ Strip-plot → three error strata

Step 2 — Write the model:

☐ List all fixed effects
☐ List all random effects (e.g.,WP error)
☐ Identify the correct denominator for each F-test ☐ ANOVA vs mixed-model approach

Step 3 — Fit in R:

☐ Factorial: lm() + car::Anova(type="III")
☐ Split-plot: aov(Error(block/A)) or lmer((1|block/A))
☐ Strip-plot: aov(Error(block + block:H + block:V))
☐ Mixed: lmer() or asreml() → wald()

Step 4 — Interpret output:

☐ Check interaction first
☐ If interaction significant → simple effects only
☐ Report variance components for random effects
☐ Adjusted means + Tukey for significant effects

If in doubt: trace back to randomization. How were treatments assigned? At what level? → That determines the error term.