> bacdata <- read.table("http://www.stat.ufl.edu/~winner/data/breathalyzer.dat",header=F,
+                       col.names=c("machine","subject","repnum","bac"))
> 
> attach(bacdata)
The following object(s) are masked from 'bacdata (position 3)':

    bac, machine, repnum, subject
The following object(s) are masked from 'bacdata (position 8)':

    bac, machine, repnum, subject
The following object(s) are masked from 'bacdata (position 9)':

    bac, machine, repnum, subject
The following object(s) are masked from 'bacdata (position 10)':

    bac, machine, repnum, subject
> 
> machine <- factor(machine)
> subject <- factor(subject)
> 
> bac.aov1 <- aov(bac ~ machine*subject)
> anova(bac.aov1)
Analysis of Variance Table

Response: bac
                 Df    Sum Sq   Mean Sq   F value    Pr(>F)    
machine           5 0.0019006 0.0003801   68.5430 < 2.2e-16 ***
subject           2 0.0219044 0.0109522 1974.9028 < 2.2e-16 ***
machine:subject  10 0.0003083 0.0000308    5.5595 4.177e-07 ***
Residuals       162 0.0008984 0.0000055                        
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> 
> # install.packages("nlme")
> library(nlme)
> 
> ms <- machine:subject
> 
> bac.aov2 <- lme(fixed = bac~machine, random = ~1|subject/ms)
> 
> summary(bac.aov2)
Linear mixed-effects model fit by REML
 Data: NULL 
        AIC       BIC   logLik
  -1541.304 -1512.872 779.6519

Random effects:
 Formula: ~1 | subject
        (Intercept)
StdDev:  0.01349148

 Formula: ~1 | ms %in% subject
        (Intercept)    Residual
StdDev: 0.001590143 0.002354929

Fixed effects: bac ~ machine 
                  Value   Std.Error  DF   t-value p-value
(Intercept)  0.07834000 0.007855002 162  9.973263  0.0000
machine2    -0.00479333 0.001433672  10 -3.343396  0.0074
machine3    -0.00209333 0.001433672  10 -1.460120  0.1749
machine4    -0.00998667 0.001433672  10 -6.965796  0.0000
machine5    -0.00239667 0.001433672  10 -1.671698  0.1255
machine6    -0.00625000 0.001433672  10 -4.359435  0.0014
 Correlation: 
         (Intr) machn2 machn3 machn4 machn5
machine2 -0.091                            
machine3 -0.091  0.500                     
machine4 -0.091  0.500  0.500              
machine5 -0.091  0.500  0.500  0.500       
machine6 -0.091  0.500  0.500  0.500  0.500

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.70630219 -0.59751478  0.06545501  0.53365314  3.53592318 

Number of Observations: 180
Number of Groups: 
        subject ms %in% subject 
              3              18 
> 
> library(lme4)
> 
> bac.aov3 <- lmer(bac~machine+(1|subject)+(1|ms))
> 
> summary(bac.aov3)
Linear mixed model fit by REML 
Formula: bac ~ machine + (1 | subject) + (1 | ms) 
   AIC   BIC logLik deviance REMLdev
 -1541 -1513  779.7    -1625   -1559
Random effects:
 Groups   Name        Variance   Std.Dev. 
 ms       (Intercept) 2.5287e-06 0.0015902
 subject  (Intercept) 1.8203e-04 0.0134917
 Residual             5.5457e-06 0.0023549
Number of obs: 180, groups: ms, 18; subject, 3

Fixed effects:
             Estimate Std. Error t value
(Intercept)  0.078340   0.007853   9.975
machine2    -0.004793   0.001434  -3.343
machine3    -0.002093   0.001434  -1.460
machine4    -0.009987   0.001434  -6.966
machine5    -0.002397   0.001434  -1.672
machine6    -0.006250   0.001434  -4.359

Correlation of Fixed Effects:
         (Intr) machn2 machn3 machn4 machn5
machine2 -0.091                            
machine3 -0.091  0.500                     
machine4 -0.091  0.500  0.500              
machine5 -0.091  0.500  0.500  0.500       
machine6 -0.091  0.500  0.500  0.500  0.500
>