Dear R-experts! I having trouble with the correct implementation of a mixed-model ANOVA in R. I my dataset subjects were tested in a cognitive performance test (numerical outcome variable 'score'). This cognitive performance test is devided into five blocks (categorical factor 'block'). All subjects were tested two times (in random order once following placebo treatment and once following drug treatment: categorical factor 'treatment'). In order to account for re-test effects I coded the testing session (categorical factor 'session' with 1 coding for the first and 2 coding for the second session). Also all subjects once underwent a blood screening before the study in which a stable blood marker was measured which is hypothesized to mediate treatment effects (covariate 'blood'). I would like to look for treatment effects on my outcome variable 'score' and whether the treatment effect is mediated by the covariate 'blood'. Also I would like to use 'subject' as a random factor, include 'block' as a within-subjects factor and control for re-test effects by including 'session'. Could you advice me regarding the proper implementation of my model? Also I am happy about any pointers to R-related literature. Please find below some dummy data and initially model formulation. Many thanks & best wishes, Jokel # some dummy data data <- structure(list(subject = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L), .Label = c("1", "2", "3", "4", "5"), class = "factor"), block = structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("1", "2", "3"), class = "factor"), treatment = structure(c(2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L), .Label = c("drug", "placebo"), class = "factor"), session = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L), .Label = c("1st", "2nd"), class = "factor"), blood = c(3.04, 3.04, 3.04, 3.04, 3.04, 3.04, 4.93, 4.93, 4.93, 4.93, 4.93, 4.93, 5.65, 5.65, 5.65, 5.65, 5.65, 5.65, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.77, 5.77, 5.77, 5.77, 5.77, 5.77), score = c(10.98, 11.66, 9.99, 9.15, 9.9, 10.3, 10.78, 9.43, 8.6, 9.67, 10.75, 9.09, 11.18, 10.33, 8.61, 10.04, 9.13, 8.71, 10.03, 9.38, 11.17, 9.82, 10.07, 8.34, 9.94, 11.38, 9.19, 8.35, 10.17, 9.04)), .Names = c("subject", "block", "treatment", "session", "blood", "score"), row.names = c(NA, -30L), class = "data.frame") # some first try which however throws a warning aov1 <- aov(score~blood*treatment*session*block+Error(subject/(treatment*session*block)), data=data) summary(aov1) [[alternative HTML version deleted]]