tmpar <- par( mar=c(3.1,2.4,0,.1), mfcol=c(1,3) )

rlabs <- c("residuals (log10 scale)","","")
xlabs <- c("(a) early period","(b) middle period","(c) late period")
names( xlabs ) <- names( rlabs ) <- levels( forage$per )
ally <- forage.mean + forage.spint
ylims <- range( forage.proj$Within[, "Residuals"] )
tmp <- seq( 14, 22 )

for ( i in levels( forage$per ) ) {
   tmpi <- forage$per==i
   mplot( jitter( ally[tmpi] ), forage.proj$Within[tmpi, "Residuals"],
      xaxt="n", ylim=ylims,
      group=forage[tmpi, "code"], xlab="", ylab="" )
   mtext( xlabs[i], 1, 2 )
   axis( 1, log10( tmp ), tmp )
   if ( i==levels( forage$per )[2] )
      axis( 1, log10( c(19.5,20.5) ), c(19.5,20.5) )
   abline( 0, 0, lty=2 )
   abline( h=c(-1,1) * std.dev( forage.fit$Within ), 0, lty=3 )
}
rm( rlabs, xlabs, ally, ylims, tmp, tmpi )

par( tmpar )