Tukey's one degree of freedom for nonadditivity?

On 3 Apr 2003 at 12:11, Spencer Graves wrote:

See the code copied below.

Kjetil Halvorsen

# From S-Plus:  Guide to Statistical & Mathemathical
# Analysis, page 11.25       slightly changed to function with R.
tukey.1  <-  function(aov.obj, data) {
	vnames <- names(aov.obj$contrasts)
	if(length(vnames) != 2)
		stop("The model must be two-way.")
	vara  <-  data[, vnames[1]]
	varb  <-  data[, vnames[2]]
	na  <-  length(levels(vara))
	nb  <-  length(levels(varb))
        where.resp <-  as.character(attr(aov.obj$terms,
	          "variables")[-1][attr(aov.obj$terms, "response" )])
	resp  <-  data[, where.resp]
	cfs  <-  coef(aov.obj)
	alpha.A  <-  aov.obj$contrasts[[vnames[1]]] %*% 
		     cfs[aov.obj$assign[aov.obj$assign==1]] 
	alpha.B  <-  aov.obj$contrasts[[vnames[2]]] %*% 
		     cfs[aov.obj$assign[aov.obj$assign==2]] 
	r.mat  <-  matrix(0, nb, na)
	r.mat[cbind(as.vector(unclass(varb)), as.vector(
		unclass(vara)))]  <-  resp
	SS.theta.num  <-  sum((alpha.B %*% t(alpha.A)) * 
		r.mat)^2
	SS.theta.den  <-  sum(alpha.A^2) * sum(alpha.B^2)
	SS.theta  <-  SS.theta.num / SS.theta.den
	SS.res    <-  sum(resid(aov.obj)^2)
	SS.res.1  <-  SS.res - SS.theta
	T.1df     <-  ((na * nb -na -nb) *  SS.theta)/SS.res.1
	p.value   <-  1 - pf(T.1df, 1, na*nb - na -nb)
	list(T.1df = T.1df, p.value = p.value)
	}

  Is there code available to decompose interactions involving at least 
one nominal factor with more than 2 levels as described, e.g., by Tukey 
or by Mandel (1971, Technometrics, 13: 1-18)?

  Tukey's model:

E(y[i,j]) = mu0 + a[i] + b[j] + c*a[i]*b[j],

estimating a, b, and c so sum(a) = sum(b)= 0.  Mandel essentially 
describes a singular value decomposition of the interaction.

Thanks,
Spencer Graves