y <- c(9.3,  9.2,  9.7,  9.4,  9.3,  9.4,  9.8,  9.5, 10.0, 10.0,  9.9, 10.2)
 tips <- as.factor(c(1,3,4, 1,2,3, 2,3,4, 1,2,4))
 coupons <- as.factor(rep(1:4, each=3))
  
 treatment.totals <- vector(length=4)
 for(i in 1:4) treatment.totals[i] <- sum(y[tips==i])
 treatment.totals
  
 block.totals <- vector(length=4)
 for(j in 1:4) block.totals[j] <- sum(y[coupons==j])
 block.totals
  
 g <- lm(y ~ coupons + tips)
 anova(g)

# For a BIBD, ANOVA depends on the order of the treatment and block factors entering the model.
# In the above, the treatment effects are adjusted for the block effects
# The following ANOVA table is NOT right.  


 h <- lm(y~ tips + coupons)
 anova(h)