DF = data.frame(
    light_level=as.factor( c( 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3) ), 
    water_level=as.factor( c( 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 2, 2, 2) ), 
    yield=c( 12, 8, 8, 13, 15, 16, 16, 12, 13, 19, 16, 17, 17, 24, 22, 24, 26, 28 )
)

# To start with, consider the classical two-factor ANOVA test (from R in Action EPage 234)
#
fit = aov( yield ~ light_level, data=DF )
print( sprintf("Classical ANOVA for light level") )
print(summary(fit))

fit = aov( yield ~ water_level, data=DF )
print( sprintf("Classical ANOVA for water level") )
print(summary(fit))

fit = aov( yield ~ light_level*water_level, data=DF )
print( sprintf("C1assical ANOVA for an interaction of light level and water level"))
print(summary(fit))

#--
# To study the main effects of sunlight we will apply a permuation test to the blocks of sunlight:
#--
source('utils.R')

F = compute_F( DF$light_level, DF$yield )

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

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

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

#--
# To study the main effects of water level we will apply the permuation test:
#--
F = compute_F( DF$water_level, DF$yield, use_F1=FALSE ) # the value of F under the original ordering

set.seed(1234)
B = 100 # 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( DF$water_level, sample(DF$yield), use_F1=FALSE )
}

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

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

#--
# Test for an interaction term:
#--
x_bar_dot_dot = mean( DF$yield ) # the grand mean
group_I = unique( DF$light_level )
I = length( group_I )
x_bar_i_dot_dot = rep( NA, I )
for( i in 1:I ){
    inds = DF$light_level == group_I[i]
    x_bar_i_dot_dot[i] = mean( DF$yield[ inds ] )
}

group_J = unique( DF$water_level )
J = length( group_J )
x_bar_dot_j_dot = rep( NA, J )
for( j in 1:J ){
    inds = DF$water_level == group_J[j]
    x_bar_dot_j_dot[j] = mean( DF$yield[ inds ] )
}

# Compute the residuals X_ijk':
#
DF_residual = DF
for( row in 1:dim(DF_residual)[1] ){
    grp_i = which( DF_residual[row,'light_level'] == group_I )
    grp_j = which( DF_residual[row, 'water_level'] == group_J )
    DF_residual[row, 'yield'] = DF_residual[row, 'yield'] - x_bar_i_dot_dot[grp_i] - x_bar_dot_j_dot[grp_j] + x_bar_dot_dot
}

# Compute the test statistics "I" using the residual data above:
#
compute_I = function( DF, group_I, group_J ){
    I_statistic = 0
    for( i in 1:length(group_I) ){
        for( j in 1:length(group_J) ){
            inds = (DF$light_level == group_I[i]) & (DF$water_level == group_J[j])
            s = sum( DF$yield[inds] )
            #print( sprintf( 'i= %d; j= %d; s= %f; s^2= %f', i, j, s, s^2 ) )
            I_statistic = I_statistic + s^2
        }
    }
    return(I_statistic)
}

I0 = compute_I( DF_residual, group_I, group_J ) # the original value of the statistic "I"

set.seed(1234)
B = 100 # the sumber of bootstraps to run
DF_boot = DF_residual
all_Is = rep( NA, B )
for( bi in 1:B ){
    DF_boot$yield = sample(DF_boot$yield)
    all_Is[bi] = compute_I( DF_boot, group_I, group_J )
}

# the bootstrap values of the F_1
#
plot( density( all_Is ), xlab="statistic", main="" )
grid()
abline(v=I0, col='red')

alpha = sum( all_Is <= I0 ) / B
print( sprintf("Interation effects (B= %d): prob. type I error: %10.6f", B, alpha) )