[R] Tukey's one degree of freedom for nonadditivity?

kjetil brinchmann halvorsen kjetil at entelnet.bo
Sat Apr 5 16:43:12 CEST 2003


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
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help



More information about the R-help mailing list