{cat("--------------- Analysis of Variance ---------------\n");T}
{
#Functions: aov, summary.aov
#Data: Sokal and Rohlf, Box 9.4 p. 220
#Reference: Sokal, R. and F. J. Rohlf. 1981. Biometry, 2nd edition.
#           W. H. Freeman and Company
#Description: 1-way balanced layout with 5 treatments, 10 replicates/trtm;
#             check df's, sum of squares and mean squares; check F value
#             and p-value using a different tolerance
tol1 <- 5e-3
tol2 <- 4e-2
y <- c(75,67,70,75,65,71,67,67,76,68,57,58,60,59,62,60,60,57,59,61,
       58,61,56,58,57,56,61,60,57,58,58,59,58,61,57,56,58,57,57,59,
       62,66,65,63,64,62,65,65,62,67)
y.treat <- factor(rep(1:5,c(10,10,10,10,10)))
y.df <- data.frame(y,y.treat)
y.aov <- aov(y ~ y.treat, data=y.df)
a.tab <- summary(y.aov)
all(c(a.tab$Df == c(4,45),
      abs(a.tab$"Sum of Sq" - c(1077.32,245.50)) < tol1,
      abs(a.tab$"Mean Sq" - c(269.33,5.46)) < tol1,
      abs(a.tab$"F Value"[1] - 49.33) < tol2,
      a.tab$"Pr(F)"[1] < 0.001))
}
{
#Function: aov
#Data: Sokal and Rohlf, Box 11.2 p. 345
#Reference: Sokal, R. and F. J. Rohlf. 1981. Biometry, 2nd edition.
#           W. H. Freeman and Company
#Description: 2-way layout without replication; check sum of squares and df's
#             and mean squares
tol1 <- 3e-3
tol2 <- 2e-2
tol3 <- 5e-2
y <- c(23.8,22.6,22.2,21.2,18.4,13.5,9.8,6.0,5.8,5.6,24.0,22.4,22.1,21.8,
       19.3,14.4,9.9,6.0,5.9,5.6,24.6,22.9,22.1,21.0,19.0,14.2,10.4,6.3,
       6.0,5.5,24.8,23.2,22.2,21.2,18.8,13.8,9.6,6.3,5.8,5.6)
y.a <- factor(rep(1:4,c(10,10,10,10)))
y.b <- factor(rep(rep(1:10,1),4))
y.anova <- aov(y ~ y.a + y.b)
a.tab <- summary(y.anova)
all(c(a.tab$Df == c(3,9,27),
      abs(a.tab$"Sum of Sq" - c(0.57,2119.66,2.23)) < tol2,
      abs(a.tab$"Mean Sq" - c(0.1900,235.5178,0.0826)) < tol1,
      abs(a.tab$"F Value"[1] - 2.30) < tol3,
      abs((a.tab$"F Value"[2]-2851.6)/a.tab$"F Value"[2]) < tol2))
}
{
#Function: aov
#Data: Sokal and Rohlf, Table 11.1 p. 325
#Reference: Sokal, R. and F. J. Rohlf. 1981. Biometry, 2nd edition.
#           W. H. Freeman and Company
#Description: 2-way layout with replicates; check sum of squares and df's
#             and mean squares
tol <- 4e-3
y.f1 <- factor(rep(1:2,c(6,6)))
y.f2 <- factor(rep(rep(1:2,c(3,3)),2))
y <- c(709,679,699,657,594,677,592,538,476,508,505,539)
y.anova <- aov(y ~ y.f2 * y.f1)
a.tab <- summary(y.anova)
all(c(a.tab$Df == c(1,1,1,8),
      abs(a.tab$"Sum of Sq" - c(3780.75,61204.08,918.75,11666.67)) < tol,
      abs(a.tab$"Mean Sq" - c(3780.75,61204.08,918.75,1458.33)) < tol))
}
{
#Function: aov
#Data: Wittkowski, Fig. 1 p. 120
#Reference: Wittkowski, K.M. 1992. Statistical Analysis of Unbalanced and 
#           Incomplete Designs---Experiences with BMDP and SAS. Computational
#           Statistics and Data Analysis. 14:119-124.
#Description: complete balanced block design; 3 treatments, 2 blocks; model
#             main effects; check sum of squares, df's, F values and
#             p-values
tol <- 1e-6
y <- c(1,1,1,2,2,2,3,3,3,3,3,3,2,2,2,1,1,1)
y.treat <- factor(c(1,1,1,2,2,2,3,3,3,1,1,1,2,2,2,3,3,3))
y.block <- factor(c(1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2))
y.anova <- aov(y ~ y.block + y.treat)
all(c(abs(summary(y.anova)$"Sum of Sq" - c(0.0000,0.0000,12.0000)) < tol,
      summary(y.anova)$"Df" == c(1,2,14),
      abs(summary(y.anova)$"F Value"[1:2] - c(0.00,0.00)) < tol,
      abs(summary(y.anova)$"Pr(F)"[1:2] - c(1.0000,1.0000)) < tol))
}
{
#Function: aov
#Data: Wittkowski, Tab. 1 p. 119
#Reference: Wittkowski, K.M. 1992. Statistical Analysis of Unbalanced and 
#           Incomplete Designs---Experiences with BMDP and SAS. Computational
#           Statistics and Data Analysis. 14:119-124.
#Description: complete unbalanced block design; 3 treatments, 2 blocks; model
#             main effects; check sum of squares, df's and
#             p-values
tol1 <- 5e-5
tol2 <- 2e-3
y <- c(1,1,2,2,2,3,3,2,1)
y.treat <- factor(c("A","A","B","B","B","C","A","B","C"))
y.block <- factor(c("1","1","1","1","1","1","2","2","2"))
y.df <- data.frame(y,y.treat,y.block)
y.anova <- aov(y ~ y.block + y.treat, data=y.df)
all(c(abs(summary(y.anova)$"Sum of Sq" - c(0.0556,0.2246,4.6087)) < tol1,
      summary(y.anova)$"Df" == c(1,2,5),
      abs(summary(y.anova)$"Pr(F)"[1:2] - c(0.8158,0.8878)) < tol1,
      abs(summary(y.anova)$"F Value"[1:2] - c(0.06,0.12)) < tol2))
}
{
#Function: aov
#Data: Wittkowski, Tab. 1 p. 119
#Reference: Wittkowski, K.M. 1992. Statistical Analysis of Unbalanced and
#           Incomplete Designs---Experiences with BMDP and SAS. Computational
#           Statistics and Data Analysis. 14:119-124.
#Description: complete unbalanced block design continued from previous
#             2 tests; model main effects with interaction; check sum of 
#             squares and df's
tol1 <- 5e-5
tol2 <- 1.01
y.anova <- aov(y ~ y.block * y.treat, data=y.df)
all(c(abs(summary(y.anova)$"Sum of Sq" - c(0.0556,0.2246,4.6087,0.0000))
           < tol1,
      summary(y.anova)$"Df" == c(1,2,2,3),
      abs((summary(y.anova)$"F Value"[1:3] - c(99999,99999,99999))/
           summary(y.anova)$"F Value"[1:3]) < tol2))
}
{
#Function: aov
#Data: Sokal and Rohlf, Box 12.1 p. 375
#Reference: Sokal, R. and F. J. Rohlf. 1981. Biometry, 2nd edition.
#           W. H. Freeman and Company
#Description: 2x3x5 unreplicated factorial design; check sum of squares,
#             df's, mean squares, F values, and p-values
tol1 <- 7e-5
tol2 <- 4e-4
tol3 <- 1e-6
y.a <- factor(rep(rep(1:3,c(3,3,3)),5))
y.b <- factor(rep(1:5,c(9,9,9,9,9)))
y.c <- factor(rep(c(1,2,3),15))
y <- c(201,246,271,124,158,207, 79,129,142,150,164,170,104,111,117, 63, 54, 93,
       131,138,149, 86, 99, 81, 50, 51, 62,130,136,127, 89, 91, 87, 51, 52, 51,
        97,102, 99, 60, 74, 72, 32, 46, 52)
y.df <- data.frame(y.a,y.b,y.c,y)
y.anova <- aov(y ~ (y.a+y.b+y.c)^2, data=y.df)
a.tab <- summary(y.anova)
all(c(abs(a.tab$"Sum of Sq" - c(57116.1333,55545.4667,3758.8000,3685.8667,
                               97.0667,5264.5333,1034.9333)) < tol1,
      a.tab$"Df" == c(2,4,2,8,4,8,16),
      abs(a.tab$"Mean Sq" - c(28558.0666,13886.3667,1879.4000,460.7333,
                             24.2667,658.0667,64.6833)) < tol1,
      abs(a.tab$"F Value"[4] - 7.123) < tol2,
      a.tab$"F Value"[5] < 1,
      abs(a.tab$"F Value"[6] - 10.174) < tol2,
      abs(a.tab$"Pr(F)"[4] - 0.000460) < tol3,
      abs(a.tab$"Pr(F)"[6] - 5.5e-5) < tol3))
}
{
#Function: aov
#Data: Snedecor and Cochran, Table 12.9.1 p. 359
#Reference: Snedecor, G. and W. Cochran. 1967. Statistical Methods, 6th ed.
#           The Iowa State University Press.
#Description: 2^k design; specifically a 2^3; check sum of squares
tol <- 4e-5
y.lysine <- factor(rep(rep(1:2,c(4,4)),8))
y.protein <- factor(rep(rep(1:2,c(2,2)),16))
y.sex <- factor(rep(1:2,32))
y <- c(1.11,1.03,1.52,1.48,1.22,0.87,1.38,1.09,0.97,0.97,1.45,1.22,1.13,1.00,
       1.08,1.09,1.09,0.99,1.27,1.53,1.34,1.16,1.40,1.47,0.99,0.99,1.22,1.19,
       1.41,1.29,1.21,1.43,0.85,0.99,1.67,1.16,1.34,1.00,1.46,1.24,1.21,1.21,
       1.24,1.57,1.19,1.14,1.39,1.17,1.29,1.19,1.34,1.13,1.25,1.36,1.17,1.01,
       0.96,1.24,1.32,1.43,1.32,1.32,1.21,1.13)
y.df <- data.frame(y.lysine,y.protein,y.sex,y)
y.anova <- aov(y ~ y.lysine*y.protein*y.sex,data=y.df)
all(c(abs(summary(y.anova)$"Sum of Sq"[1:7] - c(0.0032,0.4307,0.0570,0.2588,
                                                0.0375,0.0001,0.0113)) < tol))
}
{
#Functions: aov, summary.aov
#Data: Snedecor and Cochran, Table 12.12.1 p. 371
#Reference: Snedecor, G. and W. Cochran. 1967. Statistical Methods, 6th ed.
#           The Iowa State University Press.
#Description: split-plot design with 2 blocks, 3 varieties and 4 dates/variety;
#             check sum of squares, df's, and mean squares
tol <- 1.5e-4
y.block <- factor(rep(rep(1:6,c(4,4,4,4,4,4)),3))
y.date <- factor(rep(c("A","B","C","D"),18))
y.var <- factor(rep(1:3,c(24,24,24)))
y <- c(2.17,1.58,2.29,2.23,1.88,1.26,1.60,2.01,1.62,1.22,1.67,1.82,
       2.34,1.59,1.91,2.10,1.58,1.25,1.39,1.66,1.66,0.94,1.12,1.10,
       2.33,1.38,1.86,2.27,2.01,1.30,1.70,1.81,1.70,1.85,1.81,2.01,
       1.78,1.09,1.54,1.40,1.42,1.13,1.67,1.31,1.35,1.06,0.88,1.06,
       1.75,1.52,1.55,1.56,1.95,1.47,1.61,1.72,2.13,1.80,1.82,1.99,
       1.78,1.37,1.56,1.55,1.31,1.01,1.23,1.51,1.30,1.31,1.13,1.33)
y.df <- data.frame(y.block,y.date,y.var,y)
y.anova <- aov(y ~ y.var * y.block + Error(y.date*y.var), data=y.df)
a.tab <- summary(y.anova)
all(c(abs(a.tab$"Error: y.date"$"Sum of Sq" - 1.9625) < tol,
      abs(a.tab$"Error: y.var"$"Sum of Sq" - 0.1781) < tol,
      abs(a.tab$"Error: y.date:y.var"$"Sum of Sq" - 0.2105) < tol,
      abs(a.tab$"Error: Within"$"Sum of Sq" - c(4.1499,1.3622,1.2586)) < tol,
      a.tab$"Error: y.date"$"Df" == 3,
      a.tab$"Error: y.var"$"Df" == 2,
      a.tab$"Error: y.date:y.var"$"Df" == 6,
      a.tab$"Error: Within"$"Df" == c(5,10,45),
      abs(a.tab$"Error: y.date"$"Mean Sq" - 0.6542) < tol,
      abs(a.tab$"Error: y.var"$"Mean Sq" - 0.0890) < tol,
      abs(a.tab$"Error: y.date:y.var"$"Mean Sq" - 0.0351) < tol,
      abs(a.tab$"Error: Within"$"Mean Sq" - c(0.8300,0.1362,0.0280)) < tol))
}
{
#Function: aov
#Data: Zar, Example 12.5 p. 261
#Reference: Zar, J.H. Biostatistical Analysis, 3rd edition. 1996.
#           Prentice-Hall Inc.
#Description: repeated measures on 7 subjects; measurements taken once after
#             one of three drug treatments; check sum of squares, df's, 
#             mean squares, F value, and p-value 
tol1 <- 5e-3
tol2 <- 6e-2
tol3 <- 5e-5
y.drug <- factor(rep(1:3,c(7,7,7)))
y.subjects <- factor(rep(rep(c(1,2,3,4,5,6,7),1),3))
y <- c(164,202,143,210,228,173,161,152,181,136,194,219,159,157,
       178,222,132,216,245,182,165)
y.anova <- aov(y ~ y.drug + Error(y.subjects))
a.tab <- summary(y.anova)
all(c(abs(a.tab$"Error: y.subjects"$"Sum of Sq" - 18731.24) < tol1,
      abs(a.tab$"Error: Within"$"Sum of Sq" - c(1454.00,695.33)) < tol1,
      a.tab$"Error: y.subjects"$"Df" == 6,
      a.tab$"Error: Within"$"Df" == c(2,12),
      abs(a.tab$"Error: Within"$"Mean Sq" - c(727.00,57.94)) < tol1,
      abs(a.tab$"Error: Within"$"F Value"[1] - 12.6) < tol2,
      abs(a.tab$"Error: Within"$"Pr(F)"[1] - 0.0011) < tol3))
}
{
#Functions: manova, summary.manova
#Data: Tabachnick and Fidell, Table 9.1 p. 382
#Reference: Tabachnick, B.G. and L.S. Fidell. 1989. Using Multivariate
#           Statistics, 2nd edition. Harper and Row, Publishers
#Description: simple manova; check test statistics (Pillais's Trace, 
#             Wilk's Lambda, Hotelling-Lawley Trace and Roy's greatest root),
#             F value, numerator and demoninator df's, and p-value against
#             results from SAS using proc ANOVA
tol1 <- 5e-9
tol2 <- 8.7e-5
y.disable <- factor(rep(1:3,c(6,6,6)))
y.treat <- factor(rep(rep(1:2,c(3,3)),3))
y1 <- c(115,98,107,90,85,80,100,105,95,70,85,78,89,100,90,65,80,72)
y2 <- c(108,105,105,92,95,81,105,95,98,80,68,82,78,85,95,62,70,73)
y.df <- data.frame(y.disable,y.treat,y1,y2)
y.manova <- manova(cbind(y1,y2) ~ y.disable * y.treat)
y.sum <- summary(y.manova)
all(c(abs(y.sum$Stats[,1,1] - c(0.77249870,0.87030314,0.06184706)) < tol1,
      abs(y.sum$Stats[,2,1] - c(0.23022530,0.12969686,0.93827013)) < tol1,
      abs(y.sum$Stats[,3,1] - c(3.33173945,6.71028690,0.06566626)) < tol1,
      abs(y.sum$Stats[,4,1] - c(3.32818439,6.71028690,0.06370583)) < tol1,
      abs(y.sum$Stats[,1,2] - c(3.7760,36.9066,0.1915)) < tol2,
      abs(y.sum$Stats[,2,2] - c(5.9627,36.9066,0.1780)) < tol2,
      abs(y.sum$Stats[,3,2] - c(8.3293,36.9066,0.1642)) < tol2,
      abs(y.sum$Stats[,4,2] - c(19.9691,36.9066,0.3822)) < tol2,
      y.sum$Stats[,1,3] == c(4,2,4),
      y.sum$Stats[,2,3] == c(4,2,4),
      y.sum$Stats[,3,3] == c(4,2,4),
      y.sum$Stats[,4,3] == c(2,2,2),
      y.sum$Stats[,1,4] == c(24,11,24),
      y.sum$Stats[,2,4] == c(22,11,22),
      y.sum$Stats[,3,4] == c(20,11,20),
      y.sum$Stats[,4,4] == c(12,11,12),
      abs(y.sum$Stats[,1,5] - c(0.0161,0.0001,0.9405)) < tol2,
      abs(y.sum$Stats[,2,5] - c(0.0021,0.0001,0.9473)) < tol2,
      abs(y.sum$Stats[,3,5] - c(0.0004,0.0001,0.9541)) < tol2,
      abs(y.sum$Stats[,4,5] - c(0.0002,0.0001,0.6904)) < tol2))
}
{
#Functions: manova, summary.manova
#Data: Tabachnick and Fidell, Table 10.13 p. 474
#Reference: Tabachnick, B.G. and L.S. Fidell. 1989. Using Multivariate
#           Statistics, 2nd edition. Harper and Row, Publishers
#Description: repeated measures manova; test data consists of weight loss
#             over 3 months for 3 treatments groups; check test statistics
#             (Pillais's Trace, Wilk's Lambda, Hotelling-Lawley Trace and 
#             Roy's greatest root), F value, numerator and demoninator df's,
#             and p-value against results from SAS using proc GLM with 
#             repeated time option
tol1 <- 5e-9
tol2 <- 1.1e-4
y.group <- factor(rep(1:3,c(12,12,12)))
y1 <- c(4,4,4,3,5,6,6,5,5,3,4,5,6,5,7,6,3,5,4,4,6,7,4,7,8,3,7,4,9,2,3,6,6,
        9,7,8)
y2 <- c(3,4,3,2,3,5,5,4,4,3,2,2,3,4,6,4,2,5,3,2,5,6,3,4,4,6,7,7,7,4,5,5,6,
        5,9,6)
y3 <- c(3,3,1,1,2,4,4,1,1,2,2,1,2,1,3,2,1,4,1,1,3,4,2,3,2,3,4,1,3,1,1,2,3,
        2,4,1)
Y <- cbind(y1,y2,y3)
y.manova <- manova(Y %*% contr.poly(3) ~ y.group)
y.sum <- summary(y.manova,intercept=T)
all(c(abs(y.sum$Stats[,1,1] - c(0.88288872,0.57613065)) < tol1,
      abs(y.sum$Stats[,2,1] - c(0.11711128,0.43329986)) < tol1,
      abs(y.sum$Stats[,3,1] - c(7.53888748,1.28610617)) < tol1,
      abs(y.sum$Stats[,4,1] - c(7.53888748,1.26895473)) < tol1,
      abs(y.sum$Stats[,1,2] - c(120.6222,6.6763)) < tol2,
      abs(y.sum$Stats[,2,2] - c(120.6222,8.3067)) < tol2,
      abs(y.sum$Stats[,3,2] - c(120.6222,9.9673)) < tol2,
      abs(y.sum$Stats[,4,2] - c(120.6222,20.9378)) < tol2,
      y.sum$Stats[,1,3] == c(2,4),
      y.sum$Stats[,2,3] == c(2,4),
      y.sum$Stats[,3,3] == c(2,4),
      y.sum$Stats[,4,3] == c(2,2),
      y.sum$Stats[,1,4] == c(32,66),
      y.sum$Stats[,2,4] == c(32,64),
      y.sum$Stats[,3,4] == c(32,62),
      y.sum$Stats[,4,4] == c(32,33),
      abs(y.sum$Stats[,1,5] - c(0.0001,0.0001)) < tol2,
      abs(y.sum$Stats[,2,5] - c(0.0001,0.0001)) < tol2,
      abs(y.sum$Stats[,3,5] - c(0.0001,0.0001)) < tol2,
      abs(y.sum$Stats[,4,5] - c(0.0001,0.0001)) < tol2))
}
{
#Clean up after anova
cleanup <- T
rm(y, y.treat, y.block, y.aov, a.tab, y.anova, y.df)
rm(y.a, y.b, y.f1, y.f2, y.c, y.lysine, y.protein, y.sex)
rm(y.date, y.var, y.drug, y.subjects, y.disable, y1, y2, y.manova)
rm(y3, Y, y.group, y.sum, tol, tol1, tol2, tol3)
T
}