如何获得lsmeans()与自定义vcov的成对对比？

Question

我希望在提供稳健的系数-协方差矩阵(例如vcovHC)的同时，使用lsmeans()对调整后的平均值进行成对比较。通常，回归模型上的函数会提供一个vcov参数，但我似乎在lsmeans包中找不到任何这样的参数。

考虑这个虚拟示例，最初来自CAR：

require(car)
require(lmtest)
require(sandwich)
require(lsmeans)

mod.moore.2 <- lm(conformity ~ fcategory + partner.status, data=Moore)
coeftest(mod.moore.2)
## 
## t test of coefficients:
## 
##                     Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)        10.197778   1.372669  7.4292 4.111e-09 ***
## fcategorymedium    -1.176000   1.902026 -0.6183  0.539805    
## fcategoryhigh      -0.080889   1.809187 -0.0447  0.964555    
## partner.statushigh  4.606667   1.556460  2.9597  0.005098 ** 
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

coeftest(mod.moore.2, vcov.=vcovHAC)
## 
## t test of coefficients:
## 
##                     Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)        10.197778   0.980425 10.4014 4.565e-13 ***
## fcategorymedium    -1.176000   1.574682 -0.7468  0.459435    
## fcategoryhigh      -0.080889   2.146102 -0.0377  0.970117    
## partner.statushigh  4.606667   1.437955  3.2036  0.002626 ** 
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

lsmeans(mod.moore.2, list(pairwise ~ fcategory), adjust="none")[[2]]
##  contrast         estimate       SE df t.ratio p.value
##  low - medium   1.17600000 1.902026 41   0.618  0.5398
##  low - high     0.08088889 1.809187 41   0.045  0.9646
##  medium - high -1.09511111 1.844549 41  -0.594  0.5560
## 
## Results are averaged over the levels of: partner.status 

如您所见，lsmeans()使用默认的方差-协方差矩阵来估计p值。

如何使用vcovHAC方差估计获得成对对比？

Answer 1

事实证明，在lsmeans和multcomp包之间有一个非常好的无缝接口(请参阅?lsm)，而lsmeans提供了对glht()的支持。

require(multcomp)

x <- glht(mod.moore.2, lsm(pairwise ~ fcategory), vcov=vcovHAC)
## Note: df set to 41
summary(x, test=adjusted("none"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = conformity ~ fcategory + partner.status, data = Moore)
## 
## Linear Hypotheses:
##                    Estimate Std. Error t value Pr(>|t|)
## low - medium == 0   1.17600    1.57468   0.747    0.459
## low - high == 0     0.08089    2.14610   0.038    0.970
## medium - high == 0 -1.09511    1.86197  -0.588    0.560
## (Adjusted p values reported -- none method)

这至少是实现这一目标的一种方法。我仍然希望有人知道一种只使用lsmeans的方法...

另一种方法是侵入lsmeans对象，并在summary-ing对象之前手动替换方差-协方差矩阵。

mod.lsm <- lsmeans(mod.moore.2, ~ fcategory)
mod.lsm@V <- vcovHAC(mod.moore.2)  ##replace default vcov with custom vcov
pairs(mod.lsm, adjust = "none")
##  contrast         estimate       SE df t.ratio p.value
##  low - medium   1.17600000 1.574682 41   0.747  0.4594
##  low - high     0.08088889 2.146102 41   0.038  0.9701
##  medium - high -1.09511111 1.861969 41  -0.588  0.5597
## 
## Results are averaged over the levels of: partner.status 

Answer 2

或者使用update将自定义vcov矩阵注入到emmeans/emmGrid对象中。

示例：

# create an emmeans object from your fitted model
emmob <- emmeans(thismod, ~ predictor)

# generate a robust vcov matrix using a function
# from the sandwich or clubSandwich package
vcovR <- vcovHC(thismod, type="HC3")

# turn the resulting object into a (square) matrix
vcovRm <- matrix(vcovR, ncol=ncol(vcovR))

# update the V slot of the emmeans/emmGrid object
emmob <- update(emmob, V=vcovRm)