TukeyTest = function(y,A,B){
yNAME = deparse(substitute(y))
ANAME = deparse(substitute(A))
BNAME = deparse(substitute(B))
I = nlevels(A); J = nlevels(B)
a = lm(y~A+B,contrasts=list(A="contr.treatment",B="contr.treatment"))
coefA = rep(c(0,coef(a)[2:I]),each=J)
coefB = rep(c(0,coef(a)[-(1:I)]),I)
a = update(a,.~.+I(coefA*coefB))
aa = anova(a)
rownames(aa)[1:3] = c(ANAME,BNAME,"Interaction")
attr(aa,"heading") = c("Tukey's test for nonadditivity",paste("Response",yNAME))
aa
}