# This file is posted as  dogs.ssc

dogs <- read.table(
   "c:/courses/st511/snew/dogs.dat",
   col.names=c("Drug","Disease","Y"))
dogs$Drug <- as.factor(dogs$Drug)
dogs$Disease <- as.factor(dogs$Disease)


# Compute sample means 
# and make a profile plot.

means <- tapply(dogs$Y,
    list(dogs$Drug,dogs$Disease),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(dogs$Drug)
matplot(c(1,4), c(-10,50), type="n", 
     xlab="Drug",  ylab="Mean Response", 
     main= "Change in Systolic Blood Pressure") 

#  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(2.1,49.6,
   legend=c('Disease 1','Disease 2','Disease 3'),
               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~Drug*Disease,data=dogs)
summary.aov(lm.out1, ssType=1)
lm.out1$coef

lm.out2 <- lm(Y~Disease*Drug,data=dogs)
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(dogs$Disease,dogs$Drug,dogs$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)
  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))