Alfalfa data from Dr. Longnecker’s notes

• Six fields; three plots in each field; four sub-plots in each plot.
• Each plot randomly assigned a type of alfalfa.
• Each sub-plot randomly assigned a cutting date.
• Response for each subplot is yield in tons/acre in the following year.

Alfalfa

alfalfa <- data.frame(yield = c(2.17,1.88,1.62,2.34,1.58,1.66,
                                1.56,1.26,1.22,1.59,1.25,0.94,
                                2.29,1.60,1.67,1.91,1.39,1.12,
                                2.23,2.01,1.82,2.10,1.66,1.10,
                                2.33,2.01,1.70,1.78,1.42,1.35,
                                1.38,1.30,1.85,1.09,1.13,1.06,
                                1.86,1.70,1.81,1.54,1.67,0.88,
                                2.27,1.81,2.01,1.40,1.31,1.06,
                                1.75,1.95,2.13,1.78,1.31,1.30,
                                1.52,1.47,1.80,1.37,1.01,1.31,
                                1.55,1.61,1.82,1.56,1.23,1.13,
                                1.56,1.72,1.99,1.55,1.51,1.33),
                      variety = as.factor(c(rep("ladak",24),
                                            rep("cossack",24),
                                            rep("ranger",24))),
                      date = as.factor(rep(c(rep("none",6),
                                             rep("sep01",6),
                                             rep("sep20",6),
                                             rep("oct07",6)),
                                           3)),
                      field = as.factor(rep(1:6,12)))

Alfalfa

head(alfalfa,n = 26)

   yield variety  date field
1   2.17   ladak  none     1
2   1.88   ladak  none     2
3   1.62   ladak  none     3
4   2.34   ladak  none     4
5   1.58   ladak  none     5
6   1.66   ladak  none     6
7   1.56   ladak sep01     1
8   1.26   ladak sep01     2
9   1.22   ladak sep01     3
10  1.59   ladak sep01     4
11  1.25   ladak sep01     5
12  0.94   ladak sep01     6
13  2.29   ladak sep20     1
14  1.60   ladak sep20     2
15  1.67   ladak sep20     3
16  1.91   ladak sep20     4
17  1.39   ladak sep20     5
18  1.12   ladak sep20     6
19  2.23   ladak oct07     1
20  2.01   ladak oct07     2
21  1.82   ladak oct07     3
22  2.10   ladak oct07     4
23  1.66   ladak oct07     5
24  1.10   ladak oct07     6
25  2.33 cossack  none     1
26  2.01 cossack  none     2

Randomized complete block split plot design

• Each EU randomly assigned to one level of whole-plot factor.
• Each EU receives all levels of split-plot factor in random order.
• Groups of EUs over which this is replicated are treated as blocks.

Treatment effects model for RCB split-plot design

Assume \[ Y_{ijk} = \mu + \tau_i + \gamma_j + (\tau\gamma)_{ij} + C_k + (\tau C)_{ik} + \varepsilon_{ijk}, \] for \(i = 1,\dots,a\), \(j = 1,\dots,b\) and \(k= 1,\dots,c\), where

• \(Y_{ijk}\) is the response of sub-plot \(j\) from whole plot \(i\) in block \(k\).
• the \(\tau_i\) are treatment effects for the whole-plot factor.
• the \(\gamma_j\) are treatment effects for the split-plot factor.
• the \((\tau \gamma)_{ij}\) are interaction effects between the factors.
• the \(C_k\) are independent \(\text{Normal}(0,\sigma_C^2)\) block effects.
• the \((\tau C)_{ik}\) are indep \(\text{Normal}(0,\sigma_{AC}^2)\) whole-plot\(\times\)block effects.
• the \(\varepsilon_{ijk}\) are independent \(\text{Normal}(0,\sigma_{\varepsilon}^2)\) error terms.

Define the cell means as \[ \mu_{ij} = \mu + \tau_i + \gamma_j + (\tau\gamma)_{ij}, \quad i = 1,\dots,a, \quad j = 1,\dots,b. \]

Goals in the RCB split-plot design

In the RCB split-plot model we will focus on how to:

1. Decompose the variability in the \(Y_{ijk}\) into its sources.
2. Test for significance of main effects and interaction.
3. Estimate the variance components \(\sigma_C^2\), \(\sigma_{AC}^2\), and \(\sigma_\varepsilon^2\).
4. Make comparisons between treatment means.
5. Check whether the model assumptions are satisfied.

Sums of squares for the RCB split-plot design

• We have \(\operatorname{SS}_{\operatorname{Tot}}= \operatorname{SS}_{\operatorname{A}}+ \operatorname{SS}_{\operatorname{B}}+ \operatorname{SS}_{\operatorname{AB}}+ \operatorname{SS}_{\operatorname{C}}+ \operatorname{SS}_{\operatorname{AC}}\).

Expected mean squares in RCB split plot design

In the above

• \(\theta_A^2 = (a-1)^{-1}\sum_{i=1}^a(\bar \mu_{i.} - \bar \mu_{..})^2\)
• \(\theta_B^2 = (b-1)^{-1}\sum_{j=1}^b(\bar \mu_{.j} - \bar \mu_{..})^2\)
• \(\theta_{AB}^2 = [(a-1)(b-1)]^{-1}\sum_{i=1}^a\sum_{j=1}^b(\mu_{ij} - (\bar \mu_{i.} + \bar \mu_{.j} - \bar \mu_{..}))^2\)

ANOVA table for RCB split-plot design

1. Reject \(H_0\): \(\bar \mu_{1.} = \dots = \bar\mu_{a.}\) if \(F_{\operatorname{A}}> F_{a-1,\textcolor{red}{(a-1)(c-1)},\alpha}\).
2. Reject \(H_0\): \(\sigma_C^2 = 0\) if \(F_{\operatorname{C}}> F_{c-1,\textcolor{red}{(a-1)(c-1)},\alpha}\).
3. Reject \(H_0\): \(\sigma_{AC}^2 = 0\) if \(F_{\operatorname{AC}}> F_{(a-1)(c-1),a(b-1)(c-1),\alpha}\).
4. Reject \(H_0\): \(\bar \mu_{.1} = \dots = \bar \mu_{.b}\) if \(F_{\operatorname{B}}> F_{b-1,a(b-1)(c-1),\alpha}\).
5. Reject \(H_0\): \(\mu_{ij} = \bar \mu_{i.} + \bar \mu_{.j} - \bar \mu_{..} ~ \forall ij\) if \(F_{\operatorname{AB}}> F_{(a-1)(b-1),a(b-1)(c-1),\alpha}\).

Alfalfa

anova() on lm() output gives wrong p-vals for variety and field.

lm_out <- lm(yield ~ variety + date + variety:date + field + field:variety, 
             data = alfalfa)
anova(lm_out)

Analysis of Variance Table

Response: yield
              Df Sum Sq Mean Sq F value  Pr(>F)    
variety        2   0.18   0.088    3.12  0.0538 .  
date           3   1.97   0.658   23.41 2.8e-09 ***
field          5   4.14   0.828   29.47 3.8e-13 ***
variety:date   6   0.21   0.036    1.27  0.2883    
variety:field 10   1.36   0.136    4.83  0.0001 ***
Residuals     45   1.26   0.028                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Alfalfa

y <- alfalfa$yield
y... <- predict(lm(yield ~ 1,data = alfalfa))
yi.. <- predict(lm(yield ~ variety,data = alfalfa))
y.j. <- predict(lm(yield ~ date,data = alfalfa))
y..k <- predict(lm(yield ~ field,data = alfalfa))
yij. <- predict(lm(yield ~ variety + date + variety:date,data = alfalfa))
yi.k <- predict(lm(yield ~ variety + field + variety:field,data = alfalfa))

SST <- sum((y - y...)^2)
SSA <- sum((yi.. - y...)^2)
SSC <- sum((y..k - y...)^2)
SSAC <- sum((yi.k - (yi.. + y..k - y...))^2)
SSB <- sum((y.j. - y...)^2)
SSAB <- sum((yij. - (yi.. + y.j. - y...))^2)
SSE <- SST - (SSA + SSC + SSAC + SSB + SSAB)

Alfalfa

a <- 3
b <- 4
c <- 6
  
MSA <- SSA / (a-1)
MSC <- SSC / (c-1)
MSAC <- SSAC / ((c-1)*(a-1))
MSB <- SSB / (b-1)
MSAB <- SSAB / ((a-1)*(b-1))
MSE <- SSE / (a*(b-1)*(c-1))

FA_incorrect <- MSA / MSE
FC_incorrect <- MSC / MSE
FA <- MSA / MSAC
FC <- MSC / MSAC
FAC <- MSAC / MSE
FB <- MSB / MSE
FAB <- MSAB / MSE

pA_incorrect <- 1 - pf(FA_incorrect,a-1,a*(b-1)*(c-1))
pC_incorrect <- 1 - pf(FC_incorrect,c-1,a*(b-1)*(c-1))
pA <- 1 - pf(FA,a-1,(c-1)*(a-1))
pC <- 1 - pf(FC,c-1,(c-1)*(a-1))
pAC  <- 1 - pf(FAC,(c-1)*(a-1),a*(b-1)*(c-1))
pB <- 1 - pf(FB,b-1,a*(b-1)*(c-1))
pAB <- 1 - pf(FAB,(a-1)*(b-1),a*(b-1)*(c-1))

Alfalfa

Correct ANOVA table:

Alfalfa

Correct p values for fixed effects with lmerTest package:

library(lmerTest) # first time run install.packages("lmerTest")
lmer_out <- lmer(yield ~ variety + date + variety:date + (1|field) + (1|field:variety), 
                 data = alfalfa)
anova(lmer_out)

Type III Analysis of Variance Table with Satterthwaite's method
             Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)    
variety       0.036   0.018     2    10    0.65    0.54    
date          1.973   0.658     3    45   23.41 2.8e-09 ***
variety:date  0.215   0.036     6    45    1.27    0.29    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Rescales \(\operatorname{SS}_{\operatorname{A}}\) and \(\operatorname{MS}_{\operatorname{A}}\) by the factor \(\operatorname{MS}_{\operatorname{Error}}/ \operatorname{MS}_{\operatorname{AC}}\).

Note: If \(\hat \sigma_C^2 = 0\), it will compute the whole-plot effect p value differently!

Alfalfa

interaction.plot(alfalfa$date,alfalfa$variety,alfalfa$yield)

MoM for variance components in RCB split-plot design

The method of moments estimators for \(\sigma_C^2\), \(\sigma_{AC}^2\), and \(\sigma_{\varepsilon}^2\) are

• \(\displaystyle \dot \sigma_{\varepsilon}^2 = \operatorname{MS}_{\operatorname{Error}}\)
• \(\displaystyle \dot \sigma_{AC}^2 = \frac{\operatorname{MS}_{\operatorname{AC}}- \operatorname{MS}_{\operatorname{Error}}}{b}\)
• \(\displaystyle \dot \sigma_{C}^2 = \frac{\operatorname{MS}_{\operatorname{C}}- \operatorname{MS}_{\operatorname{AC}}}{ab}\)

It is possible to obtain negative values for \(\dot \sigma_{AC}^2\) and \(\dot \sigma_{C}^2\). Use REML.

Alfalfa

Obtain REML estimates of \(\sigma_C^2\), \(\sigma_{AC}^2\), and \(\sigma_\varepsilon^2\) from lmer().

summary(lmer_out)$varcor

 Groups        Name        Std.Dev.
 field:variety (Intercept) 0.164   
 field         (Intercept) 0.240   
 Residual                  0.168   

Variances of some means and difference in means

Some (unadjusted) CIs in RCB split plot design

In the above \(\displaystyle \nu^* = \frac{\operatorname{MS}_{\operatorname{AC}}+ (b-1)\operatorname{MS}_{\operatorname{Error}}}{\frac{\operatorname{MS}_{\operatorname{AC}}^2}{(a-1)(c-1)} + \frac{(b-1)^2\operatorname{MS}_{\operatorname{Error}}^2}{a(b-1)(c-1)}}\) à la Satterthwaite.

Alfalfa

Dunnett’s comparison of marginal date means with no cutting as baseline: \[ \bar Y_{.j.} - \bar Y_{.1.} \pm d_{b-1,a(b-1)(c-1),\alpha}\sqrt{\operatorname{MS}_{\operatorname{Error}}}\sqrt{\frac{2}{ac}} \] Replace, in previous slide, \(t_{a(b-1)(c-1),\alpha/2}\) with \(d_{b-1,a(b-1)(c-1),\alpha}\).

In Dunnett’s table we cannot find \(d_{4,45,0.05}\), so take \(d_{4,40,0.05} = 2.44\).

Dunnett’s table

Alfalfa

y.1.bar <- mean(alfalfa$yield[alfalfa$date=="none"])
y.2.bar <- mean(alfalfa$yield[alfalfa$date=="sep01"])
y.3.bar <- mean(alfalfa$yield[alfalfa$date=="sep20"])
y.4.bar <- mean(alfalfa$yield[alfalfa$date=="oct07"])


se <- sqrt(MSE) * sqrt(2/(a*c))
me <- 2.44 * se

dtab <- rbind(c(y.2.bar - y.1.bar - me,y.2.bar - y.1.bar + me),
              c(y.3.bar - y.1.bar - me,y.3.bar - y.1.bar + me),
              c(y.4.bar - y.1.bar - me,y.4.bar - y.1.bar + me))

colnames(dtab) <- c("lower","upper")
rownames(dtab) <- c("sep01 - none",
                    "sep20 - none",
                    "oct07 - none")
round(dtab,3)

              lower  upper
sep01 - none -0.578 -0.305
sep20 - none -0.343 -0.070
oct07 - none -0.226  0.046

Alfalfa

Unadjusted CIs with ls_means() from R package lmerTest:

ls_means(lmer_out,which="date")

Least Squares Means table:

          Estimate Std. Error  df t value lower upper Pr(>|t|)    
datenone     1.781      0.113 6.1    15.8 1.506 2.056  3.7e-06 ***
dateoct07    1.691      0.113 6.1    15.0 1.416 1.966  5.0e-06 ***
datesep01    1.339      0.113 6.1    11.9 1.065 1.614  2.0e-05 ***
datesep20    1.574      0.113 6.1    14.0 1.300 1.849  7.6e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  Confidence level: 95%
  Degrees of freedom method: Satterthwaite 

Alfalfa

ls_means(lmer_out, pairwise = TRUE, which = "date")

Least Squares Means table:

                      Estimate Std. Error df t value    lower    upper Pr(>|t|)
datenone - dateoct07   0.09000    0.05586 45    1.61 -0.02252  0.20252  0.11416
datenone - datesep01   0.44167    0.05586 45    7.91  0.32915  0.55418  4.7e-10
datenone - datesep20   0.20667    0.05586 45    3.70  0.09415  0.31918  0.00059
dateoct07 - datesep01  0.35167    0.05586 45    6.30  0.23915  0.46418  1.1e-07
dateoct07 - datesep20  0.11667    0.05586 45    2.09  0.00415  0.22918  0.04245
datesep01 - datesep20 -0.23500    0.05586 45   -4.21 -0.34752 -0.12248  0.00012
                         
datenone - dateoct07     
datenone - datesep01  ***
datenone - datesep20  ***
dateoct07 - datesep01 ***
dateoct07 - datesep20 *  
datesep01 - datesep20 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  Confidence level: 95%
  Degrees of freedom method: Satterthwaite