tmpar <- par( mar = c(3.1,3.1,0,0), pty = "s" )
on.exit( par( tmpar ) )

mplot( fitted( feed.ancova ), resid( feed.ancova ), feed$trt,
   yaxt = "n", xlab = "", ylab = "" )
axis( 2, seq( -20, 20, by = 10 ) )
axis( 2, 30, label = F )
mtext( "fitted", 1, 2 )
mtext( "residual", 2, 2 )
abline( h = 0, lty = 2 )
abline( h = c(-1,1) * std.dev( feed.ancova ), lty = 3 )