Setup

library(pacman); p_load(sandwich, lmtest, jtools, huxtable)

Rationale

Heteroskedasticity-robust inference in finite samples has been primarily approached with HC estimators. Unfortunately, in R, there is not an automatic robust regression function for HCs 4 and 5. Below, I have provided a function for HC1-5 using sandwich and lmtest.

The Function

HCReg <- function(Formula, Data, HC = "HC3"){
  require(sandwich); require(lmtest)
  BaseModel <- lm(Formula, Data)
  VC <- vcovHC(BaseModel, type = HC)
  Output <- coeftest(BaseModel, vcov = VC)
  return(Output)}

HCCheck <- function(Formula, Data){
  require(sandwich); require(lmtest); require(jtools); require(huxtable)
  BaseModel <- lm(Formula, Data)
  VC1 <- vcovHC(BaseModel, type = "HC1")
  VC2 <- vcovHC(BaseModel, type = "HC2")
  VC3 <- vcovHC(BaseModel, type = "HC3")
  VC4 <- vcovHC(BaseModel, type = "HC4")
  VC5 <- vcovHC(BaseModel, type = "HC5")
  Out1 <- coeftest(BaseModel, vcov = VC1)
  Out2 <- coeftest(BaseModel, vcov = VC2)
  Out3 <- coeftest(BaseModel, vcov = VC3)
  Out4 <- coeftest(BaseModel, vcov = VC4)
  Out5 <- coeftest(BaseModel, vcov = VC5)
  ModTable <- export_summs(Out1, Out2, Out3, Out4, Out5, model.names = c("HC1", "HC2", "HC3", "HC4", "HC5"),
             statistics = character(0),
             error_format = "{conf.low}, {conf.high}, {std.error}",
             digits = getOption("jtools-digits", 4))
  return(ModTable)}

Some Tests

summary(lm(Cold ~ Clothes_RF + Self_Obj * Clothes_RF + Temp_True, data))

## 
## Call:
## lm(formula = Cold ~ Clothes_RF + Self_Obj * Clothes_RF + Temp_True, 
##     data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.93764 -0.74787 -0.04511  0.76645  2.81914 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)   
## (Intercept)          4.50114    1.69073   2.662  0.00846 **
## Clothes_RF           0.49921    0.20794   2.401  0.01738 * 
## Self_Obj             0.33919    0.18622   1.821  0.07018 . 
## Temp_True           -0.04859    0.02764  -1.758  0.08043 . 
## Clothes_RF:Self_Obj -0.11332    0.05368  -2.111  0.03616 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.337 on 181 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.0604, Adjusted R-squared:  0.03964 
## F-statistic: 2.909 on 4 and 181 DF,  p-value: 0.02301

HCReg(Cold ~ Clothes_RF + Self_Obj * Clothes_RF + Temp_True, data, "HC3")

## 
## t test of coefficients:
## 
##                      Estimate Std. Error t value Pr(>|t|)   
## (Intercept)          4.501140   1.497128  3.0065 0.003018 **
## Clothes_RF           0.499206   0.210091  2.3761 0.018539 * 
## Self_Obj             0.339192   0.174309  1.9459 0.053213 . 
## Temp_True           -0.048591   0.026503 -1.8334 0.068379 . 
## Clothes_RF:Self_Obj -0.113318   0.055473 -2.0428 0.042525 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

HCReg(Cold ~ Clothes_RF + Self_Obj * Clothes_RF + Temp_True, data, "HC4")

## 
## t test of coefficients:
## 
##                      Estimate Std. Error t value Pr(>|t|)   
## (Intercept)          4.501140   1.488531  3.0239 0.002858 **
## Clothes_RF           0.499206   0.215168  2.3201 0.021453 * 
## Self_Obj             0.339192   0.175349  1.9344 0.054626 . 
## Temp_True           -0.048591   0.026324 -1.8459 0.066539 . 
## Clothes_RF:Self_Obj -0.113318   0.056714 -1.9981 0.047207 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

HCReg(Cold ~ Clothes_RF + Self_Obj * Clothes_RF + Temp_True, data, "HC5")

## 
## t test of coefficients:
## 
##                      Estimate Std. Error t value Pr(>|t|)   
## (Intercept)          4.501140   1.469425  3.0632 0.002524 **
## Clothes_RF           0.499206   0.206005  2.4233 0.016366 * 
## Self_Obj             0.339192   0.170746  1.9865 0.048483 * 
## Temp_True           -0.048591   0.026035 -1.8664 0.063605 . 
## Clothes_RF:Self_Obj -0.113318   0.054383 -2.0837 0.038590 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

HCCheck(Cold ~ Clothes_RF + Self_Obj * Clothes_RF + Temp_True, data)

## Registered S3 methods overwritten by 'broom':
##   method            from  
##   tidy.glht         jtools
##   tidy.summary.glht jtools

## Original model not retained as part of coeftest object. For additional model summary information (r.squared, df, etc.), consider passing `glance.coeftest()` an object where the underlying model has been saved, i.e.`lmtest::coeftest(..., save = TRUE)`.
## This message is displayed once per session.

	HC1	HC2	HC3	HC4	HC5
(Intercept)	4.5011 **	4.5011 **	4.5011 **	4.5011 **	4.5011 **
	1.5978, 7.4045, 1.4714	1.5927, 7.4096, 1.4740	1.5471, 7.4552, 1.4971	1.5640, 7.4382, 1.4885	1.6017, 7.4005, 1.4694
Clothes_RF	0.4992 *	0.4992 *	0.4992 *	0.4992 *	0.4992 *
	0.1037, 0.8947, 0.2004	0.0971, 0.9013, 0.2038	0.0847, 0.9137, 0.2101	0.0746, 0.9238, 0.2152	0.0927, 0.9057, 0.2060
Self_Obj	0.3392 *	0.3392 *	0.3392	0.3392	0.3392 *
	0.0061, 0.6723, 0.1688	0.0031, 0.6753, 0.1703	-0.0047, 0.6831, 0.1743	-0.0068, 0.6852, 0.1753	0.0023, 0.6761, 0.1707
Temp_True	-0.0486	-0.0486	-0.0486	-0.0486	-0.0486
	-0.1001, 0.0029, 0.0261	-0.1001, 0.0030, 0.0261	-0.1009, 0.0037, 0.0265	-0.1005, 0.0033, 0.0263	-0.1000, 0.0028, 0.0260
Clothes_RF:Self_Obj	-0.1133 *	-0.1133 *	-0.1133 *	-0.1133 *	-0.1133 *
	-0.2179, -0.0088, 0.0530	-0.2195, -0.0071, 0.0538	-0.2228, -0.0039, 0.0555	-0.2252, -0.0014, 0.0567	-0.2206, -0.0060, 0.0544
* p < 0.001; p < 0.01; * p < 0.05.

It is also possible to make a graphical model comparison if the desired coefficients to be plotted are specified. This wrapper could also be extended for non-linear models.

Robust Regressions

Setup

Rationale

The Function

Some Tests