tmpar <- par(mfrow=c(1,2),mar=c(4.1,4.1,0,0))

damage.m <- aov(y~level,damage)

attach(damage)

mplot(x,y,group=as.character(score),
	xlab="(a) ANOVA fit",ylab="")
for (i in 1:4){
	lines(range(x[score==i]),rep(mean(y[score==i]),2))
}
mplot(x,resid(damage.m),group=as.character(score),
	xlab="(b) displacement",ylab="")
mlines(x,predict(damage.xscore)-predict(damage.m),group=score,lty=1)
abline(h=0,lty=2)
abline(h=c(-1,1)*std.dev(damage.m),lty=3)
se.bar(0.5,-4.75,std.dev(damage.xscore),cap="SD")

detach()

par(tmpar)

#mplot(x,y,group=as.character(score),
#	xlab="(a) regression on logit",ylab="yield")
#tmp <- order(x)
#lines(x[tmp],predict(damage.x)[tmp],lty=2)
#se.bar(0.6,93,std.dev(damage.x),cap="SD")