#Chapter 11 Factorial ANOVA #Dataset = Obarow(2004); Obarrow.Story.1; Writing library(Rcmdr) library(foreign, pos=15) obarowStory1 <- read.spss(file.choose(), use.value.labels=TRUE, max.value.labels=Inf, to.data.frame=TRUE) numSummary(obarowStory1[,c("gnsc1.1", "gnsc1.2")], groups=obarowStory1$Treatment1, statistics=c("mean", "sd")) #boxplots par(mfrow=c(1,2)) Boxplot(gnsc1.1~Treatment1, data=obarowStory1, id.method="y", col="lightgrey", xlab="Immediate Gainscore", ylab="Gainscore") Boxplot(gnsc1.2~Treatment1, data=obarowStory1, id.method="y", col="lightgrey", xlab="Delayed Gainscore", ylab="Gainscore") #Back to normal par(mfrow=c(1,1)) #Fig11.3 plotMeans(obarowStory1$gnsc1.1, obarowStory1$MusicT1, obarowStory1$PicturesT1, error.bars="conf.int") #RCommander does the same as follows: with(obarowStory1, plotMeans(gnsc1.1, MusicT1, PicturesT1, error.bars="conf.int")) #Fig11.4 with(obarowStory1, plotMeans(gnsc1.1, Treatment1, error.bars="conf.int")) summary(obarowStory1) #Show the entire value obarowStory1$Treatment1 #Show the levels of the factor variable "Treatment1" levels(obarowStory1$Treatment1) #If you want to change the order of factors: obarowStory1$Treatment1 <- ordered(obarowStory1$Treatment1, levels= c("No music Yes pics", "Yes music Yes pics", "Yes music No pics", "No music No pics")) #Then plot the means again with(obarowStory1, plotMeans(gnsc1.1, Treatment1, error.bars="conf.int")) #If you want to set the cut off points "below 18" for pretest: new.obarow <- subset(obarow, subset=PRETEST1<18) obarow$MusicT1 <- recode(obarow$trtmnt1, '1="no music"; 2="no music"; 3="yes music"; 4="yes music"; ', as.factor.result=TRUE) #Perform a full-model factorial ANOVA using lm() AnovaModel.1 <- (lm(gnsc1.1 ~ gender*MusicT1*PicturesT1, data=obarowStory1)) Anova(AnovaModel.1) # means with(obarowStory1, (tapply(gnsc1.1, list(gender, MusicT1, PicturesT1), mean, na.rm=TRUE))) # std. deviations with(obarowStory1, (tapply(gnsc1.1, list(gender, MusicT1, PicturesT1), sd, na.rm=TRUE))) # counts with(obarowStory1, (tapply(gnsc1.1, list(gender, MusicT1, PicturesT1), function(x) sum(!is.na(x))))) summary(AnovaModel.1) #11.3 Factorial ANOVA #Create a full factorial model using aov() attach(obarowStory1) model=aov(gnsc1.1~gender*MusicT1*PicturesT1) summary(model) #remove the third-order interaction first: model2 <- update(model,~.-gender:MusicT1:PicturesT1) #evaluate the two model (model vs. model2): anova(model,model2) summary(model2) #Delete the interaction (gender x picture) and compare the models again: model3 <- update(model2,~.- gender:PicturesT1) anova(model2,model3) #Continue to delete unnecessary interactions and main effects model4 <- update(model3,~.- gender:MusicT1) model5 <- update(model4,~.- MusicT1:PicturesT1) model6 <- update(model5,~.- PicturesT1) #model7 only leaves "gender" as a variable: model7 <- update(model6,~.- MusicT1) #Then compare it with the null model: model8=aov(gnsc1.1~1) #Model7 will be selected because there is a statistical difference from the null model: anova(model7,model8) summary(model7) #The regression output summary.lm(model7) library(bootStepAIC) boot.stepAIC(model,data=obarowStory1) #Create the model selected by boot.stepAIC: model9=aov(gnsc1.1~gender+MusicT1+gender:MusicT1) summary.lm(model9) summary(model9) anova(model7,model9) #Examine ANOVA (regression) assumptions with diagnostic plots: par(mfrow=c(2,2)) plot(model7) detach(obarowStory1) #If you want to leave out the influential points (outliers) from the data when you create a full factorial model: model.adj=aov(gnsc1.1[-c(3,14,48)]~gender[-c(3,14,48)],data=obarow) summary.lm(model.adj) #11.4 Performing comparisons in a factorial ANOVA #Import "writing.txt" as Writing Writing <- read.table(file.choose(), header=TRUE, sep="", na.strings="NA", dec=".", strip.white=TRUE) #barplot with(Writing, barplot(table(L1, condition), beside=TRUE)) attach(Writing) tapply(score,list(L1,condition),mean) write=lm(score~L1*condition) Anova(write) summary(write) #Check how many levels are there in L1 --> Arabic is the baseline levels(L1) #Pairwise t-test pairwise.t.test(score,L1, p.adjust.method="fdr") pairwise.t.test(score,condition, p.adjust.method="fdr") pairwise.t.test(score,L1:condition, p.adjust.method="fdr") pairwise.t.test(score,L1:condition, p.adjust.method="none") pvalue<-c(7.4e-05,7.4e-11,2e-16,.016,2.4e-07,.00599) sorted.pvalue<-sort(pvalue) j.alpha<-(1:6)*(.05/6) dif<-sorted.pvalue-j.alpha neg.dif<-dif[dif<0] pos.dif<-neg.dif[length(neg.dif)] index<-dif==pos.dif p.cutoff<-sorted.pvalue[index] p.cutoff #11.6 Robust ANOVA # first: install dependent packages install.packages(c("MASS", "akima", "robustbase")) # second: install suggested packages install.packages(c("cobs", "robust", "mgcv", "scatterplot3d", "quantreg", "rrcov", "lars", "pwr", "trimcluster", "mc2d", "psych", "Rfit")) # third: install an additional package which provides some C functions install.packages("devtools") library("devtools") install_github("mrxiaohe/WRScpp") # fourth: install WRS install_github("nicebread/WRS", subdir="pkg") grp<-c(2,3,5,8,4,1,6,7) levels(obarowStory1$Treatment1) #find out exactly what names to use to subset O2way=list() O2way[[1]]=subset(obarowStory1, subset=Treatment1=="No music No pics", select=c(gnsc1.1)) O2way[[2]]=subset(obarowStory1, subset=Treatment1=="No music Yes pics", select=c(gnsc1.1)) O2way[[3]]=subset(obarowStory1, subset=Treatment1=="Yes music No pics", select=c(gnsc1.1)) O2way[[4]]=subset(obarowStory1, subset=Treatment1=="Yes music Yes pics", select=c(gnsc1.1)) O2way library(WRS) t2way(2,2,O2way,tr=.2,grp=c(1,2,3,4)) #the first two numbers indicate the number of levels in the first variable and the second variable, respectively pbad2way(2,2,O2way,est=mom, conall=T, grp=c(1,2,3,4)) obarowStory1Male <- subset(obarowStory1, subset=gender=="male") obarowStory1Female <- subset(obarowStory1, subset=gender=="female") O3way=list() O3way[[1]]=subset(obarowStory1Male, subset=Treatment1=="No music No pics", select=c(gnsc1.1)) O3way[[2]]=subset(obarowStory1Female, subset=Treatment1=="No music No pics", select=c(gnsc1.1)) O3way[[3]]=subset(obarowStory1Male, subset=Treatment1=="No music Yes pics", select=c(gnsc1.1)) O3way[[4]]=subset(obarowStory1Female, subset=Treatment1=="No music Yes pics", select=c(gnsc1.1)) O3way[[5]]=subset(obarowStory1Male, subset=Treatment1=="Yes music No pics", select=c(gnsc1.1)) O3way[[6]]=subset(obarowStory1Female, subset=Treatment1=="Yes music No pics", select=c(gnsc1.1)) O3way[[7]]=subset(obarowStory1Male, subset=Treatment1=="Yes music Yes pics", select=c(gnsc1.1)) O3way[[8]]=subset(obarowStory1Female, subset=Treatment1=="Yes music Yes pics", select=c(gnsc1.1)) t3way(2,2,2, O3way, tr=.2) O3way[[1]]=c(3,3,1,3,0,1,3,1,0,1,2,1) O3way[[2]]=c(2,1,1,1,1) O3way[[3]]=c(6,-1,0,3,2,3,1,0,2) O3way[[4]]=c(1,0,0,1,2,-2,1,2) O3way[[5]]=c(4,1,4,2,0) O3way[[6]]=c(2,2,1,1,-1,-1,0,-2,-1,2) O3way[[7]]=c(1,5,-1,0,2,0,1,-2) O3way[[8]]=c(1,2,1,0,1,-1,1) t3way(2,2,2,O3way)