1 Introduction
Inverse probability weight (IPW) estimators have been used across statistical literature in diverse forms: in survey sampling, in designing importance sampling in Monte Carlo techniques and in the context of average treatment effect estimation in causal inference. In survey sampling, the goal is often to estimate population mean ψ from a finite sample (${y_{1}},\dots ,{y_{n}})$, and a common approach is to weigh each observation ${y_{i}}$ by a weight ${w_{i}}$ inversely related to their probability of inclusion (probability proportional to selection, or PPS). For example, common PPS estimators admit the form $\hat{\psi }={\textstyle\sum _{i\in s}}{w_{i}}{y_{i}}/n$ or a ‘ratio’ estimator: $\hat{\psi }={\textstyle\sum _{i\in s}}{w_{i}}{y_{i}}/{\textstyle\sum _{i\in s}}{w_{i}}$, where s denotes the sample, and ${w_{i}}$’s denote the sampling weights. These weights or probabilities of inclusion might be known, or unknown depending on whether their source is sampling design or non-response. Two of the most popular examples of such estimators are the Horvitz–Thompson estimator (Horvitz and Thompson, 1952), which uses fixed weights to provide an unbiased estimate of the population mean, and the Hájek estimator (Hájek, 1971), which normalizes the weights to improve stability at the cost of introducing slight bias.
Similarly, in importance sampling, the goal is to estimate the evidence $\psi ={\mathbb{E}_{F}}\{l(X)\}=\textstyle\int l(x)dF(x)$ where F might be difficult to sample from and a common technique is to estimate ψ by drawing samples from a candidate distribution G (with density g) and calculate:
\[ \hat{\psi }={n^{-1}}{\sum \limits_{i=1}^{n}}l({x_{i}})f({x_{i}})/g({x_{i}})\doteq {n^{-1}}{\sum \limits_{i=1}^{n}}l({x_{i}})w({x_{i}}).\]
The candidate density $g(\cdot )$ is chosen such that the ‘weights’ $w({x_{i}})$ is nearly constant (Firth, 2011). The ‘ratio’ estimator analog in importance sampling would correspond to ${\textstyle\sum _{i=1}^{n}}l({x_{i}})w({x_{i}})/{\textstyle\sum _{i=1}^{n}}w({x_{i}})$ as developed by Hesterberg (1988, 1995), and further developed in (Firth, 2011). The unattainable ‘optimal’ choice of $w(\cdot )$ is of course $f(\cdot )$ itself, and a key insight in producing more accurate estimation is that self-normalization or biasing the sampler towards low probability regions can help. Methods such as nested sampling or vertical likelihood use a Lorenz curve re-ordering of summation to achieve this goal (Skilling, 2006; Chopin and Robert, 2010; Polson and Scott, 2014; Datta and Polson, 2025).A third important setting where inverse probability weighting plays a central role is causal inference. In the potential outcomes framework, IPW estimators form the backbone of average treatment effect (ATE) estimation by reweighting observed outcomes according to treatment assignment probabilities (Rubin, 1974; Imbens, 2004; Cunningham, 2021). We briefly discuss the ATE problem in subsection 2.3. IPW estimators provide a common design-based principle whether adjusting for unequal inclusion probabilities in surveys, stabilizing evidence estimation in Bayesian computation, or addressing confounding in observational studies. The wide applicability of IPW motivates our reexamination of their weaknesses and possible Bayesian resolutions.
In this paper, we first contrast and compare the popular inverse probability weight estimators—Horvitz–Thompson and Hájek—and discuss their relative merits and demerits in light of two weak paradoxes, viz., Basu’s circus example (Basu, 1988) and the Robins–Ritov–Wasserman example (Robins and Ritov, 1997; Wasserman, 2004). We then discuss a binning-and-smoothing estimator by Ghosh (2015) and a hierarchical Bayes estimator by Li (2010), in the light of the Wasserman (2004)’s example, and show how they can lead to a possible resolution. We also build connections between IPW estimators and evidence estimation techniques and argue that these innovative ideas from Monte Carlo methods can be exploited in designing IPW estimators to achieve a lower variance and higher stability. We argue that the issue of choice of weights – an optimal proposal $g(\cdot )$ in Monte Carlo or sampling weights – provides new insights on a long-standing controversy about a ‘weakness’ in Bayesian paradigm (e.g. Wasserman, 2004; Sims, 2010; Ghosh, 2015; Li, 2010).
There is a large, influential literature on statistical methods for improving survey-weighted estimates, particularly, in the context of complex survey designs. Our goal is not to attempt a review of these methodological advances in survey weight regularization or calibration, and we point the interested reader to the comprehensive review by Chen et al. (2017) and the references therein, e.g., Haziza and Beaumont (2017). Haziza and Beaumont (2017) provide a comprehensive review of methods for constructing survey weights, focusing on techniques to adjust for unequal probabilities of selection, nonresponse, and post-stratification to improve the representativeness and accuracy of survey estimates. Chen et al. (2017) discuss the limitations of basic design-based weights, derived from the inverse of inclusion probabilities, and propose modifications such as weight trimming, weight modeling, and incorporating weights into statistical models.
The structure of this article is as follows. In ğ2, we define the popular IPW estimators and recent developments, and their connections with importance sampling and ATE estimation. In ğ3, we discuss two popular examples of ‘weak paradoxes’ due to Basu, and Robins-Ritov-Wasserman. For the latter, we derive asymptotic properties of a Bayes estimator due to Li (2010) in ğ3.4. We compare different survey sampling estimators through a range of simulation studies in ğ4, to illustrate their relative merits and demerits. Finally, in ğ5, we discuss other areas of connection such as average treatment effect estimation in potential outcomes framework, and suggest a few possible future directions.
2 Inverse Probability Weight Estimators
Horvitz–Thompson and Hájek Estimators
We begin by defining two widely used estimators in survey sampling: the Horvitz–Thompson (HT) estimator (Horvitz and Thompson, 1952)1 and the Hájek estimator (Hájek, 1971) briefly. Consider a finite population $U=\{1,2,\dots ,N\}$ with values ${Y_{k}}$ for each unit $k\in U$. A sample $s\subset U$ is drawn according to a sampling design with known inclusion probabilities ${p_{k}}=P(k\in s)$.2 The Horvitz–Thompson estimator for the population mean $\psi ={N^{-1}}{\textstyle\sum _{k\in U}}{Y_{k}}$ is then
The Hájek estimator is an alternative procedure that estimates the population sum as well as N, i.e., it normalizes the estimated sum by the estimated total:
The HT estimator (2.2) is an unbiased estimator of ψ, while the Hájek estimator (2.3) is non-linear, and approximately unbiased estimator that does not require the knowledge of population size N, irrespective of whether N is known.
(2.2)
\[ {\hat{\psi }_{HT}}=\frac{1}{N}\sum \limits_{k\in s}\frac{{Y_{k}}}{{p_{k}}}.\hspace{1em}\text{(Horvitz--Thompson)}\](2.3)
\[ {\hat{\psi }_{\text{Hájek}}}=\frac{{\textstyle\sum _{k\in s}}{Y_{k}}/{p_{k}}}{{\textstyle\sum _{k\in s}}1/{p_{k}}}.\hspace{1em}\text{(Hájek)}\]Missing Data Framework and Adaptive Normalization
An equivalent framework arises in the context of missing data, i.e., when subjects are observed with nonuniform probabilities. We follow the notations of Khan and Ugander (2023) to describe it here. Here our goal is to estimate the mean ψ of ${Y_{1}},\dots ,{Y_{n}}$, with the observations missing at random. The Bernoulli indicators ${R_{k}},k=1,\dots ,n$ indicate whether ${Y_{k}},k=1,\dots ,n$ were observed or not. We assume ${R_{k}}\stackrel{\mathrm{ind}}{\sim }\text{Bernoulli}({p_{k}})$ for $k=1,\dots ,n$, which reflects non-response bias in sample surveys. Then, the Horvitz–Thompson and Hájek estimator of $\psi ={n^{-1}}{\textstyle\sum _{k=1}^{n}}{Y_{k}}$ are:
A notable generalization of the aforementioned estimators is the Trotter–Tukey estimator (Trotter and Tukey, 1956). Recently, Khan and Ugander (2023) rediscovered the Trotter–Tukey estimator as the adaptive normalization (AN) idea by letting the data choose the tuning parameter λ.
Furthermore, Khan and Ugander (2023) show that the ‘adaptive normalization’ idea dating back to Trotter and Tukey (1956) leads to a lower asymptotic variance than both while generalizing these estimators, in particular, $\mathbb{V}({\hat{\psi }_{AN}})$ is lower than both $\mathbb{V}({\hat{\psi }_{HT}})$ and $\mathbb{V}({\hat{\psi }_{\text{Hájek}}})$. Khan and Ugander (2023) also provide empirical evidence that their adaptive normalization scheme leads to a lower mean squared error of IPW estimators in different application areas in average treatment effect estimation and policy learning.
2.1 Horvitz-Thompson vs. Hájek Estimator
Horvitz–Thompson estimators possess many desirable properties: they are unbiased, admissible and consistent. Ramakrishnan (1973) provides a simple proof that the HT estimator is admissible in a class of all unbiased estimators of a finite population total. Isaki and Fuller (1982) proves that they achieve consistency under suitable conditions, such as, specifically requiring boundedness of population and weight sequences, bounded product of inclusion probabilities and second moments, and a variance that is $O({n_{t}^{-1}})$. Delevoye and Sävje (2020) provides conditions such as boundedness of moments for the outcome and inclusion probabilities and weak-design dependence. Delevoye and Sävje (2020) also point out that these conditions might be violated in ill-behaved setting with heavy-tailed outcomes or skewed sampling designs, common in natural experiments, where the practitioner has little control over the design.
It is worth noting here that IPW estimators, in particular the HT estimator, can be looked at as a weighted estimator, where the weights are not model-based, but design-based, as pointed out by several authors from Neyman (1934) to Little (2008). For example, in the HT estimator (2.2), the kth unit is assigned a weight proportional to the inverse of the selection probability as it ‘represents the $1/{p_{k}}$ units of the population’.
2.1.1 When Is the H’ajek Estimator Preferable?
We follow Särndal et al. (2003, ch. 5) for this discussion. First, note that the two estimators, HT and Hájek, can produce identical results under certain designs. However, if the population total N is unknown, only the Hájek estimator (2.3) can be used. On the other hand, if N is known—as is typically the case in the finite population estimation problem—the Hájek estimator “is usually the better estimator, despite estimation of an a priori known quantity” (Särndal et al., 2003). Särndal et al. (2003) supports this by identifying three situations where the weighted mean estimator (Hájek) outperforms the simple unbiased estimator (HT). In particular, Särndal et al. (2003, p. 183) shows that the ${\hat{\psi }_{\text{Hájek}}}$ estimator achieves lower variance than ${\hat{\psi }_{\text{HT}}}$ in each of these cases:
-
(i) When the values of ${y_{k}}$ are closer to the population mean ψ, a special case being fixed ${y_{k}}=c$ for $k=1,\dots ,N$.
-
(ii) When the sample size ${n_{s}}$ is variable with equal inclusion probabilities, e.g. the case where ${y_{k}}=c$ and ${p_{k}}={p_{0}}$ for $k=1,\dots ,N$, but the sample sizes vary.
-
(iii) Finally, when the inclusion probabilities ${p_{k}}$ are negatively correlated with the ${y_{k}}$ values, with large ${y_{k}}$ values corresponding to small ${p_{k}}$ values and vice versa. It is easy to see that the Hájek estimator is adaptable to high fluctuations in ${y_{k}}/{p_{k}}$ values that HT will suffer from, leading to high variability. Basu (1988) pointed this out in his famous circus example, recounted in ğ3.1.
Despite being popular choices, both the Horvitz–Thompson (HT) (Horvitz and Thompson, 1952), and the Hájek estimator (Hájek, 1971) have generated debates and controversies in the recent past. For example, the Horvitz–Thompson estimator where we encounter randomly missing observations and a very high-dimensional parameter space is presented as evidence of weakness of Bayesian method in this problem, in an example due to Wasserman (2004), based on Robins and Ritov (1997). In response, a few authors (Sims, 2010; Ghosh, 2015; Li, 2010; Harmeling and Touissant, 2007) have tried to provide a Bayesian answer by constructing Bayesian estimators that achieve a lower variance than the HT estimator at the cost of admitting a small (and in some cases, vanishing) bias: another example of the bias-variance trade-off.
Stein’s Paradox: Here, we argue that this phenomenon is another example of Stein’s paradox (James and Stein, 1961; Efron and Morris, 1977; Stigler, 1990) observed in a different context. The original Stein’s paradox showed the inadmissibility of the ordinary maximum likelihood estimator which is both unbiased and minimax for normal means in dimensions more than 3 by shrinking each component towards the origin (or a pre-fixed value), thereby borrowing strength from each other. The Stein’s paradox and the James–Stein shrinkage estimator has been truly transformative in statistics in both establishing Empirical Bayes’ success in high-dimensional inference and inspiring both frequentist shrinkage estimators and Bayesian shrinkage priors in such problems. We argue that the simple HT estimator can be improved upon by borrowing strength by moving from an independence framework to an exchangeable one, exhibiting Stein’s shrinkage phenomenon and effects of regularization.
Before we discuss the weak paradoxes, we briefly discuss how the connection between IPW estimators and Monte Carlo sampling, in particular importance sampling, and various improvements such as Riemann sums (Philippe, 1997; Philippe and Robert, 2001) or vertical likelihood integration that applies ‘binning and smoothing’ using a score-function heurism to choose the weight function (Polson and Scott, 2014; Madrid-Padilla et al., 2018). We then show how the idea of binning and smoothing also improves the HT estimator (Ghosh, 2015) in the apparent weakness example due to Wasserman (2004). We discuss these connections briefly in the next subsection and point out the similarities between estimators employed in survey sampling and Monte Carlo integration and their connections with statistical mechanics.
2.2 Importance Sampling
One way to look at this connection between sampling strategies and integral approximation is to represent the former as a missing data problem. Suppose that our goal is to estimate:
which can be thought of as a limiting value of ${\psi _{n}}={n^{-1}}{\textstyle\sum _{i=1}^{n}}{y_{i}}$ as $n\to \infty $, and can be approximated up to any degree of precision by ${\psi _{N}}={N^{-1}}{\textstyle\sum _{i=1}^{N}}{y_{i}}$ for a sufficiently large N. Given a random sample of size $n\ll N$, with sample probabilities π attached to ${y_{i}}$, we can view estimating θ as a problem of estimating the population quantity ${\psi _{N}}$ by ${\psi _{n}}$. The usual importance sampling estimator in this case is akin to using ${\mathbb{E}_{\pi }}[{\textstyle\sum _{i=1}^{n}}{y_{i}}/{\pi _{i}}]=\psi $, the usual HT estimator. Just like the HT estimator, the variance of importance sampling could blow up for poor choices of $g(\cdot )$, while the bias may shrink to zero. As we show below, these connections have been exploited by several authors to propose alternative importance sampling strategies.
Given a sample $({y_{1}},\dots ,{y_{n}})$, from either the density f corresponding to F itself or a suitable proposal density g, the usual importance sampling estimates ψ by the empirical average:
If the underlying probability measure is easy to sample from, i.e., if we can afford $f=g$, the above will reduce to the empirical average $\hat{\psi }={n^{-1}}{\textstyle\sum _{i=1}^{n}}l({y_{i}})$, which by the Law of Large Numbers, converge to the true value of ψ at $O({n^{-1}})$ rate. It is worth pointing out that replacing the empirical average for the naïve importance sampling estimator in (2.8) by a Riemann sum estimator provides a remarkable improvement in stability and convergence as demonstrated in (Philippe, 1997; Philippe and Robert, 2001) or (Yakowitz et al., 1978). Recently, Datta and Polson (2025) showed that combining Riemann sum estimators with nested sampling ideas can yield sharper convergence rates (up to $O({n^{-4}})$) for marginal likelihood estimation in Bayesian problems. We refer the reader to Datta and Polson (2025) for a more detailed discussion of vertical likelihood, and how it connects importance sampling with nested inference. This framework also offers a unifying perspective on several sampling-based approaches that aim to improve the stability of likelihood-weighted estimators.
(2.8)
\[ {\hat{\psi }_{IS}}=\frac{1}{n}{\sum \limits_{i=1}^{n}}\frac{l({y_{i}})f({y_{i}})}{g({y_{i}})}.\]2.2.1 A Bayesian Perspective
We conclude the discussion of evidence estimation with a pertinent and profound point raised by Diaconis (1988) in support of Bayesian approaches for probabilistic numerical integration. Diaconis (1988) starts with a algebraic function:
\[ f(x)=\exp \left\{\cosh \left(\frac{x+2{x^{2}}+\cos x}{3+\sin {x^{3}}}\right)\right\},\hspace{1em}f:[0,1]\mapsto \mathbb{R},\]
for which we want ${\textstyle\int _{0}^{1}}f(x)dx$ and asks “What does it mean to ‘know’ a function?.” In such a situation, where we might know a few properties of f and not others, a Bayesian approach seems natural, where one starts with a prior on $\mathcal{C}[0,1]$, the space of continuous functions, and estimate the integral using Bayes’ rule. Diaconis (1988) shows that Brownian motion, the ‘easiest’ prior on $\mathcal{C}[0,1]$, yields a linear interpolation leading to the trapezoid rule. Diaconis (1988) then shows that a host of well-known methods can be recovered as Bayesian estimates. This key question is revisited in (Owen, 2019) where he asks what does it mean to know the error $\Delta =|\hat{\psi }-\psi |$ and discusses advantages of Bayesian approaches over classical methods. Owen (2019) argues in favour of Bayes methods in difficult problems where extreme cost of function evaluation or skewness or heavy-tailed properties of $l(x)$ or unavailability of central limit theorem exposes the weakness of classical methods.2.3 Average Treatment Effect Estimation
A natural application of IPW estimator is average treatment effect (ATE) estimation in potential outcomes causal inference framework as pointed out in (Delevoye and Sävje, 2020; Khan and Ugander, 2023). We refer the readers to the excellent references in (Cunningham, 2021; Rubin, 1974; Imbens, 2004) for in-depth coverage and historical backgrounds. Here we measure the difference in potential outcomes observed over time points ${t_{2}}\gt {t_{1}}$, where only one of the two potential outcomes are observed for any unit. Using Khan and Ugander (2023)’s notation: we have the triplets $({Y_{k}}(0),{Y_{k}}(1),{p_{k}})$ for $k=1,\dots ,n$, where we observe ${Y_{k}}({I_{k}})$ for ${I_{k}}\sim \text{Bernoulli}({p_{k}})$. and the parameter of interest is $ATE=\tau =\mathbb{E}[{Y_{k}}(1)-{Y_{k}}(0)]$ from the observations available to us. The IPW estimators are employed here to estimate the two population means ${\psi _{1}}=\mathbb{E}[{Y_{k}}(1)]$ and ${\psi _{0}}=\mathbb{E}[{Y_{k}}(0)]$ separately and estimating τ by $\hat{\tau }=\hat{{\psi _{1}}}-\hat{{\psi _{0}}}$. Estimating the population means ${\psi _{j}}$, $j=1,2$ turns ATE into a survey sampling problem and one can use all the estimators: HT, Hájek or the various improvements such as the adaptive IPW estimator by Khan and Ugander (2023) based on the Trotter–Tukey idea (Trotter and Tukey, 1956) here. Khan and Ugander (2023) further develop an adaptive estimator by minimizing the variance of between-group differences and show that both the separate AIPW and the joint AIPW estimators achieve lower mean squared errors compared to the usual HT and Hájek.
3 Weak Paradoxes in Ratio Estimation Problem
3.1 Basu (1988)’s Example
Example 1.
Basu’s circus example: The famous circus example, due to (Basu, 1988), shows that the HT estimator (2.2) could lead to absurd estimates in some situations. Here, we imagine a circus owner, trying to estimate the total weight of his 50 elephants (say Y), picks a representative elephant from his herd (Sambo) and multiply his weight by 50. But, then a circus statistician, appalled by this estimate, devices a plan where Sambo is picked with probability $99/100$ and each of the remaining with $1/100$. Unfortunately, now the HT estimate is $\text{Sambo's weight}\times 100/99$, a serious underestimate, and moreover, if the owner picks Jumbo, the biggest in the herd, the HT estimate is $\text{Jumbo's weight}\times 4900$, an absurdity.
Discussing Basu (1988), Hájek (1971) points out that the HT estimator’s “usefulness is increased in connection with ratio estimation” and proposed an estimator in presence of auxiliary information ${A_{k}}$, related to ${Y_{k}}$ and with known total:
which would not be affected in this example like the standard HT estimator. The Hájek IPW estimate for the population mean in (2.3) can be derived from the special case ${A_{k}}\equiv 1,\forall k$.
(3.1)
\[ {\hat{Y}_{\text{Hájek}}}={\sum \limits_{k=1}^{n}}{A_{k}}\times \frac{{\textstyle\textstyle\sum _{k=1}^{n}}\frac{{Y_{k}}}{{p_{k}}}}{{\textstyle\textstyle\sum _{k=1}^{n}}\frac{{A_{k}}}{{p_{k}}}},\hspace{0.2778em}\text{where}\hspace{2.5pt}{Y_{k}}\propto {A_{k}},k=1,\dots ,n,\]Although the circus example was intended to be a pathological example, it provides at least two useful insights. First, weighted estimators can lead to nonsensical answers despite having nice large sample properties (Little, 2008). Second, problems of similar nature occur in importance sampling or Monte Carlo estimation of the marginal likelihood where empirical averages or unbiased estimator could have high or infinite variance (Li, 2010; Raftery et al., 2006). As pointed out in (Datta and Polson, 2025), strategies like Riemann sum (Philippe, 1997; Philippe and Robert, 2001), or careful choice of weight functions like vertical likelihood (Polson and Scott, 2014) or nested sampling (Skilling, 2006), or using ratio estimators (Firth, 2011), or adaptive normalization (Khan and Ugander, 2023) can resolve these issues. In particular, the same trick of ordering the draws and applying a Riemann-sum type approach is the key to avoid falling into examples like Basu’s circus.
3.2 Wasserman (2004)’s Example
We first re-state the example which is itself a simplification of an example from (Robins and Ritov, 1997), in a similar spirit. We consider IID samples $({Y_{i}},{X_{i}},{R_{i}})$, $i=1,\dots ,B$ such that ${Y_{i}}$’s are generated as a mixture of Bernoulli distributions with individual parameters indexed by the component label ${X_{i}}$, and the ‘missingness’ indicator ${R_{i}}$ denoting whether ${Y_{i}}$ was observed or not. Let the ‘success’ probabilities associated to ${R_{i}}$ be known constants ${p_{{X_{i}}}}$ satisfying:
Note that this strong condition on ${p_{j}}$’s in (3.2) ensures that the HT estimator will not lead to absurd answers like Basu’s paradox stated earlier. The hierarchical model for each draw $({Y_{i}},{X_{i}},{R_{i}})$, $i=1,\dots ,n$, is:
The parameter of interest is the average $\psi =(1/B){\textstyle\sum _{b=1}^{B}}{\theta _{b}}$. Wasserman (2004) argues that since the likelihood has little information on most ${\theta _{j}}$’s and the known constants B and ${p_{j}}$’s drop from the likelihood, Bayes’ estimates for ψ are going to be poor. On the other hand, the HT estimate:
is easily seen to be unbiased, given $\mathbb{E}({R_{i}}{Y_{i}}/{p_{{X_{i}}}})=\mathbb{E}\{\mathbb{E}({R_{i}}\mid {Y_{i}},{X_{i}})\cdot \mathbb{E}({Y_{i}}\mid {X_{i}})/{p_{{X_{i}}}}\}=\mathbb{E}(\mathbb{E}({Y_{i}}\mid {X_{i}}))=\psi $ since $\mathbb{E}({R_{i}}\mid {Y_{i}},{X_{i}})={p_{{X_{i}}}}$ by construction. The HT estimate will also satisfy $\mathbb{V}({\hat{\psi }_{HT}})\le 1/n{\delta ^{2}}$, using Hoeffding’s inequality. We would like to note here that the assumption (3.2) that the selection probabilities are between δ and $1-\delta $ is exploited explicitly in this asymptotic result, i.e., the HT estimator might have infinite variance as δ goes to zero, making the assumptions crucial.
(3.4)
\[\begin{aligned}{}[{R_{i}}& \mid {X_{i}}={x_{i}}]\sim \text{Bernoulli}({p_{{x_{i}}}})\end{aligned}\](3.5)
\[\begin{aligned}{}[{Y_{i}}& \mid {R_{i}},{X_{i}}={x_{i}}]\sim \left\{\begin{array}{l}\text{unobserved}\hspace{0.2778em}\text{if}\hspace{0.2778em}{R_{i}}=0\hspace{1em}\\ {} \text{Bernoulli}({\theta _{{x_{i}}}})\hspace{0.2778em}\text{if}\hspace{0.2778em}{R_{i}}=1.\hspace{1em}\end{array}\right.\end{aligned}\](3.6)
\[ {\hat{\psi }_{HT}}=\frac{1}{n}{\sum \limits_{i=1}^{n}}\frac{{R_{i}}{Y_{i}}}{{p_{{X_{i}}}}}\]3.3 Bayesian Solutions for Wasserman (2004)’s Problem
We review the existing Bayesian resolutions for the Robins-Ritov-Wasserman problem using Bayesian ideas and argue that, in some special cases, the improvement in accuracy of Bayes’ estimate is attained via exploiting the ‘borrowing strength’ phenomenon in Stein’s shrinkage.
3.3.1 Full Bayes Solution (Li, 2010)
Li (2010) provides a simple Bayes estimator by assuming that ${\theta _{1}},\dots ,{\theta _{B}}$ are exchangeable, not independent (see Fig. 1). Li (2010) estimates ψ by augmenting (3.5) by with Beta priors for ${\theta _{1}},\dots ,{\theta _{B}}$ and a further hyperprior on the mean of these Beta priors, as follows:
where, ψ is the mean of θ, and the shape parameters ${\alpha _{T}}$ and ${\alpha _{F}}$ control the width of range of $\boldsymbol{\theta }$ and the concentration of ψ.
Through extensive numerical simulation, Li (2010) shows that this estimator achieves a smaller mean squared error compared to the HT estimator in Wasserman (2004). Li (2010) also argues that the variance of HT estimator, given by:
will lead to inflated variance for small ${p_{b}}$ values, necessitating the bounds $\delta \le {p_{j}}\le 1-\delta $ in (3.2), but the variance for Li’s Bayes estimator remains unaffected and does not need this extra restriction. On the other hand, Li’s estimator (3.9) is prone to bias if the missingness mechanism ${p_{X}}$ and the parameter for observed outcomes ${\theta _{X}}$ are correlated, i.e., if the missingness is missing-at-random (MAR) and not missing-completely-at-random (MCAR). We shall formalize this via the derived approximation for the mean squared error of the Bayes estimator in our theoretical properties section. Linero (2024) argues that the ‘seemingly innocuous’ priors like the one used by Li (2010) encode what is effectively a-priori knowledge that the amount of selection bias is minimal.
(3.10)
\[ \mathbb{V}({\hat{\psi }_{HT}})=\frac{1}{n}\left\{\frac{1}{B}{\sum \limits_{b=1}^{B}}\frac{{\theta _{b}}}{{p_{b}}}-{\psi ^{2}}\right\},\]3.3.2 Binning and Smoothing (Ghosh, 2015)
Ghosh (2015) provides another estimator by reducing the dimension of $({p_{1}},\dots ,{p_{B}})$ by clubbing them into k ($k\ll B$) groups by utilizing the boundedness assumption (3.2). Based on this idea, we define the general binned-smoothed estimator as follows.
Definition 2.
Given a fixed $\delta \in (0,1)$ satisfying (3.2), we divide the whole range $(\delta ,1-\delta )\ni {p_{b}},b=1,\dots ,B$ into k sub-intervals $\{{\delta _{0}}=\delta ,{\delta _{1}},\dots ,{\delta _{k}}=1-\delta \}$ such that ${\delta _{i+1}}-{\delta _{i}}=(1-2\delta )/(k-1)$. We order the observed ${p_{{X_{i}}}}$’s into increasing order and define ${\tilde{p}_{j}}$ to be the mean of the ${p_{{X_{i}}}}$ values in the jth partition consisting of ${n_{j}}$ points, i.e.
where ${n_{j}}$ is the number of ${p_{{X_{i}}}}$’s falling into the jth class. Note that, if the sample size $n\gg k$, the number of partitions, the number of ${p_{{X_{i}}}}$’s in any interval, i.e. ${n_{j}}$ values would typically be non-zero, and in the unlikely case ${n_{j}}=0$ for any $j=1,\dots ,k$, the corresponding term will not contribute to the numerator in (3.11).
\[ {\tilde{p}_{j}}=\left\{\begin{array}{l@{\hskip10.0pt}l}\frac{1}{{n_{j}}}{\textstyle\sum _{i:{p_{{X_{i}}}}\in ({\delta _{j}},{\delta _{j+1}})}}{p_{{X_{i}}}},\hspace{1em}& \text{if}\hspace{0.2778em}{n_{j}}\gt 0\\ {} ({\delta _{j}}+{\delta _{j+1}})/2\hspace{1em}& \text{otherwise}\end{array}\right.\]
Then, the binned-smoothed estimator, based on (Ghosh, 2015) is:
(3.11)
\[ {\hat{\psi }_{BS:HT}}={\sum \limits_{j=1}^{k}}\frac{{n_{j}}}{n}\frac{1}{{n_{j}}}\frac{({\textstyle\textstyle\sum _{i=1}^{{n_{j}}}}{R_{i}}{Y_{i}})}{{\tilde{p}_{j}}}.\hspace{0.2778em}\]
Ghosh (2015)’s estimator performs at least as well as the HT estimator in small simulation studies. The idea of ‘binning and smoothing’ can be applied to other inverse probability weighted estimators such as Hájek (2.3). The binned-smoothed Hájek estimator would take the form:
where ${n_{j}}$ and ${\tilde{p}_{j}}$’s are as used in (3.11) before.
(3.12)
\[ {\hat{\psi }_{BS:Hajek}}=\frac{{\textstyle\textstyle\sum _{j=1}^{k}}\{({\textstyle\textstyle\sum _{i=1}^{{n_{j}}}}{R_{i}}{Y_{i}})/{\tilde{p}_{j}}\}}{{\textstyle\textstyle\sum _{j=1}^{k}}\{({\textstyle\textstyle\sum _{i=1}^{{n_{j}}}}{R_{i}})/{\tilde{p}_{j}}\}},\hspace{0.2778em}\]3.3.3 Bayesian Sieve (Sims, 2012)
We note that the idea of grouping the probabilities attached to the missingness indicators ${R_{i}}$’s, i.e., clubbing ${p_{{X_{i}}}}$’s into ${p_{j}}$’s was also proposed in (Sims, 2012). Sims (2012) argues that estimators better than the HT estimator (3.6) can be constructed using Bayesian approach if one acknowledges the existence of infinite dimensional unknown parameter in the model, e.g. assuming $p(b):1,\dots ,B\mapsto (0,1)$ is known but not the conditional distribution of $[\theta \mid p]$. Sims suggested to break up the range of the known $p(\cdot )$ into k segments, such that $\theta (b)$ and $p(b)$ are independent in each segment, and estimate the unknown $\theta (b)$ via a step function, with a constant for each segment. To complete Sim’s Bayesian sieve, one needs to estimate the probability distribution on ${p_{j}}$’s induced by the partition, which will converge to a Dirichlet distribution for a large sample, making analytical calculation for the posterior mean of ψ feasible.
3.3.4 Gaussian Likelihood (Harmeling and Touissant, 2007)
Harmeling and Touissant (2007) consider a slightly modified version of the Wasserman’s model (3.5). Instead of a Bernoulli, they consider:
\[ {X_{i}}\sim \mathcal{U}\{1,\dots ,B\},\hspace{0.2778em}{Y_{i}}\mid {\Theta _{{X_{i}}}}\sim \mathcal{N}({\theta _{{X_{i}}}},1),\hspace{0.2778em}1\le i\le n,\]
and $\mathcal{N}(\mu ,1)$ prior for each ${\theta _{b}}$, $b=1,\dots ,B$, and a ${\mathcal{N}_{\mu }}(0,\sigma )$ hyper-prior on μ. The maximum likelihood estimate and the posterior mean Bayes’ estimates for $B\to \infty $ are given by
\[\begin{aligned}{}{\hat{\psi }_{\text{MLE}}}& ={\bigg({\sum \limits_{i=1}^{n}}{R_{i}}\bigg)^{-1}}\bigg({\sum \limits_{i=1}^{n}}{R_{i}}{Y_{i}}\bigg),\hspace{1em}\text{and}\\ {} {\hat{\psi }_{\text{Bayes}}}& ={\bigg(2/\sigma +{\sum \limits_{i=1}^{n}}{R_{i}}\bigg)^{-1}}\bigg({\sum \limits_{i=1}^{n}}{R_{i}}{Y_{i}}\bigg),\end{aligned}\]
respectively. Although these estimators are not directly comparable to (3.9), Harmeling and Touissant (2007)’s likelihood-based estimators also achieve a lower MSE compared to the HT estimator, as expected.We also note that the idea of binning and smoothing ${p_{j}}$’s is also connected to importance sampling ideas in subsection 2.2, in particular, the trapezoidal rule (also called weighted Monte Carlo) by Yakowitz et al. (1978), or its generalizations, called the Riemann summation approach Philippe (1997); Philippe and Robert (2001). These trapezoidal rules based on ordered samples reduces the variance drastically and achieves faster convergence and better stability. Along these lines, a later development, called the nested sampling approach by Skilling (2006) also relies on “dividing the unit prior mass into tiny elements, and sorting them by likelihood.” For a comprehensive discussion of the different strategies for variance reduction for the evidence estimation problem, we refer the readers to (Polson and Scott, 2014; Datta and Polson, 2025).
3.4 Properties of Li’s Bayesian Estimator
While a closed form analytical expression for the variance of the Li’s estimator (3.9) is difficult, we can use the linearization strategy aka the delta theorem to derive an approximation for the mean and variance of the Bayes estimator (3.9). Doing so, we provide a proof and a closer look at an assertion in Ritov et al. (2014) that Li’s estimator is consistent only if $\mathbb{E}(Y\mid R=1)=\mathbb{E}(Y)$, and illustrate this phenomenon via simulation studies.
Recall that for a pair of random variables X, Y with mean ${\mu _{X}}$ and ${\mu _{Y}}$ and variances $\mathbb{V}(X)$, $\mathbb{V}(Y)$ and covariance $\mathrm{Cov}(X,Y)$, we have the following approximations:
To calculate the approximate mean and variance for the Li’s posterior mean estimator (3.9), we first calculate the mean and variances for the numerator and denominator separately as follows. Recall once again that the posterior mean estimator is given by ${\hat{\psi }_{Li}}=({\textstyle\sum _{i}}{R_{i}}{Y_{i}}+1)/({\textstyle\sum _{i}}{R_{i}}+2)$ assuming ${\alpha _{F}}=1$, i.e., a $\text{Beta}(1,1)$ or Uniform prior on ψ but the exact choice of shape parameter ${\alpha _{F}}$ does not influence the asymptotic nature of the variance. For notational convenience, we define the following quantities:
Hence, the expectation for the Bayes’ estimator can be approximated by taking the ratio of right hand sides from (3.15) and (3.16) as:
Similarly,
Finally, the covariance term is:
\[\begin{array}{l}\displaystyle \bar{\theta }=\frac{1}{B}{\sum \limits_{b=1}^{B}}{\theta _{b}}\doteq \psi ,\hspace{1em}\bar{p}=\frac{1}{B}{\sum \limits_{b=1}^{B}}{p_{b}},\hspace{1em}\bar{\theta .p}=\frac{1}{B}{\sum \limits_{b=1}^{B}}{\theta _{b}}{p_{b}},\\ {} \displaystyle {\sigma _{\theta ,p}}=\bar{\theta .p}-\bar{p}\bar{\theta }=\bar{\theta .p}-\bar{p}\psi .\end{array}\]
First, the expectation for the numerator and denominator follows from applying iterated expectations as: (3.15)
\[\begin{aligned}{}\mathbb{E}\bigg(\sum \limits_{i}{R_{i}}{Y_{i}}+1\bigg)& ={\mathbb{E}_{X}}\bigg\{{\sum \limits_{i=1}^{n}}\mathbb{E}({R_{i}}\mid {X_{i}})\mathbb{E}({Y_{i}}\mid {X_{i}})\bigg\}+1\\ {} & ={\mathbb{E}_{X}}\bigg\{{\sum \limits_{i=1}^{n}}{\theta _{{X_{i}}}}{p_{{X_{i}}}}\bigg\}+1\\ {} & =n\bar{\theta .p}+1=n\psi \bar{p}+n{\sigma _{\theta ,p}}+1\hspace{0.2778em}\end{aligned}\]
\[ \mathbb{E}({\hat{\psi }_{Li}})\approx \frac{n\psi \bar{p}+n{\sigma _{\theta ,p}}+1}{n\bar{p}+2}\to \psi +{\sigma _{\theta ,p}}/\bar{p},\hspace{0.2778em}\text{as}\hspace{0.2778em}n\to \infty ,\]
showing the ${\hat{\psi }_{Li}}$ is asymptotically unbiased if ${\sigma _{\theta ,p}}=0$, i.e., if ${\theta _{X}}\perp {p_{X}}\mid X$, or equivalently if $\mathbb{E}(Y\mid R=1)=\mathbb{E}(Y)$, as claimed by (Ritov et al., 2014). Now, the variances for the numerator and denominator can be calculated using the conditional variance identity. The expressions are given below.
(3.17)
\[\begin{aligned}{}\mathbb{V}\bigg(\hspace{-0.1667em}\sum \limits_{i}{R_{i}}\hspace{0.1667em}+\hspace{0.1667em}2\bigg)& \hspace{0.1667em}=\hspace{0.1667em}\sum \limits_{i}\mathbb{V}({R_{i}})\\ {} & \hspace{0.1667em}=\hspace{0.1667em}\sum \limits_{i}\mathbb{E}[\mathbb{V}({R_{i}}\mid {X_{i}}))]+\mathbb{V}[\mathbb{E}({R_{i}}\mid {X_{i}})]\\ {} & \hspace{0.1667em}=\hspace{0.1667em}\sum \limits_{i}[\mathbb{E}({p_{{x_{i}}}}(1-{p_{{X_{i}}}}))\hspace{0.1667em}+\hspace{0.1667em}\mathbb{V}({p_{{X_{i}}}})]\hspace{0.1667em}=\hspace{0.1667em}n(\bar{p}\hspace{0.1667em}-\hspace{0.1667em}{\bar{p}^{2}}).\end{aligned}\](3.18)
\[\begin{aligned}{}\mathbb{V}\bigg(\sum \limits_{i}{R_{i}}{Y_{i}}+1\bigg)& =\sum \limits_{i}\mathbb{V}({R_{i}}{Y_{i}})=n(\bar{\theta .p}-{\bar{\theta .p}^{2}}),\\ {} \hspace{1em}\text{where}\hspace{0.2778em}\bar{\theta .p}& =\frac{1}{B}{\sum \limits_{b=1}^{B}}{\theta _{b}}{p_{b}}.\end{aligned}\](3.19)
\[\begin{aligned}{}\mathrm{Cov}\left(\sum \limits_{i}{R_{i}}{Y_{i}}+1,\hspace{0.2778em}\sum \limits_{i}{R_{i}}+2\right)& =\mathrm{Cov}\left(\sum \limits_{i}{R_{i}}{Y_{i}},\hspace{0.2778em}\sum \limits_{i}{R_{i}}\right)\\ {} & =\sum \limits_{i}\mathrm{Cov}({R_{i}}{Y_{i}},\hspace{0.2778em}{R_{i}})\\ {} & =n\bar{\theta .p}\hspace{2.5pt}(1-\bar{p}).\end{aligned}\]Putting everything together from (3.15), (3.16), (3.18), (3.17) and (3.19), we get the following formula for approximate variance for the Bayes estimator as:
(3.20)
\[\begin{array}{cc}& \displaystyle \mathbb{V}({\hat{\psi }_{Li}})\approx {\left(\frac{n\bar{\theta .p}+1}{n\bar{p}+2}\right)^{2}}\times \\ {} & \displaystyle \left[\frac{n(\bar{\theta .p}-{\bar{\theta .p}^{2}})}{{(n\bar{\theta .p}+1)^{2}}}+\frac{n(\bar{p}-{\bar{p}^{2}})}{{(n\bar{p}+2)^{2}}}-\frac{2n\bar{\theta .p}(1-\bar{p})}{(n\bar{p}+2)(n\bar{\theta .p}+1)}\right].\end{array}\]A couple of immediate implications of the formula in (3.20) are as follows.
-
(i) Li’s Bayesian estimator is consistent if ${\theta _{{X_{i}}}}$ and ${p_{{X_{i}}}}$ are uncorrelated, as $\mathbb{E}({\hat{\psi }_{Li}})\to \psi $ and $\mathbb{V}({\hat{\psi }_{Li}})=O(1/n)\to 0$ as the sample size $n\to \infty $.
-
(ii) Unlike the variance of the HT estimator, given in (3.10), the variance of the Li’s Bayes estimator does not have the individual ${p_{b}}$ terms in the denominator and will not inflate for very small ${p_{b}}$ values. In other words, the restriction (3.2) is not needed for the Bayesian solution.
4 Simulation Examples
4.1 Comparing IPW Estimators for Wasserman’s Example
4.1.1 Missing Completely at Random
We compare the estimation performance for four different candidate estimators for the Robins-Ritov-Wasserman’s problem. The candidates are: the original Horvitz–Thompson (2.2), Hájek (2.3), the Li’s estimator, i.e., the Bayes posterior mean under a Beta hyperprior (3.9) and Ghosh (2015)’s binning and smoothing idea applied to the Hájek estimator. We choose the bounds for ${p_{i}}$’s $\delta =0.01$, and parameter space dimension $B=1000$, i.e., the ${p_{j}}$’s in our experiment are B equidistant grid-points in $[\delta ,1-\delta ]$. We vary the support of generative distribution for $\theta \in \mathcal{U}[a,b]$ to four different values, viz. $[a,b]\in \{[0.1,0.9],[0.1,0.4],[0.35,0.65],[0.6,0.9]\}$. Finally, we take sample size $n=100$, as was originally intended in the RRW example to make it analogous to a high-dimensional problem with $B\gg n$.
Table 1 and Fig. 2 shows the mean squared errors calculated over 100 replicates, and shows the following: (1) the Li’s estimator leads to the lowest mean squared error over replicates, (2) the Horvitz–Thompson estimator performs the worst across all situations considered and finally, (3) the Hájek estimator and the Binning-Smoothing estimator (Ghosh, 2015) achieves very similar performance and generally occupies a middle ground between the Bayes’ and the Horvitz–Thompson in terms of achieved mean-squared error.
Table 1
Mean squared error comparison between the four candiate estimators: Horvitz–Thompson (2.2), Hájek (2.3), the Li’s Bayesian estimator (3.9) and Ghosh (2015)’s estimator for different values of $[a,b]$, where $\theta \in [a,b]$. The numbers on the table are reported after multiplying the MSE by ${10^{2}}$ for better comparison and the column winner are denoted by boldfaced entries.
| [0.6, 0.9] | [0.1, 0.9] | |||
| method | mean | sd | mean | sd |
| Bayes’ (Li’s) | 0.37198 | 0.05093 | 0.47541 | 0.05980 |
| Binning-smoothing | 0.63991 | 0.08583 | 0.85225 | 0.10204 |
| Hájek | 0.71787 | 0.11643 | 0.95824 | 0.12818 |
| HT | 3.09007 | 0.83348 | 2.27093 | 0.72187 |
| [0.1, 0.4] | [0.35, 0.65] | |||
| method | mean | sd | mean | sd |
| Bayes’ (Li’s) | 0.35743 | 0.04758 | 0.47329 | 0.06272 |
| Binning-smoothing | 0.64379 | 0.08114 | 0.86459 | 0.10523 |
| Hájek | 0.73526 | 0.13494 | 0.96710 | 0.13523 |
| HT | 1.17543 | 0.51219 | 2.22213 | 0.69317 |
Figure 2
Mean square error comparison for Horvitz-Thompson, Hájek and the Li’s estimator (3.9) for different generative distributions of ${\theta _{1}},\dots ,{\theta _{B}}$, and fixed values of n, B and δ. In each case, Bayes’ estimator has the lowest MSE, followed by Hájek, and HT has the largest MSE.
For a single replication with $1,000$ evaluations, Fig. 3 shows the bias and variance for the four candidate estimators. Fig. 3 demonstrates the nature of variance reduction by Bayes’ estimator and the Binning-smoothing idea in estimating $\psi =\mathbb{E}(\boldsymbol{\theta })$ for different ranges of θ, without any significant increase in bias.
4.1.2 Missing at Random
As discussed earlier, Li’s estimator can be biased when there is a non-zero correlation between ${\theta _{X}}$ and ${p_{X}}$, e.g., when the missingness might arise due to confounding. We show an example here to illustrate the effect of this situation, and compare the four candidates. We take $\delta =0.1$, i.e., ${p_{b}}$’s to be equidistant points in $[0.1,0.9]$, and instead of ${\theta _{b}}$’s to be uniformly distributed in $[a,b]$, we take:
The remaining set-up is same as before. Clearly, in this case the true value of $\psi =\mathbb{E}(\theta )=0.5$. Figure 4a shows that the Li’s estimator has an upward bias, as expected, because of the positive correlation between θ and p, but the Hájek and Binned-Smoothed estimators are unaffected. On the other hand, the Binned-Smoothed is significantly better than all the other candidates (HT, Hájek and Li’s) in terms of MSE (Fig. 4b). This shows that the Li’s posterior mean estimator can be biased when there is confounding and one needs to be careful when using it. The binned-smoothed method, on the other hand, performs well irrespective of the correlation between ${\theta _{X}}$ and ${p_{X}}$, i.e., for both MAR and MCAR-type missingness.
5 Discussion
Our main goal in this paper was to shed light on various aspects of the family of inverse probability weight estimators, discuss its weakness and strengths and highlight the connections between survey sampling and Monte Carlo integration via these popular tool pervasive in Statistical literature. For survey sampling, we consider some of the ‘weak paradoxes’ (Ghosh, 2015) that arise from using inverse probability estimators using the well-known examples from Wasserman (2004) and Basu (1988), and review the merits and demerits of popular IPW estimators. In particular, we show that the hierarchical Bayes’ estimator due to (Li, 2010) leads to robustness against pathological situations but admits bias in presence of confounding. We provide sufficient conditions for consistency of the Li’s posterior mean estimate in Wasserman’s example and show that, under certain conditions, both the Li’s estimator and the binning-smoothing idea (Ghosh, 2015) achieves lower variation in mean squared errors. We then highlight the analogous tools in Monte Carlo integration, revisiting a few earlier works (e.g., Hesterberg, 1988), where the Horvitz–Thompson estimator is likened to the naïve importance sampling, and much like survey sampling, self-normalized weights or using control variates lead to better estimators. We conclude with a brief discussion of applications in the context of causal inference, a key use of IPW estimators, and a few possible directions for future work.
A possible future direction, borrowing from (Kong et al., 2003) is to use the general semiparametric models and the associated estimators or computational algorithm for $k\gt 1$ in the context of survey sampling that would allow one to combine data from one or more surveys that use different but known inclusion probabilities. We conjecture that this might lead to improved estimators for average treatment effect estimation, which we have briefly described in section 2.3.
In the context of causal inference, a second possible direction for future research is comparing a suitable modification of the Bayes estimator in (3.9) for the ATE problem with that of the adaptive estimators developed by Khan and Ugander (2023) and both empirically and theoretically to investigate consistency properties. It will be also worthwhile to consider the semiparametric model in (Kong et al., 2003) for ATE and compare the results of simultaneous estimation using MLE with these candidates.
Finally, the problems presented in (Wasserman, 2004) and discussed in (Sims, 2007, 2010) are inherently high-dimensional in nature, where sparsity is pervasive. There is a large and growing literature on promoting sparsity in causal inference, e.g., in presence of a large number of confounders or baseline variables while estimating the effect of an exposure on an outcome. We refer the readers to Shortreed and Ertefaie (2017); Wang et al. (2012, 2015); Kim et al. (2022) for recent advances in this area. In our simple focal example, the high-dimensional parameter $\boldsymbol{\theta }=({\theta _{1}},\dots ,{\theta _{B}})$ or $\mathbf{p}=({p_{1}},\dots ,{p_{B}})$ could exhibit sparsity with a small ${\ell _{0}}$ norm or other notions of sparsity. To handle sparsity in higher dimensions while maintaining tail-robustness and accuracy, the state-of-the-art Bayesian solution would be augmenting (3.5) with global-local shrinkage priors (Polson and Scott, 2010, 2012; Bhadra et al., 2019), i.e.,
\[\begin{aligned}{}{\theta _{b}}\mid {\alpha _{T}},\psi & \stackrel{\mathrm{ind}}{\sim }{f_{\text{HS}}}(\theta ;{\alpha _{T}}\psi ,\hspace{0.2778em}{\alpha _{T}}(1-\psi ),\hspace{0.2778em}\tau )\\ {} {f_{\text{HS}}}(\theta ;a,b,\tau )& ={\theta ^{a-1}}{(1-\theta )^{b-1}}{\{1-(1-\theta ){\tau ^{2}}\}^{-(a+b)}},\\ {} (\psi ,{\alpha _{T}}\mid {\alpha _{F}})& \sim \text{Beta}({\alpha _{F}},{\alpha _{F}})\times \pi ({\alpha _{T}}),\\ {} & \hspace{0.2778em}\theta ,\psi \in [0,1].\end{aligned}\]
Here, as before, ψ is the mean of θ, ${\alpha _{T}}$ and ${\alpha _{F}}$ are the and the shape parameters for $\boldsymbol{\theta }$ and ψ respectively and τ is a global shrinkage parameter adjusting to sparsity. Such Bayesian regularization priors have been used successfully for treatment effect estimation with more control variates than observations (Hahn et al., 2018). Designing built-in sparsity priors that can also incorporate selection mechanism for confounder selection is an interesting problem that we plan to address in a future endeavor.
