################################################################# ## Code for Book: Applied Statistical Analysis in R ## A Graduate course for Psychology, Human factors, and Data Science ## (c) 2018 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 18: ## Factorial ANOVA models with interactions. ## ##For the tooth growth data: we see that there is an interaction between ##supp and len: tapply(ToothGrowth$len,list(ToothGrowth$supp,ToothGrowth$dose),mean) tg1 <- aov(len~supp+as.factor(dose),data=ToothGrowth) tg2 <- aov(len~supp*as.factor(dose),data=ToothGrowth) tg2 <- aov(len~supp + as.factor(dose) + supp:as.factor(dose),data=ToothGrowth) tg3 <- aov(len~supp:as.factor(dose),data=ToothGrowth) summary(tg1) summary(tg2) summary(tg3) anova(tg2,tg1) ##how do the coefficients change as we add different interaction models? summary(lm(len~supp*as.factor(dose),data=ToothGrowth)) summary(lm(len~supp:as.factor(dose),data=ToothGrowth)) ##we can obtain complete tables of means. This gives the relative effects: model.tables(tg2) model.tables(tg3) ##this gives the absolute mean values model.tables(tg2,type="means") model.tables(tg3,type="means") set.seed(100) ultim1 <- data.frame(cond =rep(rep(c("N","F","M"),2),100), gender= rep(rep(c("F","M"),each=3),100), offer = rep(c(5.5, 6.5, 6.5, 5, 6.1, 6.0),100)+rnorm(600)) table(ultim1$cond,ultim1$gender) #pdf("c18-ultim1.pdf",width=8,height=5) matplot(tapply(ultim1$offer,list(ultim1$cond,ultim1$gender),mean),cex=2.5, type="b",xaxt="n",pch=16,ylab="Mean offer in ultimatum game", xlab="Gender of opponent") matplot(tapply(ultim1$offer,list(ultim1$cond,ultim1$gender),mean),cex=.8, col="white", type="p", pch=c("F","M"), add=T) axis(1,1:3,c("Female","Male","Unspecified")) legend(1.5,5.5,c("Male","Female"),pch=c("F","M"),col=1:2,lty=1:2) #dev.off() library(car) model1 <- aov(offer~cond*gender,data=ultim1) model1b <- aov(offer~cond+gender,data=ultim1) TukeyHSD(model1) TukeyHSD(model1b) Anova(model1) ############## Anova(model1,type="III") ##incorrect contrasts(ultim1$cond) <- contr.poly(3) contrasts(ultim1$gender) <- contr.poly(2) model1c <- aov(offer~cond*gender,data=ultim1) Anova(model1c,type="III") ##correct ##reset to treatment contrasts contrasts(ultim1$cond) <- contr.treatment(levels(ultim1$cond)) contrasts(ultim1$gender) <- contr.treatment(levels(ultim1$gender)) ultim1b <-ultim1[sample(1:nrow(ultim1),size=nrow(ultim1),replace=T),] table(ultim1b$cond,ultim1b$gender) model1b <- aov(offer~cond*gender,data=ultim1b) Anova(model1b,type="II") Anova(model1b,type="III") ##not right contrasts(ultim1b$cond) <- contr.poly(3) contrasts(ultim1b$gender) <- contr.poly(2) model1c <- aov(offer~cond*gender,data=ultim1b) Anova(model1c,type="III") ###There was no interaciton, so we can look at just he main effects: model1d <- aov(offer~cond+gender,data=ultim1b) Anova(model1d) TukeyHSD(model1d) ##this looks at just one effect at a time pairwise.t.test(ultim1$offer,ultim1$cond) pairwise.t.test(ultim1$offer,ultim1$gender) TukeyHSD(model1d) ##ultimatum problem #2 set.seed(99) ultim.tmp <- data.frame(cond =rep(rep(c("N","F","M"),2),100), gender= rep(rep(c("F","M"),each=3),100), offer = rep(c(5, 6.2, 6.3, 5, 6.4,5.1),100)+rnorm(600)) ultim2 <- ultim.tmp[sample(1:nrow(ultim.tmp),size=nrow(ultim.tmp),replace=T),] pdf("c18-ultim2.pdf",width=8,height=5) matplot(tapply(ultim2$offer,list(ultim2$cond,ultim2$gender),mean),cex=2.5, type="b",xaxt="n",pch=16,ylab="Mean offer in ultimatum game", xlab="Gender of opponent") matplot(tapply(ultim2$offer,list(ultim2$cond,ultim2$gender),mean),cex=.8, col="white", type="p", pch=c("F","M"), add=T) axis(1,1:3,c("Female","Male","Unspecified")) legend(1.5,6.,c("Male","Female"),pch=c("F","M"),col=1:2,lty=1:2) dev.off() model2 <- aov(offer~cond*gender,data=ultim2) Anova(model2) contrasts(ultim2$cond) <- contr.poly(3) contrasts(ultim2$gender) <- contr.poly(2) model2a <- aov(offer~cond*gender,data=ultim2) Anova(model2a,type="II") ##type II: test effect of interaction after main effects considered Anova(model2a,type="III") ##type III: test effect of main effects after interaction considered ##Approach 2: TukeyHSD(model2a) ##Approach 3: Anova(aov(offer~gender,data=ultim2[ultim2$cond=="N",])) Anova(aov(offer~gender,data=ultim2[ultim2$cond=="M",])) Anova(aov(offer~gender,data=ultim2[ultim2$cond=="F",])) ## Approach 4 Anova(aov(offer~gender*cond,data=ultim2[ultim2$cond!="N",])) model.tables(aov(offer~gender*cond,data=ultim2[ultim2$cond!="N",])) ##transform to a one-way: cond3 <- as.factor( paste(ultim2$cond,ultim2$gender,sep="-")) contrasts(cond3) <- cbind(a=c(0,0,0,0,1,-1), b=c(1,-1,-1,1,0,0), c=c(2,0,-1,0,-1,0), d=rep(0,6), e=rep(0,6), f=rep(0,6)) summary(lm(ultim2$offer~cond3)) ## This data does not exist here ##x <- aov(rt~target+foil) ################################### ##Effect sizes ##http://stats.stackexchange.com/questions/2962/omega-squared-for-measure-of-effect-in-r omega_2 <- function(aov_in, neg2zero=T){ aovtab <- summary(aov_in)[[1]] n_terms <- length(aovtab[["Sum Sq"]]) - 1 output <- rep(-1, n_terms) SSr <- aovtab[["Sum Sq"]][n_terms + 1] MSr <- aovtab[["Mean Sq"]][n_terms + 1] SSt <- sum(aovtab[["Sum Sq"]]) for(i in 1:n_terms){ SSm <- aovtab[["Sum Sq"]][i] DFm <- aovtab[["Df"]][i] output[i] <- (SSm-DFm*MSr)/(SSt+MSr) if(neg2zero & output[i] < 0){output[i] <- 0} } names(output) <- rownames(aovtab)[1:n_terms] return(output) } ##Or library(sjstats) ToothGrowth$dose2 <- as.factor(ToothGrowth$dose) aov.tootha <- aov(len~supp+dose2,data=ToothGrowth) anova_stats(aov.tootha) set.seed(99) ultim.tmp <- data.frame(cond =rep(rep(c("N","F","M"),2),100), gender= rep(rep(c("F","M"),each=3),100), offer = rep(c(5, 6.2, 6.3, 5, 6.4,5.1),100)+rnorm(600)) ultim2 <- ultim.tmp[sample(1:nrow(ultim.tmp),size=nrow(ultim.tmp),replace=T),] model1 <- aov(offer~cond*gender,data=ultim2) model1 <- lm(offer~cond*gender,data=ultim2) omega_sq(model1) anova_stats(model1)