[R] TukeyHSD and glht differ for models with a covariate

[Nicholas Negovetich]

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