Inverse probability weighting (IPW): standard errors after weighting observations

Question

When using propensity scores for inverse probability weighting (IPW) the standard errors for the parameters in the regression model may be affected. I have seen several examples of people using different types of standard errors (classical, robust, bootstrap) and am unsure which ones are correct to use and why. Classical weighting would use weights to indicate the precision of individual observation - this is not the case for IPW, where weighting indicates the importance of observations (but not their precision).

If you want to add references to R packages, that would be appreciated, but I am primarily interested in the methods and why they should or should not be used.

Answer 1

To provide a somewhat orthogonal answer to Noah's, let us see what's technically going on in the weighted least squares (WLS) case. We aim to solve $$ \theta_{WLS} = \arg\min_\theta\sum_{k=1}^N w_k (y_k - f(x_k;\theta))^2. $$ Computing the gradient, we note that the solution fulfills the weighted normal equation $$ (X^T W X) \theta_{WLS} = X^T W y, $$ and thus we have $$ \theta_{WLS} = (X^T W X)^{-1} X^T W y. $$ Recalling that $$\mathrm{Var}(Ax) = A \mathrm{Var}(x) A^T,$$ we observe that $$ \mathrm{Var}(\theta_{WLS}) = (X^T W X)^{-1} X^T W \mathrm{Var}(y) W^T X (X^T W^T X)^{-1}. $$

Now we can distinguish a few different cases.

• Basic ordinary least squares (OLS): Assume $W=I$ and $\mathrm{Var}(y)=\sigma_{\varepsilon}^2 I$. In this case, the formula above simplifies to $$\mathrm{Var}(\theta_{WLS}) = \mathrm{Var}(\theta_{OLS}) = \sigma_{\varepsilon}^2 (X^T X)^{-1}.$$
• Classic weighted least squares (WLS): Assume $W=\mathrm{Var}(y)^{-1}$, where $\mathrm{Var}(y)$ is assumed diagonal but heteroscedastic, i.e., the diagonal entries vary. In this case, we obtain $$\mathrm{Var}(\theta_{WLS}) = (X^T W X)^{-1}.$$
• Weighted least squares, but $W\neq\mathrm{Var}(y)^{-1}$: this includes the case of propensity weighting which you alluded to. In this situation, we simply have to content ourselves with the general variance equation given above: $$ \mathrm{Var}(\theta_{WLS}) = (X^T W X)^{-1} X^T W \mathrm{Var}(y) W^T X (X^T W^T X)^{-1}. $$

Answer 2

Lunceford and Davidian (2004) derive the asymptotic standard errors for IPW estimators. These rely on an M-estimation approach and assume the propensity scores are estimated using a method that can be represented as a system of estimating equations (e.g., logistic regression, but not random forests). Their proof also indicates that IPW estimators are smooth and asymptotically normal, making them amenable to bootstrapping. They also find that excluding the propensity score estimation from the estimating equations and treating the weights as fixed yields conservative estimates of the standard error (though this only applies to weights for the ATE). This is equivalent to using a robust standard error in the outcome regression model.

This leads to three ways to validly estimate the standard error in IPW:

1. Using M-estimation with the propensity scores and outcome model included together. This can be manually programmed using geex in R, and some R packages like WeightIt and PSweight can also compute them. In SAS, PROC CAUSALTRT automatically computes the correct standard errors, and in Stata, teffects ipw uses the same approach. These are asymptotically correct and can be used with any generalized linear model for the outcome.

2. Using a robust standard error for the outcome model. This will generally be conservative and is the simplest and most flexible approach because it can be used with weight-estimation methods that are not implemented in those packages or can't be represented as systems of estimating equations, like generalized boosted modeling, which is a somewhat popular method. To do this in R, you would use survey::vcovHC() after a glm() or lm() call with the outcome model, survey::svyglm(), or geepack::geeglm() as recommended by Hernán and Robins (2020). In SAS, you would use PROC SURVEYREG, and in Stata you would use supply the weights to the aweights argument in any regression model, which automatically requests robust standard errors.

3. Using the bootstrap. The bootstrap, where you include the propensity score estimation and effect estimation within each replication, is a very effective method because it does not rely on asymptotic arguments, can be used with weight-estimation methods that can't be or aren't implemented in the packages named above, and can be used for any estimand, regardless of whether analytical standard errors have been derived for them (e.g., for the rate ratio in a negative binomial outcome model). The difficulty is that one needs to know how to program a bootstrap and one needs to be prepared to wait a potentially long time when the estimation procedure takes a while, e.g., for some machine learning methods. Also, the bootstrap will tend to yield different estimates each time, adding an additional layer of uncertainty into the estimation. Using the bootstrap is easiest in R with the boot package.

I am most inclined to use M-estimation when available or robust standard errors because of their flexibility and ease of use. For a serious project where conservative standard errors could be a liability and I had a lot of time, I would use the bootstrap, which are typically the most accurate.