#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
#
#-----

source('utils.R')

# Ex 11.27
#
I = 3
J = 3
K = 3
L = 2

SST = 270024.33
SSA = 14144.44 # the main effects SS
SSB = 5111.27
SSC = 244696.39
SSAB = 1069.62 # the two-factor interaction SS
SSAC = 62.67
SSBC = 331.67
SSE = 3127.50
SSABC = SST - SSA - SSB - SSC - SSAB - SSAC - SSBC - SSE # compute SSABC from the others

# Compute the mean squares:
#
MST = SST / ( I * J * K * L - 1 )
MSA = SSA / ( I - 1 )
MSB = SSB / ( J - 1 )
MSC = SSC / ( K - 1 )
MSAB = SSAB / ( (I-1)*(J-1) )
MSAC = SSAC / ( (I-1)*(K-1) )
MSBC = SSBC / ( (J-1)*(K-1) )
MSABC = SSABC / ( (I-1)*(J-1)*(K-1) )
MSE = SSE / ( I*J*K*(L-1) )

# Compute the f values:
#
f_A = MSA / MSE
f_B = MSB / MSE
f_C = MSC / MSE
f_AB = MSAB / MSE
f_AC = MSAC / MSE
f_BC = MSBC / MSE
f_ABC = MSABC / MSE

# Compute the P-values from these f-values:
#
p_value_A = 1 - pf( f_A, I-1, I*J*K*(L-1) )
p_value_B = 1 - pf( f_B, J-1, I*J*K*(L-1) )
p_value_C = 1 - pf( f_C, K-1, I*J*K*(L-1) )
p_value_AB = 1 - pf( f_AB, (I-1)*(J-1), I*J*K*(L-1) )
p_value_AC = 1 - pf( f_AC, (I-1)*(K-1), I*J*K*(L-1) )
p_value_BC = 1 - pf( f_BC, (J-1)*(K-1), I*J*K*(L-1) )
p_value_ABC = 1 - pf( f_ABC, (I-1)*(J-1)*(K-1), I*J*K*(L-1) )

print( c( p_value_A, p_value_B, p_value_C ) )
print( c( p_value_AB, p_value_AC, p_value_BC ) )
print( p_value_ABC )


x_dot_dot_k = c( 8242, 9732, 11210 )
xbar_dot_dot_k = x_dot_dot_k / ( I * J * L )
q_crit = qtukey( 1-0.05, K, I*J*K*(L-1) )
ss = sqrt( MSE/(I*J*L) )
ws = q_crit * ss

diff( sort( xbar_dot_dot_k ) )


# Ex 11.28
#
I = 4
J = 3
K = 2
L = 2
SST = 2983164.81
SSA = 19149.73
SSB = 2589047.62
SSC = 157437.52
SSAB = 53238.21
SSAC = 9033.73
SSBC = 91880.04
SSE = 56819.50
SSABC = SST - SSA - SSB - SSC - SSAB - SSAC - SSBC - SSE # compute SSABC from the others

# Compute the mean squares:
#
MST = SST / ( I * J * K * L -  1 ) 
MSA = SSA / ( I - 1 )
MSB = SSB / ( J - 1 )
MSC = SSC / ( K - 1 )
MSAB = SSAB / ( (I-1)*(J-1) )
MSAC = SSAC / ( (I-1)*(K-1) )
MSBC = SSBC / ( (J-1)*(K-1) )
MSABC = SSABC / ( (I-1)*(J-1)*(K-1) )
MSE = SSE / ( I*J*K*(L-1) )

# Compute the f values:
#
f_A = MSA / MSE
f_B = MSB / MSE
f_C = MSC / MSE
f_AB = MSAB / MSE
f_AC = MSAC / MSE
f_BC = MSBC / MSE
f_ABC = MSABC / MSE

# Compute the P-values from these f-values:
#
p_value_A = 1 - pf( f_A, I-1, I*J*K*(L-1) )
p_value_B = 1 - pf( f_B, J-1, I*J*K*(L-1) )
p_value_C = 1 - pf( f_C, K-1, I*J*K*(L-1) )
p_value_AB = 1 - pf( f_AB, (I-1)*(J-1), I*J*K*(L-1) )
p_value_AC = 1 - pf( f_AC, (I-1)*(K-1), I*J*K*(L-1) )
p_value_BC = 1 - pf( f_BC, (J-1)*(K-1), I*J*K*(L-1) )
p_value_ABC = 1 - pf( f_ABC, (I-1)*(J-1)*(K-1), I*J*K*(L-1) )

print( c( p_value_A, p_value_B, p_value_C ) )
print( c( p_value_AB, p_value_AC, p_value_BC ) )
print( p_value_ABC )


# Ex 11.29
#
I = 3
J = 2
K = 4
L = 4
SST = 1037.833
SSAB = 1.646
SSAC = 71.021
SSBC = 1.542
SSE = 447.5

data = c( 6, 9, 7, 9, 1, 2, 6, 6, 
          1, 3, 5, 5, 0, 4, 7, 3, 
          6, 3, 8, 7, 3, 2, 7, 9, 
          1, -1, 4, 8, 1, 0, 11, 6, 
          5, 4, 10, 11, -1, 2, 10, 5,
          9, 6, 6, 4, 6, 1, 4, 8,
          4, 6, 6, 5, -1, 0, 4, 5,
          0, 1, 3, 4, 0, 1, 5, 4,
          3, 1, 6, 4, 2, 0, 9, 4,
          1, -2, 1, 3, -1, 1, 6, 3,
          6, 0, 8, 7, 0, -2, 4, 3,
          3, 7, 10, 0, 4, -4, 7, 0 )

DF = data.frame( length=data )

# Lay down the factors:
A_factor = rep( 1:I, each=K*L, times=J )
B_factor = rep( 1:J, each=I*K*L )
C_factor = rep( 1:K, each=2, times=2*I*J )

DF$A_factor = A_factor
DF$B_factor = B_factor
DF$C_factor = C_factor

# Compute SSA, SSB, and SSC from the above data:
#
res = THREE_FACTOR_ANOVA( DF$length, DF$A_factor, DF$B_factor, DF$C_factor )

print( c(res$p_value_A, res$p_value_B, res$p_value_C) )
print( c(res$p_value_AB, res$p_value_AC, res$p_value_BC) )
print( res$p_value_ABC )

q_crit = qtukey( 1-0.05, res$K, res$I*res$J*res$K*(res$L-1) )
ss = sqrt( res$MSE/(res$I*res$J*res$L) )
ws = q_crit * ss

diff( sort( res$xbar_dot_dot_k_dot ) )



# Ex 11.30:
#
I = 4
J = 2
K = 2
L = 1
SST = 0.2384
SSA = 0.22625
SSB = 0.000025
SSC = 0.0036
SSAB = 0.004325
SSAC = 0.00065
SSBC = 0.000625
SSE = 0
SSABC = SST - SSA - SSB - SSC - SSAB - SSAC - SSBC - SSE

# Compute the mean squares:
#
MST = SST / ( I * J * K * L - 1 )
MSA = SSA / ( I - 1 )
MSB = SSB / ( J - 1 )
MSC = SSC / ( K - 1 )
MSAB = SSAB / ( (I-1)*(J-1) )
MSAC = SSAC / ( (I-1)*(K-1) )
MSBC = SSBC / ( (J-1)*(K-1) )
MSABC = SSABC / ( (I-1)*(J-1)*(K-1) )

# Compute the f values:
#
f_A = MSA / MSABC
f_B = MSB / MSABC
f_C = MSC / MSABC
f_AB = MSAB / MSABC
f_AC = MSAC / MSABC
f_BC = MSBC / MSABC

# Compute the P-values from these f-values:
#
d_of_freedom_SSABC = (I-1)*(J-1)*(K-1)
p_value_A = 1 - pf( f_A, I-1, d_of_freedom_SSABC )
p_value_B = 1 - pf( f_B, J-1, d_of_freedom_SSABC )
p_value_C = 1 - pf( f_C, K-1, d_of_freedom_SSABC )
p_value_AB = 1 - pf( f_AB, (I-1)*(J-1), d_of_freedom_SSABC )
p_value_AC = 1 - pf( f_AC, (I-1)*(K-1), d_of_freedom_SSABC )
p_value_BC = 1 - pf( f_BC, (J-1)*(K-1), d_of_freedom_SSABC )

print( c( p_value_A, p_value_B, p_value_C ) )
print( c( p_value_AB, p_value_AC, p_value_BC ) )


xbar_i_dot_dot = c( 1.12, 1.3025, 1.3875, 1.43 )
q_crit = qtukey( 1-0.05, I, d_of_freedom_SSABC )
ss = sqrt( MSABC/(J*K) )
ws = q_crit * ss

diff( sort( xbar_i_dot_dot ) )


# Ex 11.32:
#
I = 2
J = 4
K = 3
L = 2

SSA = 14318.24
SSB = 9656.4
SSC = 2270.22
SSAB = 3408.93
SSAC = 1442.72
SSBC = 3096.21
SSABC = 2832.72
SSE = 8655.60
SST = SSA + SSB + SSC + SSAB + SSAC + SSBC + SSABC + SSE

# Compute the mean squares:
#
MST = SST / ( I * J * K * L - 1 )
MSA = SSA / ( I - 1 )
MSB = SSB / ( J - 1 )
MSC = SSC / ( K - 1 )
MSAB = SSAB / ( (I-1)*(J-1) )
MSAC = SSAC / ( (I-1)*(K-1) )
MSBC = SSBC / ( (J-1)*(K-1) )
MSABC = SSABC / ( (I-1)*(J-1)*(K-1) )
MSE = SSE / ( I*J*K*(L-1) )

# Compute the f values:
#
f_A = MSA / MSAC
f_B = MSB / MSBC
f_C = MSC / MSE

f_AB = MSAB / MSABC
f_AC = MSAC / MSE
f_BC = MSBC / MSE

f_ABC = MSABC / MSE

# Compute the P-values from these f-values:
#
p_value_A = 1 - pf( f_A, I-1, (I-1)*(K-1) )
p_value_B = 1 - pf( f_B, J-1, (J-1)*(K-1) )
p_value_C = 1 - pf( f_C, K-1, I*J*K*(L-1) )
p_value_AB = 1 - pf( f_AB, (I-1)*(J-1), (I-1)*(J-1)*(K-1) )
p_value_AC = 1 - pf( f_AC, (I-1)*(K-1), I*J*K*(L-1) )
p_value_BC = 1 - pf( f_BC, (J-1)*(K-1), I*J*K*(L-1) )
p_value_ABC = 1 - pf( f_ABC, (J-1)*(K-1), I*J*K*(L-1) )

print( c( p_value_A, p_value_B, p_value_C ) )
print( c( p_value_AB, p_value_AC, p_value_BC ) )
print( p_value_ABC )


# Ex 11.33:
#
N = 7
x_dot_dot_dot = 3815.8
sum_xbar_i_dot_dot_2 = 297216.9
sum_xbar_dot_j_dot_2 = 297200.64
sum_xbar_dot_dot_k_2 = 297155.01
sum_xbar_i_j_k_2 = 297317.65

xbar_dot_dot_dot = x_dot_dot_dot / N^2
SST = sum_xbar_i_j_k_2 - N^2 * xbar_dot_dot_dot^2
SSA = sum_xbar_i_dot_dot_2 - x_dot_dot_dot^2 / N^2
SSB = sum_xbar_dot_j_dot_2 - x_dot_dot_dot^2 / N^2
SSC = sum_xbar_dot_dot_k_2 - x_dot_dot_dot^2 / N^2
SSE = SST - SSA - SSB - SSC 

# Compute the mean squares:
#
MST = SST / ( N^2 - 1 ) 
MSA = SSA / ( N - 1 )
MSB = SSB / ( N - 1 )
MSC = SSC / ( N - 1 )
MSE = SSE / ( (N-1)*(N-2) )

# Compute the f values:
#
f_A = MSA / MSE
f_B = MSB / MSE
f_C = MSC / MSE

# Compute the P-values from these f-values:
#
p_value_A = 1 - pf( f_A, N-1, (N-1)*(N-2) )
p_value_B = 1 - pf( f_B, N-1, (N-1)*(N-2) )
p_value_C = 1 - pf( f_C, N-1, (N-1)*(N-2) )

print( c( p_value_A, p_value_B, p_value_C ) )


# Ex 11.34:
#
N = 6
data = c( 27, 14, 18, 34, 31, 34, 39, 67, 31, 40, 57, 39, 15, 15, 
          11, 16, 15, 14, 35, 28, 22, 46, 37, 23, 49, 38, 48, 70,
          37, 50, 9, 18, 17, 12, 19, 22)

A_factor = rep( 1:N, each=3, times=2 )
B_factor = c( rep( 1:3, each=1, times=6 ), rep( 4:6, each=1, times=6 ) )
C_factor = c( 5, 4, 3, 6, 5, 4, 2, 6, 5, 3, 1, 2, 4, 3, 1, 1, 2, 6,
              1, 6, 2, 3, 2, 1, 4, 1, 3, 6, 4, 5, 2, 5, 6, 5, 3, 4 ) 
DF = data.frame( sales=data, A_factor=A_factor, B_factor=B_factor, C_factor=C_factor )

res = LATIN_SQUARE( DF$sales, DF$A_factor, DF$B_factor, DF$C_factor ) 

print( c( res$p_value_A, res$p_value_B, res$p_value_C ) ) 


# Ex 11.35:
#
N = 5
data = c( 6.67, 7.15, 8.29, 5.40, 4.77, 5.40, 7.32, 8.53, 8.50, 
          4.92, 5.00, 7.29, 4.88, 6.16, 7.83, 8.95, 9.62, 7.54, 
          6.93, 9.99, 9.68, 7.85, 7.08, 5.83, 8.51 )

A_factor = c( rep( 1:5, each=3, times=1 ), rep( 1:5, each=2, times=1 ) )
B_factor = c( rep( 1:3, each=1, times=5 ), rep( 4:5, each=1, times=5 ) )
C_factor = c( 5, 4, 1, 2, 5, 4, 3, 2, 5, 1, 3, 2, 4, 1, 3, 3, 2, 1, 3, 4, 1, 5, 4, 2, 5 ) 

DF = data.frame( moisture=data, A_factor=A_factor, B_factor=B_factor, C_factor=C_factor )

res = LATIN_SQUARE( DF$moisture, DF$A_factor, DF$B_factor, DF$C_factor ) 

print( c( res$p_value_A, res$p_value_B, res$p_value_C ) ) 


# Ex 11.36:
#
I = 4
J = 3
K = 6
L = 3
SSA = 39.171
SSB = 0.665
SSC = 21.508
SSAB = 1.432
SSAC = 15.953
SSBC = 1.382
SSABC = 9.016
SSE = 115.82
SST = SSA + SSB + SSC + SSAB + SSAC + SSBC + SSABC # compute SST from the others

# Compute the mean squares:
#
MST = SST / ( I * J * K * L - 1 )
MSA = SSA / ( I - 1 )
MSB = SSB / ( J - 1 )
MSC = SSC / ( K - 1 )
MSAB = SSAB / ( (I-1)*(J-1) )
MSAC = SSAC / ( (I-1)*(K-1) )
MSBC = SSBC / ( (J-1)*(K-1) )
MSABC = SSABC / ( (I-1)*(J-1)*(K-1) )
MSE = SSE / ( I*J*K*(L-1) )

# Compute the f values:
#
f_A = MSA / MSE
f_B = MSB / MSE
f_C = MSC / MSE
f_AB = MSAB / MSE
f_AC = MSAC / MSE
f_BC = MSBC / MSE
f_ABC = MSABC / MSE

# Compute the P-values from these f-values:
#
p_value_A = 1 - pf( f_A, I-1, I*J*K*(L-1) )
p_value_B = 1 - pf( f_B, J-1, I*J*K*(L-1) )
p_value_C = 1 - pf( f_C, K-1, I*J*K*(L-1) )
p_value_AB = 1 - pf( f_AB, (I-1)*(J-1), I*J*K*(L-1) )
p_value_AC = 1 - pf( f_AC, (I-1)*(K-1), I*J*K*(L-1) )
p_value_BC = 1 - pf( f_BC, (J-1)*(K-1), I*J*K*(L-1) )
p_value_ABC = 1 - pf( f_ABC, (I-1)*(J-1)*(K-1), I*J*K*(L-1) )

print( c( p_value_A, p_value_B, p_value_C ) )
print( c( p_value_AB, p_value_AC, p_value_BC ) )
print( p_value_ABC )


# Ex 11.37:
# 
# WARNING: for some reason my error degree of freedom calcuation does not seem to be correct.
#
I = 2
J = 2
K = 3
L = 2
MSA = 2207.329
MSB = 47.255
MSC = 491.783
MSD = 0.044
MSAB = 15.303
MSAC = 275.446
MSAD = 0.470
MSBC = 2.141
MSBD = 0.273
MSCD = 0.247
MSABC = 3.714
MSABD = 4.072
MSACD = 0.767
MSBCD = 0.280
MSE = 0.977
MST = 93.621

# Compute the f values:
#
f_A = MSA / MSE
f_B = MSB / MSE
f_C = MSC / MSE
f_D = MSD / MSE
f_AB = MSAB / MSE
f_AC = MSAC / MSE
f_AD = MSAD / MSE
f_BC = MSBC / MSE
f_CD = MSCD / MSE
f_ABC = MSABC / MSE
f_ABD = MSABD / MSE
f_ACD = MSACD / MSE
f_BCD = MSBCD / MSE

# Compute the P-values from these f-values:
#
# IJKL-1 = all DOF             first order DOF                     second order DOF                                                          third order DOF
error_df = ( I*J*K*L - 1 ) - ( (I-1) + (J-1) + (K-1) + (L-1) ) - ( (I-1)*(J-1) + (I-1)*(K-1) + (I-1)*(L-1) + (J-1)*(K-1) + (K-1)*(L-1) ) - ( (I-1)*(J-1)*(K-1) + (I-1)*(J-1)*(L-1) + (I-1)*(K-1)*(L-1) + (J-1)*(K-1)*(L-1) )
p_value_A = 1 - pf( f_A, I-1, error_df )
p_value_B = 1 - pf( f_B, J-1, error_df )
p_value_C = 1 - pf( f_C, K-1, error_df )
p_value_D = 1 - pf( f_D, L-1, error_df )
p_value_AB = 1 - pf( f_AB, (I-1)*(J-1), error_df )
p_value_AC = 1 - pf( f_AC, (I-1)*(K-1), error_df )
p_value_AD = 1 - pf( f_AC, (I-1)*(L-1), error_df )
p_value_BC = 1 - pf( f_BC, (J-1)*(K-1), error_df )
p_value_CD = 1 - pf( f_CD, (K-1)*(L-1), error_df )
p_value_ABC = 1 - pf( f_ABC, (I-1)*(J-1)*(K-1), error_df )
p_value_ABD = 1 - pf( f_ABD, (I-1)*(J-1)*(L-1), error_df )
p_value_ACD = 1 - pf( f_ACD, (I-1)*(K-1)*(L-1), error_df )
p_value_BCD = 1 - pf( f_BCD, (J-1)*(K-1)*(L-1), error_df )