### multiple comparison procedures
### set up a one-way ANOVA
data(warpbreaks)
amod <- aov(breaks ~ tension, data = warpbreaks)
### specify all pair-wise comparisons among levels of variable "tension"
tuk <- glht(amod, linfct = mcp(tension = "Tukey"))
### extract information
tuk.cld <- cld(tuk)
### use sufficiently large upper margin
old.par <- par(mai=c(1,1,1.25,1), no.readonly=TRUE)
### plot
plot(tuk.cld)
par(old.par)
### now using covariates
amod2 <- aov(breaks ~ tension + wool, data = warpbreaks)
tuk2 <- glht(amod2, linfct = mcp(tension = "Tukey"))
tuk.cld2 <- cld(tuk2)
old.par <- par(mai=c(1,1,1.25,1), no.readonly=TRUE)
### use different colors for boxes
plot(tuk.cld2, col=c("green", "red", "blue"))
par(old.par)
### get confidence intervals
ci.glht <- confint(tuk)
### plot them
plot(ci.glht)
old.par <- par(mai=c(1,1,1.25,1), no.readonly=TRUE)
### use 'confint.glht' object to plot all pair-wise comparisons
plot(cld(ci.glht), col=c("white", "blue", "green"))
par(old.par)
### set up all pair-wise comparisons for count data
data(Titanic)
mod <- glm(Survived ~ Class, data = as.data.frame(Titanic),
weights = Freq, family = binomial())
### specify all pair-wise comparisons among levels of variable "Class"
glht.mod <- glht(mod, mcp(Class = "Tukey"))
### extract information
mod.cld <- cld(glht.mod)
### use sufficiently large upper margin
old.par <- par(mai=c(1,1,1.5,1), no.readonly=TRUE)
### plot
plot(mod.cld)
par(old.par)
### set up all pair-wise comparisons of a Cox-model
if (require("survival") && require("MASS")) {
### construct 4 classes of age
Melanoma$Cage <- factor(sapply(Melanoma$age, function(x){
if( x <= 25 ) return(1)
if( x > 25 & x <= 50 ) return(2)
if( x > 50 & x <= 75 ) return(3)
if( x > 75 & x <= 100) return(4) }
))
### fit Cox-model
cm <- coxph(Surv(time, status == 1) ~ Cage, data = Melanoma)
### specify all pair-wise comparisons among levels of "Cage"
cm.glht <- glht(cm, mcp(Cage = "Tukey"))
# extract information & plot
old.par <- par(no.readonly=TRUE)
### use mono font family
if (dev.interactive())
old.par <- par(family = "mono")
plot(cld(cm.glht), col=c("black", "red", "blue", "green"))
par(old.par)
}
if (require("nlme") && require("lme4")) {
data("ergoStool", package = "nlme")
stool.lmer <- lmer(effort ~ Type + (1 | Subject),
data = ergoStool)
glme41 <- glht(stool.lmer, mcp(Type = "Tukey"))
old.par <- par(mai=c(1,1,1.5,1), no.readonly=TRUE)
plot(cld(glme41))
par(old.par)
}
Run the code above in your browser using DataLab