prodres <- product
prodres.fit <- lm( y ~ B, product )
prodres$y <- resid( prodres.fit ) + lsmean( prodres.fit, factors=NULL )$pred
prodres.fit  <- aov( y ~ A + B + D + E + A:E, prodres )

tmpoff <- .025
tmpy <- c(3,7)
tmpfun <- function( x ) paste( as.character( unique( x ) ), collapse = "" )

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

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

factors <- c("A","D")
lsm <- lsmean( prodres.fit, prodres, factors )

lsd.plot( prodres.fit, factors = factors, lsm = lsm,
   xlab = "", ylab = "", xpos = 2.1, ypos = 4.8, ylim=c(4.4,6.1) )
mtext("(a) A by D (LH) with E (XYZ)",1,2)

text( codes( lsm[[factors[1]]]), lsm$pred + tmpoff * diff( tmpy ),
   c( tapply( prodres$E, prodres[,rev(factors)], tmpfun ) ) )

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

factors <- c("A","E")
lsm <- lsmean( prodres.fit, prodres, factors)

lsd.plot( prodres.fit, factors = factors, lsm = lsm,
   xlab = "", ylab = "", xpos = 2.5, ypos = 4.5, ylim=c(3,7) )
mtext("(b) A by E (XYZ) with D (LH)",1,2)
#mtext("response",2,2)

tmpoff <- .06
text( codes( lsm[[factors[1]]]), lsm$pred + tmpoff * diff( tmpy ),
   c( tapply( prodres$D, prodres[,rev(factors)], tmpfun ) ) )

par( tmpar )
rm( tmpy, tmpoff, tmpfun, tmpar )