"mmc.multicomp" object from the formula and
other arguments. The constructed object must be explicitly plotted
with the mmcplot function.mmc(model, ...) ## R
## S3 method for class 'glht':
mmc(model, ...)
## S3 method for class 'default':
mmc(model, ## lm object
linfct=NULL,
focus=
if (is.null(linfct))
{
if (length(model$contrasts)==1) names(model$contrasts)
else stop("focus or linfct must be specified.")
}
else
{
if (is.null(names(linfct)))
stop("focus must be specified.")
else names(linfct)
},
focus.lmat,
ylabel=deparse(terms(model)[[2]]),
lmat=if (missing(focus.lmat)) {
t(linfct)
} else {
lmatContrast(t(none.glht$linfct), focus.lmat)
},
lmat.rows=lmatRows(model, focus),
lmat.scale.abs2=TRUE,
estimate.sign=1,
order.contrasts=TRUE,
level=.95,
calpha=NULL,
alternative = c("two.sided", "less", "greater"),
...
)
multicomp.mmc(x, ## S-Plus
focus=dimnames(attr(x$terms,"factors"))[[2]][1],
comparisons="mca",
lmat,
lmat.rows=lmatRows(x, focus),
lmat.scale.abs2=TRUE,
ry,
plot=TRUE,
crit.point,
iso.name=TRUE,
estimate.sign=1,
x.offset=0,
order.contrasts=TRUE,
main,
main2,
focus.lmat,
...)
## S3 method for class 'mmc.multicomp':
[(x, ..., drop = TRUE)"aov" object in "lm" method.multicomp.
The convention for lmat in R is to use
the transpose of the linfct component produced by
glht. Required for user-specified contrasts.lmat for the focus factor.glht.
#endif
#ifdef S-Plus
multicomp.
#endiflmat from the none
component to create the lmat for tglht.
#endif
#ifdef S-Plus
multicomp.
#endifalternative and
base are frequently used with glht.multicomplmat to make the sum of the absolute values of each column equal 2.0, leave contrasts in the
default lexicographic direction. If positive, force all contrasts to positive,
reversing their names if needed (if contrast A-B is negative, reverse it
to B-A). If negative, the force almca, none,
lmat) components by height on the MMC plot. This will place the
contrasts in the multicomp plots in the same order as in the MMC plot.glht
#endif
#ifdef S-Plus
glht
#endif
in R. S-Plus multicomp
uses the argument boundsmulticomp method is used for the user-specified
contrasts when lmat is specified. This argument is called
crit.point with multicompTRUE.plot.mmc.multicomp."["."mmc.multicomp" object contains either the first two or all
three of the "multicomp" components mca, none,
lmat described here. Each "multicomp" component in
R also contains a "glht" object.lmat argument."[.mmc.multicomp" is a subscript method.lmat is not specified, we plot the isomeans grid
and the pairwise comparisons for the focus factor. By default,
we plot the specified contrasts if the lmat is specified.
Each contrast is plotted at a height which is the weighted average of
the means being compared. The weights are scaled to the sum of their
absolute values equals 2.
We get the right contrasts automatically if the aov is oneway. If we
specify an lmat for oneway it must have a leading row of 0.
For any more complex design, we must study the lmat from the mca
component of the result to see how to construct the lmat (with the
extra rows as needed) and how to specify the lmat.rows
corresponding to the rows for the focus factor.
mmc in R works from either an "glht" object or an
"aov" object. multicomp.mmc in S-Plus works from an
"aov" object.as.multicomp, plot.mmc.multicomp## Use mmc with R.
## Use multicomp.mmc with S-Plus.
## data and ANOVA
## catalystm example
data(catalystm)
if.R(r=
bwplot(concent ~ catalyst, data=catalystm,
scales=list(cex=1.5),
ylab=list("concentration", cex=1.5),
xlab=list("catalyst",cex=1.5))
,s=
t(bwplot(catalyst ~ concent, data=catalystm,
scales=list(cex=1.5),
xlab=list("concentration", cex=1.5),
ylab=list("catalyst",cex=1.5)))
)
catalystm1.aov <- aov(concent ~ catalyst, data=catalystm)
summary(catalystm1.aov)
catalystm.mca <-
if.R(r=glht(catalystm1.aov, linfct = mcp(catalyst = "Tukey")),
s=multicomp(catalystm1.aov, plot=FALSE))
## plot(catalystm.mca)
if.R(s=catalystm.mca,
r=confint(catalystm.mca))
## pairwise comparisons
catalystm.mmc <-
if.R(r=mmc(catalystm1.aov, linfct = mcp(catalyst = "Tukey")),
s=multicomp.mmc(catalystm1.aov, plot=FALSE))
catalystm.mmc
if.R(s={old.omd <- par(omd=c(0,.95,0,1))
plot(catalystm.mmc, x.offset=1)
plotMatchMMC(catalystm.mmc$mca, xlabel.print=FALSE)},
r=mmcplot(catalystm.mmc, style="both")
)
## user-specified contrasts
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.mmc <-
if.R(r=mmc(catalystm1.aov,
linfct = mcp(catalyst = "Tukey"),
focus.lmat=catalystm.lmat),
s=multicomp.mmc(catalystm1.aov,
focus.lmat=catalystm.lmat, plot=FALSE))
catalystm.mmc
if.R(s={plot(catalystm.mmc, x.offset=1)
plotMatchMMC(catalystm.mmc$lmat, xlabel.print=FALSE, col.signif='blue')
par(old.omd)},
r=mmcplot(catalystm.mmc, style="both", type="lmat"))
## Dunnett's test
## weightloss example
data(weightloss)
if.R(r=
bwplot(loss ~ group, data=weightloss,
scales=list(cex=1.5),
ylab=list("Weight Loss", cex=1.5),
xlab=list("group",cex=1.5))
,s=
t(bwplot(group ~ loss, data=weightloss,
scales=list(cex=1.5),
xlab=list("Weight Loss", cex=1.5),
ylab=list("group",cex=1.5)))
)
weightloss.aov <- aov(loss ~ group, data=weightloss)
summary(weightloss.aov)
if.R(r={
group.count <- table(weightloss$group)
},s={})
if.R(r={
tmp.dunnett <-
glht(weightloss.aov,
linfct=mcp(group=contrMat(group.count, base=4)),
alternative="greater")
mmcplot(tmp.dunnett)
},s={
tmp.dunnett <-
multicomp(weightloss.aov,
method="dunnett", comparisons="mcc",
bounds="lower", control=4,
valid.check=FALSE)
plot(tmp.dunnett)
})
if.R(r={
tmp.dunnett.mmc <-
mmc(weightloss.aov,
linfct=mcp(group=contrMat(group.count, base=4)),
alternative="greater")
mmcplot(tmp.dunnett.mmc)
},s={
tmp.dunnett.mmc <-
multicomp.mmc(weightloss.aov,
method="dunnett", comparisons="mcc",
bounds="lower", control=4,
valid.check=FALSE, plot=FALSE)
plot(tmp.dunnett.mmc)
})
tmp.dunnett.mmc
## two-way ANOVA
## display example
data(display)
if.R(r=
interaction2wt(time ~ emergenc * panel.ordered, data=display)
,s=
interaction2wt(time ~ emergenc * panel.ordered, data=display,
xlim=c(.5,4.5), key.in=list(x=-1.8))
)
displayf.aov <- aov(time ~ emergenc * panel, data=display)
anova(displayf.aov)
## multiple comparisons
## MMC plot
if.R(r={
displayf.mmc <-
mmc(displayf.aov,
linfct=mcp(panel="Tukey", `interaction_average`=TRUE, `covariate_average`=TRUE))
mmcplot(displayf.mmc)
},
s={
displayf.mmc <-
multicomp.mmc(displayf.aov, "panel", plot=FALSE)
plot(displayf.mmc)
})
panel.lmat <- cbind("3-12"=c(-1,-1,2),
"1-2"=c( 1,-1,0))
dimnames(panel.lmat)[[1]] <- levels(display$panel)
if.R(r={
displayf.mmc <-
mmc(displayf.aov,
linfct=mcp(panel="Tukey", `interaction_average`=TRUE, `covariate_average`=TRUE),
focus.lmat=panel.lmat)
mmcplot(displayf.mmc, type="lmat")
},s={
displayf.mmc <-
multicomp.mmc(displayf.aov, "panel",
focus.lmat=panel.lmat, plot=FALSE)
plot(displayf.mmc)
}
)
## 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)
if.R(s={old.omd <- par(omd=c(.1,1,.05,1))},
r={})
interaction2wt(yield ~ hibrido+nitrogeno+bloque, data=maiz,
key.in=list(x=-5), ## ignored by R
par.strip.text=list(cex=.7))
interaction2wt(yield ~ hibrido+nitrogeno, data=maiz)
if.R(s={par(old.omd)},
r={})
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)))
## multicomp(maiz.aov, focus="hibrido") ## can't use 'aovlist' objects
## glht(maiz.aov, linfct=mcp(hibrido="Tukey")) ## can't use 'aovlist' objects
sapply(maiz[-1], contrasts)
if.R(r={
## R mmc requires treatment contrasts
contrasts(maiz$nitrogeno) <- "contr.treatment"
contrasts(maiz$bloque) <- "contr.treatment"
sapply(maiz[-1], contrasts)
},
s={})
## Both R glht() and S-Plus multicomp() require 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.
if.R(s={
maiz2.mmc <- multicomp.mmc(maiz2.aov, focus="hibrido", plot=FALSE)
old.omd <- par(omd=c(.05,.85,0,1))
plot(maiz2.mmc, ry=c(145,170), x.offset=4)
par(omd=c(.05,.85,0,1))
plotMatchMMC(maiz2.mmc$mca)
par(old.omd)
},r={
maiz2.mmc <- mmc(maiz2.aov,
linfct=mcp(hibrido="Tukey", interaction_average=TRUE))
mmcplot(maiz2.mmc, style="both")
})
## 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 <-
if.R(s=multicomp.mmc(maiz2.aov, focus="hibrido",
focus.lmat=hibrido.lmat,
plot=FALSE),
r=mmc(maiz2.aov, linfct=mcp(hibrido="Tukey",
`interaction_average`=TRUE), focus.lmat=hibrido.lmat)
)
if.R(s={
old.omd <- par(omd=c(.05,.85,0,1))
plot(maiz2.mmc, ry=c(145,170), x.offset=4)
par(omd=c(.05,.85,0,1))
plotMatchMMC(maiz2.mmc$lmat, col.signif='blue')
par(old.omd)
},r={
mmcplot(maiz2.mmc, style="both", type="lmat")
})Run the code above in your browser using DataLab