# Exercise 5.1 library(phia) library(multcomp) Boik obj <- aov(edr ~ medication, data = Boik) summary(obj) pairwise <- glht(obj, linfct = mcp(medication = "Tukey")) summary(pairwise) ctr <- rbind(complex = c(-1, 1/3, 1/3, 1/3)) custom <- glht(obj, linfct = mcp(medication = ctr)) summary(custom) dunnett <- glht(obj, linfct = mcp(medication = "Dunnett")) boxplot(edr ~ medication, data = Boik) describe(Boik) library(QuantPsyc) plotNormX(Boik$edr) qqnorm(Boik$edr) qqline(Boik$edr) library(Rcmdr) leveneTest(edr ~ medication, data = Boik, center = mean) bartlett.test(edr ~ medication, data = Boik) # Exercise 5.2 obj2 <- aov(edr ~ medication*therapy, data = Boik) summary(obj2) library(car) Anova(obj2, type = 2) Anova(obj2, type = 3) library(phia) testInteractions(obj2, fixed = "medication", across = "therapy", adjustment = "bonferroni") testInteractions(obj2, fixed = "therapy", across = "medication", adjustment = "bonferroni") testInteractions(obj2, fixed = "medication", pairwise = "therapy", adjustment = "bonferroni") testInteractions(obj2, fixed = "therapy", pairwise = "medication", adjustment = "bonferroni") ctr <- cbind(c(-1, 1/3, 1/3, 1/3)) total.ctr <- list(medication = ctr) testInteractions(obj2, custom = total.ctr, fixed = "therapy") aggregate(edr ~ medication*therapy, data = Boik, FUN = mean) aggregate(edr ~ medication*therapy, data = Boik, FUN = sd) plot(interactionMeans(obj2)) # Assignment 5 x11 <- c(51, 54, 45, 55, 52, 64, 45, 47, 64, 59) x12 <- c(66, 79, 74, 79, 70, 71, 67, 63, 70, 79) x13 <- c(62, 67, 75, 72, 65, 66, 73, 60, 68, 66) x21 <- c(44, 35, 41, 37, 39, 31, 43, 46, 37, 41) x22 <- c(52, 50, 64, 52, 47, 60, 49, 60, 63, 52) x23 <- c(72, 72, 72, 65, 73, 64, 60, 63, 71, 69) x31 <- c(46, 41, 36, 45, 39, 31, 41, 49, 31, 44) x32 <- c(35, 44, 43, 34, 45, 49, 44, 49, 40, 46) x33 <- c(45, 49, 40, 34, 45, 44, 40, 38, 34, 31) part <- factor(rep(c(6, 9, 12), each = 30)) hrs <- factor(rep(rep(c(10, 15, 20), each = 10), 3)) score <- c(x11, x12, x13, x21, x22, x23, x31, x32, x33) dat5 <- data.frame(part, hrs, score) oneway1 <- aov(score ~ part, data = dat5) summary(oneway1) library(multcomp) pairwise <- glht(oneway1, linfct = mcp(part = "Tukey")) summary(pairwise) ctr <- t(poly(c(6, 9, 12), 2)) custom1 <- glht(oneway1, linfct = mcp(part = ctr)) summary(custom1) oneway2 <- aov(score ~ hrs, data = dat5) summary(oneway2) custom2 <- glht(oneway2, linfct = mcp(hrs = ctr)) summary(custom2) ctr3 <- rbind(c(-1, 0.5, 0.5), c(0, -1, 1)) custom3 <- glht(oneway1, linfct = mcp(part = ctr3)) summary(custom3, test = adjusted(type = "bonferroni")) twoway <- aov(score ~ part*hrs, data = dat5) library(car) Anova(twoway, type = 3) aggregate(score ~ part*hrs, data = dat5, FUN = mean) aggregate(score ~ part*hrs, data = dat5, FUN = sd) library(phia) plot(interactionMeans(twoway)) testInteractions(twoway, fixed = "hrs", across = "part", adjustment = "bonferroni") testInteractions(twoway, fixed = "part", across = "hrs", adjustment = "bonferroni") testInteractions(twoway, fixed = "hrs", pairwise = "part", adjustment = "bonferroni") testInteractions(twoway, fixed = "part", pairwise = "hrs", adjustment = "bonferroni") ctr4 <- t(ctr) total.ctr1 <- list(part = ctr4) testInteractions(twoway, fixed = "hrs", custom = total.ctr1) total.ctr2 <- list(hrs = ctr4) testInteractions(twoway, fixed = "part", custom = total.ctr2)