# This file is posted as macid.ssc

#  Enter the data into a data frame and 
#  define factors

macid <- read.table("c:\\courses\\st511\\sas\\macid.dat", 
   col.names=c("plant","leaf","method","y"))
macid$plant <- as.factor(macid$plant)
macid$leaf <- as.factor(macid$leaf)
macid$method <- as.factor(macid$method)


# Plot meean responses for the levels of each factor

plot.design(macid,fun=mean)


# Use the lme function to fit a model
# with random plant and leaf effects. 

options(contrasts=c(factor="contr.treatment",
                    ordered="contr.poly"))
macid.lme <- lme(fixed=y ~ method , data=macid,
    random = ~1|plant/leaf)

# Compute an ANOVA table for fixed effects

anova(macid.lme)

# Print estimates of variance components
VarCorr(macid.lme)

# Print estimates of fixed effects
fixef(macid.lme)

# Print predictions of random effects
ranef(macid.lme)

# Print approximate confidence intervals
intervals(macid.lme)

# Print the covariance matrix for estimates of fixed effects
macid.lme$varFix


# Obtain estimates of mean responses for the macidoline 
# methods and their standard errors

cc <- matrix(c(1, 0, 0, 1, 1, 0, 1, 0, 1 ),
        nrow=3,byrow=T)
amean <- cc%*%fixef(macid.lme)
vmean<- cc%*%macid.lme$varFix%*%t(cc)
std.amean <- as.matrix(sqrt(diag(vmean)))
amean<-data.frame(amean, std.amean) 
dimnames(amean)[[2]]<-c("method mean","Standard error")
amean
 
# Construct simultaneous confidence intervals for 
# differences in method means using the Tukey 
# honest significant difference (HSD)

tt <- qtukey(.95, k=3, df=22)/sqrt(2)

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

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

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(i, j, c1.est.lower, c1.est, c1.est.upper))
c1 <- matrix(0,1,3)

  } }

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