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

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

data(landrent)
attach(landrent)

postscript("../../WriteUp/Graphics/Chapter9/prob_10_scatter_plot_matrix.eps", onefile=FALSE, horizontal=FALSE)
pairs( ~X1+X2+X3+X4+Y )
dev.off()

cor(landrent)

# result of this matrix:
#
# X1 most correlated, then X3, then X2 (which is not very correlated with X1)

mAll <- lm( Y ~ X1+X2+X3+X4 )

library(car)
postscript("../../WriteUp/Graphics/Chapter9/prob_10_avp.eps", onefile=FALSE, horizontal=FALSE)
avp( mAll, ask=FALSE, one.page=TRUE, identify.points=FALSE)
dev.off()


## # to see if X4 should be a factor lets look at scatterplot matrices of the two possible factors
## #
## X1s <- X1[X4==0]
## X2s <- X2[X4==0]
## X3s <- X3[X4==0]
## Ys  <- Y[X4==0] 
## pairs( ~X1s+X2s+X3s+Ys )
## scatterplot.matrix( ~X1s+X2s+X3s+Ys )

## X1s <- X1[X4==1]
## X2s <- X2[X4==1]
## X3s <- X3[X4==1]
## Ys  <- Y[X4==1] 
## pairs( ~X1s+X2s+X3s+Ys )
## scatterplot.matrix( ~X1s+X2s+X3s+Ys )

# remove X4:
# 
m0 <- lm( Y ~ X1 + X2 + X3 )
summary(m0) 

# compute the most general model: 
x4Factor <- factor( X4, ordered=FALSE )
m1       <- lm( Y ~ X1 + X2 + X3 + x4Factor + X1:x4Factor + X2:x4Factor + X3:x4Factor ) 
summary(m1) # none of the parameters are very well estimated ...

# lets compare the two models ... 
anova(m1,m0) 

# a model with only a different "mean" beta_{0j} 
m2  <- lm( Y ~ X1 + X2 + X3 + x4Factor ) 
summary(m2)

# is it worth including a different mean ... looks like NO
anova(m2,m0) 


# we drop X4 to consider 
m0 <- lm( Y ~ X1 + X2 + X3 )

postscript("../../WriteUp/Graphics/Chapter9/prob_10_avp_no_X4.eps", onefile=FALSE, horizontal=FALSE)
avp( m0, ask=FALSE, one.page=TRUE, identify.points=FALSE)
dev.off()

mprime <- lm( Y ~ X1 + X2 )
summary(mprime)

anova(mprime,m0) # we conclude that X3 does NOT contribute significantly to the reduction in RSS ... we will use this model.

ri <- rstandard( m0 )
ti <- ri * sqrt( ( m0$df - 1 ) / (m0$df-ri^2) ) 
spot <- match( max(abs(ti)), abs(ti) )
landrent[spot,]

n <- length(X1)
pprime <- length(coefficients(m0))
n * 2 * ( 1 - pt(ti[spot], n-pprime-1) )

#
postscript("../../WriteUp/Graphics/Chapter9/prob_10_avp_no_X3_X4.eps", onefile=FALSE, horizontal=FALSE)
avp( mprime, ask=FALSE, one.page=TRUE, identify.points=FALSE)
dev.off()


# lets look for transformations of the independent variables 
bct <- bctrans( Y ~ X1 + X2 )
summary(bct)   # <- seems to indicate no transformation needed

# look for a transformation of the dependent variable ... not done.