# Data from a split plot experiment with completely randomized whole plot structure # from Roger Kirk 'Experimental Design: Procedures for the Behavioral Sciences' 3rd edition page 516 # Each observation has the block, factor a level, factor b level, then response listed data <- matrix(scan(),ncol=4,byrow=TRUE) 1 1 1 3 1 1 2 4 1 1 3 7 1 1 4 7 2 1 1 6 2 1 2 5 2 1 3 8 2 1 4 8 3 1 1 3 3 1 2 4 3 1 3 7 3 1 4 9 4 1 1 3 4 1 2 3 4 1 3 6 4 1 4 8 5 2 1 1 5 2 2 2 5 2 3 5 5 2 4 10 6 2 1 2 6 2 2 3 6 2 3 6 6 2 4 10 7 2 1 2 7 2 2 4 7 2 3 5 7 2 4 9 8 2 1 2 8 2 2 3 8 2 3 6 8 2 4 11 data KirkCh12SPF <- data.frame(data) colnames(KirkCh12SPF) <- c("block","a","b","y") KirkCh12SPF$block <- as.factor(KirkCh12SPF$block) KirkCh12SPF$a <- as.factor(KirkCh12SPF$a) KirkCh12SPF$b <- as.factor(KirkCh12SPF$b) rm(data) library(afex) library(emmeans) KCh12SPF.afex1 <- mixed(y ~ a + (1|a:block) +b +a:b, data = KirkCh12SPF) summary(KCh12SPF.afex1) nice(KCh12SPF.afex1) lsmeans(KCh12SPF.afex1, ~a) lsmeans(KCh12SPF.afex1, ~b) lsmeans(KCh12SPF.afex1, ~a:b) ## -------------- two different approaches for Tukeys method ------------------- # multiple comparisons with Tukeys method # (note that they are not of interest here due to interaction!!) emmeans(KCh12SPF.afex1, list(pairwise ~ a), adjust = "tukey") emmeans(KCh12SPF.afex1, list(pairwise ~ b), adjust = "tukey") # multiple comparisons with the Tukey method # (note that they are not of interest here due to interaction!!) refa <- lsmeans(KCh12SPF.afex1,specs = c("a")) refb <- lsmeans(KCh12SPF.afex1,specs = c("b")) contrast(refa,method="pairwise",adjust="tukey") contrast(refb,method="pairwise",adjust="tukey") ## ----------------------------------------------------------------------------- # examine possible interactions lsmip(KCh12SPF.afex1, a~b) lsmip(KCh12SPF.afex1, b~a) # examples of simple effect (slice) tests to examine interactions # slice by b refsliceb <- lsmeans(KCh12SPF.afex1, ~ a | b) csourceslb <- contrast(refsliceb, "consec") test(csourceslb,joint=TRUE,by="b") # slice by a refslicea <- lsmeans(KCh12SPF.afex1, ~ b | a) csourcesla <- contrast(refslicea, "consec") test(csourcesla,joint=TRUE,by="a") # checking assumptions # apparently the result from the mixed function in afex does not # recognize the plot function that we have used to check assumptions in other # analyses. I will continue to look for ways to check assumptions, but for now # I have two approaches below: # for a normal plot, use the qqnorm and qqline functions applied as below: qqnorm(summary(KCh12SPF.afex1)$residuals) qqline(summary(KCh12SPF.afex1)$residuals) # I have not found an easy way to get predicted values, so here we can look at the # variability of the original variable across the factor combinations: library(lattice) bwplot(y ~ a | b, data = KirkCh12SPF)