tmpar <- par( mar = c(3.1,4.1,0,0), mfrow = c(1,2) )
on.exit( { 
   par( tmpar )
   rm( tmplevel, tmpx, tmpg, tmpy, tmpw, ylims, tmp )
} )

# Cow Averages over Periods

tmp <- forage
tmp$ally <- forage.mean + forage.wpint
lsm <- list( se = std.dev( forage.fit$cow ) / sqrt( sample.size( forage.fit$cow,
   forage, c("trt","hc") ) ) )

attach( tmp )
ci.plot( trt, ally, hc, ci = ci.width( forage.fit$cow, lsm = lsm ),
   yaxt = "n", xlab = "", ylab = "" )
detach( )

mtext( "(a) forage treatment", 1, 2 )
mtext( "dry matter intake (dmi)", 2, 2 )
tmp <- seq( 14, 22, 2 )
axis( 2, log10( tmp ), tmp )
abline( h = forage.mean, lty = 3 )

#############################################################################

mplot( jitter( forage.mean+forage.wpint ),
   jitter( forage.proj$cow[, "Residuals"] ), group = forage$code,
   xlab = "", ylab = "", xaxt = "n" )
mtext( "(b) predicted", 1, 2 )
mtext( "residuals (log10 scale)", 2, 2 )
axis( 1, log10( tmp ), tmp )
axis( 1, log10( tmp + 1 ), F )
abline( h = 0, lty = 2 )
abline( h = c(-1,1) * std.dev( forage.fit$cow ), lty = 3 )