Numerical validation of CANE against independent reference analyses

Software that produces F-tests should be able to show that those F-tests are right. This vignette checks every design supported by CANE against reference analyses specified independently of the package: sums of squares computed by hand from cell means, the stratified aov() fit with explicit Error() terms, REML variance components from nlme, and the base-R distribution functions underlying the three multiple comparison procedures.

library(CANE)
chk <- function(label, a, b, tol = 1e-6) {
  ok <- isTRUE(all.equal(as.vector(a), as.vector(b), tolerance = tol))
  cat(sprintf("%-58s %s\n", label, if (ok) "PASS" else "FAIL"))
  invisible(ok)
}
results <- logical(0)

1. Split plot: every stratum against aov() with Error()

aov() with an Error() specification is the canonical stratified fit in base R and knows nothing of CANE’s internals.

data(spd_example)
o <- PooledSPD(spd_example, "Yield", "Location", "Replication", "MainPlot",
               "SubPlot", 0.05, 1)
tab <- o$Pooled_ANOVA
g <- function(s, c) tab[[c]][tab$Source == s]

ref <- suppressWarnings(summary(aov(
  Yield ~ Location * MainPlot * SubPlot + Error(Location:Replication / MainPlot),
  data = spd_example)))
s1 <- ref[["Error: Location:Replication"]][[1]]
s2 <- ref[["Error: Location:Replication:MainPlot"]][[1]]
s3 <- ref[["Error: Within"]][[1]]
for (z in c("s1", "s2", "s3")) {
  tmp <- get(z); rownames(tmp) <- trimws(rownames(tmp)); assign(z, tmp)
}

results <- c(results,
  chk("Environment SS",            g("Environment", "Sum Sq"),     s1["Location", "Sum Sq"]),
  chk("Rep within Env SS",         g("Rep within Env", "Sum Sq"),  s1["Residuals", "Sum Sq"]),
  chk("Environment F",             g("Environment", "F value"),    s1["Location", "F value"]),
  chk("Main plot SS",              g("Main plot", "Sum Sq"),       s2["MainPlot", "Sum Sq"]),
  chk("Main plot x Env SS",        g("Main plot x Env", "Sum Sq"), s2["Location:MainPlot", "Sum Sq"]),
  chk("Main plot error SS",        g("Main plot error", "Sum Sq"), s2["Residuals", "Sum Sq"]),
  chk("Main plot F vs MP error",   g("Main plot", "F value"),      s2["MainPlot", "F value"]),
  chk("Main plot P",               g("Main plot", "Pr(>F)"),       s2["MainPlot", "Pr(>F)"]),
  chk("Subplot SS",                g("Subplot", "Sum Sq"),         s3["SubPlot", "Sum Sq"]),
  chk("Subplot F vs SP error",     g("Subplot", "F value"),        s3["SubPlot", "F value"]),
  chk("Subplot x Env SS",          g("Subplot x Env", "Sum Sq"),   s3["Location:SubPlot", "Sum Sq"]),
  chk("Main plot x SP SS",         g("Main plot x SP", "Sum Sq"),  s3["MainPlot:SubPlot", "Sum Sq"]),
  chk("MP x SP x Env SS",          g("MP x SP x Env", "Sum Sq"),   s3["Location:MainPlot:SubPlot", "Sum Sq"]),
  chk("Subplot error SS",          g("Subplot error", "Sum Sq"),   s3["Residuals", "Sum Sq"]))
#> Environment SS                                             PASS
#> Rep within Env SS                                          PASS
#> Environment F                                              PASS
#> Main plot SS                                               PASS
#> Main plot x Env SS                                         PASS
#> Main plot error SS                                         PASS
#> Main plot F vs MP error                                    PASS
#> Main plot P                                                PASS
#> Subplot SS                                                 PASS
#> Subplot F vs SP error                                      PASS
#> Subplot x Env SS                                           PASS
#> Main plot x SP SS                                          PASS
#> MP x SP x Env SS                                           PASS
#> Subplot error SS                                           PASS

The main plot F-ratio is the check that matters. Testing the main plot effect against the subplot error instead, as a single-stratum fit does, gives a different and incorrect answer:

F_right <- g("Main plot", "Mean Sq") / g("Main plot error", "Mean Sq")
F_wrong <- g("Main plot", "Mean Sq") / g("Subplot error", "Mean Sq")
data.frame(
  denominator = c("main plot error (correct)", "subplot error (incorrect)"),
  F = c(F_right, F_wrong),
  df2 = c(g("Main plot error", "Df"), g("Subplot error", "Df")),
  P = c(pf(F_right, g("Main plot", "Df"), g("Main plot error", "Df"), lower.tail = FALSE),
        pf(F_wrong, g("Main plot", "Df"), g("Subplot error", "Df"), lower.tail = FALSE)))
#>                 denominator        F df2            P
#> 1 main plot error (correct)  39.4833   6 7.559601e-04
#> 2 subplot error (incorrect) 341.3372  24 1.074257e-15

2. Strip plot: three strata against aov()

data(stpd_example)
o2 <- PooledStPD(stpd_example, "Yield", "Location", "Replication", "RowFactor",
                 "ColumnFactor", 0.05, 1)
t2 <- o2$Pooled_ANOVA
gg <- function(s, c) t2[[c]][t2$Source == s]
r2 <- suppressWarnings(summary(aov(
  Yield ~ Location * RowFactor * ColumnFactor +
    Error(Location:Replication / (RowFactor + ColumnFactor)), data = stpd_example)))
nm <- names(r2)
B <- r2[[grep("Replication:RowFactor", nm)]][[1]]
C <- r2[[grep("Replication:ColumnFactor", nm)]][[1]]
D <- r2[["Error: Within"]][[1]]
for (z in c("B", "C", "D")) { tmp <- get(z); rownames(tmp) <- trimws(rownames(tmp)); assign(z, tmp) }

results <- c(results,
  chk("Row factor SS",         gg("Row factor", "Sum Sq"),         B["RowFactor", "Sum Sq"]),
  chk("Row factor F",          gg("Row factor", "F value"),        B["RowFactor", "F value"]),
  chk("Row strip error SS",    gg("Row strip error", "Sum Sq"),    B["Residuals", "Sum Sq"]),
  chk("Column factor SS",      gg("Column factor", "Sum Sq"),      C["ColumnFactor", "Sum Sq"]),
  chk("Column factor F",       gg("Column factor", "F value"),     C["ColumnFactor", "F value"]),
  chk("Column strip error SS", gg("Column strip error", "Sum Sq"), C["Residuals", "Sum Sq"]),
  chk("Row x Column SS",       gg("Row x Column", "Sum Sq"),       D["RowFactor:ColumnFactor", "Sum Sq"]),
  chk("Row x Col x Env SS",    gg("Row x Col x Env", "Sum Sq"),    D["Location:RowFactor:ColumnFactor", "Sum Sq"]),
  chk("Intersection error SS", gg("Intersection error", "Sum Sq"), D["Residuals", "Sum Sq"]))
#> Row factor SS                                              PASS
#> Row factor F                                               PASS
#> Row strip error SS                                         PASS
#> Column factor SS                                           PASS
#> Column factor F                                            PASS
#> Column strip error SS                                      PASS
#> Row x Column SS                                            PASS
#> Row x Col x Env SS                                         PASS
#> Intersection error SS                                      PASS

3. Sums of squares computed by hand from cell means

The strongest check involves no modelling function at all: each sum of squares is accumulated directly from the deviations of the relevant cell means.

d <- spd_example
y <- d$Yield; E <- d$Location; P <- d$Replication; A <- d$MainPlot; B2 <- d$SubPlot
l <- nlevels(E); r <- nlevels(P); a <- nlevels(A); b <- nlevels(B2); G <- mean(y)
m <- function(...) tapply(y, list(...), mean)
mE <- m(E); mA <- m(A); mB <- m(B2); mEP <- m(E, P); mEA <- m(E, A)
mEB <- m(E, B2); mAB <- m(A, B2); mEPA <- m(E, P, A); mEAB <- m(E, A, B2)
acc <- function(n, f, ...) { gr <- expand.grid(...); n * sum(apply(gr, 1, function(z) f(z)^2)) }

hand <- c(
  "Environment"     = acc(r*a*b, function(z) mE[z[1]] - G, j = 1:l),
  "Rep within Env"  = acc(a*b,   function(z) mEP[z[1], z[2]] - mE[z[1]], j = 1:l, k = 1:r),
  "Main plot"       = acc(l*r*b, function(z) mA[z[1]] - G, i = 1:a),
  "Main plot x Env" = acc(r*b,   function(z) mEA[z[1], z[2]] - mE[z[1]] - mA[z[2]] + G,
                          j = 1:l, i = 1:a),
  "Main plot error" = acc(b, function(z) mEPA[z[1], z[2], z[3]] - mEP[z[1], z[2]] -
                            mEA[z[1], z[3]] + mE[z[1]], j = 1:l, k = 1:r, i = 1:a),
  "Subplot"         = acc(l*r*a, function(z) mB[z[1]] - G, q = 1:b),
  "Subplot x Env"   = acc(r*a,   function(z) mEB[z[1], z[2]] - mE[z[1]] - mB[z[2]] + G,
                          j = 1:l, q = 1:b),
  "Main plot x SP"  = acc(l*r,   function(z) mAB[z[1], z[2]] - mA[z[1]] - mB[z[2]] + G,
                          i = 1:a, q = 1:b),
  "MP x SP x Env"   = acc(r, function(z) mEAB[z[1], z[2], z[3]] - mEA[z[1], z[2]] -
                            mEB[z[1], z[3]] - mAB[z[2], z[3]] + mE[z[1]] + mA[z[2]] +
                            mB[z[3]] - G, j = 1:l, i = 1:a, q = 1:b))
for (s in names(hand)) results <- c(results, chk(paste0(s, " SS (by hand)"), g(s, "Sum Sq"), hand[[s]]))
#> Environment SS (by hand)                                   PASS
#> Rep within Env SS (by hand)                                PASS
#> Main plot SS (by hand)                                     PASS
#> Main plot x Env SS (by hand)                               PASS
#> Main plot error SS (by hand)                               PASS
#> Subplot SS (by hand)                                       PASS
#> Subplot x Env SS (by hand)                                 PASS
#> Main plot x SP SS (by hand)                                PASS
#> MP x SP x Env SS (by hand)                                 PASS
results <- c(results, chk("Subplot error SS (by subtraction)",
  g("Subplot error", "Sum Sq"), sum((y - G)^2) - sum(hand)))
#> Subplot error SS (by subtraction)                          PASS

4. Degrees of freedom against design theory

With environments, replications, main plot levels and subplot levels, theory gives for replication-within-environment, for the main plot error and for the subplot error; a strip plot has , and for its three errors; and a Latin square series has each for rows and columns and for the error.

results <- c(results,
  chk("SPD: rep within env df",      g("Rep within Env", "Df"),  l*(r-1)),
  chk("SPD: main plot error df",     g("Main plot error", "Df"), l*(r-1)*(a-1)),
  chk("SPD: subplot error df",       g("Subplot error", "Df"),   l*a*(r-1)*(b-1)),
  chk("SPD: df partition complete",  sum(tab$Df[tab$Source != "Total"]), nrow(spd_example)-1),
  chk("StPD: row strip error df",    gg("Row strip error", "Df"),    2*(4-1)*(2-1)),
  chk("StPD: column strip error df", gg("Column strip error", "Df"), 2*(4-1)*(3-1)),
  chk("StPD: intersection error df", gg("Intersection error", "Df"), 2*(4-1)*(2-1)*(3-1)),
  chk("StPD: df partition complete", sum(t2$Df[t2$Source != "Total"]), nrow(stpd_example)-1))
#> SPD: rep within env df                                     PASS
#> SPD: main plot error df                                    PASS
#> SPD: subplot error df                                      PASS
#> SPD: df partition complete                                 PASS
#> StPD: row strip error df                                   PASS
#> StPD: column strip error df                                PASS
#> StPD: intersection error df                                PASS
#> StPD: df partition complete                                PASS

data(lsd_example)
tl <- PooledLSD(lsd_example, "Yield", "Location", "Treatment", "Row", "Column",
                0.05, 1)$Pooled_ANOVA
results <- c(results,
  chk("LSD: row within env df",      tl$Df[tl$Source == "Row within Env"], 3*(4-1)),
  chk("LSD: column within env df",   tl$Df[tl$Source == "Col within Env"], 3*(4-1)),
  chk("LSD: error df",               tl$Df[tl$Source == "Pooled error"],   3*(4-1)*(4-2)),
  chk("LSD: df partition complete",  sum(tl$Df[tl$Source != "Total"]), nrow(lsd_example)-1))
#> LSD: row within env df                                     PASS
#> LSD: column within env df                                  PASS
#> LSD: error df                                              PASS
#> LSD: df partition complete                                 PASS

data(sugarcane)
tr <- PooledRBD(sugarcane, "ccsy", "env", "gen", "rep", 0.05, 1)$Pooled_ANOVA
results <- c(results,
  chk("RBD: rep within env df",      tr$Df[tr$Source == "Rep within Env"],  6*(2-1)),
  chk("RBD: treatment x env df",     tr$Df[tr$Source == "Treatment x Env"], 5*17),
  chk("RBD: error df",               tr$Df[tr$Source == "Pooled error"],    6*(2-1)*17),
  chk("RBD: df partition complete",  sum(tr$Df[tr$Source != "Total"]), nrow(sugarcane)-1))
#> RBD: rep within env df                                     PASS
#> RBD: treatment x env df                                    PASS
#> RBD: error df                                              PASS
#> RBD: df partition complete                                 PASS

5. REML variance components reproduce the split plot strata

Expected mean squares link the ANOVA strata to the variance components: the subplot error mean square estimates , and the main plot error mean square estimates . Fitting the same model by REML in nlme, which uses an entirely different estimation algorithm, should recover both.

s <- spd_example
s$EP  <- interaction(s$Location, s$Replication, drop = TRUE)
s$EPA <- interaction(s$Location, s$Replication, s$MainPlot, drop = TRUE)
fm <- nlme::lme(Yield ~ Location * MainPlot * SubPlot, random = ~ 1 | EP / EPA,
                data = s, method = "REML")
vc <- suppressWarnings(as.numeric(nlme::VarCorr(fm)[, "Variance"]))
vc <- vc[!is.na(vc)]
sig_rep <- vc[1]; sig_mp <- vc[2]; sig_e <- vc[length(vc)]
cat(sprintf("REML: sigma^2_rep = %.4f   sigma^2_mainplot = %.4f   sigma^2_e = %.4f\n",
            sig_rep, sig_mp, sig_e))
#> REML: sigma^2_rep = 1969.2820   sigma^2_mainplot = 3750.9440   sigma^2_e = 1471.9070
results <- c(results,
  chk("subplot error MS == REML residual variance", g("Subplot error", "Mean Sq"),
      sig_e, tol = 1e-4),
  chk("main plot error MS == sigma^2_e + b sigma^2_mp", g("Main plot error", "Mean Sq"),
      sig_e + b * sig_mp, tol = 1e-3))
#> subplot error MS == REML residual variance                 PASS
#> main plot error MS == sigma^2_e + b sigma^2_mp             PASS

6. Bartlett’s statistic computed by hand

CANE computes Bartlett’s statistic from the environment-wise error mean squares and their error degrees of freedom. Applying stats::bartlett.test() to the stacked residuals instead, as version 0.1.1 did, estimates each variance with the wrong divisor and the wrong degrees of freedom, and gives a materially different statistic.

ob <- PooledRBD(sugarcane, "ccsy", "env", "gen", "rep", 0.05, 1)
ems <- ob$Transformation$error_mean_squares
dfe <- ob$Transformation$error_df
k <- length(ems); N <- sum(dfe); sp2 <- sum(dfe * ems) / N
chi <- N * log(sp2) - sum(dfe * log(ems))
Cc  <- 1 + (sum(1/dfe) - 1/N) / (3 * (k - 1))
results <- c(results,
  chk("Bartlett K-squared (by hand)", ob$Bartlett_Test$statistic, chi / Cc),
  chk("Bartlett p-value (by hand)",   ob$Bartlett_Test$p.value,
      pchisq(chi / Cc, k - 1, lower.tail = FALSE)),
  chk("pooled error variance",        ob$Bartlett_Test$pooled_error_variance, sp2))
#> Bartlett K-squared (by hand)                               PASS
#> Bartlett p-value (by hand)                                 PASS
#> pooled error variance                                      PASS

res_list <- lapply(levels(sugarcane$env), function(e) {
  s <- droplevels(sugarcane[sugarcane$env == e, ])
  residuals(lm(ccsy ~ rep + gen, data = s))
})
rv <- unlist(res_list)
gl <- rep(levels(sugarcane$env), times = vapply(res_list, length, 1L))
data.frame(
  method = c("CANE 0.2.0: on error mean squares, error df",
             "CANE 0.1.1: bartlett.test() on stacked residuals"),
  statistic = c(unname(ob$Bartlett_Test$statistic),
                unname(bartlett.test(rv, gl)$statistic)),
  df_per_env = c(unique(dfe), length(res_list[[1]]) - 1))
#>                                             method statistic df_per_env
#> 1      CANE 0.2.0: on error mean squares, error df  81.17245         17
#> 2 CANE 0.1.1: bartlett.test() on stacked residuals 169.06426         35

7. The Aitken identity

Dividing each observation by the square root of the error mean square of its own environment makes each environment contribute an error sum of squares equal to its own error degrees of freedom. The pooled error mean square of a single-error design must therefore be exactly 1 after the transformation. This is a sharp test of the per-row lookup: any misalignment of observations to environments breaks it.

pe <- ob$Pooled_ANOVA
results <- c(results,
  chk("pooled error MS == 1",  pe$`Mean Sq`[pe$Source == "Pooled error"], 1),
  chk("pooled error SS == df", pe$`Sum Sq`[pe$Source == "Pooled error"],
                               pe$Df[pe$Source == "Pooled error"]))
#> pooled error MS == 1                                       PASS
#> pooled error SS == df                                      PASS

8. Multiple comparison procedures against base R

Tukey’s HSD is checked against stats::TukeyHSD(); Fisher’s LSD against an unadjusted pooled-variance t test; and Duncan’s procedure against the identity that its two-mean critical range coincides with Fisher’s LSD, since the protection level for a span of two means is and .

set.seed(7)
w <- data.frame(E = factor(rep(c("A", "B"), each = 40)),
                T = factor(rep(paste0("T", 1:4), 20)), y = rnorm(80, 10, 1))
oc <- PooledCRD(w, "y", "E", "T", 0.05, 1, transform = "never")
tc <- oc$Treatments_Comparison
MSE <- tc$statistics$MSerror; DFE <- tc$statistics$Df
th <- TukeyHSD(aov(y ~ E + T + E:T, data = w), "T")$T
results <- c(results,
  chk("Tukey critical difference", tc$statistics$Critical_difference,
      qtukey(0.95, 4, DFE) * sqrt(MSE / 20)),
  chk("Tukey p-values == stats::TukeyHSD",
      sort(as.vector(tc$comparison$p.value)), sort(as.vector(th[, "p adj"]))))
#> Tukey critical difference                                  PASS
#> Tukey p-values == stats::TukeyHSD                          PASS

ol <- suppressWarnings(PooledCRD(w, "y", "E", "T", 0.05, 3, transform = "never"))
mu <- tapply(w$y, w$T, mean); pr <- combn(names(mu), 2)
tstat <- abs(mu[pr[1, ]] - mu[pr[2, ]]) / sqrt(2 * MSE / 20)
ord <- match(paste(pr[1, ], "-", pr[2, ]), ol$Treatments_Comparison$comparison$Comparison)
results <- c(results,
  chk("LSD critical difference", ol$Treatments_Comparison$statistics$Critical_difference,
      qt(0.975, DFE) * sqrt(2 * MSE / 20)),
  chk("LSD p-values == unadjusted t test",
      ol$Treatments_Comparison$comparison$p.value[ord],
      2 * pt(tstat, DFE, lower.tail = FALSE)),
  chk("Duncan R_2 == Fisher LSD", qtukey(0.95, 2, DFE) * sqrt(MSE / 20),
      qt(0.975, DFE) * sqrt(2 * MSE / 20)))
#> LSD critical difference                                    PASS
#> LSD p-values == unadjusted t test                          PASS
#> Duncan R_2 == Fisher LSD                                   PASS

The letter-based grouping is checked against multcompView::multcompLetters(), by comparing the relation “shares a letter with” that each labelling induces.

tcs <- ob$Treatments_Comparison
muS <- setNames(tcs$means$Mean, tcs$means$Treatment)
prS <- combn(names(muS), 2)
pvS <- ptukey(abs(muS[prS[1, ]] - muS[prS[2, ]]) /
              sqrt(tcs$statistics$MSerror / tcs$parameters$r),
              tcs$parameters$ntr, tcs$statistics$Df, lower.tail = FALSE)
names(pvS) <- paste(prS[1, ], prS[2, ], sep = "-")
refL <- multcompView::multcompLetters(pvS, threshold = 0.05)$Letters
mine <- setNames(tcs$groups$Groups, tcs$groups$Treatment)
rel <- function(lab) {
  n <- names(lab)
  outer(n, n, Vectorize(function(i, j)
    any(strsplit(lab[[i]], "")[[1]] %in% strsplit(lab[[j]], "")[[1]])))
}
results <- c(results,
  chk("grouping relation == multcompView", sum(rel(as.list(mine)) !=
      rel(as.list(refL[names(mine)]))), 0))
#> grouping relation == multcompView                          PASS

Summary

cat(sprintf("%d of %d checks passed.\n", sum(results), length(results)))
#> 62 of 62 checks passed.
if (!all(results)) {
  cat("FAILED:\n"); print(which(!results))
}
stopifnot(all(results))

Every pooled sum of squares, every F-ratio denominator, every degree of freedom, Bartlett’s statistic, the Aitken identity and all three multiple comparison procedures agree with reference analyses specified independently of the package.