Learn R
HH (version 2.2-17)

mmc: MMC (mean--mean multiple comparisons) plots.

Description

Constructs a "mmc.multicomp" object from the formula and other arguments. The object must be explicitly plotted.

Usage

glht.mmc(model, ...)  ## R

## S3 method for class 'glht':
glht.mmc(model, ...)


## S3 method for class 'lm':
glht.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)

Arguments

Value

An "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.mcaObject containing the pairwise comparisons.noneObject comparing each mean to 0.lmatObject for the contrasts specified in the lmat argument."[.mmc.multicomp" is a subscript method.

Details

By default, if 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. glht.mmc in R works from either an "glht" object or an "aov" object. multicomp.mmc in S-Plus works from an "aov" object.

References

Heiberger, Richard M. and Holland, Burt (2004b). Statistical Analysis and Data Display: An Intermediate Course with Examples in S-Plus, R, and SAS. Springer Texts in Statistics. Springer. ISBN 0-387-40270-5. Heiberger, R.~M. and Holland, B. (2006). "Mean--mean multiple comparison displays for families of linear contrasts." Journal of Computational and Graphical Statistics, 15:937--955. Hsu, J. and Peruggia, M. (1994). "Graphical representations of Tukey's multiple comparison method." Journal of Computational and Graphical Statistics, 3:143--161.

See Also

as.multicomp, plot.mmc.multicomp

Examples

Run this code
## Use glht.mmc with R.
## Use multicomp.mmc with S-Plus.

## data and ANOVA
## catalystm example
catalystm <- read.table(hh("datasets/catalystm.dat"), header=FALSE,
                       col.names=c("catalyst","concent"))
catalystm$catalyst <- factor(catalystm$catalyst, labels=c("A","B","C","D"))

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
old.omd <- par(omd=c(0,.95,0,1))
catalystm.mmc <-
if.R(r=glht.mmc(catalystm1.aov, linfct = mcp(catalyst = "Tukey")),
     s=multicomp.mmc(catalystm1.aov, plot=FALSE))
catalystm.mmc
if.R(s=plot(catalystm.mmc, x.offset=1),
     r=plot(catalystm.mmc, ry=c(50,58), x.offset=1.8))
## tiebreaker
plot.matchMMC(catalystm.mmc$mca, xlabel.print=FALSE)

## 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=glht.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),
     r=plot(catalystm.mmc, ry=c(50,58), x.offset=1.8))
## tiebreaker
plot.matchMMC(catalystm.mmc$lmat, xlabel.print=FALSE, col.signif='blue')
par(old.omd)


## Dunnett's test
## weightloss example
weightloss <- read.table(hh("datasets/weightloss.dat"), header=TRUE)
weightloss <- data.frame(loss=unlist(weightloss),
                         group=rep(names(weightloss), rep(10,5)))
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={})

tmp.dunnett <- 
if.R(r=
glht(weightloss.aov,
     linfct=mcp(group=contrMat(group.count, base=4)),
     alternative="greater")
,s=
multicomp(weightloss.aov,
          method="dunnett", comparisons="mcc",
          bounds="lower", control=4,
          valid.check=FALSE)
)
plot(tmp.dunnett)

tmp.dunnett.mmc <-
if.R(r=
   glht.mmc(weightloss.aov,
            linfct=mcp(group=contrMat(group.count, base=4)),
            alternative="greater")
,s=
   multicomp.mmc(weightloss.aov,
                 method="dunnett", comparisons="mcc",
                 bounds="lower", control=4,
                 valid.check=FALSE, plot=FALSE)
)

tmp.dunnett.mmc
plot(tmp.dunnett.mmc)



## two-way ANOVA
## display example
display <- read.table(hh("datasets/display.dat"), header=TRUE)
display$panel <- factor(display$panel)
position(display$panel) <- (1:3)+.5
if.R(r={
        contrasts(display$panel) <- "contr.treatment"
       }
    ,s={})
display$emergenc <- factor(display$emergenc)

if.R(r=
interaction2wt(time ~ emergenc * panel, data=display)
     ,s=
interaction2wt(time ~ emergenc * panel, 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
displayf.mmc <-
if.R(r={
        glht.mmc(displayf.aov, linfct=mcp(panel="Tukey"),
                 interaction_average=TRUE)
},
     s=multicomp.mmc(displayf.aov, "panel", plot=FALSE))

if.R(s=
       plot(displayf.mmc)
    ,r=
       plot(displayf.mmc, x.offset=1.7, ry=c(17,26.5))
)


panel.lmat <- cbind("3-12"=c(-1,-1,2),
                     "1-2"=c( 1,-1,0))
dimnames(panel.lmat)[[1]] <- levels(display$panel)

displayf.mmc <-
if.R(r={
        glht.mmc(displayf.aov, linfct=mcp(panel="Tukey"),
                 interaction_average=TRUE, focus.lmat=panel.lmat)
},
     s=multicomp.mmc(displayf.aov, "panel",
       focus.lmat=panel.lmat, plot=FALSE))

if.R(s=
       plot(displayf.mmc)
    ,r=
       plot(displayf.mmc, x.offset=1.7, ry=c(17,26.5))
)


## 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.



maiz <- read.table(hh("datasets/maiz.dat"), header=TRUE)
maiz$hibrido <- factor(maiz$hibrido,
                       levels=c("P3747","P3732","Mol17","A632","LH74"))
maiz$nitrogeno <- factor(maiz$nitrogeno)
position(maiz$nitrogeno) <- c(1.2, 2.4, 3.6, 4.8) ## inherits(maiz$nitrogeno, "ordered")
position(maiz$bloque) <- c(2.25, 3.75)            ## inherits(maiz$bloque, "ordered")

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 glht.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))
  plot.matchMMC(maiz2.mmc$mca)
  par(old.omd)
},r={
  maiz2.mmc <- glht.mmc(maiz2.aov, linfct=mcp(hibrido="Tukey"),
                    interaction_average=TRUE)
  old.omd <- par(omd=c(.05,.85,.35,1)) ## x1 x2 y1 y2
  plot(maiz2.mmc)
  par(omd=c(.05,.85,0,.5), new=TRUE)
  plot.matchMMC(maiz2.mmc$mca, cex.axis=.7)
  par(old.omd)
})


## 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=glht.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))
  plot.matchMMC(maiz2.mmc$lmat, col.signif='blue')
  par(old.omd)
},r={
  old.omd <- par(omd=c(.05,.85,.35,1)) ## x1 x2 y1 y2
  plot(maiz2.mmc)
  par(omd=c(.05,.85,0,.45), new=TRUE)
  plot.matchMMC(maiz2.mmc$lmat, cex.axis=.7, col.signif='blue')
  par(old.omd)
})

Run the code above in your browser using DataCamp Workspace