The code below is intended to analyse a textbook example of a balanced incomplete block design: # # Data taken from pp. 219-230 in # Cox, D.R. (1958) Planning of Experiments. John Wiley and Son, Inc. New York. 308 pp. # day <- factor(rep(1:10, each = 3)) T <- factor(c("T4","T5","T1","T4","T2","T5","T2","T4","T1","T5", "T3","T1","T3","T4","T5","T2","T3","T1","T3","T1", "T4","T3","T5","T2","T2","T3","T4","T5","T1","T2")) response <- c(4.43,3.16,1.40,5.09,1.81,4.54,3.91,6.02,3.32,4.66, 3.09,3.56,3.66,2.81,4.66,1.60,2.13,1.31,4.26,3.86, 5.87,2.57,3.06,3.45,3.31,5.10,5.42,5.53,4.46,3.94) incomplt.blk.modelaov <- aov(response ~ T + Error(day)) summary(incomplt.blk.modelaov) model.tables(incomplt.blk.modelaov, type = "means", se = TRUE) It gives the correct anova, but the attempt to extract the treatment means using function model.tables() produces the following error: Error in FUN(X[[1]], ...) : eff.aovlist: non-orthogonal contrasts would give an incorrect answer Is there a way to obtain these means, either with or without recovery of inter-block information? (N.B. The means with recovery of inter-block information can be obtained from lme() after loading package nlme, but I'd like to do this using an anova approach rather than a mixed modelling approach if possible.) Thank you in anticipation. Nick Galwey