#
# 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(BGSall)
attach(BGSall)

indsM <- Sex==0
indsF <- Sex==1

postscript("../../WriteUp/Graphics/Chapter6/prob_6_sp.eps", onefile=FALSE, horizontal=FALSE)
plot( HT9[indsM], HT18[indsM], pch=c("m"), col=c("red"),
     xlim=range(HT9[indsM],HT9[indsF]),
     ylim=range(HT18[indsM],HT18[indsF]),
     xlab="HT9", ylab="HT18")
par(new=T)
plot( HT9[indsF], HT18[indsF], pch=c("f"), col=c("black"),
     xlim=range(HT9[indsM],HT9[indsF]),
     ylim=range(HT18[indsM],HT18[indsF]),
     xlab="HT9", ylab="HT18" )
par(new=F)
dev.off()

# create a factor based on sex:
BGSall$DF <- factor( BGSall$Sex, ordered=FALSE )

# set this code up to do the model specifcation in general:
#
Xvar    = BGSall$HT9
Fvar    = BGSall$DF 
Yvar    = BGSall$HT18

#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 )


# Make a figure with lattice fitting all four models.
#
library(lattice) # this is needed for the following lattice plot ... 

NC <- length(unique(Fvar)) # the number of classes
n  <- length(Fvar) # the total number of observations 

# create a data frame with the needed information for the plot:
# 
tmpDataFrame      <- BGSall
tmpDataFrame$D    <- Fvar
tmpDataFrame$Xvar <- Xvar
tmpDataFrame$Yvar <- Yvar 

# compute various linear models depending:
# 
n1 <- lm(Yvar~ -1+D+D:Xvar, tmpDataFrame,na.action=na.omit) # model #1: this is the most general model
n2 <- lm(Yvar~ -1+D+Xvar,   tmpDataFrame,na.action=na.omit) # model #2: 
n3 <- lm(Yvar~ +1+D:Xvar,   tmpDataFrame,na.action=na.omit) # model #3
n4 <- lm(Yvar~ Xvar,        tmpDataFrame,na.action=na.omit) # model #4

postscript("../../WriteUp/Graphics/Chapter6/prob_6_model_comparison.eps", onefile=FALSE, horizontal=FALSE)

if(is.null(version$language) == FALSE) {
  
  d      <- rbind(tmpDataFrame,tmpDataFrame,tmpDataFrame,tmpDataFrame)
  labels <- c("(a) General", "(b) Parallel","(c) Common intercept", "(d) Common regression")
  d$f    <- rep(c(1,2,3,4), rep(n,4))
  d$fig  <- rep(labels, rep(n,4))
  xyplot(Yvar~Xvar|fig,group=D, d, subscripts=TRUE,ID=d$f,
         as.table=TRUE, key=simpleKey(as.character(1:NC),FALSE,FALSE,T),
         xlab=expression(Xvar), cex=.75, between=list(x=.32,y=.32),
         panel = function(x,y,groups,subscripts,ID){
           panel.superpose(x,y,groups=groups,subscripts=subscripts,
                           pch=as.character(1:NC))
           g <- unique(ID[subscripts]) # which panel
           if(g==1) {panel.abline(coef(n1)[c(1,3)],lty=1)
                     panel.abline(coef(n1)[c(2,4)],lty=2)}
           if(g==3) {panel.abline(coef(n3)[c(1,2)],lty=1)
                     panel.abline(coef(n3)[c(1,3)],lty=2)}
           if(g==2) {panel.abline(coef(n2)[c(1,3)],lty=1)
                     panel.abline(coef(n2)[c(2,3)],lty=2)}
           if(g==4) {panel.abline(coef(n4)[c(1,2)],lty=1)}
         })
  
}

dev.off()