# This file is posted as carrots.ssc carrot <- read.table( "c:/courses/st511/snew/carrots.dat", col.names=c("Soil","Variety","Days")) carrot\$Soil <- as.factor(carrot\$Soil) carrot\$Variety <- as.factor(carrot\$Variety) # Compute sample means of germination # times for all combinations of soil # type and varieties of carrot seeds # and make a profile plot. # At this point UNIX users should open # a graphics window with the motif( ) # function means <- tapply(carrot\$Days, list(carrot\$Variety,carrot\$Soil),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(carrot\$Variety) matplot(c(1,3,1), c(0,40,10), type="n", xlab="Variety", ylab="Mean Time", main= "Average Time to Carrot Seed Germination") # Add a profile for each soil type matlines(x.axis,means,type='l',lty=c(1,3),lwd=3) # Plot points for the individual observations matpoints(x.axis,means, pch=c(1,16)) # Add a legend to the plot legend(2.2,38.6, legend=c('Soil Type 1','Soil Type 2'), lty=c(1,3),bty='n') # Fit a model with main effects and interaction # Compute both sets of Type I sums of squares # The default contrast matrices can be changed by # resetting the contrasts options. The contr.sum # option restricts parameters to sum to zero # across the levels of any single factor. This # affects the values of parameter estimates # and their interpretation, but does not affect # the ANOVA tables. options(contrasts=c('contr.sum','contr.ploy')) # Produce ANOVA's with Type I sums of squares lm.out1 <- lm(Days~Soil*Variety,data=carrot) anova(lm.out1) lm.out2 <- lm(Days~Variety*Soil,data=carrot) anova(lm.out2) # print the model matrix for the second model model.matrix(lm.out2) # Create a data frame containing the original # data and the residuals and estimated means data.frame(carrot\$Soil,carrot\$Variety,carrot\$Days, Pred=lm.out1\$fitted, Resid=round(lm.out1\$resid,3)) # Create residual plots frame( ) par(fin=c(7,7),cex=1.2,lwd=3,pch=2,mex=1.5) plot(lm.out1\$fitted, lm.out1\$resid, xlab="Estimated Means", ylab="Residuals") abline(h=0, lty=2, lwd=3) qqnorm(lm.out1\$resid) 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), main="Studentized Residuals") qqline(studres(lm.out1)) # Compute Type III sums of squares and F-tests. summary.aov(lm.out1,ssType=3)