How to obtain lsmeans() pairwise contrasts with custom vcov?

Question

I would like to get pairwise comparisons of adjusted means using lsmeans(), while supplying a robust coefficient-covariance matrix (e.g. vcovHC). Usually functions on regression models provide a vcov argument, but I can't seem to find any such argument in the lsmeans package.

Consider this dummy example, originally from 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

As you can see, lsmeans() estimates p-values using the default variance-covariance matrix.

How can I obtain pairwise contrasts using the vcovHAC variance estimate?

landroni · Answer 1 · 2015-07-09T07:05:56.987

It turns out that there is a wonderful and seamless interface between lsmeans and multcomp packages (see ?lsm), whereas lsmeans provides support for 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)

This is at least one way to achieve this. I'm still hoping someone knows of an approach using lsmeans only...

Another way to approach this is to hack into the lsmeans object, and manually replace the variance-covariance matrix prior to summary-ing the object.

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

JWilliman · Answer 2 · 2021-07-11T10:51:15.943

I'm not sure if this was possible using the 'lsmeans' package but it is using the updated emmeans package.


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

^{Created on 2021-07-08 by the reprex package (v0.3.0)}

Note, you can't just write the following

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

due to conflict with the sandwich::vcovHAC command which also has an adjust option. (I had incorrectly thought this was a bug).

You can use the argument `options = list(adjust = "none")` to insulate it from being passed elsewhere. And BTW **lsmeans** is now just a front end to **emmeans**. Even the `lsmeans()` function is in **emmeans** — Russ Lenth, Dec 09 '21 at 15:20

score 0 · Answer 3 · answered Dec 09 '21 at 09:22

OR use update to inject a custom vcov matrix into your emmeans/emmGrid object.

Example:

# 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)

How to obtain lsmeans() pairwise contrasts with custom vcov?

3 Answers3