[R] Means from balanced incomplete block design

[N.W.Galwey]

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