#
# Epage 
#
# 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.
# 
#-----

if(is.null(version$language) == FALSE){ 
  require(alr3)
  require(car)
}else{
  library(alr3)
  library(car)
}

data(cloud)
attach(cloud)

postscript("../../WriteUp/Graphics/Chapter9/prob_11_scatter_plot_matrix.eps", onefile=FALSE, horizontal=FALSE)
pairs( ~S+E+C+P+Rain )
dev.off()

cor(cloud) # as a result of this matrix we expect that S, E, C, P are most correlated with Rain 

aF <- factor( A, ordered=FALSE ) # a factor for whether we seed or not 

m0 <- lm( Rain ~ -1 + aF + aF:S + aF:E + aF:C + aF:P )
summary(m0) 

#postscript("../../WriteUp/Graphics/Chapter9/prob_11_avp.eps", onefile=FALSE, horizontal=FALSE)
avp( m0, ask=FALSE, one.page=TRUE, identify.points=FALSE)
#dev.off()

m1 <- lm( Rain ~ -1 + aF + aF:S + E + aF:C + P ) # drop the factor splitting of P and E 
summary(m1)
anova(m1,m0) # the added complexity of factors for E and P is not needed

avp( m1, ask=FALSE, one.page=TRUE, identify.points=FALSE)

m2 <- lm( Rain ~ -1 + aF + S + E + aF:C + P ) # attempt to drop the factor splitting of S 
summary(m2)
anova(m2,m1) # we cant drop this factor splitting seem significant at the 5% level 

m3 <- lm( Rain ~ -1 + aF + aF:S + E + C + P ) # attempt to drop the factor splitting of C 
summary(m3)
anova(m3,m1) # we cant drop this factor splitting seem significant at the < 0.01% level 

m4 <- lm( Rain ~ -1 + aF + aF:S + E + aF:C ) # attempt to drop P all together 
summary(m4)
anova(m1,m4) # we can drop this term and model 4 is our current "best" 

m5 <- lm( Rain ~ -1 + aF + aF:S + aF:C ) # attempt to drop E from model the expression E 
summary(m5)
anova(m4,m5) # 10% chance reduction is due to chance ... lets assume we can drop this term ... model 5 is our current best 

postscript("../../WriteUp/Graphics/Chapter9/prob_11_m5_avp.eps", onefile=FALSE, horizontal=FALSE)
avp( m5, ask=FALSE, one.page=TRUE, identify.points=FALSE)
dev.off()