# 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