#
# 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(lathe1)
lathe1$LLife <- log(lathe1$Life,2)

attach(lathe1)

jLLife = jitter(LLife)
jSpeed = jitter(Speed)
jFeed = jitter(Feed)
jLife = jitter(Life)

# a scatterplot matrix
#
postscript("../../WriteUp/Graphics/Chapter6/prob_2.eps", onefile=FALSE, horizontal=FALSE)
pairs(jLLife~jSpeed+jFeed+jLife)
dev.off()

# a model with an interaction term Speed:Feed:
#
m1 <- lm(LLife ~ Speed + Feed + I(Speed^2) + I(Feed^2) + Speed:Feed, data=lathe1)
s1 <- summary( m1 )
s1 

m0 <- lm(LLife ~ Speed + Feed + I(Speed^2) + I(Feed^2), data=lathe1)
summary(m0)

anova(m0,m1)

X1 <- Speed
X2 <- Feed
Y  <- LLife

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

# duplicate Fig. 6.3 (a model with an interaction term)  "../EPS-figs/cake2.eps"
#
oldpar<-par(mfrow=c(1,2),mar=c(4,3,0,.5)+.1,mgp=c(2,1,0))
# part (a)
#
if(is.null(version$language) == FALSE)
   plot(X1,Y,type="n",xlab=expression(paste("(a)  ",X[1]))) else
   plot(X1,Y,type="n",xlab="(a) X1")
X1new <- seq(-1.414,1.414,len=50)
lines(X1new,predict(m1,newdata=data.frame(Speed=X1new,Feed=rep(-1.414,50))))
lines(X1new,predict(m1,newdata=data.frame(Speed=X1new,Feed=rep(0.0,50))))
lines(X1new,predict(m1,newdata=data.frame(Speed=X1new,Feed=rep(1.414,50))))
text(-1.414,1.0,"Feed=-1.414",adj=0,cex=0.7)
text(0.0,2.0,"Feed=0.0",adj=0,cex=0.7)
text(0.5,4.0,"Feed=+1.414",adj=0,cex=0.7)
# part (b)
if(is.null(version$language) == FALSE)
   plot(X2,Y,type="n",xlab=expression(paste("(b)  ",X[2]))) else
   plot(X2,Y,type="n",xlab="(b) X2")
X2new <- seq(-1.414,1.414,len=50)
lines(X2new,predict(m1,newdata=data.frame(Speed=rep(-1.414,50),Feed=X2new)))
lines(X2new,predict(m1,newdata=data.frame(Speed=rep(0,50),Feed=X2new)))
lines(X2new,predict(m1,newdata=data.frame(Speed=rep(+1.414,50),Feed=X2new)))
text(-1.0,0,  "Speed=33",adj=0,cex=0.7)
text(0,2.,"Speed=35",adj=0,cex=0.7)
text(.45,5,"Speed=37",adj=0,cex=0.7)
par(oldpar)

dev.off()


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

# duplicate Fig. 6.4 (a model with NO interaction term)  "../EPS-figs/cake2.eps"
#
oldpar<-par(mfrow=c(1,2),mar=c(4,3,0,.5)+.1,mgp=c(2,1,0))
# part (a)

if(is.null(version$language) == FALSE)
   plot(X1,Y,type="n",xlab=expression(paste("(a)  ",X[1]))) else
   plot(X1,Y,type="n",xlab="(a) X1")
X1new <- seq(-1.414,1.414,len=50)
lines(X1new,predict(m0,newdata=data.frame(Speed=X1new,Feed=rep(-1.414,50))))
lines(X1new,predict(m0,newdata=data.frame(Speed=X1new,Feed=rep(0.0,50))))
lines(X1new,predict(m0,newdata=data.frame(Speed=X1new,Feed=rep(1.414,50))))
text(-1.414,1.0,"Feed=-1.414",adj=0,cex=0.7)
text(0.0,2.0,"Feed=0.0",adj=0,cex=0.7)
text(0.5,4.0,"Feed=+1.414",adj=0,cex=0.7)
# part (b)
if(is.null(version$language) == FALSE)
   plot(X2,Y,type="n",xlab=expression(paste("(b)  ",X[2]))) else
   plot(X2,Y,type="n",xlab="(b) X2")
X2new <- seq(-1.414,1.414,len=50)
lines(X2new,predict(m0,newdata=data.frame(Speed=rep(-1.414,50),Feed=X2new)))
lines(X2new,predict(m0,newdata=data.frame(Speed=rep(0,50),Feed=X2new)))
lines(X2new,predict(m0,newdata=data.frame(Speed=rep(+1.414,50),Feed=X2new)))
text(-1.0,0,  "Speed=33",adj=0,cex=0.7)
text(0,2.,"Speed=35",adj=0,cex=0.7)
text(.45,5,"Speed=37",adj=0,cex=0.7)
par(oldpar)

dev.off()



#
# Note the code below is not correct for this problem but could be modified to be correct ... 
# 


#--
# the delta method for the estimates of tilde{X}_1 and tilde{X}_2
#-- 

# Extract the estimated coefficients:
#
b0  <- coef(m1)[1]
b1  <- coef(m1)[2]
b2  <- coef(m1)[3]
b11 <- coef(m1)[4]
b22 <- coef(m1)[5]
b12 <- coef(m1)[6]

# Extract the covariance matrix for the beta_hat estimates.  This is \sigma^2 (X' X)^{-1}
#
Var <- vcov(m1)

# Compute the needed analytic derivatives (the gradient) of \tidle{X}_1:
# 
x1m <- "-(2*b1*b22+b2*b12)/(4*b11*b22 - b12^2)"
x1m.expr <- parse(text=x1m)

# this will evaluate x1m at the estimated above beta's:
eval(x1m.expr) 

# this will compute the analytic derivatives:
derivs <- c(D(x1m.expr,"b0"),D(x1m.expr,"b1"),D(x1m.expr,"b2"),D(x1m.expr,"b11"),D(x1m.expr,"b22"),D(x1m.expr,"b12"))
derivs

# this will compute the numeric gradient at our estimated beta_hat: 
eval.derivs<-c(eval(D(x1m.expr,"b0")),eval(D(x1m.expr,"b1")),eval(D(x1m.expr,"b2")),eval(D(x1m.expr,"b11")),eval(D(x1m.expr,"b22")),eval(D(x1m.expr,"b12")))

# now we can evalaute the standard error of \tilde{X}_1:  
sqrt(t(eval.derivs) %*% Var %*% eval.derivs)

# or we can use the ALR3 package "delta.method" (but we need to specify the parameter names correctly):
delta.method(m1,"-(2*b1*b4+b2*b5)/(4*b3*b4 - b5^2)")