I have a question about TukeyHSD and the glht function because I'm
getting different answers when a covariate is included in model for
ANCOVA. I'm using the cabbages dataset in the 'MASS' package for
repeatability. If I include HeadWt as a covariate, then I get different
answers when performing multiple comparisons using TukeyHSD and the glht
function. The difference appears related to the predicted means used in
the Tukey posthoc test. The glht function uses the marginal means, which
I can obtain using the effect function in the 'effects' package.
TukeyHSD generates the means using the model.tables function. Per the
help file, "the implementation is incomplete, and only simpler cases
have been tested thoroughly". Does this mean that TukeyHSD shouldn't be
used when covariates are included in the model? Could anyone elaborate
on why the model.tables function is generating means that differ from
the effect function? Thanks...
#--------------------------------
# Load libraries and data
library(multcomp)
library(effects)
data(cabbages, package='MASS')
#create the model
mod1 <- lm(VitC~HeadWt+Cult+Date, data=cabbages)
# Using TukeyHSD
TukeyHSD(aov(mod1), which='Date')
# Tukey multiple comparisons of means
# 95% family-wise confidence level
#
#Fit: aov(formula = mod1)
#
#$Date
# diff lwr upr p adj
#d20-d16 -0.9216847 -5.5216345 3.678265 0.8797985
#d21-d16 3.4237706 -1.1761792 8.023720 0.1814431
#d21-d20 4.3454553 -0.2544945 8.945405 0.0678038
# Tukey contrasts in glht should generate the same difference in means,
but it does not
summary(glht(mod1, linfct=mcp(Date='Tukey')))
#
# Simultaneous Tests for General Linear Hypotheses
#
#Multiple Comparisons of Means: Tukey Contrasts
#
#
#Fit: lm(formula = VitC ~ HeadWt + Cult + Date, data = cabbages)
#
#Linear Hypotheses:
# Estimate Std. Error t value Pr(>|t|)
#d20 - d16 == 0 -1.213 1.926 -0.630 0.8042
#d21 - d16 == 0 4.186 2.018 2.074 0.1044
#d21 - d20 == 0 5.400 2.112 2.556 0.0351 *
#---
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#(Adjusted p values reported -- single-step method)
# Differences used for glht could be obtained with the effect function
effect('Date', mod1)
#
# Date effect
#Date
# d16 d20 d21
#56.95889 55.74578 61.14533
# model.tables() is used to generate the means for TukeyHSD
model.tables(aov(mod1), type='means')$tables$Date
#Date
# d16 d20 d21
#57.11597 56.19429 60.53974
#Warning message:
#In replications(paste("~", xx), data = mf) : non-factors ignored: HeadWt
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
