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