data <- matrix(scan(),ncol=3,byrow=TRUE)
   1 2 40.89  1 5 37.99  1 4 37.18  1 1 34.98  1 6 34.89  1 3 42.07 
   2 1 41.22  2 3 49.42  2 4 45.85  2 6 50.15  2 5 41.99  2 2 46.69 
   3 6 44.57  3 3 52.68  3 5 37.61  3 1 36.94  3 2 46.65  3 4 40.23 
   4 2 41.90  4 4 39.2   4 6 43.29  4 5 40.45  4 3 42.91  4 1 39.97 

data
disp81            <-  data.frame(data)
colnames(disp81)  <-  c("block","treatment","nit")
rm(data)

disp81$block      <-  as.factor(disp81$block)
disp81$treatment      <-  as.factor(disp81$treatment)

 # blocks as fixed effects

rcb81  <-  lm(nit ~ block + treatment, data = disp81)
anova(rcb81)

 # check assumptions 
par(mfrow=c(2,2))  
plot(rcb81)
par(mfrow=c(1,1))  





 # blocks as random effects

library(nlme)

  # regular parametrization for ANOVA result, but not means and std errors

rcb81aov  <- lme(nit ~ treatment, data = disp81, random = ~ 1 | block)
anova(rcb81aov)

  # means model parametrization for means and std errors, but ANOVA is incorrect

rcb81stderr  <- lme(nit ~ treatment -1, data = disp81, random = ~ 1 | block)
summary(rcb81stderr)