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

if(is.null(version$language) == FALSE) data(turk0)
attach(turk0)

postscript("../../WriteUp/Graphics/Chapter6/prob_15_sp.eps", onefile=FALSE, horizontal=FALSE)
plot(A,Gain, xlab="Amount (percent of diet)", ylab="Weight gain (g)")
dev.off()

# fit the data without any weights:
#
m1 <- lm( Gain ~ A )

# do a quadratic OLS regression on yBar:
#
m2 <- update( m1, ~.+I(A^2) )

# plot everything on one graph:
#
postscript("../../WriteUp/Graphics/Chapter6/prob_15_sp.eps", onefile=FALSE, horizontal=FALSE)
plot(A,Gain, xlab="Amount (percent of diet)", ylab="Weight gain (g)")
abline(m1)
rA <- range(A)
aseq <- seq( rA[1], rA[2],length=50)
lines(aseq,predict(m2,data.frame(A=aseq)),lty=2)
dev.off()

# get information needed to do lack-of-fit tests:
#
nIn = c()
aIn = c()
yBar = c()
SD = c()
for ( x in unique(A) ){
  inds <- A == x
  nIn  <- c(nIn,sum(inds))
  aIn  <- c(aIn,x)
  yBar <- c(yBar,mean(turk0$Gain[inds]))
  SD   <- c(SD,sqrt(var(turk0$Gain[inds])))
}

SD2 <- SD^2

# compute the mean square for pure error:
#
SSPE <- sum( (nIn-1) * SD2 ) 
dfPE <- sum( (nIn-1) )

s2PE <- SSPE/dfPE

n <- length(A)

# for each model compute a F lack-of-fit statistic:
#
#linear:
pprime <- 1 + 1
RSS    <- sum( residuals(m1)^2 )

F     <- ( (RSS - SSPE)/(n-pprime-dfPE) ) / ( SSPE / dfPE )
alpha <- 1 - pf( F, n-pprime-dfPE, dfPE )
print(alpha)

#quadratic: 
pprime <- 2 + 1
RSS    <- sum( residuals(m2)^2 )

F     <- ( (RSS - SSPE)/(n-pprime-dfPE) ) / ( SSPE / dfPE )
alpha <- 1 - pf( F, n-pprime-dfPE, dfPE )
print(alpha)

# if we just wish to compute a $F$-statistic threshold we can use: 

#qf( 1 - 0.05, n-pprime-dfPE, dfPE )

#qf( 1 - 0.05, 2, 6 )