DF_first_month = data.frame(
    laboratory=as.factor( rep( c( 'A', 'B', 'C', 'D', 'E', 'F' ), ) ),
    month=as.factor( c( rep( '1', 6 ), rep( '1', 6 ), rep( '1', 6 ) ) ), 
    value=c(
        221.1, 208.8, 211.1, 208.3, 221.1, 224.2,
        224.2, 206.9, 198.4, 214.1, 208.8, 206.9,
        217.8, 205.9, 213.0, 209.1, 211.1, 198.4
    )
)
DF = DF_first_month

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

# Next lets apply a permutation test and verify that we make the same conclusion:
#
compute_F = function( group_names, data ){
    #
    # Compute the F_1 metric
    #
    ugroups = unique( group_names )
    ng = length(ugroups)

    x_bar_dot_dot = mean( data ) # the grand mean
    x_bar_i_dot = rep( NA, ng )
    for( gi in 1:ng ){
        d = data[ group_names == ugroups[gi] ]
        x_bar_i_dot[gi] = mean(d)
    }
    F = sum( abs( x_bar_i_dot - x_bar_dot_dot ) )
    return(F)
}

F = compute_F( DF$laboratory, DF$value ) # the value of F under the original ordering
set.seed(1234)
B = 200 # the sumber of bootstraps to run
all_Fs = rep( NA, B ) # the bootstrap values of the F_1
for( bi in 1:B ){
    all_Fs[bi] = compute_F( sample(DF$laboratory), DF$value ) # reorder the labels
}

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

alpha = sum( all_Fs >= F ) / B
print( sprintf("Permutation Bootstrap (laboratory ME) (B: %d): prob. type I error: %10.6f", B, alpha) )



DF_second_month = data.frame(
    laboratory=as.factor( rep( c( 'A', 'B', 'C', 'D', 'E', 'F' ), ) ),
    month=as.factor( c( rep( '2', 6 ), rep( '2', 6 ), rep( '2', 6 ) ) ), 
    value=c(
        208.8, 211.4, 208.9, 207.7, 208.3, 214.1,
        212.6, 205.8, 206.0, 216.2, 208.8, 212.6,
        213.3, 202.5, 209.8, 203.7, 211.4, 205.8
    )
)

# Combine these two dataframes into one:
# 
DF = rbind( DF_first_month, DF_second_month )

# Do one-factor AVNOVA with the month a the dependent variable:
#
fit = aov( value ~ month, data=DF )
sfit = summary(fit)
print( sprintf("Classical ANOVA (month ME): prob. type I error: %10.6f", sfit[[1]][["Pr(>F)"]][1]) )

# Next lets apply a permutation test and verify that we make the same conclusion:
#
F = compute_F( DF$laboratory, DF$value ) # the value of F under original ordering

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

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

alpha = sum( all_Fs >= F ) / B
print( sprintf("Permutation Bootstrap (month ME) (B: %d) prob type I error: %10.6f", B, alpha) )