# This file is posted as batteries.ssc

#  First create a data frame

battery <- read.table(
   "c:/courses/st511/sas/batteries.dat",
   col.names=c("Material","TempF","Y"))
battery$Material <- as.factor(battery$Material)
battery$Temp <- as.factor(battery$TempF)


# Compute sample means 
# and make a profile plot.

means <- tapply(battery$Y,
    list(battery$Temp,battery$Material),mean)
means


# Set up the axes and title of the 
# profile plot.  

par(fin=c(7,7),cex=1.2,lwd=3,mex=1.5)
x.axis <- unique(battery$TempF)
matplot(c(15,125), c(0,200), type="n", 
     xlab="Temperature (F)",  ylab="Mean Life (hrs)", 
     main= "Battery Life") 

#  Add a profile for each soil type

matlines(x.axis,means,type='l',lty=c(1,3,7),lwd=3)     

# Plot points for the individual observations

matpoints(x.axis,means, pch=c(1,16,18))    

# Add a legend to the plot

legend(20,65,
   legend=c('Material A','Material B','Material C'),
               lty=c(1,3,7),bty='n')


# The default contrast matrices can be changed by
# resetting the contrasts options.  The contr.treatment
# restricts some parameters to be zero.

options(contrasts=c('contr.treatment','contr.ploy'))


# Fit a model with main effects and interaction
# Compute both sets of Type I sums of squares

lm.out1 <- lm(Y~Material*Temp,data=battery)
summary.aov(lm.out1, ssType=1)
lm.out1$coef

lm.out2 <- lm(Y~Temp*Material,data=battery)
summary.aov(lm.out2, ssType=1)
lm.out2$coef


# Compute Type III sums of squares and F-tests.

summary.aov(lm.out1, ssType=3)


# Create a data frame containing the original 
# data and the residuals and estimated means

data.frame(battery$Temp,battery$Material,battery$Y,
              Pred=lm.out1$fitted,
             Resid=round(lm.out1$resid,3))


# Create residual plots

frame( )
par(cex=1.4,mex=1.0,lwd=3,pch=2,mkh=0.1)

  plot(lm.out1$fitted, lm.out1$resid,
       xlab="Estimated Means",
       ylab="Residuals",
        main="Residual Plot")
  abline(h=0, lty=2, lwd=3)

    qqnorm(lm.out1$resid, 
        xlab="Standard normal quantiles",
        ylab="Residuals", 
        main="Normal Probability Plot")
  qqline(lm.out1$resid)


# Create plots for studentized residulas
# You must attach the MASS library
# to have access to the studres( )
# function that computes studentized
# residuals in the following code

library(MASS)

frame( )
par(cex=1.0,mex=1.0,lwd=3,pch=2,
     mkh=0.1,fin=c(6.5,6.5))
  plot(lm.out1$fitted, studres(lm.out1),
       xlab="Estimated Means",
       ylab="Studentized Residuals",
       main="Studentized Residual Plot")
  abline(h=0, lty=2, lwd=3)

  qqnorm(studres(lm.out1),
        xlab="Standard normal quantiles",
        ylab="Residuals", 
       main="Studentized Residuals")
  qqline(studres(lm.out1))