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

Question

我想使用 lsmeans() 对调整后的均值进行成对比较，同时提供稳健的系数协方差矩阵(例如 vcovHC)。通常回归模型上的函数提供一个 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 包之间有一个美妙而无缝的接口(interface)(参见 ?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 对象之前手动替换方差-协方差矩阵。

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