# Section 2; Question 1:
#
source('chap_13_sect_2_question_1_data.R')

fit = aov( n_correct ~ state + student, data=DF_melt )
print( summary( fit ) )

# Section 2; Question 2:
#
source('chap_13_sect_2_question_2_data.R')

fit = aov( share ~ network + city, data=DF_melt )
print( summary( fit ) )

# Section 2; Question 3:
#
source('chap_13_sect_2_question_3_data.R')

fit = aov( sumens ~ additive + batch, data=DF_melt )
print( summary( fit ) )

# Section 2; Question 4:
#
source('chap_13_sect_2_question_4_data.R')
fit = aov( percent_increase ~ region + year, data=DF_melt )
print( summary( fit ) )

# Section 2; Question 5:
#
# Following the notes here: https://www.r~bloggers.com/anova-and-tukeys-test-on-r/
#
source('case_study_13_2_3_data.R')
fit = aov( admission_rate ~ time + month, data=DF_melt )
print( summary( fit ) )
print( TukeyHSD(fit, 'time') )

# Section 2; Question 6:
#
source('chap_13_sect_2_question_6_data.R')
fit = aov( return ~ quarter + year, data=DF_melt )
print( summary( fit ) )

# Section 2; Question 7:
#
source('chap_13_sect_2_question_2_data.R')

fit = aov( share ~ network + city, data=DF_melt )
print( summary( fit ) )
print( TukeyHSD(fit, 'network') )

# Section 2; Question 8:
#
source('chap_13_sect_2_question_8_data.R')

fit = aov( time ~ system + subject, data=DF_melt )
print( summary( fit ) )
print( TukeyHSD(fit, 'system') )

# Section 2; Question 9:
#
source('chap_13_sect_2_question_9_data.R')

fit = aov( heart_rate ~ sleep_stage + shrew, data=DF_melt )
print( summary( fit ) )

# Use some routines written for the previous chapter to extract summary statistics:
#
source('../Chapter12/utils.R')
res = wANOVA( DF_melt, 'sleep_stage', 'heart_rate' )

c1 = c( -1, 1/2, 1/2 ) # ordered to match according to the res$factor field in the res=wANOVA() output
c1_hat = sum( c1 * res$Y_bar_dot_j )
SSc1 = c1_hat^2 / sum( c1^2 / res$n_j )

c2 = c( 0, 1/2, -1/2 ) # an orthogonal contrast
c2_hat = sum( c2 * res$Y_bar_dot_j )
SSc2 = c2_hat^2 / sum( c2^2 / res$n_j )

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= %6.3f f= %6.3f P-value= %6.3f", SSc1, f_1, p_1 ) )
print( sprintf( "Contrast 2: SS= %6.3f f= %6.3f P-value= %6.3f", SSc2, f_2, p_2 ) )


# Section 2; Question 10:
#
source('case_study_13_2_2_data.R')

# Use some routines written for the previous chapter to extract summary statistics:
#
source('../Chapter12/utils.R')
res = wANOVA( DF_melt, 'type', 'acceptance_percentage' )

c1 = c( 0, 1, -1, 0 ) # ordered to match according to the res$factor field in the res=wANOVA() output
c1_hat = sum( c1 * res$Y_bar_dot_j )
SSc1 = c1_hat^2 / sum( c1^2 / res$n_j )

c2 = c( -1, 0, 0, 1 ) # the first orthogonal contrast
c2_hat = sum( c2 * res$Y_bar_dot_j )
SSc2 = c2_hat^2 / sum( c2^2 / res$n_j )

c3 = c( 1/2, -1/2, -1/2, 1/2 ) # the second orthogonal contrast
c3_hat = sum( c3 * res$Y_bar_dot_j )
SSc3 = c3_hat^2 / sum( c3^2 / res$n_j )

f_1 = ( SSc1 / 1 ) / ( res$SSE / (res$n - res$k) )
f_2 = ( SSc2 / 1 ) / ( res$SSE / (res$n - res$k) )
f_3 = ( SSc3 / 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 )
p_3 = 1 - pf( f_3, res$k - 1, res$n - res$k )

# Print the ANOVA table for these three contrasts:
#
print( sprintf( "Contrast 1: SS= %6.3f f= %6.3f P-value= %6.3f", SSc1, f_1, p_1 ) )
print( sprintf( "Contrast 2: SS= %6.3f f= %6.3f P-value= %6.3f", SSc2, f_2, p_2 ) )
print( sprintf( "Contrast 3: SS= %6.3f f= %6.3f P-value= %6.3f", SSc3, f_3, p_3 ) )


# Section 2; Question 12:
#

M = matrix(data=c(
               1, 0, 0, 1/2, 1/2,
               0, 1, 0, 1/2, 1/2,
               0, 0, 1, 1/2, 1/2,
               1/3, 1/3, 1/3, 1, 0,
               1/3, 1/3, 1/3, 0, 1),
           nrow=5, ncol=5, byrow=TRUE )

library(Matrix)
print( rankMatrix(M) )