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.gos...@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.fra...@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
>
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
