Inverse probability weighting

Inverse probability weighting is a statistical technique for estimating quantities related to a population other than the one from which the data was collected. Study designs with a disparate sampling population and population of target inference (target population) are common in application.^[1] There may be prohibitive factors barring researchers from directly sampling from the target population such as cost, time, or ethical concerns.^[2] A solution to this problem is to use an alternate design strategy, e.g. stratified sampling. Weighting, when correctly applied, can potentially improve the efficiency and reduce the bias of unweighted estimators.

One very early weighted estimator is the Horvitz–Thompson estimator of the mean.^[3] When the sampling probability is known, from which the sampling population is drawn from the target population, then the inverse of this probability is used to weight the observations. This approach has been generalized to many aspects of statistics under various frameworks. In particular, there are weighted likelihoods, weighted estimating equations, and weighted probability densities from which a majority of statistics are derived. These applications codified the theory of other statistics and estimators such as marginal structural models, the standardized mortality ratio, and the EM algorithm for coarsened or aggregate data.

Inverse probability weighting is also used to account for missing data when subjects with missing data cannot be included in the primary analysis.^[4] With an estimate of the sampling probability, or the probability that the factor would be measured in another measurement, inverse probability weighting can be used to inflate the weight for subjects who are under-represented due to a large degree of missing data.

The inverse probability weighting estimator can be used to demonstrate causality when the researcher cannot conduct a controlled experiment but has observed data to model. Because it is assumed that the treatment is not randomly assigned, the goal is to estimate the counterfactual or potential outcome if all subjects in population were assigned either treatment.

We consider random variables ${\displaystyle (X,A,Y)\in \mathbb {R} ^{p}\times \{0,1\}\times \mathbb {R} }$ jointly distributed according to a law ${\displaystyle P}$ where

• ${\displaystyle X\in \mathbb {R} ^{p}}$ are the covariates
• ${\displaystyle A\in \{0,1\}}$ are the two possible treatments
• ${\displaystyle Y\in \mathbb {R} }$ is the response
• No assumptions such as random assignment of treatment are made.

Following Rubin's potential outcomes framework we also stipulate the existence of random variables ${\displaystyle Y^{\ast }(a)\in \mathbb {R} }$ for each ${\displaystyle a=0,1}$. Semantically, ${\displaystyle Y^{*}{\bigl (}a{\bigr )}}$ denotes the potential outcome that would be observed if the subject were assigned treatment ${\displaystyle a}$. Technically speaking, we actually work with the full joint distribution ${\displaystyle P^{\ast }}$ of ${\displaystyle (X,A,Y,Y^{\ast }(0),Y^{\ast }(1))}$; in that case ${\displaystyle P}$ is the marginal distribution for only the observed components of ${\displaystyle P^{\ast }}$. Special assumptions are needed to infer properties about ${\displaystyle P^{\ast }}$ using ${\displaystyle P}$, which will be detailed below.

Now suppose we have observations ${\displaystyle \{(X_{i},A_{i},Y_{i})\}_{i=1}^{n}}$ distributed identically and independently according to ${\displaystyle P}$. The goal is to use the observed data to estimate properties of the potential outcome ${\displaystyle Y^{\ast }(a)}$. For instance, we may wish to compare the mean outcome if all patients in the population were assigned either treatment: ${\displaystyle \mu _{a}=\mathbb {E} Y^{*}(a)}$. We want to estimate ${\displaystyle \mu _{a}}$ using observed data ${\displaystyle \{{\bigl (}X_{i},A_{i},Y_{i}{\bigr )}\}_{i=1}^{n}}$.

${\displaystyle {\hat {\mu }}_{a,n}^{IPWE}={\frac {1}{n}}\sum _{i=1}^{n}Y_{i}{\frac {\mathbf {1} _{A_{i}=a}}{{\hat {p}}_{n}(A_{i}|X_{i})}}}$

1. ${\displaystyle \mu _{a}=\mathbb {E} \left[{\frac {\mathbf {1} _{A=a}Y}{p(A|X)}}\right]}$ where ${\displaystyle p(a|x)={\frac {P(A=a,X=x)}{P(X=x)}}}$
2. construct ${\displaystyle {\hat {p}}_{n}(a|x)}$ or ${\displaystyle p(a|x)}$ using any propensity model (often a logistic regression model)
3. ${\displaystyle {\hat {\mu }}_{a,n}^{IPWE}=\sum _{i=1}^{n}{\frac {Y_{i}1_{A_{i}=a}}{n{\hat {p}}_{n}(A_{i}|X_{i})}}}$

With the mean of each treatment group computed, a statistical t-test or ANOVA test can be used to judge difference between group means and determine statistical significance of treatment effect.

Recall the full joint probability model ${\displaystyle (X,A,Y,Y^{\ast }(0),Y^{\ast }(1))\sim P^{\ast }}$ for the covariate ${\displaystyle X}$, action ${\displaystyle A}$, response ${\displaystyle Y}$, and potential outcomes ${\displaystyle Y^{\ast }(0),Y^{\ast }(1)}$. Recall also that ${\displaystyle P}$ is the marginal distribution of the observed data ${\displaystyle (X,A,Y)}$.

We make the following assumptions on ${\displaystyle P^{\ast }}$ relating the potential outcomes to the observed data. These allow us to infer properties of ${\displaystyle P^{\ast }}$ via ${\displaystyle P}$.

• (A1) Consistency: ${\displaystyle Y=\mathbf {1} \{A=0\}Y^{\ast }(0)+\mathbf {1} \{A=1\}Y^{\ast }(1)}$. So ${\displaystyle \mathbf {1} \{A=a\}Y=\mathbf {1} \{A=a\}Y^{\ast }(a)}$ for any ${\displaystyle a}$.
• (A2) No unmeasured confounders: ${\displaystyle \{Y^{*}(0),Y^{*}(1)\}\perp A\,|\,X}$. Formally, for any bounded, Borel-measurable functions ${\displaystyle f}$ and ${\displaystyle g}$, ${\displaystyle {\begin{aligned}\mathbb {E} _{(A,Y^{*}(a))}\left[f(Y^{*}(a))g(A)\,|\,X\right]=\mathbb {E} _{Y^{*}(a)}\left[f(Y^{*}(a))\,|\,X\right]\mathbb {E} _{A}\left[g(A)\,|\,X\right]\end{aligned}}}$for any ${\displaystyle a}$. This means that treatment assignment is based solely on covariate data and independent of potential outcomes.
• (A3) Positivity: ${\displaystyle P(A=a|X=x)>0}$ for all ${\displaystyle a}$ and ${\displaystyle x}$.

Under the assumptions (A1)-(A3), we will derive the following identities ^[5]

${\displaystyle {\begin{aligned}\mathbb {E} _{Y^{*}(a)}\left[Y^{*}(a)\right]=\mathbb {E} _{X}\left[\mathbb {E} _{Y}\left[Y\,|\,A=a,X\right]\right]=\mathbb {E} _{(X,A,Y)}\left[{\frac {\mathbf {1} \{A=a\}Y}{P(A|X)}}\right].\qquad \cdots \cdots (*)\end{aligned}}}$

The first equality follows from iterated expectation, (A3), (A2) and (A1), applied in that order:

${\displaystyle {\begin{aligned}\mathbb {E} _{Y^{*}(a)}\left[Y^{\ast }(a)\right]&=\mathbb {E} _{X}\left[\mathbb {E} _{Y^{*}(a)}[Y^{\ast }(a)\,|\,X]\right]\qquad \cdots \cdots {\text{ iterated expectation }}\\&=\mathbb {E} _{X}\left[\mathbb {E} _{Y^{*}(a)}[Y^{\ast }(a)\,|\,A=a,X]\right]\qquad \cdots \cdots {\text{ (A3) and (A2) }}\\&=\mathbb {E} _{X}\left[\mathbb {E} _{Y}\left[Y\,|\,A=a,X\right]\right]\qquad \cdots \cdots {\text{ (A1) }}\end{aligned}}}$

For the second equality, first note from the proof above that

${\displaystyle {\begin{aligned}\mathbb {E} _{X}\left[\mathbb {E} _{Y}\left[Y\,|\,A=a,X\right]\right]=\mathbb {E} _{X}\left[\mathbb {E} _{Y^{*}(a)}[Y^{\ast }(a)\,|\,X]\right].\end{aligned}}}$

Now by (A3), ${\displaystyle P(A\,|\,X)>0}$ almost surely. Furthermore, note that

${\displaystyle \mathbb {E} _{A}{\Big [}{\frac {\mathbf {1} \{A=a\}}{P(A\,|\,X)}}\,{\Big |}\,X{\Big ]}=\sum _{\alpha =0,1}{\frac {\mathbf {1} \{\alpha =a\}}{P(A=\alpha |X)}}P(A=\alpha |X)=1}$. Hence we can write

${\displaystyle {\begin{aligned}\mathbb {E} _{X}\left[\mathbb {E} _{Y^{*}(a)}[Y^{\ast }(a)\,|\,X]\right]&=\mathbb {E} _{X}{\Big [}\mathbb {E} _{A}{\Big [}{\frac {\mathbf {1} \{A=a\}}{P(A\,|\,X)}}\,{\Big |}\,X{\Big ]}\,\mathbb {E} _{Y^{*}(a)}[Y^{\ast }(a)\,|\,X]{\Big ]}\\&=\mathbb {E} _{X}{\Big [}\mathbb {E} _{(A,Y^{\ast }(a))}{\Big [}{\frac {\mathbf {1} \{A=a\}Y^{\ast }(a)}{P(A\,|\,X)}}\,{\Big |}\,X{\Big ]}{\Big ]}\qquad \cdots \cdots {\text{ (A2) }}\\&=\mathbb {E} _{X}{\Big [}\mathbb {E} _{(A,Y)}{\Big [}{\frac {\mathbf {1} \{A=a\}Y}{P(A\,|\,X)}}\,{\Big |}\,X{\Big ]}{\Big ]}\qquad \cdots \cdots {\text{ (A1) }}\\&=\mathbb {E} _{(X,A,Y)}\left[{\frac {\mathbf {1} \{A=a\}Y}{P(A|X)}}\right]\end{aligned}}}$

The Inverse Probability Weighted Estimator (IPWE) is known to be unstable if some estimated propensities are too close to 0 or 1. In such instances, the IPWE can be dominated by a small number of subjects with large weights. To address this issue, a smoothed IPW estimator using Rao-Blackwellization has been proposed, which reduces the variance of IPWE by up to 7-fold and helps protect the estimator from model misspecification. ^[6]

An alternative estimator is the augmented inverse probability weighted estimator (AIPWE) combines both the properties of the regression based estimator and the inverse probability weighted estimator. It is therefore a 'doubly robust' method in that it only requires either the propensity or outcome model to be correctly specified but not both. This method augments the IPWE to reduce variability and improve estimate efficiency. This model holds the same assumptions as the Inverse Probability Weighted Estimator (IPWE).^[7]

${\displaystyle {\begin{aligned}{\hat {\mu }}_{a,n}^{AIPWE}&={\frac {1}{n}}\sum _{i=1}^{n}{\Biggl (}{\frac {Y_{i}1_{A_{i}=a}}{{\hat {p}}_{n}(A_{i}|X_{i})}}-{\frac {1_{A_{i}=a}-{\hat {p}}_{n}(A_{i}|X_{i})}{{\hat {p}}_{n}(A_{i}|X_{i})}}{\hat {Q}}_{n}(X_{i},a){\Biggr )}\\&={\frac {1}{n}}\sum _{i=1}^{n}{\Biggl (}{\frac {1_{A_{i}=a}}{{\hat {p}}_{n}(A_{i}|X_{i})}}Y_{i}+(1-{\frac {1_{A_{i}=a}}{{\hat {p}}_{n}(A_{i}|X_{i})}}){\hat {Q}}_{n}(X_{i},a){\Biggr )}\\&={\frac {1}{n}}\sum _{i=1}^{n}{\Biggl (}{\hat {Q}}_{n}(X_{i},a){\Biggr )}+{\frac {1}{n}}\sum _{i=1}^{n}{\frac {1_{A_{i}=a}}{{\hat {p}}_{n}(A_{i}|X_{i})}}{\Biggl (}Y_{i}-{\hat {Q}}_{n}(X_{i},a){\Biggr )}\end{aligned}}}$

With the following notations:

1. ${\displaystyle 1_{A_{i}=a}}$ is an indicator function if subject i is part of treatment group a (or not).
2. Construct regression estimator ${\displaystyle {\hat {Q}}_{n}(x,a)}$ to predict outcome ${\displaystyle Y}$ based on covariates ${\displaystyle X}$ and treatment ${\displaystyle A}$, for some subject i. For example, using ordinary least squares regression.
3. Construct propensity (probability) estimate ${\displaystyle {\hat {p}}_{n}(A_{i}|X_{i})}$. For example, using logistic regression.
4. Combine in AIPWE to obtain ${\displaystyle {\hat {\mu }}_{a,n}^{AIPWE}}$

The later rearrangement of the formula helps reveal the underlying idea: our estimator is based on the average predicted outcome using the model (i.e.: ${\displaystyle {\frac {1}{n}}\sum _{i=1}^{n}{\Biggl (}{\hat {Q}}_{n}(X_{i},a){\Biggr )}}$). However, if the model is biased, then the residuals of the model will not be (in the full treatment group a) around 0. We can correct this potential bias by adding the extra term of the average residuals of the model (Q) from the true value of the outcome (Y) (i.e.: ${\displaystyle {\frac {1}{n}}\sum _{i=1}^{n}{\frac {1_{A_{i}=a}}{{\hat {p}}_{n}(A_{i}|X_{i})}}{\Biggl (}Y_{i}-{\hat {Q}}_{n}(X_{i},a){\Biggr )}}$). Because we have missing values of Y, we give weights to inflate the relative importance of each residual (these weights are based on the inverse propensity, a.k.a. probability, of seeing each subject observations) (see page 10 in ^[8]).

The "doubly robust" benefit of such an estimator comes from the fact that it's sufficient for one of the two models to be correctly specified, for the estimator to be unbiased (either ${\displaystyle {\hat {Q}}_{n}(X_{i},a)}$ or ${\displaystyle {\hat {p}}_{n}(A_{i}|X_{i})}$, or both). This is because if the outcome model is well specified then its residuals will be around 0 (regardless of the weights each residual will get). While if the model is biased, but the weighting model is well specified, then the bias will be well estimated (And corrected for) by the weighted average residuals.^[8][9][10]

The bias of the doubly robust estimators is called a second-order bias, and it depends on the product of the difference ${\displaystyle {\frac {1}{{\hat {p}}_{n}(A_{i}|X_{i})}}-{\frac {1}{{p}_{n}(A_{i}|X_{i})}}}$ and the difference ${\displaystyle {\hat {Q}}_{n}(X_{i},a)-Q_{n}(X_{i},a)}$. This property allows us, when having a "large enough" sample size, to lower the overall bias of doubly robust estimators by using machine learning estimators (instead of parametric models).^[11]

• Propensity score matching