# 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