## Solutions

### Problem 1 ###

# Analysis of the original data
peak = read.table("peak.txt",header=TRUE)
y = peak$y; method = as.factor(peak$method)
obj1 = aov(y ~ method)
anova(obj1)
# Compare with Table 3-8
plot(resid(obj1) ~ fitted(obj1), pch=16)
abline(h=0, lty=2)
# Compare with Figure 3-7
bartlett.test(peak$y, peak$method)

# Finding the variance-stabilising transformation
m = tapply(y, method, mean)
S = tapply(y, method, sd)
plot(log(S) ~ log(m), pch=16)
lm(log(S) ~ log(m))
abline(.Last.value)
# The slope is 0.4465, but 0.5 is close enough, indicating the square
# root transformation according to Table 3-9.

# Analysis of the transformed data
obj2 = aov(sqrt(y) ~ method)
anova(obj2) 
# Hey! The degrees of freedom deviates from Table 3-10.
# Montgomery subtracts one degree of freedom, because the exponent
# (1/2 for square root) is estimated
MSE = 2.688/19; MSE
F0 = 10.895/MSE; F0
pf(F0, 3, 19, lower.tail=FALSE)
bartlett.test(sqrt(y), method)
plot(resid(obj2) ~ fitted(obj2), pch=16)
abline(h=0, lty=2)
# Compare with Figure 3.9

# Cleaning
rm(list=ls())


### Problem 2 ###

y = scan("pr0322.txt")
type = gl(3, 5)
type  # Checking the factor
## Question (a):
av = aov(y ~ type)
M = anova(av); M
# Three-starred significance (two stars were asked for)
## Question (b):
MSE = M[2,3]; MSE
tukeyspan = qtukey(0.95, 3, 12) * sqrt(MSE/5); tukeyspan
m = tapply(y, type, mean)
sort(m)
# Or, simpler
TukeyHSD(av, "type")
# Types 1 and 3 cannot be distinguished, but Type 2 is different
## Question (d):
# Notice that there is a reason for the choice of contrasts
# Hence the Scheffe protection is not required
c1 = c(1, 0, -1)
c2 = c(-1, 2, -1)
C1 = sum(c1*m); C2 = sum(c2*m)
c(C1, C2)
SS1 = C1^2/(sum(c1^2)/5); SS2 = C2^2/(sum(c2^2)/5)
c(SS1, SS2)
# Their sum is SS_Treatment; SS2 takes the main part
F1 = SS1/MSE; F2 = SS2/MSE
c(F1, F2)
1 - pf(c(F1,F2), 1, 12)  
# As expected, the first contrast is not significant
## Question (e):
# Type 1 and Type 3 are candidates for minimal response time
# More experiments are desirable, but forced to make a decision based
# on these data, Type 3 would be chosen -- with the risk of criticism
# if Type 1 actually shows to be better
# Question (f):
res = resid(av)
yhat = fitted(av)
plot(res ~ yhat)
abline(h=0, lty=2)
# A nice residual plot
# No discrepancy in variance, as confirmed by
bartlett.test(y, type)  # p=0.567
# Cleaning
rm(list=ls())