# default plots from S
plot(hardy.fit)

# diagnostic plot of resid vs predict with 0 line and +/- 1 SD
mplot(jitter(predict(hardy.fit)),resid(hardy.fit),group=hardy$code)
abline(h=0,lty=2)
abline(h=c(-1,1)*std.dev(hardy.fit),lty=3)

# interaction plot with LSD bars (but play with this!)
lsd.plot(hardy.fit,hardy,c("potato","regime"))
# the following give a slightly different plot -- why?
attach(hardy)
lsd.plot(potato,leak,regime)
detach()

# side-by-side interaction plots for the subsets of t=1,2
tmpar <- par(mfrow=c(1,2))
lsd.plot(hardy.fit,hardy,c("potato","regime"),hardy.lsm[hardy.lsm$temp==1,])
lsd.plot(hardy.fit,hardy,c("potato","regime"),hardy.lsm[hardy.lsm$temp==2,])
par(tmpar)