## 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
## Not run:
# ## 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
# ## End(Not run)
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\nof the two levels in the contrasts")
tmp.dunnett.mmc
## Not run:
# ## 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")
# ## End(Not run)
## Not run:
# ## 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")
# ## End(Not run)
Run the code above in your browser using DataLab