# This file is posted as  dogs.ssc

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


#  Make a new data frame with the missing 
#  values deleted
dogs2<-dogs[!is.na(dogs$Y),]


# Compute sample means 
# and make a profile plot.

means <- tapply(dogs2$Y, 
    list(dogs2$Drug,dogs2$Disease),FUN=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(dogs2$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')


# Plot mean responses for the levels of each factor

plot.design(dogs2,fun=mean)



# 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

dogs.out1 <- lme(fixed=Y~Drug*Disease,data=dogs2,
    random = ~1| Clinic)

# Print F-tests for fixed effects.  These
# are Type I sums of squares.
anova(dogs.out1)

# Print estimates of variance components
VarCorr(dogs.out1)

# Print estimates of fixed effects
fixef(dogs.out1)

# Print predictions of random effects
ranef(dogs.out1)

# Print approximate confidence intervals
intervals(dogs.out1)

# Print the covariance matrix for estimates of fixed effects
dogs.out1$varFix


#  Now fit Drug effects after adjusting for 
#  Disease effects

dogs.out2 <- lme(fixed=Y~Disease*Drug,data=dogs2,
    random = ~1 | Clinic)

# Print F-tests for fixed effects.  These
# are Type I sums of squares.
anova(dogs.out2)


# Create residual plots

frame( )
par(cex=1.4,mex=1.5,lwd=3,pch=2,mkh=0.1)
  plot(dogs.out1$fitted, dogs.out1$resid,
       xlab="Estimated Means",
        ylab="Residuals",
        main="Residual Plot")
  abline(h=0, lty=2, lwd=3)

  qqnorm(dogs.out1$resid,ylab="Residuals",
        main="Normal Plot: Residuals")
  qqline(dogs.out1$resid)




# Obtain estimates of mean responses for 
# conbunations of drugs and diseases
# and their standard errors

# Use the contrast matrices in the object
# created by lme( ) to obtain the appropriate
# coefficient matrix.

da <- dogs.out1$contrast$Drug
db <- dogs.out1$contrast$Disease
nrowa<- nrow(da)
nrowb<- nrow(db)
nab<-nrowa*nrowb

oneab<- matrix(1,nrow=nab,ncol=1)
onea<- matrix(1,nrow=nrowa,ncol=1)
oneb<-matrix(1,nrow=nrowb,ncol=1)
cc<-cbind(oneab, kronecker(oneb,da), kronecker(db,onea),
      kronecker(db,da))
 
amean <- cc%*%fixef(dogs.out1)
vmean<- cc%*%dogs.out1$varFix%*%t(cc)
std.amean <- as.matrix(sqrt(diag(vmean)))
amean<-data.frame(amean, std.amean) 
dimnames(amean)[[2]]<-c("Estimated mean","Standard error")
labels <- 10*kronecker(oneb,1:nrowa)+kronecker(1:nrowb,onea)
dimnames(amean)[[1]]<-labels
amean
 
# Construct simultaneous confidence intervals for 
# differences in mean blood pressure increase using 
# the Tukey honest significant difference (HSD)

nc <- nab*(nab-1)/2
dferror <- nrow(dogs2)-nab-nrow(ranef(dogs.out1))
tt <- qtukey(.95, k=nab, df=dferror)/sqrt(2)

c1<-matrix(0,1,nab)
at<-matrix(0,1,5) 

for(i in 1:(nab-1)){ ip1<-i+1
	for(j in ip1:nab){

c1[1,i]<- 1
c1[1,j]<- -1
c1.est <- c1%*%amean[ ,1]
c1.std <- sqrt(c1%*%vmean%*%t(c1)) 
c1.est.lower <- c1.est - tt*c1.std
c1.est.upper <- c1.est + tt*c1.std

at<-rbind(at, cbind(labels[i], -labels[j], c1.est.lower, c1.est, c1.est.upper))
c1 <- matrix(0,1,nab)

  } }

at<-at[-1, ]
dimnames(at)[[2]]<-c("Level", "-Level", "Lower limit",
    "Difference","Upper limit")
at