### R code from vignette source '/home/jimrc/classes/stat505/notes/ARMchapter3-2.Rnw'

###################################################
### code chunk number 1: setup
###################################################
kidiq <- read.csv("../data/arm_KidsIQ.csv")
options(width=60, digits=2)
require(xtable)


###################################################
### code chunk number 2: fit1
###################################################
summary(warp.fit <- lm(breaks ~ tension, warpbreaks))$coef
coef(warp.fit2 <- lm(breaks ~ tension -1, warpbreaks))
Lambda <- matrix(c( -1,1,0, -1,0,1, 0,1,-1),byrow=TRUE,3,3)
as.numeric(Lambda %*% coef(warp.fit2))
sqrt(diag( Lambda %*% summary(warp.fit2)$cov.unscaled %*%
           t(Lambda))) * summary(warp.fit2)$sigma


###################################################
### code chunk number 3: cov1
###################################################
xtable(summary(warp.fit)$cov.unscaled*summary(warp.fit)$sigma^2)


###################################################
### code chunk number 4: cov2
###################################################
xtable(summary(warp.fit2)$cov.unscaled*summary(warp.fit2)$sigma^2)


###################################################
### code chunk number 5: line1
###################################################
kidfit.2 <- lm (kid.score ~ mom.iq, kidiq)
plot ( kid.score~mom.iq,data=kidiq, xlab="Mother IQ score", 
      ylab="Child test score",  col= mom.hs+2, pch=20)
curve (coef(kidfit.2)[1] + coef(kidfit.2)[2]*x, add=TRUE, col="blue")
fit.3 <- lm (kid.score ~ mom.hs + mom.iq, kidiq)
curve (cbind (1, 1, x) %*% coef(fit.3), add=TRUE, col="darkgreen",lwd=2)
curve (cbind (1, 0, x) %*% coef(fit.3), add=TRUE, col="red",lwd=2)


###################################################
### code chunk number 6: line2
###################################################
kidfit.4 <- lm (kid.score ~  mom.hs * mom.iq, kidiq)
plot(kid.score ~ mom.iq, data=kidiq, xlab="Moms IQ", 
      ylab="Child test score",  col= mom.hs+2, pch=20)
curve (cbind (1, 1, x, 1*x) %*% coef(kidfit.4), add=TRUE, 
       col="darkgreen",lwd=2)
curve (cbind (1, 0, x, 0*x) %*% coef(kidfit.4), add=TRUE, 
       col="red",lwd=2)


###################################################
### code chunk number 7: ggplot
###################################################
 require(ggplot2)
 myplot <- qplot( mom.iq, kid.score, data=kidiq, shape= factor(mom.hs),
       colour = factor(mom.hs) ) +theme_bw()
 print( myplot + geom_smooth(method = "lm"))


###################################################
### code chunk number 8: predError
###################################################
 predictDF <- data.frame(mom.iq=rep(7:14*10,2), mom.hs=rep(0:1,each=8))
 predictDF <- cbind(predictDF, predict(kidfit.4, newdata=predictDF, interval="pred"))
  names(predictDF)[3] <- "kid.score"
  print( myplot +  geom_smooth(aes(ymin = lwr, ymax = upr), data = predictDF, stat="identity"))


###################################################
### code chunk number 9: diagnostics
###################################################
par(mfrow=c(1,4))
plot(kidfit.4)