## This example is from Hsu and Peruggia
## This is the R version
## See ?aov.sufficient for S-Plus
if.R(s={},
r={
pulmonary <- read.table(hh("datasets/pulmonary.dat"), header=TRUE,
row.names=NULL)
names(pulmonary)[3] <- "FVC"
names(pulmonary)[1] <- "smoker"
pulmonary$smoker <- factor(pulmonary$smoker, levels=pulmonary$smoker)
row.names(pulmonary) <- pulmonary$smoker
pulmonary
pulmonary.aov <- aov.sufficient(FVC ~ smoker,
data=pulmonary)
summary(pulmonary.aov)
pulmonary.mca <- glht(pulmonary.aov,
linfct=mcp(smoker="Tukey"),
df=pulmonary.aov$df.residual,
vcov.=vcov.sufficient)
plot(pulmonary.mca)
pulm.lmat <- cbind("npnl-mh"=c( 1, 1, 1, 1,-2,-2), ## not.much vs lots
"n-pnl" =c( 3,-1,-1,-1, 0, 0), ## none vs light
"p-nl" =c( 0, 2,-1,-1, 0, 0), ## {} arbitrary 2 df
"n-l" =c( 0, 0, 1,-1, 0, 0), ## {} for 3 types of light
"m-h" =c( 0, 0, 0, 0, 1,-1)) ## moderate vs heavy
dimnames(pulm.lmat)[[1]] <- row.names(pulmonary)
if.R(r=pulm.lmat <- rbind(Int=0, pulm.lmat[-1,]),
s={})
pulm.lmat
pulmonary.mmc <- glht.mmc(pulmonary.aov,
linfct=mcp(smoker="Tukey"),
df=pulmonary.aov$df.residual,
vcov.=vcov.sufficient,
lmat=pulm.lmat,
calpha=attr(confint(pulmonary.mca)$confint,"calpha"))
plot(pulmonary.mmc, print.mca=TRUE, print.lmat=FALSE)
## tiebreaker plot, with contrasts ordered to match MMC plot,
## with all contrasts forced positive and with names also reversed,
## and with matched x-scale.
plot(confint(as.glht(pulmonary.mmc$mca)),
xlim=par()$usr[1:2], xaxs="i",
main="", xlab="")
## orthogonal contrasts
plot(pulmonary.mmc, print.lmat=TRUE, col.lmat.signif='blue', col.iso='gray')
## pairwise and orthogonal contrasts on the same plot
plot(pulmonary.mmc, print.mca=TRUE, print.lmat=TRUE,
col.mca.signif='red', col.lmat.signif='blue', col.iso='gray',
lty.lmat.not.signif=2)
})
Run the code above in your browser using DataLab