#
# Epage 191
#
# 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(pipeline)
attach(pipeline)

# 8.3.1:
# 
postscript("../../WriteUp/Graphics/Chapter8/prob_3_orig_scatter_plot.eps", onefile=FALSE, horizontal=FALSE)

plot( Field, Lab )

dev.off() 

# 8.3.2: the score test for non constant variance i.e. does $\sigma^2$ depend on Field:
#
postscript("../../WriteUp/Graphics/Chapter8/prob_3_orig_scatter_plot.eps", onefile=FALSE, horizontal=FALSE)
plot( Field, Lab )
m0 <- lm( Lab ~ Field ) 
abline(m0)
dev.off()

# this is one of the residual plots
postscript("../../WriteUp/Graphics/Chapter8/prob_3_residual_plot.eps", onefile=FALSE, horizontal=FALSE)
plot( predict(m0), residuals(m0,type="pearson") )
dev.off()
#postscript("../../WriteUp/Graphics/Chapter8/prob_3_residual_plots.eps", onefile=FALSE, horizontal=FALSE)
#residual.plots( m0 )
#dev.off()

# we should plot the residuals vs. the independent variable Field to suggest what variance model to take

if(is.null(version$language) == FALSE) library(xtable)

# perform a score test for constant variance:
# 
round(coef(m0),2)
sig2 <- sum(residuals(m0,type="pearson")^2)/length(pipeline$Lab)
U <- residuals(m0,type="pearson")^2/sig2
m2 <- lm(U~pipeline$Field)
anova(m2)
xtable(anova(m2))

S = 0.5 * 59.17 
1 - pchisq( S, 1 )

# use the car package (we don't have a constant variance \sigma) ... gives the same result as above:
library(car)
ncv.test(m) 

# 8.3.3:
#
m1 <- lm( Lab ~ Field, weights=Field ) 

# make a residual plot (which still looks like it has a great deal of variance): 
plot( predict(m1), residuals(m1,type="pearson") ) 

# compute the new score test:
# 
round(coef(m1),2)
sig2 <- sum(residuals(m1,type="pearson")^2)/length(pipeline$Lab)
U <- residuals(m1,type="pearson")^2/sig2
m2 <- lm(U~ I(1/Field))
anova(m2)
xtable(anova(m2))

S = 0.5 * 49.97
1 - pchisq( S, 1 )

ncv.test(m1, ~ I(1/Field))


# factor addition skipped ... 
#bFactor <- factor( Batch, ordered=FALSE ) 



# 8.3.4:
#
m2 <- lm( Lab ~ Field, weights=Field^2 ) 

# make a residual plot (which still looks like it has a great deal of variance): 
plot( predict(m2), residuals(m2,type="pearson") ) 

# compute the new score test:
# 
round(coef(m2),2)
sig2 <- sum(residuals(m2,type="pearson")^2)/length(pipeline$Lab)
U <- residuals(m2,type="pearson")^2/sig2
m3 <- lm(U~I(1/Field^2))
anova(m3)
xtable(anova(m3))

S = 0.5 * 22.23
1 - pchisq( S, 1 )

ncv.test(m2, ~ I(1/Field^2))