# Variance Estimation for Least Squares with Probability Weights

I'm running a simulation study and finding that the nominal SEs of the estimated coefficients when using weights in lm in R are an underestimate of the simulation SE. I have confirmed that $$\hat{\beta}$$ is estimated correctly.

I believe that this is a well known phenomenon, and that there are many ways to make the SEs robust. One such method is called "HC3" and is estimated using a sandwich estimator that incorporates the residuals and the leverages. The weights that I am using are "Average Treatment Effect" weights that are derived from a propensity score model, which is estimated in each iteration of my simulation (because new data is generated each iteration). So it's not surprising that lm doesn't account for the fact that the weights aren't fixed.

My question: Given that I am using propensity score weights in OLS regression, does anyone know about a package in R that will help give me more conservative (hopefully similar to the simulation variance) SE estimates? I see there's a survey package but I am not sure if it applies.

A little more info, in case it helps; based on my findings, I don't think that what they are doing in figure 10 on page 6 of this paper produce robust SEs. I am guessing that these SEs are an underestimate (possibly a gross underestimate). https://pareonline.net/getvn.asp?v=20&n=13

Thanks!

You can use the survey package, but if you want more flexibility in how the standard errors are computed, you need to use the sandwich package, which implements all the HC SE estimators. If you use vcovHC() on the output of a call to lm(), it will produce an HC variance-covariance matrix of the coefficients, and you can choose the type you want (I believe HC3 is the default). To get the SEs you can take the square root of the diagonal elements.
The jtools package makes this much easier. Its primary function, summ() is a substitute for the base R function summary() for lm() and other regression objects. You can set robust = "HC3" in the call to summ() to produce output similar to summary() but with robust standard errors automatically computed.