# choccake.R
library(lsmeans)
library(lme4)
library(car)
#library(nlme)
#library(Matrix)
#library(pbkrtest)
# get the data
cake<-read.table(file="choccake.dat",header=F,
colClasses=c("factor","factor","factor","numeric"))
names(cake)<-c("recipe","replicat","baketemp","angle")
# the name "break" has problems and cannot be used in R
# create a replicat factor (called allreplicat below)
# that has a unique level for each replicat within each recipe.
cake <- within(cake, batch <- factor(replicat:recipe))
head(cake)
is.numeric(cake$recipe)
is.numeric(cake$angle)
is.factor(cake$batch)
# the aov function gives the traditoinal type III anova F tests. It actually
# fits two separate models, one at each level of experimental unit
m0 <-aov(angle ~ recipe*baketemp + Error(batch/baketemp), data=cake)
summary(m0)
# the lmer will fit the model as a single mixed effect model and the anova
# function will give the standard F tests.
m1 <- lmer(angle ~ recipe*baketemp + (1|batch), data=cake)
anova(m1)
# Now for contrasts. Construct the contrasts of interest
c1<-c(1,1,-2) # recipe contrast 1: normal vs extra sugar
c2<-c(1,-1,0) # recipe contrast 2: choc temp
c3<-c(-5,-3,-1,1,3,5) # linear in baketemp
c41<-c(5,-1,-4,-4,-1,5) # quadratic in baketemp
c42<-c(-5,7,4,-4,-7,5) # cubic in baketemp
c43<-c(1,-3,2,2,-3,1) #quartic in baketemp
c44<-c(-1,5,-10,10,-5,1) #quintic in baketemp
c5 <- c(1,0,0,-1,rep(0,14)) # mu11-mu12
c6 <- c(1,-1,rep(0,16)) # mu11-m21
lsmeans(m1,specs=list(lsm1~recipe,lsm2~baketemp,lsm3~recipe:baketemp),
contr=list(lsm1=list(normal_vs_extra.sugar=c1,
chocolate_temp=c2),
lsm2=list(linear_baketemp=c3),
lsm3=list(mu11_mu12=c5,
mu11_mu21=c6)
))
K1<-rbind("nonlinear baketemp I"=c41,"nonlinear baketemp II"=c42,
"nonlinear baketemp III"=c43,"nonlinear baketemp IV"=c44)
# It seems like the following line (commented out) should work to get the
# test of nonlinearity in baking temperature, but it does not
#summary(glht(m1,linfct=mcp(baketemp=K1)),test=Ftest())
# Instead, I am going to use the linearHypothesis function. To make
# this easier, I will refit the model with a cell means parameterization
# for the treatments. Then create the contrast matrix that will test
# nonlinearity in baking temperature (averaged over recipe) if applied to the
# treatment means. That matrix is K1a below.
m1a <- lmer(angle ~ -1+recipe:baketemp + (1|batch), data=cake)
m1a
(K1a <- K1%x%matrix(1,nrow=1,ncol=3))
linearHypothesis(m1a,K1a,test="F") # gives nonlinearity in bakiing temp test
# Finally, get the profile plot
lsmip(m1,recipe~baketemp)