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