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?

Answer 1

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)

Answer 2

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

Answer 3

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