Is there a method for estimating heteroskedasticity-robust standard errors for a Heckman/heckit model in R?

Question

I am trying to figure out how to estimate a Type II Tobit/Heckman/heckit model in R and extract heteroskedasticity robust standard errors. It would be nice if the standard errors were robust for both stages of estimation, but the final one is my primary concern.

The standard package for estimating Heckman/heckit models in R seems to be sampleSelection, which I use below. The usual packages, sandwich and lmtest, do not have support for this type of model.

Here is the minimum reproducible example to estimate the model:

    library(sampleSelection)
    library(wooldridge)

    d <- wooldridge::mroz

    heck_fit <- heckit(
        selection = inlf ~ educ + nwifeinc + age + kidslt6, 
        outcome = lwage ~ educ + nwifeinc + age, 
        data = d, 
        method = "ml")

    summary(heck_fit)

I can extract the constant-variance standard errors for the model with vcov(heck_fit) and only for the outcome regression with vcov(heck_fit, part = "outcome"), but I do not know how to get heteroskedasticity robustness. I have played around with vcovHC() and coeftest() and they both fail.

Heteroskedasticity-robustness is never mentioned in the documentation for sampleSelection either. broom and many other regression analysis packages do not seem to support these models.

I'll note that heck_fit is equivalent to the following in Stata:

    heckman lwage educ nwifeinc age, select(educ nwifeinc age kidslt6)

The standard errors are equal to the 5th decimal place.

However, I would like to be estimating the R equivalent of this in Stata:

    heckman lwage educ nwifeinc age, select(educ nwifeinc age kidslt6) vce(robust)

If you are testing, here are the standard errors for the outcome equation.

The normal standard errors per Stata and R are educ = .017882, nwifeinc = .003657, and age = .004631.

Whereas, per Stata, the heteroskedasticity robust standard errors are educ = .0198926, nwifeinc = .0037428, age = .0049529.

Thank you for your help!

Answer 1

The sampleSelection package is built on the maxLik infrastructure which in turn supports the object-oriented approach of the sandwich package. Thus, you can use functions like sandwich() for "robust" cross-section sandwich covariances and vcovCL() for clustered sandwich covariances.

coeftest(heck_fit, vcov = sandwich)
## z test of coefficients:
## 
##               Estimate Std. Error z value  Pr(>|z|)    
## (Intercept)  0.1074434  0.4169059  0.2577 0.7966259    
## educ         0.1459250  0.0230264  6.3373 2.338e-10 ***
## nwifeinc    -0.0201104  0.0046496 -4.3252 1.524e-05 ***
## age         -0.0269064  0.0072367 -3.7180 0.0002008 ***
## kidslt6     -0.7045530  0.1468855 -4.7966 1.614e-06 ***
## (Intercept)  0.1386951  0.3282874  0.4225 0.6726742    
## educ         0.0604852  0.0198794  3.0426 0.0023454 ** 
## nwifeinc     0.0096951  0.0037403  2.5921 0.0095399 ** 
## age          0.0108875  0.0049496  2.1997 0.0278298 *  
## sigma        0.7897007  0.0760488 10.3841 < 2.2e-16 ***
## rho         -0.7345077  0.1283075 -5.7246 1.037e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This is very close to the Stata results. The differences are just do to slightly different degrees-of-freedom adjustments, see: https://stackoverflow.com/questions/27367974/different-robust-standard-errors-of-logit-regression-in-stata-and-r/27368468#27368468

If you use vcov = vcovCL (without indicating any actual clustering) instead of vcov = sandwich, then the results coincidentally exactly align between Stata and R.

The reason that vcovHC() cannot be applied is that vcovHC() only works for models with a single linear predictor and associated working residuals. See the discussion of Equations 7 and 8 at the end of Section 2 in vignette("sandwich-OOP", package = "sandwich").