#
# Epage 139
#
# 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(allshoots)
attach(allshoots)

# compute the mean square for pure error for both types of shoots:
#
shortType <- Type == 0
longType  <- Type == 1

SSPEShort <- sum( ( n[shortType] - 1 ) * SD[shortType]^2 )
SSPELong  <- sum( ( n[longType] - 1 ) * SD[longType]^2 )

dfShort <- sum( ( n[shortType] - 1 ) )
dfLong <- sum( ( n[longType] - 1 ) )

sigma2PEShort = SSPEShort/dfShort
sigma2PELong = SSPELong/dfLong

F = sigma2PELong / sigma2PEShort 

# compute the significance of this statistic:
qf( 1-0.05, dfLong, dfShort )

# the probability under the null hypothesis we would get this value or larger:
1 - pf( F, dfLong, dfShort )

# compute the pooled estimate:
#
sigma2Pooled = sigma2PEShort + (1/2) * sigma2PELong 

# 6.9.2:
#
postscript("../../WriteUp/Graphics/Chapter6/prob_9_sp.eps", onefile=FALSE, horizontal=FALSE)
plot( Day[shortType], ybar[shortType], pch=c("s"), col=c("red"),
     ylim=range(ybar[shortType],ybar[longType]),
     xlim=range(Day[shortType],Day[longType]),
     xlab="Day", ylab="ybar")
par(new=T)
plot( Day[longType], ybar[longType], pch=c("l"), col=c("black"),
     ylim=range(ybar[shortType],ybar[longType]),
     xlim=range(Day[shortType],Day[longType]),
     xlab="Day", ylab="ybar" )
par(new=F)
dev.off()


# 6.9.3:
#

# create a factor based on the long/short:
allshoots$D <- factor( allshoots$Type, ordered=FALSE )

# set this code up to do the model specifcation in general:
#
Xvar    = allshoots$Day
Fvar    = allshoots$D
Yvar    = allshoots$ybar

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

# compute weights for each of these sample points (weight = n_i)
# 
# 
theWeights           <- allshoots$n
theWeights[longType] <- theWeights[longType]/2 # normalize the variance of the long shoots
theWeights           <- theWeights / sigma2Pooled 

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

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

print( c(m1$df.residual, sum( residuals(m1,type="pearson")^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      <- allshoots
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, weights=theWeights) # model #1: this is the most general model
n2 <- lm(Yvar~ -1+D+Xvar,   tmpDataFrame,na.action=na.omit, weights=theWeights) # model #2: 
n3 <- lm(Yvar~ +1+D:Xvar,   tmpDataFrame,na.action=na.omit, weights=theWeights) # model #3
n4 <- lm(Yvar~ Xvar,        tmpDataFrame,na.action=na.omit, weights=theWeights) # model #4

postscript("../../WriteUp/Graphics/Chapter6/prob_9_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()