## Use glht.mmc with R.
## Use multicomp.mmc with S-Plus.
## data and ANOVA
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)
catalystm.mca <-
if.R(r=glht(catalystm1.aov, linfct = mcp(catalyst = "Tukey")),
s=multicomp(catalystm1.aov, plot=FALSE))
plot(catalystm.mca)
catalystm.mca
## pairwise comparisons
catalystm.mmc <-
if.R(r=glht.mmc(catalystm1.aov, linfct = mcp(catalyst = "Tukey")),
s=multicomp.mmc(catalystm1.aov, plot=FALSE))
catalystm.mmc
plot(catalystm.mmc)
plot(catalystm.mmc$mca)
plot(catalystm.mmc$none)
## user-specified contrasts
catalystm.lmat <- cbind("AB-D" =c(0, 1, 1, 0,-2),
"A-B" =c(0, 1,-1, 0, 0),
"ABD-C"=c(0, 1, 1,-3, 1))
if.R(r=catalystm.lmat <- catalystm.lmat[-2,],
s={})
dimnames(catalystm.lmat)[[1]] <- dimnames(catalystm.mmc$mca$lmat)[[1]]
zapsmall(catalystm.lmat)
if.R(s=dimnames(catalystm.mca$lmat)[[1]],
r=dimnames(catalystm.mca$linfct)[[2]])
catalystm.mmc <-
if.R(r=glht.mmc(catalystm1.aov, linfct = mcp(catalyst = "Tukey"),
lmat=catalystm.lmat)
,s=multicomp.mmc(catalystm1.aov, lmat=catalystm.lmat,
plot=FALSE)
)
catalystm.mmc
plot(catalystm.mmc)
plot(catalystm.mmc$mca)
plot(catalystm.mmc$none)
plot(catalystm.mmc$lmat)
## Dunnett's test
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 <- read.table(hh("datasets/display.dat"), header=TRUE)
display$panel <- factor(display$panel) ## display$panel <- positioned(display$panel, value=(1:3)+.5)
display$emergenc <- factor(display$emergenc)
displayf.aov <- aov(time ~ emergenc * panel, data=display)
anova(displayf.aov)
## multiple comparisons
tmp <- if.R(
r=glht(displayf.aov, linfct=mcp(panel="Tukey")),
s=multicomp(displayf.aov, "panel", plot=FALSE))
zapsmall(
if.R(r=t(tmp$linfct),
s=tmp$lmat)
)
## MMC plot
displayf.mmc <-
if.R(r=glht.mmc(displayf.aov, linfct=mcp(panel="Tukey"), focus="panel", lmat.rows=5:6),
s=multicomp.mmc(displayf.aov, "panel", lmat.rows=6:8, plot=FALSE))
plot(displayf.mmc)
## orthogonal contrasts
zapsmall(mca.lmat <- displayf.mmc$mca$lmat)
panel.lmat <- cbind("3-12"=mca.lmat[,1] + mca.lmat[,2],
"1-2"=mca.lmat[,3])
displayf.mmc <-
if.R(r=glht.mmc(displayf.aov, linfct=mcp(panel="Tukey"), focus="panel",
lmat.rows=5:6, lmat=panel.lmat),
s=multicomp.mmc(displayf.aov, "panel", lmat.rows=6:8,
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.
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.5, 4, 5.5) ## forces class="ordered"
interaction2wt(yield ~ hibrido+nitrogeno+bloque, data=maiz)
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)))
## 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"
sapply(maiz[-1], contrasts)
},
s={})
## Both R and S-Plus require aov, not aovlist
maiz2.aov <- aov(terms(yield ~ bloque*nitrogeno + hibrido/nitrogeno,
keep.order=TRUE), data=maiz)
summary(maiz2.aov)
if.R(s={
maiz2.mca <- multicomp(maiz2.aov, focus="hibrido")
## plot(maiz2.mca)
dimnames(maiz2.mca$lmat)[[1]]
maiz2.mmc <- multicomp.mmc(maiz2.aov, focus="hibrido",
lmat.rows=16:20, plot=FALSE)
old.mar <- par(mar=c(15,4,4,7)+.1)
plot(maiz2.mmc)
par(mar=c(2,4,28,7)+.1, new=TRUE)
old.cex <- par(cex=.8)
plot(maiz2.mmc$mca, col.signif=8, lty.signif=1, xlabel.print=FALSE,
xaxs="d", plt=par()$plt+c(0,0,-.25,.05), xrange.include=c(-30,40))
par(old.cex)
par(old.mar)
},r={
maiz2.mca <- glht(maiz2.aov, linfct=mcp(hibrido="Tukey"))
dimnames(maiz2.mca$linfct)[[2]]
maiz2.mmc <- glht.mmc(maiz2.aov, linfct=mcp(hibrido="Tukey"), lmat.rows=9:12)
old.oma <- par(oma=c(12,3,0,4))
plot(maiz2.mmc)
par(oma=c(0,3,22,4), new=TRUE)
plot(maiz2.mmc$mca,
xlim=par()$usr[1:2], xaxs="i",
main="", xlab="", cex.axis=.7)
par(old.oma)
})
Run the code above in your browser using DataCamp Workspace