#
# 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.1
#
I = 5
J = 4
SSA = 30.6
SSB = 44.1
SSE = 59.2

# Compute the mean squares:
MSA = SSA / (I-1)
MSB = SSB / (J-1)
MSE = SSE / ( (I-1)*(J-1) )

# Compute the f-values
#
f_A = MSA / MSE
f_A_crit = qf( 1 - 0.05, I-1, (I-1)*(J-1) )
1 - pf( f_A, I-1, (I-1)*(J-1) )

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


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

res = TWO_FACTOR_ANOVA( DF$Response, DF$Coating, DF$Soil.Type )

print( c(res$f_A, res$f_B) )
print( c(res$p_value_A, res$p_value_B) )


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

res = TWO_FACTOR_ANOVA( DF$Response, DF$Gas.Rate..A., DF$Liquid.Rate..B. )

mc_res_1eve1_A = mean_diff_CI_utils( res$I, res$J, res$MSE, res$xbar_i_dot, q_alpha=0.01 )

mc_res_leve1_B = mean_diff_CI_utils( res$J, res$I, res$MSE, res$xbar_dot_j, q_alpha=0.01 )



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

res = TWO_FACTOR_ANOVA( DF$Response, DF$Paint.Brand, DF$Roller.Brand )

mc_res_1eve1_A = mean_diff_CI_utils( res$I, res$J, res$MSE, res$xbar_i_dot, q_alpha=0.05 )


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

res = TWO_FACTOR_ANOVA( DF$Force, DF$Angle..A., DF$Connector..B. )


# Ex 11.6
#
I = 3
J = 5
SSA = 11.7
SSB = 113.5
SSE = 25.6

# Compute the mean squares:
MSA = SSA / (I-1)
MSB = SSB / (J-1)
MSE = SSE / ( (I-1)*(J-1) )

# Compute the f-values
#
f_A = MSA / MSE
f_A_crit = qf( 1 - 0.05, I-1, (I-1)*(J-1) )
1 - pf( f_A, I-1, (I-1)*(J-1) )

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


# Ex 11.7
#
I = 3
J = 18

SST = 3476.00
SSTr = 28.78
SSBl = 2977.67
SSE = SST - SSTr - SSBl

# Compute the mean squares:
MSTr = SSTr / (I-1)
MSBl = SSBl / (J-1)
MSE = SSE / ( (I-1)*(J-1) )

# Compute the f-values
#
f_Tr = MSTr / MSE
f_Tr_crit = qf( 1 - 0.05, I-1, (I-1)*(J-1) )
1 - pf( f_Tr, I-1, (I-1)*(J-1) )

f_B1 = MSBl / MSE
f_B1_crit = qf( 1 - 0.05, J-1, (I-1)*(J-1) )
1 - pf( f_B1, J-1, (I-1)*(J-1) )


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

res = TWO_FACTOR_ANOVA( DF$loss..kg., DF$Subject, DF$temperatures )

mc_res_level_B = mean_diff_CI_utils( res$J, res$I, res$MSE, res$xbar_dot_j, q_alpha=0.01 )

# A qq-plot of the residuals:
#
qqnorm( res$residuals )
qqline( res$residuals )

# residuals vs. the fitted values:
plot( res$fitted_values, res$residuals, type='p', col='black', xlab='fitted value', ylab='residual' )
abline(h=0, col='green')
grid()


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

res = TWO_FACTOR_ANOVA( DF$Response, DF$Type, DF$Subject )


# Ex 11.10 (this looks like the right data but the headers in this dataframe look strange:
# (they look like they belong to the previous problem):
#
DF = read.csv( "../../Data/CH11/ex11-10.txt", header=TRUE, quote="'" )
res = TWO_FACTOR_ANOVA( DF$Effor, DF$Subject, DF$Stool.Type )


# Ex 11.11
#
DF = read.csv( "../../Data/CH11/exp11-01.txt", header=TRUE, quote="'" )
res = TWO_FACTOR_ANOVA( DF$Response, DF$Brand, DF$Treatment )

# A qq-plot of the residuals:
qqnorm( res$residuals )
qqline( res$residuals )

# residuals vs. the fitted values:
plot( res$fitted_values, res$residuals, type='p', col='black', xlab='fitted value', ylab='residual' )
abline(h=0, col='green')
grid()


# Ex 11.15
#
I = 4
J = 3
s = 4
phi2 = (J/I) * ( 1/s^2 ) * ( 4^2 + 2*4 )
phi = sqrt(phi2) #=1.06066

phi2 = (J/I) * ( 1/s^2 ) * ( 6^2 + 2*9 )
phi = sqrt(phi2) #=1.59099

I = 4
J = 3
s = .3
phi2 = (I/J) * ( 1/s^2 ) * ( 0.3^2 + 3*(0.1)^2 )
phi = sqrt(phi2) #=1