lines = '"Trt" "Initial" "Final"
1 27.2 32.6
1 31.0 35.6
1 32.0 35.6
1 27.8 30.8
2 29.5 33.7
2 27.8 31.3
2 26.3 30.4
2 27.0 31.0
3 28.9 33.8
3 23.1 29.2
3 24.4 27.6
3 25.0 30.8
4 29.3 34.8
4 30.2 36.5
4 25.5 30.8
4 22.7 25.9'
con = textConnection(lines)
DF = read.table( con, header=TRUE )
close(con)

DF$Trt = as.factor( DF$Trt )
DF$WeightGain = DF$Final - DF$Initial


# The classical ANOVA test (from R in Action EPage 234)
#
fit = aov( WeightGain ~ Trt, data=DF )
sfit = summary(fit)
print( sprintf("Classical ANOVA: prob. type I error: %10.6f", sfit[[1]][["Pr(>F)"]][1]) )

#
# To study the main effects of Treatment plant we will apply a permuation test:
#
source('utils.R')

F = compute_F( DF$Trt, DF$WeightGain )

set.seed(1234)
B = 200 # the sumber of bootstraps to run
all_Fs = rep( NA, B ) # the bootstrap values of the F
for( bi in 1:B ){
    all_Fs[bi] = compute_F( DF$Trt, sample(DF$WeightGain) )
}

plot( density( all_Fs ), xlab="statistic", main="" )
grid()
abline(v=F, col='red')

alpha = sum( all_Fs >= F ) / B
print( sprintf("Main effects for Trt (B= %d): prob. type I error: %10.6f", B, alpha) )