tmpar <- par( mar = c(3,3.1,0,0.4), mfrow = c(1,2) )
on.exit( { par( tmpar ) ; rm( tmpar ) } )

mplot( wasp$G1Wa, wasp$HW, wasp$caste, xlab ="", ylab = "" )
mlines( wasp$G1Wa, predict( wasp.cov ), wasp$caste, lty = 2 )
mlines( wasp$G1Wa, predict( wasp.par ), wasp$caste, lty = 1 )
mtext("head width (HW)",2,2)
mtext("(a) gonadium width (G1Wa)",1,2)

tmpdata <- wasp[,c("Hres","Gres","G1Wa","caste","TW")]
tmpdata$TW <- rep( mean( tmpdata$TW ), nrow( tmpdata ) )

mplot( wasp$G1Wa, wasp$Hres, wasp$caste, xlab ="", ylab = "" )
mlines( wasp$G1Wa, predict( wasp.covr, tmpdata ), wasp$caste, lty = 1 )
mlines( wasp$G1Wa, predict( wasp.covwr ), wasp$caste, lty = 2 )
mtext("head width residuals on TW",2,2)
mtext("(b) gonadium width (G1Wa)",1,2)