# full model: you may want to consider various reduced models */
> hardy.fit <- aov(leak ~ potato * regime * temp, hardy)

Call:
   aov(formula = leak ~ potato * regime * temp, data = hardy)

Terms:
                  potato   regime     temp potato:regime potato:temp 
 Sum of Squares 58.24120 31.70007 27.93869      30.41464     0.01684
Deg. of Freedom        1        1        1             1           1
                regime:temp potato:regime:temp Residuals 
 Sum of Squares     2.65151            0.91904  73.73749
Deg. of Freedom           1                  1        67

Residual standard error: 1.049076 
Estimated effects may be unbalanced

# Type I sums of squares
> summary(hardy.fit)

                   Df Sum of Sq  Mean Sq  F Value     Pr(F) 
potato              1  58.24120 58.24120 52.91963 0.0000000
regime              1  31.70007 31.70007 28.80360 0.0000011
temp                1  27.93869 27.93869 25.38590 0.0000038
potato:regime       1  30.41464 30.41464 27.63562 0.0000016
potato:temp         1   0.01684  0.01684  0.01530 0.9019128
regime:temp         1   2.65151  2.65151  2.40924 0.1253323
potato:regime:temp  1   0.91904  0.91904  0.83507 0.3640873
Residuals          67  73.73749  1.10056                   

# Type III sums of squares
> drop1(hardy.fit,formula(hardy.fit))

Single term deletions

Model:
leak ~ potato * regime * temp
                   Df Sum of Sq      RSS  F Value     Pr(F) 
            <none>               73.7375                   
            potato  1  32.45736 106.1948 29.49168 0.0000008
            regime  1  24.83339  98.5709 22.56433 0.0000111
              temp  1  19.23303  92.9705 17.47569 0.0000863
     potato:regime  1  30.22579 103.9633 27.46402 0.0000018
       potato:temp  1   0.32667  74.0642  0.29682 0.5876926
       regime:temp  1   3.00040  76.7379  2.72625 0.1033910
potato:regime:temp  1   0.91904  74.6565  0.83507 0.3640873

# Least Squares means
> hardy.lsm <- lsmean(hardy.fit)

      potato regime temp     pred        se 
1.0.1      1      0    1 2.208690 0.4691609
1.0.2      1      0    2 2.489456 0.4691609
1.1.1      1      1    1 1.664653 0.3028420
1.1.2      1      1    2 3.290713 0.2909612
2.0.1      2      0    1 4.624648 0.2909612
2.0.2      2      0    2 5.670452 0.2909612
2.1.1      2      1    1 1.811254 0.3965133
2.1.2      2      1    2 3.243771 0.3965133

> lsmean(hardy.fit,factors=c("potato","regime"))

    potato regime     pred        se 
1.0      1      0 2.349073 0.3317468
1.1      1      1 2.477683 0.2099832
2.0      2      0 5.147550 0.2057406
2.1      2      1 2.527513 0.2803772