# 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)