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

ylims <- range( c( resid( forage.fit$cow ), resid( forage.fit$Within ) ) )

plot( forage.mean + fitted( forage.fit$cow ), resid( forage.fit$cow ),
   ylim=ylims, xlab="", ylab="", xaxt="n", pch = "o" )
mtext( "(a) cow predicteds", 1, 2 )
mtext( "cow residuals", 2, 2 )
axis( 1, log10( c(10,15,20,25) ), c(10,15,20,25) )
abline( h=0, lty=2 )
abline( h=c(-1,1) * std.dev( forage.fit$cow ), lty=3 )

plot( forage.mean + fitted( forage.fit$Within ), resid( forage.fit$Within ),
   ylim=ylims, xlab="", ylab="", xaxt="n", pch = "o" )
mtext( "(b) period predicteds", 1, 2 )
mtext( "period residuals", 2, 2 )
axis( 1, log10( 16:21 ), 16:21 )
abline( h=0, lty=2 )
abline( h=c(-1,1) * std.dev( forage.fit$Within ), lty=3 )

rm( ylims )

par( tmpar )