3

我想使用 对调整后的均值进行成对比较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方差估计获得成对对比？

4

3 回答 3

4

事实证明，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。

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
于 2015-07-08T19:07:23.837 回答
0

我不确定这是否可以使用“lsmeans”包，但它正在使用更新的emmeans包。


Moore <- within(carData::Moore, {
  
  partner.status <- factor(partner.status, c("low", "high"))
  fcategory      <- factor(fcategory, c("low", "medium", "high"))
  
})

mod.moore.2 <- lm(conformity ~ fcategory + partner.status, data=Moore)
lmtest::coeftest(mod.moore.2, vcov.= sandwich::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
emmeans::emmeans(
  mod.moore.2, trt.vs.ctrl ~ fcategory, 
  vcov = sandwich::vcovHAC(mod.moore.2),
  adjust = "none")$contrasts
#>  contrast     estimate   SE df t.ratio p.value
#>  medium - low  -1.1760 1.57 41 -0.747  0.4594 
#>  high - low    -0.0809 2.15 41 -0.038  0.9701 
#> 
#> Results are averaged over the levels of: partner.status

由reprex 包（v0.3.0）于 2021-07-08 创建

请注意，您不能只写以下内容

emmeans::emmeans(
  mod.moore.2, trt.vs.ctrl ~ fcategory, 
  vcov = sandwich::vcovHAC,
  adjust = "none")$contrasts

adjust由于与也有一个选项的 Sandwich::vcovHAC 命令冲突。（我错误地认为这是一个错误）。

于 2021-07-08T05:22:49.043 回答
0

或用于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)
于 2021-12-09T09:22:31.717 回答