# mpg.R
library(lsmeans)
library(car)
library(plyr)
#library(phia)
#library(doBy)
# get the data
mpg <- read.table(file="mpg.dat",header=T,
colClasses=c("factor","factor","numeric"))
head(mpg)
is.factor(mpg$car)
is.factor(mpg$mileage)
# the intra-block analysis of a BIBD using ANOVA
m1<-aov(mileage~add+car,data=mpg)
summary(m1)
# note that the summary function gives F tests based on type I
# or sequential sums of squares. To get F tests based on the type III
# SS, we need the Anova() function from the car package.
Anova(m1,type=3)
# get the contrasts of interest
c1<-c(2,2,2,-3,-3)
c2<-c(-1,0,1,0,0)
c3<-c(1,-2,1,0,0)
c4<-c(0,0,0,1,-1)
lsmeans(m1,specs=lsm~add,contr=list(lsm=list(chemicals=c1,
linear_chem1=c2,
nonlinear_chem1=c3,
doses_chem2=c4)))
# use package plyr to get means instead of lsmeans for each level of add
# and car. These are not the right estimators of these means, but are produced
# here for comparison with the lsmeans, which should be used instead.
ddply(mpg, .(add), summarize, mean_mileage=mean(mileage))
ddply(mpg, .(car), summarize, mean_mileage=mean(mileage))
# Plot of estimated means for each additive, by car. Note that this is an
# additive model, so the mean profiles are necessarily parallel (the model
# constrains them to be so). Allowing them to be otherwise would require a
# replicated experiment and a more complex model (e.g., one with a car*add
# interaction)
lsmip(m1,add~car)
# plot of add means to show the "dose" effect for each chemical
# create variables "dose" and "chemical"
addmeans<-lsmeans(m1,specs=~add)[[1]]
addmeans$dose<-c(1,2,3,1,2)
addmeans$chemical<-c(1,1,1,2,2)
head(addmeans)
plot(addmeans$dose[1:3],addmeans$lsmean[1:3],ylim=c(10,15),
type="b",lty=1,pch=1,col=1,
main="Estimated Treatment Means, MPG Example",ylab="lsmeans of add",
xlab="dose")
lines(addmeans$dose[4:5],addmeans$lsmean[4:5],type="b",lty=2,pch=2,col=2)
legend(2,15,legend=c("chemical 1","chemical 2"),lty=c(1,2),pch=c(1,2),col=c(1,2))