hardy <- read.table("hardy.dat",header=T)
hardy$potato <- factor(hardy$potato)
hardy$regime <- factor(hardy$regime)
hardy$temp <- factor(hardy$temp)

# full model: you may want to consider various reduced models */
hardy.fit <- aov(leak ~ potato * regime * temp, hardy)
# Type I sums of squares
summary(hardy.fit)
# Type III sums of squares
drop1(hardy.fit,formula(hardy.fit))

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