# 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))