#### Sample Code for Two-Way ANOVA
rm(list=ls(all=TRUE))
#### Battery Data. Life time per dollar of batteries from Dean and Voss (1999).
battery <- read.table("http://www.unc.edu/~zhangk/STOR665Spring15/battery.dat", header=T)
battery ### Balanced design. Four observations in each Duty/Brand cell.
## Duty Brand Cost
## 1 Alkaline Name 611
## 2 Alkaline Name 537
## 3 Alkaline Name 542
## 4 Alkaline Name 593
## 5 Alkaline Store 923
## 6 Alkaline Store 794
## 7 Alkaline Store 827
## 8 Alkaline Store 898
## 9 Heavy Name 445
## 10 Heavy Name 490
## 11 Heavy Name 384
## 12 Heavy Name 413
## 13 Heavy Store 476
## 14 Heavy Store 569
## 15 Heavy Store 480
## 16 Heavy Store 460
par(mfrow=c(1,1),mar=c(3,3,2,2),mgp=c(1.8,0.5,0))
interaction.plot(battery$Brand,battery$Duty,battery$Cost)
title("Interaction Plot", sub="Battery")
options()$contrasts ### Current contrasts. The first for unordered and the second for ordered.
## unordered ordered
## "contr.treatment" "contr.poly"
options(contrasts=c("contr.sum", "contr.sum")) ### Setting the sum of alphas, betas, and gammas to be 0.
summary(aov(Cost~Brand*Duty, data=battery)) ### Complete model
## Df Sum Sq Mean Sq F value Pr(>F)
## Brand 1 124609 124609 52.63 1.01e-05 ***
## Duty 1 252004 252004 106.43 2.55e-07 ***
## Brand:Duty 1 51302 51302 21.67 0.000556 ***
## Residuals 12 28412 2368
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
X=model.matrix(aov(Cost~Duty*Brand,data=battery)) ### In balanced designs, the matrix for the two variables is orthogonal.
Y=battery$Cost
XtX= t(X)%*%X
XtXinv = solve(XtX)
beta_hat =XtXinv %*% t(X) %*% Y
beta_hat
## [,1]
## (Intercept) 590.125
## Duty1 125.500
## Brand1 -88.250
## Duty1:Brand1 -56.625
summary(lm(Cost~Brand*Duty, data=battery)) ### Estimates of mu, alphas, betas, and gammas.
##
## Call:
## lm(formula = Cost ~ Brand * Duty, data = battery)
##
## Residuals:
## Min 1Q Median 3Q Max
## -66.50 -33.56 -18.12 38.19 72.75
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 590.12 12.16 48.511 3.85e-15 ***
## Brand1 -88.25 12.16 -7.255 1.01e-05 ***
## Duty1 125.50 12.16 10.317 2.55e-07 ***
## Brand1:Duty1 -56.63 12.16 -4.655 0.000556 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 48.66 on 12 degrees of freedom
## Multiple R-squared: 0.9377, Adjusted R-squared: 0.9222
## F-statistic: 60.24 on 3 and 12 DF, p-value: 1.662e-07
#main effect model
model.matrix(aov(Cost~Duty+Brand,data=battery)) ### In balanced designs, the matrix for the two variables is orthogonal.
## (Intercept) Duty1 Brand1
## 1 1 1 1
## 2 1 1 1
## 3 1 1 1
## 4 1 1 1
## 5 1 1 -1
## 6 1 1 -1
## 7 1 1 -1
## 8 1 1 -1
## 9 1 -1 1
## 10 1 -1 1
## 11 1 -1 1
## 12 1 -1 1
## 13 1 -1 -1
## 14 1 -1 -1
## 15 1 -1 -1
## 16 1 -1 -1
## attr(,"assign")
## [1] 0 1 2
## attr(,"contrasts")
## attr(,"contrasts")$Duty
## [1] "contr.sum"
##
## attr(,"contrasts")$Brand
## [1] "contr.sum"
summary(aov(Cost~Brand+Duty, data=battery)) ### So the ANOVA tables are the same.
## Df Sum Sq Mean Sq F value Pr(>F)
## Brand 1 124609 124609 20.32 0.000589 ***
## Duty 1 252004 252004 41.10 2.3e-05 ***
## Residuals 13 79715 6132
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(Cost~Duty+Brand, data=battery))
## Df Sum Sq Mean Sq F value Pr(>F)
## Duty 1 252004 252004 41.10 2.3e-05 ***
## Brand 1 124609 124609 20.32 0.000589 ***
## Residuals 13 79715 6132
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lm(Cost~Brand+Duty, data=battery)) ### Estimates of mu, alphas, and betas
##
## Call:
## lm(formula = Cost ~ Brand + Duty, data = battery)
##
## Residuals:
## Min 1Q Median 3Q Max
## -92.875 -73.875 -1.125 44.625 119.125
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 590.12 19.58 30.144 2.04e-13 ***
## Brand1 -88.25 19.58 -4.508 0.000589 ***
## Duty1 125.50 19.58 6.411 2.30e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 78.31 on 13 degrees of freedom
## Multiple R-squared: 0.8253, Adjusted R-squared: 0.7984
## F-statistic: 30.71 on 2 and 13 DF, p-value: 1.188e-05
options(contrasts=c("contr.treatment", "contr.treatment")) ### Set the first level (alpha1, beta1) of each factor to be 0
model.matrix(aov(Cost~Duty+Brand,data=battery))
## (Intercept) DutyHeavy BrandStore
## 1 1 0 0
## 2 1 0 0
## 3 1 0 0
## 4 1 0 0
## 5 1 0 1
## 6 1 0 1
## 7 1 0 1
## 8 1 0 1
## 9 1 1 0
## 10 1 1 0
## 11 1 1 0
## 12 1 1 0
## 13 1 1 1
## 14 1 1 1
## 15 1 1 1
## 16 1 1 1
## attr(,"assign")
## [1] 0 1 2
## attr(,"contrasts")
## attr(,"contrasts")$Duty
## [1] "contr.treatment"
##
## attr(,"contrasts")$Brand
## [1] "contr.treatment"
summary(aov(Cost~Duty+Brand, data=battery)) ### The ANOVA table remains the same
## Df Sum Sq Mean Sq F value Pr(>F)
## Duty 1 252004 252004 41.10 2.3e-05 ***
## Brand 1 124609 124609 20.32 0.000589 ***
## Residuals 13 79715 6132
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lm(Cost~Duty+Brand, data=battery)) ### Coefficient estimates.
##
## Call:
## lm(formula = Cost ~ Duty + Brand, data = battery)
##
## Residuals:
## Min 1Q Median 3Q Max
## -92.875 -73.875 -1.125 44.625 119.125
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 627.38 33.91 18.502 1.01e-10 ***
## DutyHeavy -251.00 39.15 -6.411 2.30e-05 ***
## BrandStore 176.50 39.15 4.508 0.000589 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 78.31 on 13 degrees of freedom
## Multiple R-squared: 0.8253, Adjusted R-squared: 0.7984
## F-statistic: 30.71 on 2 and 13 DF, p-value: 1.188e-05
#对比两种contrasts
options(contrasts=c("contr.sum", "contr.sum"));options()$contrasts
## [1] "contr.sum" "contr.sum"
summary(aov(Cost~Brand*Duty, data=battery))
## Df Sum Sq Mean Sq F value Pr(>F)
## Brand 1 124609 124609 52.63 1.01e-05 ***
## Duty 1 252004 252004 106.43 2.55e-07 ***
## Brand:Duty 1 51302 51302 21.67 0.000556 ***
## Residuals 12 28412 2368
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
options(contrasts=c("contr.treatment", "contr.treatment"));options()$contrasts
## [1] "contr.treatment" "contr.treatment"
summary(aov(Cost~Brand*Duty, data=battery))
## Df Sum Sq Mean Sq F value Pr(>F)
## Brand 1 124609 124609 52.63 1.01e-05 ***
## Duty 1 252004 252004 106.43 2.55e-07 ***
## Brand:Duty 1 51302 51302 21.67 0.000556 ***
## Residuals 12 28412 2368
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
options(contrasts=c("contr.sum", "contr.sum"))
summary(lm(Cost~Brand*Duty, data=battery))
##
## Call:
## lm(formula = Cost ~ Brand * Duty, data = battery)
##
## Residuals:
## Min 1Q Median 3Q Max
## -66.50 -33.56 -18.12 38.19 72.75
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 590.12 12.16 48.511 3.85e-15 ***
## Brand1 -88.25 12.16 -7.255 1.01e-05 ***
## Duty1 125.50 12.16 10.317 2.55e-07 ***
## Brand1:Duty1 -56.63 12.16 -4.655 0.000556 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 48.66 on 12 degrees of freedom
## Multiple R-squared: 0.9377, Adjusted R-squared: 0.9222
## F-statistic: 60.24 on 3 and 12 DF, p-value: 1.662e-07
options(contrasts=c("contr.treatment", "contr.treatment"))
summary(lm(Cost~Brand*Duty, data=battery))
##
## Call:
## lm(formula = Cost ~ Brand * Duty, data = battery)
##
## Residuals:
## Min 1Q Median 3Q Max
## -66.50 -33.56 -18.12 38.19 72.75
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 570.75 24.33 23.459 2.15e-11 ***
## BrandStore 289.75 34.41 8.421 2.21e-06 ***
## DutyHeavy -137.75 34.41 -4.004 0.001751 **
## BrandStore:DutyHeavy -226.50 48.66 -4.655 0.000556 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 48.66 on 12 degrees of freedom
## Multiple R-squared: 0.9377, Adjusted R-squared: 0.9222
## F-statistic: 60.24 on 3 and 12 DF, p-value: 1.662e-07
#对比两种contrasts
library(faraway)
data(pvc)
options(contrasts=c("contr.sum", "contr.sum"));options()$contrasts
## [1] "contr.sum" "contr.sum"
summary(aov(psize~operator*resin, data=pvc))
## Df Sum Sq Mean Sq F value Pr(>F)
## operator 2 20.72 10.36 7.007 0.00401 **
## resin 7 283.95 40.56 27.439 5.66e-10 ***
## operator:resin 14 14.34 1.02 0.693 0.75987
## Residuals 24 35.48 1.48
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
options(contrasts=c("contr.treatment", "contr.treatment"));options()$contrasts
## [1] "contr.treatment" "contr.treatment"
summary(aov(psize~operator*resin, data=pvc))
## Df Sum Sq Mean Sq F value Pr(>F)
## operator 2 20.72 10.36 7.007 0.00401 **
## resin 7 283.95 40.56 27.439 5.66e-10 ***
## operator:resin 14 14.34 1.02 0.693 0.75987
## Residuals 24 35.48 1.48
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
options(contrasts=c("contr.sum", "contr.sum"))
summary(aov(psize~operator*resin, data=pvc))
## Df Sum Sq Mean Sq F value Pr(>F)
## operator 2 20.72 10.36 7.007 0.00401 **
## resin 7 283.95 40.56 27.439 5.66e-10 ***
## operator:resin 14 14.34 1.02 0.693 0.75987
## Residuals 24 35.48 1.48
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
options(contrasts=c("contr.treatment", "contr.treatment"))
summary(aov(psize~operator*resin, data=pvc))
## Df Sum Sq Mean Sq F value Pr(>F)
## operator 2 20.72 10.36 7.007 0.00401 **
## resin 7 283.95 40.56 27.439 5.66e-10 ***
## operator:resin 14 14.34 1.02 0.693 0.75987
## Residuals 24 35.48 1.48
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#### Female Labor Supply Data (unbalanced design)
rm(list=ls(all=TRUE))
labor <- read.table("http://www.unc.edu/~zhangk/STOR665Spring15/labor.dat", header=T)
earning <- labor[,3]
edu <- labor[,5]
eduf <- rep("A", length(edu))
eduf[edu<16] <- "B"
eduf[edu<13] <- "C"
eduf[edu<12] <- "D"
eduf <- factor(eduf)
child <- factor(labor[,8])
table(child,eduf) ### Unbalanced design
## eduf
## child A B C D
## 0 18 67 51 95
## 1 52 140 149 35
interaction.plot(eduf,child,earning)
title("Interaction Plot", sub="Female Labor Supply")
options(contrasts=c("contr.sum", "contr.sum"))
summary(aov(earning~child+eduf)) ### Different SS
## Df Sum Sq Mean Sq F value Pr(>F)
## child 1 4 4.1 0.256 0.613
## eduf 3 1936 645.5 40.690 <2e-16 ***
## Residuals 602 9550 15.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(earning~eduf+child))
## Df Sum Sq Mean Sq F value Pr(>F)
## eduf 3 1885 628.3 39.605 <2e-16 ***
## child 1 56 55.7 3.511 0.0614 .
## Residuals 602 9550 15.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### Adding interaction
summary(aov(earning~child*eduf)) ### Different SS for child and eduf. same interaction
## Df Sum Sq Mean Sq F value Pr(>F)
## child 1 4 4.1 0.257 0.613
## eduf 3 1936 645.5 40.852 <2e-16 ***
## child:eduf 3 85 28.4 1.798 0.146
## Residuals 599 9465 15.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(earning~eduf*child))
## Df Sum Sq Mean Sq F value Pr(>F)
## eduf 3 1885 628.3 39.762 <2e-16 ***
## child 1 56 55.7 3.525 0.0609 .
## eduf:child 3 85 28.4 1.798 0.1464
## Residuals 599 9465 15.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### Coeffcient estimates are the same
summary(lm(earning~child+eduf))
##
## Call:
## lm(formula = earning ~ child + eduf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.866 -2.160 -0.118 1.738 58.143
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.5905 0.1818 58.245 < 2e-16 ***
## child1 0.3375 0.1801 1.874 0.0614 .
## eduf1 3.6473 0.3837 9.506 < 2e-16 ***
## eduf2 0.2741 0.2654 1.032 0.3023
## eduf3 -1.3106 0.2713 -4.830 1.73e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.983 on 602 degrees of freedom
## Multiple R-squared: 0.1689, Adjusted R-squared: 0.1634
## F-statistic: 30.58 on 4 and 602 DF, p-value: < 2.2e-16
summary(lm(earning~eduf+child))
##
## Call:
## lm(formula = earning ~ eduf + child)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.866 -2.160 -0.118 1.738 58.143
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.5905 0.1818 58.245 < 2e-16 ***
## eduf1 3.6473 0.3837 9.506 < 2e-16 ***
## eduf2 0.2741 0.2654 1.032 0.3023
## eduf3 -1.3106 0.2713 -4.830 1.73e-06 ***
## child1 0.3375 0.1801 1.874 0.0614 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.983 on 602 degrees of freedom
## Multiple R-squared: 0.1689, Adjusted R-squared: 0.1634
## F-statistic: 30.58 on 4 and 602 DF, p-value: < 2.2e-16
summary(lm(earning~child*eduf))
##
## Call:
## lm(formula = earning ~ child * eduf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.926 -2.199 -0.210 1.792 57.428
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.5656 0.2002 52.786 < 2e-16 ***
## child1 0.3610 0.2002 1.804 0.071773 .
## eduf1 3.6161 0.4333 8.345 4.93e-16 ***
## eduf2 0.1244 0.2892 0.430 0.667188
## eduf3 -1.0509 0.3034 -3.464 0.000571 ***
## child1:eduf1 -0.1389 0.4333 -0.321 0.748626
## child1:eduf2 -0.5185 0.2892 -1.793 0.073542 .
## child1:eduf3 0.4558 0.3034 1.502 0.133556
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.975 on 599 degrees of freedom
## Multiple R-squared: 0.1763, Adjusted R-squared: 0.1667
## F-statistic: 18.32 on 7 and 599 DF, p-value: < 2.2e-16
summary(lm(earning~eduf*child))
##
## Call:
## lm(formula = earning ~ eduf * child)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.926 -2.199 -0.210 1.792 57.428
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.5656 0.2002 52.786 < 2e-16 ***
## eduf1 3.6161 0.4333 8.345 4.93e-16 ***
## eduf2 0.1244 0.2892 0.430 0.667188
## eduf3 -1.0509 0.3034 -3.464 0.000571 ***
## child1 0.3610 0.2002 1.804 0.071773 .
## eduf1:child1 -0.1389 0.4333 -0.321 0.748626
## eduf2:child1 -0.5185 0.2892 -1.793 0.073542 .
## eduf3:child1 0.4558 0.3034 1.502 0.133556
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.975 on 599 degrees of freedom
## Multiple R-squared: 0.1763, Adjusted R-squared: 0.1667
## F-statistic: 18.32 on 7 and 599 DF, p-value: < 2.2e-16
#### alloy data
alloy <- matrix(scan(file="http://www.unc.edu/~zhangk/STOR665Spring15/alloy.txt"), ncol=7, byrow=T)
Strength <- as.vector(alloy[,2:7])
Lab <- factor(c(rep(1,4), rep(2,4), rep(3,4), rep(4,4), rep(5,4), rep(6,4)))
Bar <- factor(rep(1:4,6))
alloy <- data.frame(Strength, Bar, Lab)
alloy ### One observation per Lab/Bar cell
## Strength Bar Lab
## 1 3.06 1 1
## 2 2.59 2 1
## 3 3.32 3 1
## 4 3.08 4 1
## 5 3.01 1 2
## 6 2.75 2 2
## 7 2.93 3 2
## 8 2.87 4 2
## 9 2.84 1 3
## 10 3.13 2 3
## 11 3.27 3 3
## 12 2.76 4 3
## 13 3.48 1 4
## 14 3.58 2 4
## 15 3.80 3 4
## 16 2.94 4 4
## 17 3.24 1 5
## 18 2.35 2 5
## 19 3.55 3 5
## 20 2.24 4 5
## 21 2.42 1 6
## 22 2.58 2 6
## 23 2.77 3 6
## 24 2.58 4 6
interaction.plot(Bar,Lab,Strength)
title("Interaction Plot", sub="Alloy")
summary(aov(Strength~Bar+Lab, alloy))
## Df Sum Sq Mean Sq F value Pr(>F)
## Bar 3 0.9814 0.3271 3.978 0.0286 *
## Lab 5 1.6049 0.3210 3.904 0.0182 *
## Residuals 15 1.2334 0.0822
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(Strength~Bar*Lab, alloy))
## Df Sum Sq Mean Sq
## Bar 3 0.9814 0.3271
## Lab 5 1.6049 0.3210
## Bar:Lab 15 1.2334 0.0822
###
Lab.mean <- c(mean(alloy$Strength[alloy$Lab==1]), mean(alloy$Strength[alloy$Lab==2]),mean(alloy$Strength[alloy$Lab==3]),mean(alloy$Strength[alloy$Lab==4]),mean(alloy$Strength[alloy$Lab==5]),mean(alloy$Strength[alloy$Lab==6]))
Bar.mean <- c(mean(alloy$Strength[alloy$Bar==1]), mean(alloy$Strength[alloy$Bar==2]),mean(alloy$Strength[alloy$Bar==3]),mean(alloy$Strength[alloy$Bar==4]))
tukey.1df <- function(y, factorA, factorB) {
factorA <- as.factor(factorA)
factorB <- as.factor(factorB)
a <- length(unique(factorA))
b <- length(unique(factorB))
Y <- matrix(nrow=a, ncol=b)
for(i in 1:a) {for (j in 1:b) Y[i,j] <- y[factorA == unique(factorA)[i] & factorB == unique(factorB)[j]]}
Ameans <- apply(Y, 1, mean)
Bmeans <- apply(Y, 2, mean)
term1 <- t(Ameans)%*%Y%*%Bmeans
y.. <- mean(Y)
SSA <- b*sum((Ameans-y..)^2)
SSB <- a*sum((Bmeans-y..)^2)
term2 <- SSA + SSB + a*b*y..^2
SSN <- a*b*((term1 - y..*term2)^2)/(SSA*SSB)
SST <- sum((Y-y..)^2)
SSE <- SST-SSN-SSA-SSB
SS <- c(SSA, SSB, SSN, SSE, SST)
df <- round(c(a-1, b-1, 1, (a-1)*(b-1)-1, a*b-1))
MS <- SS[1:4]/df[1:4]
F0 <- c(MS[1]/MS[4], MS[2]/MS[4], MS[3]/MS[4])
p <- vector(length=3)
for(i in 1:3) p[i] <-1-pf(F0[i], df[i], df[4])
SS <- round(SS,3)
MS <- round(MS,3)
F0 <- round(F0,3)
p <- round(p,4)
MS <- c(MS, " ")
F0 <- c(F0, " ", " ")
p <- c(p, " ", " ")
out <- cbind(SS, df, MS, F0, p)
dimnames(out) <- list(c("A","B","N","Err","Tot"),c("SS","df","MS", "F0", "p"))
#print.matrix(out, quote=F)
list(out = out)
}
tukey.1df(alloy$Strength, alloy$Bar, alloy$Lab)
## $out
## SS df MS F0 p
## A "0.981" "3" "0.327" "3.772" "0.0356"
## B "1.605" "5" "0.321" "3.701" "0.0241"
## N "0.019" "1" "0.019" "0.221" "0.6452"
## Err "1.214" "14" "0.087" " " " "
## Tot "3.82" "23" " " " " " "