[R] multcomp two-way anova with interactions within and between

Schreiber, Stefan Stefan.Schreiber at ales.ualberta.ca
Thu Jan 12 01:20:11 CET 2012


Hi all,

I'd like to compare all levels of my interaction with each other. I read
the pdf 'Additional multcomp Examples' but even though there is an
example with an interaction it doesn't work for me when I want to
compare within and between groups.

Here is an example:

####
d.fr<-data.frame(id=rep(1:16,3),treat1=rep(as.factor(LETTERS[1:3]),each=
16),treat2=rep(as.factor(letters[4:7]),each=4),response=rnorm(48))

require(multcomp)
require(lme4)

fit1<-lmer(response~treat1*treat2+(1|id),data=d.fr)

temp<-expand.grid(treat1=unique(d.fr$treat1),treat2=unique(d.fr$treat2))
X<-model.matrix(~treat1*treat2,data=temp)

# X gives me a matrix with the dimensions 12 by 12.

glht(fit1, linfct = X)

Tukey<-contrMat(table(d.fr$treat1),'Tukey')

K1<-cbind(Tukey,matrix(0,nrow=nrow(Tukey),ncol=ncol(Tukey)))
rownames(K1) <- paste(levels(d.fr$treat2)[1],rownames(K1), sep = ":")

K2 <- cbind(matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)), Tukey)
rownames(K2) <- paste(levels(d.fr$treat2)[2],rownames(K1), sep = ":")

K3 <- cbind(matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)), Tukey)
rownames(K3) <- paste(levels(d.fr$treat2)[3],rownames(K1), sep = ":")

K4 <- cbind(matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)), Tukey)
rownames(K4) <- paste(levels(d.fr$treat2)[4],rownames(K1), sep = ":")

K <- rbind(K1, K2, K3,K4)
colnames(K) <- c(colnames(Tukey), colnames(Tukey))

#K gives me a matrix with the dimension 12 by 6 and this will obviously
not work in the next step

summary(glht(fit2, linfct = K %*% X))

####

So, my questions is how would a 12 by 12 matrix for K look like, in
order to do K %*% X correctly for all levels within and between Tukey
adjusted.

But since I am not in interested in every *single* comparison within and
between, I was wondering to how to come up with a potential maximal but
empty matrix and code the comparisons myself.

Thanks for any hints!

Stefan



More information about the R-help mailing list