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

# version$language == "R" for R version$language == NULL for SPlus

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

data(wm2)
attach(wm2)

# create a factor based on the bin:
D <- factor( wm2$Bin, ordered=FALSE )

# set this code up to do the model specifcation in general:
#
Xvar    = wm2$RSpd
Fvar    = D
Yvar    = wm2$CSpd

#postscript("../../WriteUp/Graphics/Chapter6/prob_5_boxplot.eps", onefile=FALSE, horizontal=FALSE)
plot(Fvar,Yvar,xlab="factor index" ,ylab="Response")
#dev.off()

# fit the various four cases of models to this dateset:
#
m1 <- lm(Yvar~ Xvar + Fvar + Fvar:Xvar, na.action=na.omit)  # this is model 1 the most general 
m2 <- lm(Yvar~ Xvar + Fvar            , na.action=na.omit)  # this is model 2 parallel
m3 <- lm(Yvar~ Xvar + Fvar:Xvar       , na.action=na.omit)  # this is model 3 common intercept
m4 <- lm(Yvar~ Xvar                   , na.action=na.omit)  # this is model 4 the least general (all the same) 

# get some statistics on m1 (residual degrees of freedom) and residual sum of squares:
#
dgts <- 4

print( c(m1$df.residual, sum( m1$residuals^2 )), digits=dgts ) 

# compare all models to the most general model (model 1):
#
a2<-anova(m2,m1)
print( c( a2$Res.Df[1], a2$RSS[1], a2$F[2], a2$'Pr(>F)'[2] ), digits=dgts ) 

a3<-anova(m3,m1)
print( c( a3$Res.Df[1], a3$RSS[1], a3$F[2], a3$'Pr(>F)'[2] ), digits=dgts )

a4<-anova(m4,m1)
print( c( a4$Res.Df[1], a4$RSS[1], a4$F[2], a4$'Pr(>F)'[2] ), digits=dgts )


# 6.11.2: 
#
data(wm3)

m  = sum( wm3$bin.count ) 
mI = wm3$bin.count / m