# Analyze the penicillin data from Box,
# Hunter, and Hunter.  This code is 
# posted as  penclln.ssc

# Enter the data into a data frame and 
# change the Batch and Process variables
# into factors 

penclln <- read.table("c:/courses/st511/snew/penclln.dat", 
     col.names=c("Batch","Process","Yield"))
penclln$Batch <- as.factor(penclln$Batch)
penclln$Process <- as.factor(penclln$Process)
penclln


# Construct a profile plot.  UNIX users
# should use the  motif( ) command to open 
# a graphics window
 
 attach(penclln)
 means <- tapply(Yield,list(Process,Batch),mean)
 
 par(fin=c(6,7),cex=1.2,lwd=3,mex=1.5)
 x.axis <- unique(Process)
 matplot(c(1,4), c(75,100), type="n", xaxt="n",
      xlab="Process",  ylab="Yield", 
      main= "Penicillin Production Results") 
 axis(1, at=(1:4)*1, labels=c("A", "B", "C", "D"))
 matlines(x.axis,means,type='l',lty=1:5,lwd=3)         
 
 legend(4.2,95, legend=c('Batch 1','Batch 2',
     'Batch 3','Batch 4','Batch 5'), lty=1:5,bty='n') 
 detach( )


#  Use the  lme( )  function to fit a model
#  with additive batch (random) and process 
#  (fixed) effects and create diagnostic plots.  

options(contrasts=c("contr.treatment",
    "contr.poly"))

penclln.lme <- lme(Yield ~ Process, 
    random= ~ 1|Batch, data=penclln,
    method=c("REML"))

summary(penclln.lme)
names(penclln.lme)

#  Contruct ANOVA table for fixed effects
anova(penclln.lme)

#  Estimated parameters for fixed effects
coef(penclln.lme)

# BLUP's for random effects
ranef(penclln.lme)

# Confidence intervals for fixed effects
# and estimated standard deviations
intervals(penclln.lme)

# Create a listing of the original data 
# residuals and predicted values

data.frame(penclln$Process,penclln$Batch,
             penclln$Yield,
             Pred=penclln.lme$fitted,
             Resid=round(penclln.lme$resid,3))


# Create residual plots
frame( )
par(fin=c(7,7),cex=1.2,lwd=3,mex=1.5)
  plot(penclln.lme$fitted, penclln.lme$resid,
       xlab="Estimated Means",
       ylab="Residuals",
       main="Residual Plot")
  abline(h=0, lty=2, lwd=3)


  qqnorm(penclln.lme$resid)
  qqline(penclln.lme$resid)