The fracture load data

I = 3
J = 4

x1 = c(2.62, 2.99, 3.39, 2.86)
x2 = c(3.47, 3.85, 3.77, 3.63)
x3 = c(4.78, 4.41, 4.91, 5.06)

dat = data.frame(load = c(x1, x2, x3), level=rep(c("42mm", "36mm","31.2mm"), each=4))
print(dat)
##    load  level
## 1  2.62   42mm
## 2  2.99   42mm
## 3  3.39   42mm
## 4  2.86   42mm
## 5  3.47   36mm
## 6  3.85   36mm
## 7  3.77   36mm
## 8  3.63   36mm
## 9  4.78 31.2mm
## 10 4.41 31.2mm
## 11 4.91 31.2mm
## 12 5.06 31.2mm

The treament means / variance and grand mean can be computed

treat.mean = aggregate(dat$load, by=list(dat$level), FUN=mean)
treat.var = aggregate(dat$load, by=list(dat$level), FUN=var)
grand.mean = mean(dat$load)

The sum of squares and mean squares can be computed

sse = (J-1) * sum(treat.var$x) 
sstr = J * (I-1) * var(treat.mean$x)
mse = sse / (I*J-I)
mstr = sstr / (I-1)

The F-statistic and its p-value:

f.stat = mstr/mse 
p.value = 1 - pf(f.stat, I-1, I*J-I)
print(c(f.stat, p.value))
## [1] 4.857779e+01 1.504349e-05

It can be done by the aov function.

res.anova = aov(load~level, data=dat)
print(summary(res.anova))
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## level        2  6.765   3.383   48.58 1.5e-05 ***
## Residuals    9  0.627   0.070                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The Tukey’s test can be done by TukeyHSD function.

res.tukey = TukeyHSD(res.anova)
print(res.tukey)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = load ~ level, data = dat)
## 
## $level
##               diff       lwr        upr     p adj
## 36mm-31.2mm -1.110 -1.630967 -0.5890334 0.0005650
## 42mm-31.2mm -1.825 -2.345967 -1.3040334 0.0000115
## 42mm-36mm   -0.715 -1.235967 -0.1940334 0.0100950