library(psych) library(car) library(phia) library(multcomp) ### ANCOVA sat.act$education <- factor(sat.act$education, labels = paste0("L", 0:5)) sat.act$gender <- factor(sat.act$gender, labels = c("male", "female")) ancova1 <- aov(ACT ~ SATV + SATQ + education, data = sat.act) ancova2 <- aov(ACT ~ education + SATV + SATQ, data = sat.act) summary(ancova1) # Type 1 SS summary(ancova2) # Type 1 SS Anova(ancova1, type = 3) Anova(ancova2, type = 3) sat.act$SATVC <- scale(sat.act$SATV, scale = FALSE) sat.act$SATQC <- scale(sat.act$SATQ, scale = FALSE) ancova1 <- aov(ACT ~ SATVC + SATQC + education, data = sat.act) ancova2 <- aov(ACT ~ education + SATV + SATQ, data = sat.act) Anova(ancova1, type = 3) Anova(ancova2, type = 3) #Unadjusted means vs. Adjusted means aggregate(ACT ~ education, data = sat.act, FUN = mean) interactionMeans(ancova1, factors = "education") # Pairwise comparison in ANCOVA pairwisec <- glht(ancova1, linfct = mcp(education = "Tukey")) summary(pairwisec, test = adjusted(type = "bonferroni")) confint(pairwisec) plot(confint(pairwisec)) # A priori contrast in ANCOVA linear <- c(-5, -3, -1, 1, 3, 5) quadratic <- c(5, -1, -4, -4, -1, 5) # Linear change ctr <- rbind(linear) linearc <- glht(ancova1, linfct = mcp(education = ctr)) summary(linearc) # Both linear and quadratic changes mctr <- rbind(linear, quadratic) polynomialc <- glht(ancova1, linfct = mcp(education = mctr)) summary(polynomialc, test = adjusted(type="bonferroni")) # Pairwise comparison testInteractions(ancova1, pairwise = "education") # Both linear and quadratic changes ctr <- cbind(linear, quadratic) total.ctr <- list(education = ctr) testInteractions(ancova1, custom = total.ctr) ######## Two-way ANCOVA twoancova <- aov(ACT ~ SATV + SATQ + education*gender, data = sat.act) summary(twoancova) Anova(twoancova, type = 2) Anova(twoancova, type = 3) # Visualize interaction based on adjusted means interactionMeans(twoancova) plot(interactionMeans(twoancova)) # Simple main effect testInteractions(twoancova, fixed = "education", across = "gender") testInteractions(twoancova, fixed = "gender", across = "education") testInteractions(twoancova, pairwise = "education", fixed = "gender") testInteractions(twoancova, custom = total.ctr, fixed = "gender") ################################ ## ANOVA for repeated-measures designs head(OBrienKaiser) # Traditional repeated-measures ANOVA # 1. Rearrange data to long format fac <- cbind(id = 1:16, OBrienKaiser[,1:2]) newdata <- rbind(fac, fac, fac) timefac <- rep(c("pre", "post", "fup"), each = 16) score <- c(OBrienKaiser$pre.1, OBrienKaiser$post.1, OBrienKaiser$fup.1) newdata <- data.frame(newdata, time = timefac, score = score) # Repeated-measures ANOVA ranova <- aov(score ~ time + Error(id/(time)), data = newdata) summary(ranova) # Traditional mixed-design ANOVA mixanova <- aov(score ~ treatment*gender*time + Error(id/(time)), data = newdata) summary(mixanova) # Repeated-measures ANOVA library(ez) ezANOVA(data = newdata, dv = score, wid = id, within = .(time)) ezANOVA(data = newdata, dv = score, wid = id, within = .(time), between = c("treatment", "gender")) ### MANOVA approach for repeated-measures design # Require wide-format data # manova function is not designed for repeated-measures design obj1 <- manova(cbind(pre.1, post.1, fup.1) ~ treatment * gender, data = OBrienKaiser) summary(obj1, type = "Wilks") # Use lm function with idata to obtain MANOVA for repeated-measures design mod <- lm(cbind(pre.1, post.1, fup.1) ~ treatment * gender, data = OBrienKaiser) idata <- data.frame(phase = c("pre", "post", "fup")) Anova(mod, idata = idata, idesign = ~ phase, type = 3) # NOTE: for more than one within factor # cbind(y1, y2, y3, y4) # idata <- data.frame(gender = c("male", "female", "male", "female"), university = c("cu", "cu", "tu", "tu")) testInteractions(mod, across = c("phase", "treatment"), fixed = "gender", idata = idata) # Simple interaction effect testing interaction between phase and treatment within each gender testInteractions(mod, across = c("phase"), fixed = c("treatment", "gender"), idata = idata) # test simple main effect of phase within each combination of treatment and gender testInteractions(mod, pairwise = c("phase"), fixed = c("treatment", "gender"), idata = idata)