How to estimate robust sandwich standard errors when estimating parameters using optim() in R?

Currently I am using numerical optimization in R via the optim() function to estimate some parameters in a complicated likelihood function. I know that optim can return the "Hessian" matrix which can be used to calculate model-based standard errors for my estimated parameters, however, I am wondering if it is possible to calculate robust-sandwich standard errors instead for my estimated parameters by purely using numerical optimization?

I am interested in using numerical optimization for estimation and inference. While an analytic solution to the robust sandwich error estimator may be available, it may not be possible to derive in all cases where I am using optim() for estimation. Is there some general way I can get robust standard errors in R by numerical optimization of a likelihood function? Or perhaps this is not possible and hence why optim only returns the Hessian.

It is in principle possible but potentially burdensome, especially if you have many observations. The reason is the following: For the full sandwich you need the "bread" for which the observed Hessian can be used (either computed analytically or numerically via optim() or numDeriv::hessian() etc.). And you also need the "meat" which in the simplest case is the outer product of gradients. See vignette("sandwich-OOP", package = "sandwich") for more details.
Thus, for the meat you not only need the complete gradient (which should be zero if the model is estimated well) but observation-wise contributions to the gradient. The extractor function estfun() hence returns an $n \times k$ matrix with the contributions of each of the $n$ observations to the gradient of each of the $k$ parameters. I think you should be able to set this up numerically, e.g., using numDeriv::grad(), but you would have to call it separately for each of observation, i.e., $n$ times.
With small sample sizes, the jackknife error estimator would be a better alternative. It too shares most of the robustness properties of the bootstrap (and sandwich), but rather than using Monte-Carlo methods (or exact analytic solutions), the jackknife requires only $n$ iterations, one for each omitted observation.
• Good point, thanks. In the sandwich package there is also the vcovBS() function which provides an object-oriented bootstrap covariance implementation. This can be easily applied but, of course, is computationally intensive. Commented Mar 9, 2018 at 20:29