# This file is posted as gas.additive.ssc

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

gas <- read.table("c:\\courses\\st511\\sas\\gas.additive.dat", 
   col.names=c("driver","auto","additive","y"))
gas$driver <- as.factor(gas$driver)
gas$auto <- as.factor(gas$auto)
gas$additive <- as.factor(gas$additive)


# Plot meean responses for the levels of each factor

plot.design(gas,fun=mean)


# Use the lme function to fit a model
# with random block effects. The pdBlocked
# function is used to specify a positive definite
# block diagonal matrix for the Covariance matrix
# of the joint set of 4 driver random effects and 4 
# auto random effects.  We called this matrix G in
# the course notes.  Results match with PROC MIXED
# in SAS with the except of the denominator df in 
# the F-tests and t-tests for fixed effects.   

options(contrasts=c(factor="contr.treatment",
                    ordered="contr.poly"))
gas.lme <- lme(fixed=y ~ additive , data=gas,
    random = pdBlocked(list(pdIdent(~driver-1),
               pdIdent(~auto-1))))

# Print F-tests for fixed effects.  The denominator
# degrees of freedom are too large.  Degrees of 
# for the block effects were not removed.

anova(gas.lme)

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

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

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

# Print approximate confidence intervals
intervals(gas.lme)

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


# Obtain estimates of mean responses for the gasoline 
# additives and their standard errors

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

tt <- qtukey(.95, k=4, df=6)/sqrt(2)

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

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

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,4)

  } }

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