[R] tukey.1

Raphael Fraser raphael.fraser at gmail.com
Tue Nov 9 19:49:36 CET 2010


aov.pen = aov(y ~ blend + treatment)
summary(aov.pen)


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))
resp=data[,as.character(attr(aov.obj$terms,
     "variables")[attr(aov.obj$terms, "response"
     )])]
cfs=coef(aov.obj)
alpha.A=aov.obj$contrasts[[vnames[1]]] %*% cfs[
     aov.obj$assign[[vnames[1]]]]
alpha.B=aov.obj$contrasts[[vnames[2]]] %*% cfs[
     aov.obj$assign[[vnames[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)
}




On Tue, Nov 9, 2010 at 10:22 AM, Sarah Goslee <sarah.goslee at gmail.com> wrote:
> Including pen.df is a good start, but
> we likely can't help without more info:
> Where'd you get tukey.1.r and what's in it?
> What's aov.pen?
>
> On Tue, Nov 9, 2010 at 10:13 AM, Raphael Fraser
> <raphael.fraser at gmail.com> wrote:
>> I have been trying to do tukey's test to no avail. This is what I have
>> and the error produced:
>>
>> pen.df = data.frame(blend, treatment, y)
>> source("tukey.1.r")
>> tukey.1(aov.pen, pen.df)
>> Error in tukey.1(aov.pen, pen.df) : the model must be two-way
>>
>> This is a two-way design. Therefore I am confused. Can anyone help?
>>
>>
>> Raphael
>>
>>> pen.df
>>   blend treatment  y
>> 1      1         A 89
>> 2      2         A 84
>> 3      3         A 81
>> 4      4         A 87
>> 5      5         A 79
>> 6      1         B 88
>> 7      2         B 77
>> 8      3         B 87
>> 9      4         B 92
>> 10     5         B 81
>> 11     1         C 97
>> 12     2         C 92
>> 13     3         C 87
>> 14     4         C 89
>> 15     5         C 80
>> 16     1         D 94
>> 17     2         D 79
>> 18     3         D 85
>> 19     4         D 84
>> 20     5         D 88
>>
>
>
>
> --
> Sarah Goslee
> http://www.functionaldiversity.org
>



More information about the R-help mailing list