## Use mmc with R.
## Use multicomp.mmc with S-Plus.
## data and ANOVA
## catalystm example
data(catalystm)
bwplot(concent ~ catalyst, data=catalystm,
scales=list(cex=1.5),
ylab=list("concentration", cex=1.5),
xlab=list("catalyst",cex=1.5))
catalystm1.aov <- aov(concent ~ catalyst, data=catalystm)
summary(catalystm1.aov)
catalystm.mca <-
glht(catalystm1.aov, linfct = mcp(catalyst = "Tukey"))
confint(catalystm.mca)
plot(catalystm.mca) ## multcomp plot
mmcplot(catalystm.mca, focus="catalyst") ## HH plot
## pairwise comparisons
catalystm.mmc <-
mmc(catalystm1.aov, focus="catalyst")
catalystm.mmc
## these three statements are identical for a one-way aov
mmc(catalystm1.aov) ## simplest
mmc(catalystm1.aov, focus="catalyst") ## generalizes to higher-order designs
mmc(catalystm1.aov, linfct = mcp(catalyst = "Tukey")) ## glht arguments
mmcplot(catalystm.mmc, style="both")
## User-Specified Contrasts
## Row names must include all levels of the factor.
## Column names are the names the user assigns to the contrasts.
## Each column must sum to zero.
catalystm.lmat <- cbind("AB-D" =c( 1, 1, 0,-2),
"A-B" =c( 1,-1, 0, 0),
"ABD-C"=c( 1, 1,-3, 1))
dimnames(catalystm.lmat)[[1]] <- levels(catalystm$catalyst)
catalystm.lmat
catalystm.mmc <-
mmc(catalystm1.aov,
linfct = mcp(catalyst = "Tukey"),
focus.lmat=catalystm.lmat)
catalystm.mmc
mmcplot(catalystm.mmc, style="both", type="lmat")
## Dunnett's test
## weightloss example
data(weightloss)
bwplot(loss ~ group, data=weightloss,
scales=list(cex=1.5),
ylab=list("Weight Loss", cex=1.5),
xlab=list("group",cex=1.5))
weightloss.aov <- aov(loss ~ group, data=weightloss)
summary(weightloss.aov)
group.count <- table(weightloss$group)
tmp.dunnett <-
glht(weightloss.aov,
linfct=mcp(group=contrMat(group.count, base=4)),
alternative="greater")
mmcplot(tmp.dunnett, main="contrasts in alphabetical order", focus="group")
tmp.dunnett.mmc <-
mmc(weightloss.aov,
linfct=mcp(group=contrMat(group.count, base=4)),
alternative="greater")
mmcplot(tmp.dunnett.mmc,
main="contrasts ordered by average value of the means
of the two levels in the contrasts")
tmp.dunnett.mmc
## two-way ANOVA
## display example
data(display)
interaction2wt(time ~ emergenc * panel.ordered, data=display)
displayf.aov <- aov(time ~ emergenc * panel, data=display)
anova(displayf.aov)
## multiple comparisons
## MMC plot
displayf.mmc <- mmc(displayf.aov, focus="panel")
displayf.mmc
## same thing using glht argument list
displayf.mmc <-
mmc(displayf.aov,
linfct=mcp(panel="Tukey", `interaction_average`=TRUE, `covariate_average`=TRUE))
mmcplot(displayf.mmc)
panel.lmat <- cbind("3-12"=c(-1,-1,2),
"1-2"=c( 1,-1,0))
dimnames(panel.lmat)[[1]] <- levels(display$panel)
panel.lmat
displayf.mmc <-
mmc(displayf.aov, focus="panel", focus.lmat=panel.lmat)
## same thing using glht argument list
displayf.mmc <-
mmc(displayf.aov,
linfct=mcp(panel="Tukey", `interaction_average`=TRUE, `covariate_average`=TRUE),
focus.lmat=panel.lmat)
mmcplot(displayf.mmc, type="lmat")
## split plot design with tiebreaker plot
##
## This example is based on the query by Tomas Goicoa to R-news
## http://article.gmane.org/gmane.comp.lang.r.general/76275/match=goicoa
## It is a split plot similar to the one in HH Section 14.2 based on
## Yates 1937 example. I am using the Goicoa example here because its
## MMC plot requires a tiebreaker plot.
data(maiz)
interaction2wt(yield ~ hibrido+nitrogeno+bloque, data=maiz,
par.strip.text=list(cex=.7))
interaction2wt(yield ~ hibrido+nitrogeno, data=maiz)
maiz.aov <- aov(yield ~ nitrogeno*hibrido + Error(bloque/nitrogeno), data=maiz)
summary(maiz.aov)
summary(maiz.aov,
split=list(hibrido=list(P3732=1, Mol17=2, A632=3, LH74=4)))
try(glht(maiz.aov, linfct=mcp(hibrido="Tukey"))) ## can't use 'aovlist' objects in glht
## R glht() requires aov, not aovlist
maiz2.aov <- aov(terms(yield ~ bloque*nitrogeno + hibrido/nitrogeno,
keep.order=TRUE),
data=maiz)
summary(maiz2.aov)
## There are many ties in the group means.
## These are easily seen in the MMC plot, where the two clusters
## c("P3747", "P3732", "LH74") and c("Mol17", "A632")
## are evident from the top three contrasts including zero and the
## bottom contrast including zero. The significant contrasts are the
## ones comparing hybrids in the top group of three to ones in the
## bottom group of two.
## We have two graphical responses to the ties.
## 1. We constructed the tiebreaker plot.
## 2. We construct a set of orthogonal contrasts to illustrate
## the clusters.
## pairwise contrasts with tiebreakers.
maiz2.mmc <- mmc(maiz2.aov,
linfct=mcp(hibrido="Tukey", interaction_average=TRUE))
mmcplot(maiz2.mmc, style="both") ## MMC and Tiebreaker
## orthogonal contrasts
## user-specified contrasts
hibrido.lmat <- cbind("PPL-MA" =c(2, 2,-3,-3, 2),
"PP-L" =c(1, 1, 0, 0,-2),
"P47-P32"=c(1,-1, 0, 0, 0),
"M-A" =c(0, 0, 1,-1, 0))
dimnames(hibrido.lmat)[[1]] <- levels(maiz$hibrido)
hibrido.lmat
maiz2.mmc <-
mmc(maiz2.aov, focus="hibrido", focus.lmat=hibrido.lmat)
maiz2.mmc
## same thing using glht argument list
maiz2.mmc <-
mmc(maiz2.aov, linfct=mcp(hibrido="Tukey",
`interaction_average`=TRUE), focus.lmat=hibrido.lmat)
mmcplot(maiz2.mmc, style="both", type="lmat")
Run the code above in your browser using DataLab