lines = 'Subject Prestress Five_Minute Ten_Minute Treatment
A 9.3 11.7 10.5 Brand_A
B 8.4 10.0 10.5 Brand_A
C 7.8 10.4 9.0 Brand_A
D 7.5 9.2 9.0 New_Drug
E 8.9 9.5 10.2 New_Drug
F 8.3 9.5 9.5 New_Drug'
con = textConnection( lines )
DF = read.table( con, header=TRUE )
close(con)

DF$Treatment = as.factor( DF$Treatment )

# Normalize the stress observed by the prestress baseline:
DF$Five_Minute = ( DF$Five_Minute - DF$Prestress )
DF$Ten_Minute = ( DF$Ten_Minute - DF$Prestress )

# Do one-factor AVNOVA with the Treatment a the dependent variable:
#
fit = aov( Five_Minute ~ Treatment, data=DF )
sfit = summary(fit)
print( sprintf("Main Effects (Five Minute): prob. type I error: %10.6f", sfit[[1]][["Pr(>F)"]][1]) )

fit = aov( Ten_Minute ~ Treatment, data=DF )
sfit = summary(fit)
print( sprintf("Main Effects (Ten Minute): prob. type I error: %10.6f", sfit[[1]][["Pr(>F)"]][1]) )

#
# Use a permutation test based on either T^2 or Wald and Wolfowitz's T' statistic (which is easier to compute):
#

# Get the pooled feature means U:
U = colMeans( DF[,3:4] )

# Compute the matrix C:
#
Five_Min_D = DF$Five_Minute - U[1] # this is X_i - U_i
Ten_Min_D = DF$Ten_Minute - U[2]

XmU = data.frame( Treatment=DF$Treatment, Five_Min_D=Five_Min_D, Ten_Min_D=Ten_Min_D ) # X m(inus) U

c = matrix( 0, nrow=2, ncol=2 )
for( ii in 1:2 ){
    for( jj in 1:2 ){
        c[ii,jj] = mean( XmU[,ii+1] * XmU[,jj+1] )
    }
}

# Get the inverse of the matrix C (the same for every resampling):
#
c_inv = solve(c)

# Lets now compute bootstrapped replications of The {T'}^2 statistic:
#
compute_Tprime = function( DF, c_inv ){
    # Get the means in each class under the given ordering:
    #
    mask = DF$Treatment == 'Brand_A'
    u_N = colMeans( DF[mask,3:4] )
    mask = DF$Treatment == 'New_Drug'
    u_Y = colMeans( DF[mask,3:4] )

    TP2 = t( u_N - u_Y ) %*% c_inv %*% (u_N - u_Y)

    return(as.double(TP2))
}

T0 = compute_Tprime( DF, c_inv )

# Do bootstraps:
#
set.seed(1234)
B = 500 # the number of bootstraps to run
all_T2s = rep( NA, B ) # the bootstrap values of the T^2
for( bi in 1:B ){
    DF$Treatment = sample(DF$Treatment) # reorder the labels
    all_T2s[bi] = compute_Tprime( DF, c_inv )
}


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

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