product <- data(product)
product$A <- ordered(product$A)
product$B <- ordered(product$B)
product$C <- ordered(product$C)
product$D <- ordered(product$D)
product$E <- factor(product$E)

product.fit <- aov(y ~ .,product)
product.afit <- aov(y~A+B+D+E,product)
product.ifit <- aov(y~B+A*E,product)
tmpar <- par(mar=c(4,4,1,1),mfrow=c(1,2))

attach(product)

tmp <- resid(lm(y~B+E))+mean(y)
margin.plot(A,D,tmp,margin=F,sd=1.85877/4,rdf=11,
   ylim=c(3,7),xlab="(a) A by D by E (-ox) after B+E",ylab="response")
mn <- tapply(resid(lm(y~B+E))+mean(y),list(A,D),mean)
ma <- tapply(A,list(A,D),unique)
me <- tapply(E,list(A,D),function(x)paste(as.character(unique(x)),collapse=""))
text(as.numeric(levels(A)[ma]),mn+.1,me)

tmp <- resid(lm(y~B))+mean(y)
margin.plot(A,E,tmp,margin=F,sd=1.25531/4,rdf=6,
   ylim=c(3,7),xlab="(b) A by E by D (-/+) after B",ylab="")
mn <- tapply(resid(lm(y~B))+mean(y),list(A,E),mean)
ma <- tapply(A,list(A,E),unique)
md <- tapply(D,list(A,E),function(x)paste(as.character(unique(x)),collapse=""))
text(as.numeric(levels(A)[ma]),mn+.1,md)

detach(product)

par(tmpar)