#
# 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.50 (EPage 442):
#
DF = read.csv( "../../Data/CH11/ex11-50.txt", header=TRUE, quote="'" )
res = TWO_FACTOR_ANOVA( DF$Response, DF$Fabric, DF$Drying )


# Ex 11.51:
#
I = 2
J = 4
K = 3
SST = 372.113
SSA = 322.667
SSB = 35.623
SSAB = 8.557
SSE = SST - SSA - SSB - SSAB

MST = SST / (I*J*K-1)
MSA = SSA / (I-1)
MSB = SSB / (J-1)
MSE = SSE / ( I*J*(K-1) )
MSAB = SSAB / ( (I-1)*(J-1) )

f_A = MSA/MSE
f_A_crit = qf( 1 - 0.05, I-1, I*J*(K-1) )
p_value_A = 1 - pf( f_A, I-1, I*J*(K-1) )

f_B = MSB / MSE
f_B_crit = qf( 1 - 0.05, J-1, I*J*(K-1) )
p_value_B = 1 - pf( f_B, J-1, I*J*(K-1) )

f_AB = MSAB / MSE 
f_AB_crit = qf( 1 - 0.05, (I-1)*(J-1), I*J*(K-1) )
p_value_AB = 1 - pf( f_AB, (I-1)*(J-1), I*J*(K-1) )

print( c( p_value_A, p_value_B, p_value_AB ) ) 


# Ex 11.52:
#
DF = read.csv( "../../Data/CH11/ex11-52.txt", header=TRUE, quote="'" )
res = TWO_FACTOR_ANOVA( DF$Response, DF$Plot, DF$Sowing.Rate )


# EX 11.53: 
#
DF = read.csv( "../../Data/CH11/ex11-53.txt", header=TRUE, quote="'" )
D = DF[,c(5,6)]
res = YATES_METHOD( D )


# Ex 11.54: 
#
DF = read.csv( "../../Data/CH11/ex11-54.txt", header=TRUE, quote="'" )


# Ex 11.55: 
#
DF = read.csv( "../../Data/CH11/ex11-54.txt", header=TRUE, quote="'" )
#res = YATES_METHOD( DF$Response )


# Ex 11.56: 
#
DF = read.csv( "../../Data/CH11/ex11-56.txt", header=TRUE, quote="'" )
res = TWO_FACTOR_KIJ_ANOVA( DF$Response, DF$Health, DF$pH )


# Ex 11.57: 
#
DF = read.csv( "../../Data/CH11/ex11-57.txt", header=TRUE, quote="'" )
res = THREE_FACTOR_ANOVA( DF$Response, DF$Pressure, DF$Temperature..Degrees., DF$Fabric )

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 )


# Ex 11.58:
#
I = 2
J = 4
K = 3
L = 2
SST = 70.82
SSA = 6.94
SSB = 5.61
SSC = 12.33
SSAB = 4.05
SSAC = 7.32
SSBC = 15.8
SSE = 14.40 
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 stats:
#
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.59:
#
DF = read.csv( "../../Data/CH11/ex11-59.txt", header=TRUE, quote="'" )
res = THREE_FACTOR_ANOVA( DF$Response, DF$Conductor.Material, DF$Adhesive, DF$Cure.Time )

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 )


# Ex 11.60:
#
# I'm not sure this data matches the summary stats given for this exercise:
# DF = read.csv( "../../Data/CH11/ex11-60.txt", header=TRUE, quote="'" )
#
I = 3
J = 3
SSA = 18138
SSB = 5182
SSAB = 1737
SST = 36348
SSE = SST - SSA - SSB - SSAB
error_df = 36 # this equals I * J * (K-1) so
K = error_df/(I*J) + 1

MST = SST / (I*J*K-1)
MSA = SSA / (I-1)
MSB = SSB / (J-1)
MSAB = SSAB / ( (I-1)*(J-1) )
MSE = SSE / ( I*J*(K-1) )

f_A = MSA / MSE
f_B = MSB / MSE
f_AB = MSAB / MSE

p_value_A = 1 - pf( f_A, I-1, I*J*(K-1) )
p_value_B = 1 - pf( f_B, J-1, I*J*(K-1) )
p_value_AB = 1 - pf( f_AB, (I-1)*(J-1), I*J*(K-1) )

print( c(f_A, f_B, f_AB) )
print( c(p_value_A, p_value_B, p_value_AB) )