data <- matrix(scan(),ncol=3,byrow=TRUE)
 1 1 42.5  1 1 43.3  1 1 42.9  1 2 42.2  1 2 41.4  1 2 41.8 
 2 3 48.0  2 3 44.6  2 3 43.7  2 4 42.0  2 4 42.8  2 4 42.8 
 3 5 41.7  3 5 43.4  3 5 42.5  3 6 40.6  3 6 41.8  3 6 41.8 

data
table77  <- data.frame(data)
colnames(table77)  <- c("day","run","glucose")
rm(data)

table77$day  <- as.factor(table77$day)
table77$run  <- as.factor(table77$run)

 # Gives SS and MS, but not easy to use

summary(aov(glucose ~ day + Error(day/run), data = table77))

 # or instead use this, which looks better, but has the 
 #  wrong test for day :

summary(aov(glucose ~ day + run %in% day, data = table77))


 # the lme4 package estimates variance components and leads 
 # to useful diagnostics:

library(lme4)

table77.lme1  <- lmer(glucose ~ (1|day) + (1|run:day), data = table77)

summary(table77.lme1)

plot(table77.lme1)
windows()                                 # use a different command in Mac or Linux
qqnorm(summary(table77.lme1)$residuals)
qqline(summary(table77.lme1)$residuals)


 # another way to fit the same lmer model:

table77.lme1a  <- lmer(glucose ~ (1|day/run), data = table77)

summary(table77.lme1a)

  

 # the afex package provides fixed effects tests as well as variance 
 # component estimates, but does not seem to give diagnostic information 
 # as easily

library(afex)

table77.afex1 <- mixed(glucose ~ day + (1|run:day), data = table77)

summary(table77.afex1)

nice(table77.afex1)

lsmeans(table77.afex1, ~day)