################################################################# ## Code for Book: Applied Statistical Analysis in R ## A Graduate course for Psychology, Human factors, and Data Science ## (c) 2018-2020 Shane T. Mueller, Ph.D. ## Michigan Technological University ## shanem@mtu.edu ## ## Textbook and videos available at http://pages.mtu.edu/~shanem/psy5210 ## ## ## This file contains example software code relevant for Chapter 17: Multi-way ANOVA models. ## #pdf("c17-tooth.pdf",width=8,height=5) interaction.plot(ToothGrowth$dose,ToothGrowth$supp,ToothGrowth$len,pch=15, type="b",xlab="Dosage (mg)", ylab="Tooth length (mm)", bty="n", las=1,xtick=T,leg.bg="white",leg.bty="n", ylim=c(0,30), col=c("orange","red"),trace.label="Supplement") #dev.off() ToothGrowth$dosage <- as.factor(ToothGrowth$dose) lm.tooth3 <- lm(len~supp+dosage,data=ToothGrowth) summary(lm.tooth3) anova(lm.tooth3) ##plot the model predictions: modelfit <- tapply(lm.tooth3$fit,list(dosage=ToothGrowth$dosage,supp=ToothGrowth$supp),mean) #pdf("c17-toothmodel.pdf",width=8,height=5) interaction.plot(ToothGrowth$dose,ToothGrowth$supp,ToothGrowth$len,pch=15, type="b",xlab="Dosage (mg)", ylab="Tooth length (mm)", bty="n", las=1,xtick=T,leg.bg="white",leg.bty="n", ylim=c(0,30), col=c("orange","red"),trace.label="Supplement") matplot(modelfit,type="o",pch=1,add=T,lty=1,col='black') #dev.off() #Compare to just dosage: l1 <- lm(len~supp,data=ToothGrowth) l2 <- lm(len~dosage,data=ToothGrowth) l3 <- lm(len~supp+dosage,data=ToothGrowth) l0 <- lm(len~1,data=ToothGrowth) summary(lm(len~dosage,data=ToothGrowth)) #What changed? ##again, these are the same: anova(lm.tooth3) summary(aov(len~supp+dosage,data=ToothGrowth) ) aov(len~dosage,data=ToothGrowth) summary(aov(len~supp+dosage,data=ToothGrowth) ) summary(aov(len~dosage,data=ToothGrowth) ) ##Building the ANOVA table: ##underlying models lm.int <- lm(len~1,data=ToothGrowth) lm.supp <- lm(len~supp,data=ToothGrowth) lm.dose <- lm(len~dosage,data=ToothGrowth) lm.both <- lm(len~dosage+supp,data=ToothGrowth) ##anova version: m.int <- aov(len~1,data=ToothGrowth) m.supp <- aov(len~supp,data=ToothGrowth) m.dose <- aov(len~dosage,data=ToothGrowth) m.both <- aov(len~dosage+supp,data=ToothGrowth) m.int summary(m.int) summary(m.both) summary(m.supp) anova(m.both) ##can we find where the values from here come from? anova(m.supp) anova(m.dose) anova(m.supp,m.int) anova(m.both,m.supp) anova(m.both,m.dose) anova(m.both) m.int <- aov(len~0,data=ToothGrowth) m.supp <- aov(len~supp,data=ToothGrowth) m.dose <- aov(len~dosage,data=ToothGrowth) m.both <- aov(len~dosage+supp,data=ToothGrowth) ############post-hoc tests for multiple predictor ANOVA TukeyHSD(m.both) library(agricolae) ##from agricolae package, there is a nice implementation HSD.test(m.both,"dosage", group=TRUE, main="Effect of dosage",console=T) ############# ##Multi-way Bayes factor ANOVA: library(BayesFactor) #abf2 <- anovaBF(len~dose+supp,data=ToothGrowth) abf2 <- anovaBF(len~dosage+supp,data=ToothGrowth) summary(abf2) ##model comparison: abf2[3]/abf2[2] abf2[3]/abf2[1] ############## ##Exercise} Perform an ANOVA on the insect multi-way ANOVA OrchardSprays outcome variable 'decrease'. Use rowpos, colpos, and treatment as categorical predictors. Use a Tukey test to determine the relative effects of row, column, and treatment. data(OrchardSprays) lm1 <- lm(decrease~rowpos+colpos+treatment, data=OrchardSprays) am1 <- aov(decrease~rowpos+colpos+treatment, data=OrchardSprays) summary(am1) am2 <- aov(decrease~as.factor(rowpos)+as.factor(colpos)+treatment, data=OrchardSprays) summary(am2) TukeyHSD(am2,"treatment") library(agricolae) ##from agricolae package, there is a nice implementation HSD.test(am2,"treatment", group=TRUE, main="Effect of treatment",console=T) emmeans() ################################# ##Non-orthogonal predictors ##The following has an orthogonal design: color <- c("r","r","r","r","r","r","r","r","r", "g","g","g","g","g","g","g","g","g", "b","b","b","b","b","b","b","b","b") setsize <- c(1,2,3,1,2,3,1,2,3, 1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1,2,3) stimsize <- c(1,1,1,2,2,2,3,3,3,1,1,1,2,2,2,3,3,3,1,1,1,2,2,2,3,3,3) set.seed(100) rt <- rnorm(27) table(color,setsize,stimsize) lm1 <- lm(rt~color+setsize+stimsize) lm1a <- lm(rt~setsize+stimsize+color) summary(lm1) anova(lm1) anova(lm1a) anova(lm(rt~color+setsize+stimsize),lm(rt~setsize+stimsize)) ##how should we test the effect of setsize? ##note the delta sum of sq is the same for each of these: anova(lm(rt~color+setsize+stimsize),lm(rt~color+stimsize)) anova(lm(rt~setsize+stimsize),lm(rt~stimsize)) anova(lm(rt~color+setsize),lm(rt~color)) anova(lm(rt~setsize)) #### non-orthogonal design: set.seed(100) color <- sample(c("r","g","b"),27,replace=T) setsize <- sample(1:3,27,replace=T) stimsize <- sample(1:3,27,replace=T) rt <- rnorm(27) table(color,setsize,stimsize) ##now the delta sum of squares depends on where we look in the model lattice anova(lm(rt~color+setsize+stimsize),lm(rt~color+stimsize)) anova(lm(rt~setsize+stimsize),lm(rt~stimsize)) anova(lm(rt~color+setsize),lm(rt~color)) anova(lm(rt~setsize)) ##unbalanced design: table(color,setsize) ##moreover, the F values now depend on the order in the lm: lmz1 <- lm(rt~color+setsize) lmz2 <- lm(rt~setsize+color) lmz1 anova(lmz1) lmz2 anova(lmz2) ## The direct way of making the model: ## z.anova1 <- aov(rt ~ setsize + color) z.anova2 <- aov(rt ~ color+setsize) summary(z.anova1) summary(z.anova2) ##where do the SSE values actually come from? lm0 <- lm(rt~1) lmcolor <-lm(rt~color) lmsetsize <- lm(rt~setsize) lmcolorsetsize <- lm(rt~color+setsize) lmsetsizecolor <- lm(rt~setsize+color) anova(lm0,lmcolor) anova(lmsetsizecolor,lmsetsize) anova(lm0,lmsetsize) anova(lmcolor,lmsetsizecolor) ##here are the two top models as anovas: anova(lmsetsizecolor) anova(lmcolorsetsize) ##anova() in R is a 'type-I' ANOVA, which compares up and down the chain in the order that the ##arguments appear ## Alternative is compare the full model to each next-smaller model. ## drawback is that when you are not orthonal, SS do not add up. ##Solutions to non-orthogonality ##type-II/III ANOVA anova(lm(rt~color + setsize)) anova(lm(rt~setsize+color)) #Add up Sum Sq for each model, and new table: # model 1 5.748+.2457+18.0159 #model 2 .1705+5.8233+18.0159 #New combined ANOVA: 5.823+.2457+18.0159 drop1(lm(rt~color+setsize)) library(car) anova(lmsetsizecolor) anova(lmcolorsetsize) Anova(lmsetsizecolor) Anova(lmcolorsetsize) ##Model Lattice Exercise: set.seed(1000) #Everyone should get the same numbers. ## Build the full model lattice for the set of models with three predictors a <- factor(sample(c("I","II"),replace=T,50)) b <- factor(sample(c("i","ii"),replace=T,50)) c <- factor(sample(c("a","b","c"),replace=T,50)) y <- c(.3,-5)[a] - c(50,35)[b] + c(-.35,1.5,-5)[c] + rnorm(50) ##these are comparisons of different submodels of the model lattice anova(lm(y~a+b+c)) anova(lm(y~a+c+b)) anova(lm(y~b+a+c)) anova(lm(y~b+c+a)) anova(lm(y~c+a+b)) anova(lm(y~c+b+a)) anova(lm(y~a+b+c)) anova(lm(y~a),lm(y~1)) anova(lm(y~a+b),lm(y~a)) anova(lm(y~a+b+c),lm(y~a+b)) anova(lm(y~a),lm(y~1)) Anova(lm(y~a+b+c))