Statistics for Laboratory Scientists ( 140.615 )

Analysis of Variance - Multiple Comparisons

Example from class: pea sections and sugar

The data.

x <- 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)
treat <- factor(rep(c("C","G","F","G+F","S"), rep(10,5)))
sugar <- data.frame(rsp=x, trt=treat)
head(sugar)

##   rsp trt
## 1  75   C
## 2  67   C
## 3  70   C
## 4  75   C
## 5  65   C
## 6  71   C

str(sugar)

## 'data.frame':    50 obs. of  2 variables:
##  $ rsp: num  75 67 70 75 65 71 67 67 76 68 ...
##  $ trt: Factor w/ 5 levels "C","F","G","G+F",..: 1 1 1 1 1 1 1 1 1 1 ...

summary(sugar)

##       rsp         trt    
##  Min.   :56.00   C  :10  
##  1st Qu.:58.00   F  :10  
##  Median :60.50   G  :10  
##  Mean   :61.94   G+F:10  
##  3rd Qu.:65.00   S  :10  
##  Max.   :76.00

Plot the data.

set.seed(1)
par(las=1)
stripchart(sugar$rsp ~ sugar$trt, method="jitter", pch=1, xlab="response", ylab="group")

Analysis of variance.

aov.out <- aov(rsp ~ trt, data=sugar)
summary(aov.out)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## trt          4 1077.3  269.33   49.37 6.74e-16 ***
## Residuals   45  245.5    5.46                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Tukey’s HSD

sugar.tukey <- TukeyHSD(aov.out)
par(las=1)
plot(sugar.tukey)
abline(v=0, lty=2, col="red")

sugar.tukey

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = rsp ~ trt, data = sugar)
## 
## $trt
##        diff        lwr       upr     p adj
## F-C   -11.9 -14.868072 -8.931928 0.0000000
## G-C   -10.8 -13.768072 -7.831928 0.0000000
## G+F-C -12.1 -15.068072 -9.131928 0.0000000
## S-C    -6.0  -8.968072 -3.031928 0.0000072
## G-F     1.1  -1.868072  4.068072 0.8291029
## G+F-F  -0.2  -3.168072  2.768072 0.9996878
## S-F     5.9   2.931928  8.868072 0.0000100
## G+F-G  -1.3  -4.268072  1.668072 0.7256157
## S-G     4.8   1.831928  7.768072 0.0003242
## S-G+F   6.1   3.131928  9.068072 0.0000052

Bonferroni corrected confidence intervals (a customized function), based on a pooled estimate of the within-group standard deviation.

library(devtools)

## Loading required package: usethis

devtools::install_github("bllfrg/SPH.140.615",quiet=TRUE)
library(SPH.140.615)
sugar.bonf <- ci.bonf(sugar$rsp, sugar$trt)
sugar.bonf

##          est     lower     upper
## C - F   11.9  8.816368 14.983632
## C - G   10.8  7.716368 13.883632
## C - G+F 12.1  9.016368 15.183632
## C - S    6.0  2.916368  9.083632
## F - G   -1.1 -4.183632  1.983632
## F - G+F  0.2 -2.883632  3.283632
## F - S   -5.9 -8.983632 -2.816368
## G - G+F  1.3 -1.783632  4.383632
## G - S   -4.8 -7.883632 -1.716368
## G+F - S -6.1 -9.183632 -3.016368

An alternative built-in function to get Bonferroni corrected pairwise comparisons p-values, based on a pooled estimate of the within-group standard deviation.

pairwise.t.test(sugar$rsp, sugar$trt, p.adjust="bonf")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  sugar$rsp and sugar$trt 
## 
##     C       F       G       G+F    
## F   7.5e-14 -       -       -      
## G   1.8e-12 1.00000 -       -      
## G+F 4.3e-14 1.00000 1.00000 -      
## S   7.5e-06 1.0e-05 0.00035 5.4e-06
## 
## P value adjustment method: bonferroni

The Newman-Keuls procedure, using a customized function.

newman.keuls(sugar$rsp, sugar$trt)