source('utils.R')

# Section 4; Question 1:
#
source('chap_12_sect_4_question_1_data.R')
res = wANOVA( DF, 'type', 'time' )

# Due to the order of the factors output above we need to permute them to get them back into the original order:
#
Y_bar_dot_j = res$Y_bar_dot_j[ res$increasing_order_to_original_order ] # get the means in the original factor order
n_j = res$n_j[ res$increasing_order_to_original_order ]

# Compute the sum of squares for the two contrasts:
#
contrast_coeffs_1 = c( 1, 0, -1 )
c1 = sum( contrast_coeffs_1 * Y_bar_dot_j ) # the first contrast
SSc1 = c1^2 / sum( contrast_coeffs_1^2 / n_j ) # sum of squares

contrast_coeffs_2 = c( 1/2, -1, 1/2 )
c2 = sum( contrast_coeffs_2 * Y_bar_dot_j ) # the second contrast
SSc2 = c2^2 / sum( contrast_coeffs_2^2 / n_j ) # sum of squares
print( sprintf( 'SSTR= %10.6f; SSc1+SSc2= %10.6f', res$SSTR, SSc1+SSc2 ) )

# Test each contrast for significance:
#
f_1 = ( SSc1 / 1 ) / ( res$SSE / (res$n - res$k) )
f_2 = ( SSc2 / 1 ) / ( res$SSE / (res$n - res$k) )

# The P-value for each contrast:
#
p_1 = 1 - pf( f_1, res$k - 1, res$n - res$k )
p_2 = 1 - pf( f_2, res$k - 1, res$n - res$k )

# Print the ANOVA table for these two contrasts:
#
print( sprintf( "Contrast 1: SS= %10.6f f= %10.6f P-value= %0.6f", SSc1, f_1, p_1 ) )
print( sprintf( "Contrast 2: SS= %10.6f f= %10.6f P-value= %0.6f", SSc2, f_2, p_2 ) )


# Section 4; Question 2:
#
source('chap_12_sect_2_question_4_data.R')
res = wANOVA( DF, 'sector', 'yield' )

# Due to the order of the factors output above we need to permute them to get them back into the original order:
#
Y_bar_dot_j = res$Y_bar_dot_j[ res$increasing_order_to_original_order ] # get the means in the original factor order
n_j = res$n_j[ res$increasing_order_to_original_order ]

contrast_coeffs = c( 1/3, 1/3, 1/3, -1/2, -1/2 )
c = sum( contrast_coeffs * Y_bar_dot_j ) # the contrast
SSc = c^2 / sum( contrast_coeffs^2 / n_j )

f = ( SSc / 1 ) / ( res$SSE / (res$n - res$k) )
p = 1 - pf( f, res$k - 1, res$n - res$k )
print( sprintf( "SS= %10.6f f= %10.6f P-value= %0.6f", SSc, f, p ) )


# Section 4; Question 3:
#
source('case_study_12_2_1_data.R')
res = wANOVA( DF, 'smoker_type', 'value' )

# Due to the order of the factors output above we need to permute them to get them back into the original order:
#
Y_bar_dot_j = res$Y_bar_dot_j[ res$increasing_order_to_original_order ] # get the means in the original factor order
n_j = res$n_j[ res$increasing_order_to_original_order ]

contrast_coeffs = c( 0, 1/2, 1/2, -1 )
c = sum( contrast_coeffs * Y_bar_dot_j ) # the contrast
SSc = c^2 / sum( contrast_coeffs^2 / n_j )

f = ( SSc / 1 ) / ( res$SSE / (res$n - res$k) )
p = 1 - pf( f, res$k - 1, res$n - res$k )
print( sprintf( "SS= %10.6f f= %10.6f P-value= %0.6f", SSc, f, p ) )


# Section 4; Question 4:
#
source('chap_12_sect_4_question_4_data.R')
res = wANOVA( DF, 'size', 'asserts' )

# Due to the order of the factors output above we need to permute them to get them back into the original order:
#
Y_bar_dot_j = res$Y_bar_dot_j[ res$increasing_order_to_original_order ] # get the means in the original factor order
n_j = res$n_j[ res$increasing_order_to_original_order ]

contrast_coeffs = c( 1/2, 1/2, -1 )
c = sum( contrast_coeffs * Y_bar_dot_j ) # the contrast
SSc = c^2 / sum( contrast_coeffs^2 / n_j )
f = ( SSc / 1 ) / ( res$SSE / (res$n - res$k) )

p = 1 - pf( f, res$k - 1, res$n - res$k )
print( sprintf( "SS= %10.6f f= %10.6f P-value= %0.6f", SSc, f, p ) )

# Section 4; Question 5:
#
source('case_study_12_4_1_data.R')
res = wANOVA( DF, 'group', 'age' )

# Due to the order of the factors output above we need to permute them to get them back into the original order:
#
Y_bar_dot_j = res$Y_bar_dot_j[ res$increasing_order_to_original_order ] # get the means in the original factor order
n_j = res$n_j[ res$increasing_order_to_original_order ]

# What are the contrasts:
#
contrast_1 = c( 1, -1, 0, 0 )
contrast_2 = c( 0, 0, 1, -1 )
contrast_3 = c( 11/12, 11/12, -1, -5/6 )

# Verify that the contrasts are orthogonal:
#
sum( ( contrast_1 * contrast_2 ) / n_j )
sum( ( contrast_1 * contrast_3 ) / n_j )
sum( ( contrast_2 * contrast_3 ) / n_j )

c_1 = sum( contrast_1 * Y_bar_dot_j )
SSc_1 = c_1^2 / sum( contrast_1^2 / n_j )
c_2 = sum( contrast_2 * Y_bar_dot_j )
SSc_2 = c_2^2 / sum( contrast_2^2 / n_j )
c_3 = sum( contrast_3 * Y_bar_dot_j )
SSc_3 = c_3^2 / sum( contrast_3^2 / n_j )

f = ( SSc_1 / 1 ) / ( res$SSE / (res$n - res$k) )
p = 1 - pf( f, res$k - 1, res$n - res$k )
print( sprintf( "Contrast 1: c= %.3f SS= %10.6f P-value= %0.6f", c_1, SSc_1, p ) )

f = ( SSc_2 / 1 ) / ( res$SSE / (res$n - res$k) )
p = 1 - pf( f, res$k - 1, res$n - res$k )
print( sprintf( "Contrast 2: c= %.3f SS= %10.6f P-value= %0.6f", c_2, SSc_2, p ) )

f = ( SSc_3 / 1 ) / ( res$SSE / (res$n - res$k) )
p = 1 - pf( f, res$k - 1, res$n - res$k )
print( sprintf( "Contrast 3: c= %.3f SS= %10.6f P-value= %0.6f", c_3, SSc_3, p ) )

print( sprintf("SSTR= %10.6f; sum of all contrast SSS= %10.6f", res$SSTR, SSc_1+SSc_2+SSc_3) );



# Section 4; Question 6:
#
source('chap_12_sect_4_question_6_data.R')
res = wANOVA( DF, 'type', 'amount' )

# Due to the order of the factors output above we need to permute them to get them back into the original order:
#
Y_bar_dot_j = res$Y_bar_dot_j[ res$increasing_order_to_original_order ] # get the means in the original factor order
n_j = res$n_j[ res$increasing_order_to_original_order ]

# The constrast of interest:
#
contrast_1 = c( 1, -1, 0, 0 )
contrast_2 = c( 1/2, 1/2, -1/2, -1/2 )

c_1 = sum( contrast_1 * Y_bar_dot_j )
SSc_1 = c_1^2 / sum( contrast_1^2 / n_j )
c_2 = sum( contrast_2 * Y_bar_dot_j )
SSc_2 = c_2^2 / sum( contrast_2^2 / n_j )

# What must the SS of the third contrast be:
res$SSTR - SSc_1 - SSc_2

# Find the null space of the above linear conditions:
#
library(MASS)
m = matrix( data=c( rep( 1, 4 ), contrast_1/3, contrast_2/3 ), nrow=3, ncol=4, byrow=TRUE )
Null(t(m))

#
contrast_3 = c( 0, 0, -1/2, +1/2 )
c_3 = sum( contrast_3 * Y_bar_dot_j )
SSc_3 = c_3^2 / sum( contrast_3^2 / n_j )

# Lets look at what contrasts are orthogonal:
#
f = ( SSc_1 / 1 ) / ( res$SSE / (res$n - res$k) )
p = 1 - pf( f, res$k - 1, res$n - res$k )
print( sprintf( "Contrast 1: SS= %10.6f P-value= %0.6f", SSc_1, p ) )

f = ( SSc_2 / 1 ) / ( res$SSE / (res$n - res$k) )
p = 1 - pf( f, res$k - 1, res$n - res$k )
print( sprintf( "Contrast 2: SS= %10.6f P-value= %0.6f", SSc_2, p ) )

f = ( SSc_3 / 1 ) / ( res$SSE / (res$n - res$k) )
p = 1 - pf( f, res$k - 1, res$n - res$k )
print( sprintf( "Contrast 3: SS= %10.6f P-value= %0.6f", SSc_3, p ) )