Shape-constrained partial identification of a population mean under unknown probabilities of sample selection

TL;DR

This paper introduces a method to improve bounds on population mean estimates under unknown sample selection probabilities by incorporating shape constraints like symmetry or log-concavity, demonstrated through an educational test case.

Contribution

It develops a novel approach that leverages auxiliary shape constraints to tighten bounds on population means with unknown sampling probabilities.

Findings

Tighter bounds achieved using shape constraints.

Application to estimating Aymara students' test performance.

Implementation in R package 'scbounds'.

Abstract

A prevailing challenge in the biomedical and social sciences is to estimate a population mean from a sample obtained with unknown selection probabilities. Using a well-known ratio estimator, Aronow and Lee (2013) proposed a method for partial identification of the mean by allowing the unknown selection probabilities to vary arbitrarily between two fixed extreme values. In this paper, we show how to leverage auxiliary shape constraints on the population outcome distribution, such as symmetry or log-concavity, to obtain tighter bounds on the population mean. We use this method to estimate the performance of Aymara students---an ethnic minority in the north of Chile---in a national educational standardized test. We implement this method in the new statistical software package scbounds for R.

Equations

$$\hat{\mu}^{*}=\sum_{i=1}^{n}\pi_{i}^{-1}Y_{i}\,\Big{/}\,\sum_{i=1}^{n}\pi_{i}^{-1}.$$

$${\hat{\mu }}_{AL}^{+}=\sup \left\{\sum _{i=1}^n{w}_i{Y}_i:\sum _{i=1}^n{w}_i=1,\frac{\max \left\{w_i\right\}}{\min \left\{w_i\right\}}\le \gamma \right\}$$

$${\hat{F}}_w(y):=\sum _i{w}_{i}1\left\{Y_i\le y\right\}$$

$$\widehat{F}^{*}(y)=\sum_{i=1}^{n}\pi^{-1}_{i}\,1\!\left\{i:Y_{i}\leq y\right\}\,\Big{/}\,\sum_{i=1}^{n}\pi^{-1}_{i}=\sum_{i=1}^{n}w^{*}_{i}1\!\left\{i:Y_{i}\leq y\right\}.$$

$${\hat{\mu }}^{+}=\sup \left\{\sum _{i=1}^n{w}_i{Y}_i:\sum _{i=1}^n{w}_i=1,\frac{\max \left\{w_i\right\}}{\min \left\{w_i\right\}}\le \gamma ,{\hat{F}}_w\in {\mathcal{C}}_{\gamma ,n}\right\},$$

$${\rho }_{\alpha ,n}:=\mathbb{P}\left[\sqrt{n}\underset{y\in \mathbb{R}}{\sup }\left|{\hat{F}}^{\ast }(y)-F(y)\right|\ge \sqrt{\frac{{\sigma }_{\gamma }^2\left(1+\gamma \right)\left(1+{\gamma }^{-1}\right)\log \left({\alpha }^{-1}\right)}{2}}\right],$$

$$\limsup_{\alpha\rightarrow 0}\left\{\limsup_{n\rightarrow\infty}\ \log\left(\rho_{\alpha,\,n}\right)\,\Big{/}\,\log\left(\alpha\right)\right\}\geq 1.$$

$$\begin{array}{rl}& {\mathcal{C}}_{\gamma ,n}\left(\Theta \right):={\cup }_{\theta \in \Theta }\left(H:\underset{y\in \mathbb{R}}{\sup }\left|H(y)-F_{\theta }(y)\right|\le {\delta }_{\gamma ,n}\right),\\ & {\delta }_{\gamma ,n}:=\sqrt{\frac{{\sigma }_{\gamma }^2\left(1+\gamma \right)\left(1+{\gamma }^{-1}\right)\log \left(n\right)}{4n}}.\end{array}$$

$${\hat{\mu }}_{\theta }^{+}=\sup \left\{\sum _{i=1}^n{w}_i{Y}_i:\sum _{i=1}^n{w}_i=1,\frac{\max \left\{w_i\right\}}{\min \left\{w_i\right\}}\le \gamma ,\underset{y\in \mathbb{R}}{\sup }\left|{\hat{F}}_w(y)-F_{\theta }(y)\right|\le {\delta }_{\gamma ,n}\right\},$$

$${\mathcal{C}}_{\gamma ,n}\left(\Theta \right):={\cup }_{\theta \in \Theta }\left(H:\underset{y\in \mathbb{R}}{\sup }\left|H(y)-F_{\theta }(y)\right|\le {\delta }_{\gamma ,n}+{\delta }^{\ast }\right).$$

$$\begin{array}{rl}& \underset{n\to \infty }{\lim }\mathbb{P}\left[\sqrt{n}\underset{q\in [0,1]}{\sup }\left|{\hat{F}}^{\ast }\left(F^{-1}(q)\right)+{\hat{F}}^{\ast }\left(F^{-1}(1-q)\right)-1\right|\ge {\zeta }_{\gamma ,\alpha }\right]\le \alpha ,\\ & {\zeta }_{\gamma ,\alpha }:={\Phi }^{-1}\left(1-\frac{\alpha }{4}\right)\sqrt{\frac{\left(1+\gamma \right)\left(1+{\gamma }^{-1}\right)}{4}}\end{array}$$

$$\begin{array}{rl}& {\hat{\mu }}_m^{+}=\sup \left\{\sum _{i=1}^n{w}_i{Y}_i:\sum _{i=1}^n{w}_i=1,\frac{\max \left\{w_i\right\}}{\min \left\{w_i\right\}}\le \gamma ,{\hat{F}}_w\in {\mathcal{C}}_{\gamma ,n}^{SYM}\right\}\\ & {\mathcal{C}}_{\gamma ,n}^{SYM}=\bigcup _{m\in \mathbb{R}}\left\{H:\underset{y}{\sup }\left|H(m+y)+H(m-y)-1\right|\le {\zeta }_{\gamma ,n^{-1/2}}\right\}.\end{array}$$

$$\widehat{U}_{m}(y)=\left(\gamma\widehat{S}_{KS}(y)-\left(\gamma-1\right)\widehat{S}_{KS}\left(\min\left\{y,\,m\right\}\right)\right)\,\Big{/}\,\left(\gamma-\left(\gamma-1\right)\widehat{S}_{KS}(m)\right).$$

$${\hat{L}}_m(y)=\mathrm{argmin}\left\{\int L(y)dy:\log (L(y))\text{ is concave, and }L(y)\ge {\hat{U}}_m(y)\text{ for all }y\in \mathbb{R}\right\}.$$

$$C_{\gamma ,n}^{LC+}=\bigcup _{m\in \mathbb{R}}\left\{H:H(y)\ge {\hat{L}}_m(y),y\in \mathbb{R}\right\}.$$

$${\hat{\mu }}^{+}=\underset{w}{\sup }\left\{\sum _{i=1}^n{w}_i{Y}_i:\sum _{i=1}^n{w}_i=1,\frac{\max \left\{w_i\right\}}{\min \left\{w_i\right\}}\le \gamma ,{\hat{F}}_w(y)\in C_{\gamma ,n}^{LC+}\right\}.$$

$$N\operatorname{Var}\left[\hat{H}(y)\right]=\omega^{2}y,\ \ \ \ \ \hat{H}(y):=\frac{1}{N}\sum_{\left\{i:Y_{i}\leq y\right\}}\frac{1}{\pi_{i}}\,\Big{/}\,\mathbb{E}_{obs}\left[\frac{1}{\pi_{i}}\right],$$

$$N\operatorname{Var}\left[\hat{H}(y)\right]=\mathbb{E}_{obs}\left[\frac{1\!\left\{Y_{i}\leq y\right\}}{\pi_{i}^{2}}\right]\,\Big{/}\,\mathbb{E}_{obs}\left[\frac{1}{\pi_{i}}\right]^{2}.$$

$$\sqrt{n}\left({\hat{F}}^{\ast }(\cdot )-F(\cdot )\right)\Rightarrow G(y):=\omega \left(W(y)-F(y)W(1)\right),$$

$$\underset{u\to \infty }{\lim }\frac{1}{u^2}\log \left(\mathbb{P}\left[\underset{t\in [0,1]}{\sup }\left|X(t)\right|>u\right]\right)\le \frac{-1}{2{\sigma }_{\gamma }^2}.$$

$$\limsup_{\alpha\rightarrow 0}\ \ \log\left(\mathbb{P}\left[\sup_{t\in[0,\,1]}\left\lvert X(t)\right\rvert>\omega\sigma_{\gamma}\sqrt{2\log{\alpha^{-1}}}\right]^{-1}\right)\,\bigg{/}\,\log\left(\alpha^{-1}\right)\geq 1.$$

$$\mathbb{E}_{obs}\left[{\pi_{i}}^{-1}\right]=\mathbb{E}_{pop}\left[\pi_{i}\right]^{-1},\ \ \mathbb{E}_{obs}\left[{\pi_{i}^{-2}}\right]=\mathbb{E}_{pop}\left[{\pi_{i}}^{-1}\right]\,\big{/}\,\mathbb{E}_{pop}\left[\pi_{i}\right].$$

$${\omega }^2={\mathbb{E}}_{pop}\left[{{\pi }_i}^{-1}\right]{\mathbb{E}}_{pop}\left[{\pi }_i\right]\le \underset{a\in [0,1]}{\sup }\left(1+a\left(\gamma -1\right)\right)\left(1+a\left({\gamma }^{-1}-1\right)\right),$$

$$\sqrt{n}\left({\hat{F}}^{\ast }(y)-F(y)\right)\Rightarrow \omega \left(W(y)-F(y)W(1)\right),$$

$${\omega }^{-1}\sqrt{n}$$

$$=W(F^{-1}(q))-(W(1)-W(F^{-1}(1-q)))\frac{\text{buildrel}d}{=}W(F^{-1}(q)+1-F^{-1}(1-q))),$$

$$\underset{n\to \infty }{\lim }\mathbb{P}\left[\sqrt{n}\underset{y\in [0,1]}{\sup }\left({\hat{F}}^{\ast }(F^{-1}(q))+{\hat{F}}^{\ast }(F^{-1}(1-q))-1\right)\ge \omega t\right]=\mathbb{P}\left[\underset{y\in [0,1]}{\sup }W(y)\ge t\right]=2\Phi (-t),$$

$$S(m)\,\big{/}\,(S(m)+\gamma(1-S(m)))=r,$$

$$F(y)\geq\left(S(y)+(\gamma-1)\left(S(y)-S(m)\right)_{+}\right)\,\big{/}\,\left(S(m)+\gamma(1-S(m))\right).$$

$$\underset{u\to \infty }{\lim }\frac{1}{u^2}\log \left(\mathbb{P}\left[\underset{t\in [0,1]}{\sup }\left\{\left|X(t)\right|\right\}>u\right]\right)=\frac{-1}{2{\sup }_{t\in [0,1]}\left\{\mathrm{Var}\left[X(t)\right]\right\}};$$

$$X(t)=B(t)+W(1)(t-F(t)),$$

Full text

Keywords: Partial identification; Sensitivity analysis; Survey sampling; Survey weights.

1 Introduction

A common challenge in the biomedical and social sciences is to estimate a population mean from a sample obtained with unknown probabilities of sample selection. This is often the case when drawing inferences about mobile populations, such as the homeless and hospital outpatients, as well as with hard-to-reach populations, such as injection drug users and some ethnic minorities. In general, this problem arises when a sampling frame is unavailable or unreliable, and when there is no or limited information about the sampling design.

In brief, the estimation problem can be formalized as follows. Let $\mathcal{P}$ denote a potentially infinite population, and let $F$ denote the cumulative distribution function of our outcome of interest $Y$ over $\mathcal{P}$. Our goal is to estimate the population mean $\mu=\mathbb{E}_{F}(Y)$. To do so, we have access to a random sample $\mathcal{S}=\left\{Y_{i}\right\}$ of size $n$ obtained via biased sampling. Concretely, we can imagine that $\mathcal{S}$ was generated using an accept reject scheme as follows: until we have $n$ observations, repeatedly draw independent and identically distributed pairs $(Y,\,\pi)\in\mathbb{R}\times(0,\,1]$ where $Y\sim F$, and then add $Y$ to our sample $\mathcal{S}$ with probability $\pi$. Whenever the inverse sampling probabilities $\pi_{i}^{-1}$ are correlated with $Y_{i}$, the sample mean will be an inconsistent estimator for the population mean.

If these sampling probabilities $\pi_{i}$ for our $n$ sampled observations were known, then we would have access to the following ratio estimator that is consistent for $\mu$ under weak conditions (Hájek 1958; Cochran 1977)

[TABLE]

Here, however, we are interested in the setting where the sampling weights $\pi_{i}$ are unknown. In a recent advance, Aronow and Lee (2013) showed that it is possible to obtain meaningful identification intervals for $\mu$ in the sense of, e.g., Manski (2003), even if all we have is bounds on the sampling weights $\pi_{i}$. Suppose that we know that $\max\left\{\pi_{i}\right\}/\min\left\{\pi_{i}\right\}\leq\gamma$ for some constant $\gamma<\infty$. This gives an asymptotically consistent identification interval $\mathcal{I}_{AL}:=[\hat{\mu}_{AL}^{-},\,\hat{\mu}_{AL}^{+}]$ for $\mu$, where

[TABLE]

and $\hat{\mu}_{AL}^{-}$ is the $\inf$ over the same set. They also develop an efficient algorithm for computing these bounds.

While this is a creative and fertile approach that can help us get identification intervals for $\mu$ under weak assumptions, the Aronow-Lee bounds can be unnecessarily pessimistic in many applications. To understand the root of this pessimism, it is helpful to consider their method as first estimating the true population $F$ as

[TABLE]

where $w_{i}$ are the maximizing or minimizing weights in (2), and then setting the limits of the interval as $\mathbb{E}_{\widehat{F}_{w}}(Y)$ for the two resulting extreme sets of weights. The problem is that the population distributions implied by these extreme weights are often rather implausible in practice.

Specifically, the weights $w_{i}$ induced by the optimization problem (2) correspond to a step function depending on whether or not $Y_{i}$ falls below some threshold, and so the weighted empirical distribution functions $\widehat{F}_{AL}^{+}$ and $\widehat{F}_{AL}^{-}$have a sharp change in slope at that threshold, as illustrated in Figure 1(a) below. This threshold can also be interpreted as a substantial discontinuity in an associated density, provided we are willing to posit the existence of such a density (See Figure 1(b)). Such sharp elbows in the estimated cumulative distribution function often contradict expert knowledge about what the true population distribution $F$ should look like. For example, physical measurements (such as height and weight in some populations) in the biological and medical sciences often exhibit a bell-shaped distribution, as do stock returns and other indicators in finance and mechanical error measurements in industry.

This paper studies how to use such auxiliary information about the shape of the population outcome distribution $F$ to get shorter identification intervals for $\mu$ by ruling out “implausible” weightings in order to tighten the resulting identification bounds. We allow for various types of specifications for $F$, such as parametric assumptions, and shape constraints based on symmetry or log-concavity. In general, the more we are willing to assume about $F$, the shorter the resulting identification intervals. At one extreme, if we know that $F$ is Gaussian, then we can substantially shorten identification intervals, while if we make weaker assumptions, e.g., that $F$ has a log-concave density, then we get smaller but still noticeable improvements over $\mathcal{I}_{AL}$. We focus on the situation when $F$ has real-valued support; when $F$ has categorical (or binary) support, it is less common to have access to plausible shape constraints. This paper relates to the literature on biased sampling, empirical likelihood, and exponential tilting (see, for example, Owen (1988), Efron and Tibshirani (1996), De Carvalho and Davison (2014), and Fithian and Wager (2015)). We implement these methods in the new statistical software package scbounds for R.

2 Tighter identification bounds via shape constraints

As discussed above, our goal is to use existing information about the population distribution $F$, e.g., that $F$ is symmetric or log-concave, to obtain tighter identification bounds for $\mu=\mathbb{E}_{F}(Y)$. Operationally, we seek to encode such shape information about $F$ into constraints that can be added to the optimization problem (2). Throughout the paper, we assume that $\mu$ is in fact well defined and finite.

Our analysis focuses on the weighted empirical distribution function $\widehat{F}^{*}$ induced by the oracle ratio estimator $\hat{\mu}^{*}$ (1) that has access to the true sampling probabilities $\pi_{i}$ that give corresponding oracle weights $w_{i}^{*}$:

[TABLE]

As shown in the result below, any shape constraint we make on $F$ that lets us control the behavior of $\widehat{F}^{*}$ induces an asymptotically consistent identification interval for the population mean $\mu$. The intuition here is that if we can construct sets of distribution functions that are as constrained as possible, while still containing our oracle $\hat{F}^{*}$ with high probability, then the optimization problem will contain the oracle estimate $\hat{\mu}^{*}$, giving short but still consistent bounds.

Theorem 2.1**.**

Suppose that we have access to auxiliary information on $F$ that lets us construct sets of distribution functions $\mathcal{C}_{\gamma,\,n}$ with the property that $\lim_{n\rightarrow\infty}\mathbb{P}[\widehat{F}^{*}\in\mathcal{C}_{\gamma,\,n}]=1$. Then, if we write

[TABLE]

and $\hat{\mu}^{-}$ as the infimum in the analogous optimization problem, the resulting identification interval $\mathcal{I}:=[\hat{\mu}^{-},\,\hat{\mu}^{+}]$ is asymptotically valid in the sense that $\Delta(\mu,\mathcal{I})\rightarrow_{p}0$, where $\Delta(\mu,\mathcal{I})$ is the distance between $\mu$ and the nearest point in the interval $\mathcal{I}$.

We note that the resulting intervals $\mathcal{I}$ can never be wider than the intervals of Aronow and Lee (2013), because the optimization problem (4) has strictly more constraints than the original optimization problem (2).

At a high level, this result shows that if we have any auxiliary information about $F$, then the identification bounds of Aronow and Lee (2013) are needlessly long. However, the above result is of course rather abstract, and cannot directly guide practical data analysis. First of all, it leaves open the problem of how to turn shape constraints on $F$ into plausibility sets $\mathcal{C}_{\gamma,\,n}$ that contain $\widehat{F}^{*}$ with high probability. Second, any guarantee of the above form is not useful if we cannot solve the optimization problem (4) in practice. Our next concern is to address these issues given specific side-information about $F$.

2.1 Identification bounds in parametric families

Although our main goal is to provide inference under shape constraints on $F$, we begin by considering the parametric case of $F=F_{\theta}$ for some $\theta\in\Theta$, as the parametric setting allows us to construct particularly simple plausibility sets $\mathcal{C}_{\gamma,\,n}$. Our approach is built around a Kolmogorov-Smirnov type concentration bound for ratio estimators. Our proof relies on finding the worst-case weighted distribution with $\pi_{\max}/\pi_{\min}\leq\gamma$ in terms of the characterization of Marcus and Shepp (1972) for the tails of Gaussian processes. (See Appendix.)

Lemma 2.2**.**

Suppose that we have a population and sampling scheme for which $\pi_{\min}\leq\pi_{i}\leq\pi_{\max}$ with $\pi_{\max}/\pi_{\min}\leq\gamma$. Then, defining the tail probability of $\widehat{F}^{*}(\cdot)$ deviating far from the true $F(\cdot)$ as

[TABLE]

where $\sigma_{\gamma}^{2}\leq 1$ is a constant defined as the maximum of a concave function specified in the proof, we have

[TABLE]

In other words, for large $n$, the limiting probability of $\rho_{\alpha,\,n}$ is bounded by $\alpha^{1+o(1)}$.

Setting $\alpha=1/\sqrt{n}$, the above bound on $\widehat{F}^{*}$ suggests using as our plausibility set the union of all possible sets for all candidate $F(\cdot)$:

[TABLE]

We see that, regardless of the true parameter value $\theta\in\Theta$, we have $\hat{F}^{*}\in\mathcal{C}_{\gamma,\,n}(\Theta)$ with probability tending to 1, and so we immediately have:

Corollary 2.3**.**

Suppose that, under the conditions of Lemma 2.2, we know that $F=F_{\theta}$ for some $\theta\in\Theta$, and set $\mathcal{C}_{\gamma,\,n}\left(\Theta\right)$ as in (6). Then, (4) provides an asymptotically valid identification interval for $\mu$.

Furthermore, we can also check that the resulting intervals are asymptotically sharp.

Operationally, we implement this procedure by first solving the optimization problem

[TABLE]

over a grid of candidate values $\theta\in\Theta$, and then set $\hat{\mu}^{+}=\sup_{\theta}\hat{\mu}^{+}_{\theta}$. The problem (7) is a fractional programming problem that can be solved as a linear program using standard optimization methods; see, e.g., section 4.3 of Boyd and Vandenberghe (2004).

2.2 Relaxation of sharp shape constraints

Consider the parametric family case, above. Lemma 2.2 implies that $\widehat{F}^{*}(\cdot)$ will grow arbitrarily close to some $F_{\theta}$ with probability 1. This can be a very strong assumption: under a Gaussian family, for example, this implies that, with increasing $n$, not only are our bounds sharp, but they will collapse to a single point as the $\mathcal{C}_{\gamma,\,n}\left(\Theta\right)$ shrinks due to the tightening of $\delta_{\gamma,\,n}$. If we instead impose a shape constraint of $\max_{y}\left|F(y)-F_{\delta}(y)\right|\leq\delta^{*}$, for some $\delta^{*}$ and unknown $\theta\in\Theta$, we can expand our set $\mathcal{C}_{\gamma,\,n}\left(\Theta\right)$ to

[TABLE]

Due to the triangle inequality on the KS-distances, we still have, in the limit, $\widehat{F}^{*}$ in our set with probability 1, and therefore valid bounds on our mean. This relaxation allows for restricting the shape of our unknown distribution to be near a given parametric family without imposing a strong parametric assumption. The $\delta^{*}$ is a sensitivity parameter, and practitioners could examine how the bounds respond to different choices. We, however, instead investigate methods that make general shape constraint assumptions instead of these parametric ones.

2.3 Identification bounds with symmetry

We now move back to our main topic of interest, i.e., how to leverage shape constraints on $F$ to obtain improved identification bounds for $\mu$. The difference between parametric versus shape-constrained side information about $F$ is that it is not usually practical to do a grid search over all distributions in some shape-constrained class, and so the algorithm based on (7) does not generalize. Rather, for each candidate shape-constrained class, we may need to find an ad-hoc way to avoid a full grid search.

First, we consider the case where $F$ is symmetric, i.e., there is some value $m\in\mathbb{R}$ such that $F(m+y)=1-F(m-y)$ for all $y\in\mathbb{R}$. Such symmetry constraints interface particularly nicely with our approach. In order to make use of such constraints, we begin by establishing the following analogue to Lemma 2.2. Note that, unlike Lemma 2.2, this result holds for any value of $\alpha>0$ and not only for “small” $\alpha\rightarrow 0$; from a technical perspective, this result follows directly from Donsker arguments used to prove the classical Kolmogorov-Smirnov theorem, and does not require the additional machinery of Marcus and Shepp (1972).

Lemma 2.4**.**

Suppose that we have a population and sampling scheme for which $\pi_{\min}\leq\pi_{i}\leq\pi_{\max}$ for some $\pi_{\max}/\pi_{\min}\leq\gamma$. Then, for any $\alpha>0$,

[TABLE]

To draw a connection to the previous result, we note that $\Phi^{-1}\left(1-\frac{\alpha}{4}\right)\asymp\sqrt{2\log\alpha^{-1}}$ for small values of $\alpha$.

Now, if the distribution $F(\cdot)$ is symmetric around a point $m$, then any pair $(F^{-1}(q),\,F^{-1}(1-q))$ with $q\in[0,\,1]$ can be written as $(m-y,\,m+y)$ for some $y\in\mathbb{R}$. Thus, in the case of symmetric distributions, (8) immediately provides a tail bound on the supremum of $\widehat{F}^{*}(m+y)+\widehat{F}^{*}(m-y)-1$ over $y\in\mathbb{R}$, and suggests using the following estimator:

[TABLE]

The lower bound $\hat{\mu}^{-}$ is computed analogously. This algorithm thus enables us to leverage symmetry constraints while only performing a grid search over a single parameter $m$, i.e., the center of symmetry. These identification intervals are again asymptotically valid and sharp. They can also be relaxed in a similar manner to the parametric approach described above.

2.4 Identification bounds with log-concavity

Finally, we consider the case where $F$ is known to have a log-concave density. Imposing log-concavity constraints appears to be a promising, light-weight method for encoding side information about $F$: the class of log-concave distributions is quite flexible, including most widely used parametric distributions with continuous support, while enforcing regularity properties such as unimodality and exponentially decaying tails (Walther 2009).

Unlike in the case of symmetry, there does not appear to be a simple way to turn log-concavity constraints into asymptotically sharp identification bounds for $\mu$ using only linear programming paired with a low-dimensional grid search. Below, we detail our procedure for obtaining $\hat{\mu}^{+}$; to obtain $\hat{\mu}^{-}$ we can apply the same procedure to $-Y_{i}$.

Let $\widehat{S}(y)=n^{-1}\sum_{i}1\!\left\{Y_{i}\leq y\right\}$ and $\widehat{S}_{KS}(y)=\max\{\widehat{S}(y)-D_{KS}(1-1/\sqrt{n})/\sqrt{n},\,0\}$, where $D_{KS}(\cdot)$ denotes the Kolmogorov-Smirnov cumulative distribution function. By the Kolmogorov-Smirnov theorem, we know that, with probability tending to 1, $S(y)\geq\widehat{S}_{KS}(y)$ for all $y\in\mathbb{R}$, where $S(y)$ is the distribution function of the observed (i.e., biased) sample. 2. 2.

Next, in the proof, we show that for some $m\in\mathbb{R}$, the following is a lower bound for the population distribution of interest, $F(y)$ with probability tending to 1:

[TABLE] 3. 3.

We are now ready to use the fact that our distribution is log-concave. It is well known that if $F$ has a log-concave density, then the function $\log(F(y))$ must itself also be concave (Prékopa 1973). Thanks to this fact, we know that if $F(y)\geq\widehat{U}_{m}(y)$, then also $F(y)\geq\widehat{L}_{m}(y)$ where

[TABLE] 4. 4.

Finally, we define $C_{\gamma,\,n}$ as the set of distributions satisfying at least one of these lower bounds:

[TABLE]

Given this construction, we can obtain an upper endpoint for our identification interval as usual,

[TABLE]

The following result shows that $C_{\gamma,\,n}^{LC+}$ does in fact contain the population sampling distribution with probability tending to 1, and so Theorem 2.1 establishes the validity of our identification intervals. Unlike in the parametric or symmetric cases, our log-concave identification bounds are not asymptotically sharp (i.e., they may not converge to the shortest possible identification interval given our assumptions about log-concavity and bounded sampling ratios); however, they still provide tighter intervals than the baseline results of Aronow and Lee (2013).

Lemma 2.5**.**

Suppose that we have a population for which $\pi_{\min}\leq\pi_{i}\leq\pi_{\max}$ for some $\pi_{\max}/\pi_{\min}\leq\gamma$, and that $F$ has a log-concave density. Then, $\mathbb{P}[\widehat{F}^{*}\in\mathcal{C}_{\gamma,\,n}^{LC+}]=1$.

3 Application: sampling ethnic minorities

The Aymara are an indigenous population of the Andean plateau of South America. At the present, they live predominantly in Bolivia and Peru, and only a small proportion of them live in the north of Argentina and Chile. In Chile, they constitute a minority of nearly 50,000 in a country of approximately 18 million. Across the world, it is of great importance to understand how ethnic minorities fare in order to design effective affirmative action policies. Here, we use the proposed method to bound the average performance of the Aymara students in the national standardized test held in Chile for admission to higher education. This test is called PSU (for Prueba de Selección Universitaria) and nearly 90% of enrolled high school student take it every year; however, this figure is known to be lower in vulnerable populations such as the Aymara in northern Chile.

Using the sample of 847 Aymara students that took the PSU in mathematics in 2008, we seek identification intervals for the population mean counterfactual test score had everyone taken the test. We assume that the sampling ratio is bounded by $\max\left\{\pi_{i}\right\}\,/\,\min\left\{\pi_{i}\right\}\leq\gamma=9$, and consider two inferences, one under the assumptions that the population test score distribution is symmetric and one that it is log-concave. We also consider the approach of Aronow and Lee (2013) that does not use any shape constraints.

Here, the observed test scores have a mean of 502 with a sample standard deviation of 104. Given $\gamma=9$, we obtain population identification intervals of $(426,\,578)$ assuming symmetry, $(414,\,589)$ assuming log-concavity, and $(410,\,591)$ without any constraints. Thus, in this example, assuming symmetry buys us shorter identification intervals than assuming log-concavity.

Figure 1 depicts the weighted distribution functions $\widehat{F}_{w}(\cdot)$ underlying the upper endpoints of all 3 identification intervals. The sharp threshold of the weights $w_{i}$ resulting from the unconstrained method of Aronow and Lee (2013) is readily apparent. Assuming either symmetry or log-concavity of the population sampling distribution yields more regular-looking distributions.

Appendix: Proofs

Proof of Theorem 2.1

First, mirroring the argument of Aronow and Lee (2013), we note that $\hat{F}^{*}$ is itself an estimator of the form $\hat{F}^{*}(y)=\sum_{i}w_{i}1\left(Y_{i}\leq y\right)$ with $\sum_{i}w_{i}=1$ and $\max{w_{i}}/\min{w_{i}}\leq\gamma$, for some vector of weights $w$. Further, if $\hat{F}^{*}\in\mathcal{C}_{\gamma,\,n}$ then $\mu^{*}\in\mathcal{I}$. Now, if $\hat{F}^{*}\in\mathcal{C}_{\gamma,\,n}$ with probability tending to 1, then $\Delta(\mu,\,\mathcal{I})$ is asymptotically stochastically dominated by $|\hat{\mu}^{*}-\mu|$, which itself converges in probability to 0 by the weak law of large numbers.

Proof of Lemma 2.2

In order to derive this type of uniform tail bound, we proceed in two steps. First, we verify that $\sqrt{n}\,(\widehat{F}^{*}(\cdot)-F(\cdot))$ converges in distribution to a tight Gaussian process $G(\cdot)$; then, we can bound the tail probabilities of $\sup_{y\in\mathbb{R}}\left\lvert G(y)\right\rvert$ directly. The first step is a routine application of Donsker’s theorem as presented, e.g., in Chapter 19 of Van der Vaart (1998); here, we take the existence of a limiting process $G(\cdot)$ as given. Because the supremum of $\sqrt{n}\,(\widehat{F}^{*}(\cdot)-F(\cdot))$ is invariant to monotone transformations of $Y$, and is stochastically maximized when $Y$ has a continuous density without atoms, we can without loss of generality assume that $Y\in[0,\,1]$.

Now, to move forward, it is convenient to assume that the actual sample size we observe is Poisson, $n\sim\text{Poisson}(N)$, with $N\rightarrow\infty$; this assumption will not affect our conclusions, but makes the exposition more transparent. Given this Poisson assumption, we can, again without loss of generality, assume that $y$ is scaled such that there is a $\omega>0$ such that

[TABLE]

where $\mathbb{E}_{obs}$ denotes expectations with respect to the observed (i.e., biased) sample. The above is simply rescaling $y$ so the variance of $\hat{H}(y)$, the total mass sampled below $y$, linearly increases with $Y$, allowing conception of this object as, effectively, a random walk. We note that $\hat{H}(y)$ is the un-normalized weighted estimator for $F$, and, using standard results on compound Poisson processes, we have

[TABLE]

In connecting the above two displays, note the additional variability due to the random sample size $n$ plays a critical role. Rearranging the above gives $\omega^{2}=\mathbb{E}_{obs}\left[\pi_{i}^{-2}\right]\,/\,\mathbb{E}_{obs}\left[\pi_{i}^{-1}\right]^{2}$, where $\mathbb{E}_{obs}$ denotes the observed sampling distribution.

Given these preliminaries, an examination of pairwise covariances implies that, if Gaussian limits exist—and we know that they do by Donsker’s theorem—we must have $\sqrt{N}\,(\hat{H}(y)-F(y))\Rightarrow\omega W(y),$ where $W(y)$ is a standard Wiener process on $[0,\,1]$. Further, noting that $\widehat{F}^{*}(y)=\hat{H}(y)/\hat{H}(1)$, we also have

[TABLE]

noting that, here can use $\sqrt{n}$ instead of $\sqrt{N}$ thanks to Slutsky’s theorem, since $\sqrt{n}-\sqrt{N}=\mathcal{O}_{P}(1)$.

If we had $F(y)=y$, then $G(y)$ would be a standard Brownian bridge. However, in our case, $F(y)$ can take on a wider range of values. To bound the maximum of $G(y)$ we will make use of the following technical Lemma, proved at the end of the end of this Section.

Lemma 3.1**.**

Suppose that $W(t)$ is a standard Wiener process, and that $F(t)$ is a monotone-increasing, differentiable function on $[0,\,1]$ with $F(0)=0$, $F(1)=1$ and $\sup F^{\prime}(t)/\inf F^{\prime}(t)\leq\gamma$ for some constant $\gamma>0$. Then, there exists a constant $\sigma_{\gamma}$, bounded by 1, for which the stochastic process $X(t)=W(t)-F(t)W(1)$ satisfies the following bound:

[TABLE]

Inverting the above inequality, we find that

[TABLE]

Finally, in order to bound $\omega$, we define $\mathbb{E}_{pop}$ as the expectation with respect to the underlying (unbiased) population, and can then check that

[TABLE]

This implies that

[TABLE]

where the last inequality follows form the fact that, because $1/x$ is a convex function of $x$, our expression of interest is maximized when $\pi_{i}$ only takes on two discrete values representing the endpoints of the allowed range. We can check by calculus that the above bound is maximized at $a^{*}=1/2$, resulting in the bound $\omega^{2}\leq(1+\gamma)(1+\gamma^{-1})/4$.

Proof of Corollary 2.3

By Lemma 2.2, we know that $\widehat{F}^{*}\in\mathcal{C}_{\gamma,\,n}$ with probability tending to 1, so the result follows immediately from Theorem 2.1.

Proof of Lemma 2.4

As in the proof of Lemma 2.2, we start by re-scaling our problem such that $Y\in[0,\,1]$, and that $\hat{H}(y)$ has linearly increasing variance. Given this setup, we again have that, over $y\in[0,\,1]$,

[TABLE]

where $W(y)$ is a standard Wiener process and $0<\omega^{2}\leq(1+\gamma)(1+\gamma^{-1})/4$. It follows that, over $q\in[0,\,0.5]$,

[TABLE]

where for the last statement we note that $W(F^{-1}(q))$ and $W(1)-W(F^{-1}(1-q))$ are two independent Gaussian processes over $q\in[0,\,0.5]$, because the Wiener process has independent increments. Finally, we note that, for $q\in[0,\,0.5]$, $F^{-1}(q)+1-F^{-1}(1-q)$ takes all values in $[0,\,1]$, and so we find that for any threshold $t$,

[TABLE]

where the last equality is a consequence of the well known reflection principle for Brownian motion (e.g., Mörters and Peres 2010). The desired conclusion then follows by noting our bound on $\omega$ and applying the bound to both tails.

Proof of Lemma 2.5

To prove this result, it suffices to verify that all claims made in the 4 steps leading to our construction of $\mathcal{C}_{\gamma,\,n}^{LC+}$ in fact holds. Step 1 is a direct consequence of the Kolmogorov-Smirnov theorem. Steps 3 and 4 are also immediate given the result of Prékopa (1973). The interesting fact here is that the optimization problem used to define $\widehat{L}_{m}(y)$ is computationally tractable: Finding a concave upper-bound for a function $l(y)$ is equivalent to taking the convex hull of the curve $(y,\,l(y))$.

It remains to check step 2. Define $r=\min\left\{\pi_{i}\right\}/\mathbb{E}\left[\pi_{i}\right]$; given a known value of $r$ and the fact that $\max{\pi_{i}}/\min{\pi_{i}}\leq\gamma$ we can then verify that $F$ is stochastically dominated by the following function, whose sampling weights $\pi_{i}$ jump from a low value $a$ to a high value $\gamma a$ at a threshold $m$ characterized by

[TABLE]

where $S(\cdot)$ is the observed sampling distribution. In other words, this implies that

[TABLE]

The claim made in Step 2 then follows by scanning over all $m$, and noting that $S(y)\geq\widehat{S}(y)$ with probability tending to 1 thanks to Step 1.

Proof of Lemma A3.1

We first recall the classic result of Marcus and Shepp (1972) which, in our situation, implies that the Gaussian process $X(t)$ satisfies

[TABLE]

it thus remains to use the shape restrictions on $F(t)$ to lower-bound the right-hand expression. To do so, it is helpful to decompose the stochastic process $X(t)$ into two independent parts:

[TABLE]

where $B(t)=W(t)-t$ is a standard Brownian bridge (here, we use the well known fact that $B(t)$ and $W(1)$ are independent of each other). Given our assumptions on $F(t)$, we can moreover show that

[TABLE]

for $0\leq t\leq 0.5$, and $A_{\gamma}(t)=A_{\gamma}(1-t)$ for all $t\in[0\,1]$. Thus, noting that $\operatorname{Var}\left[B(t)\right]=t(1-t)$ and $\operatorname{Var}\left[W(1)\right]=1$, we see that

[TABLE]

Finally, the above function is concave on $[0,\,0.5]$ (in particular, $A^{\prime\prime}_{\gamma}(t)=-2\gamma(\gamma-1)\,/\,(1-t+\gamma t)^{3}$), and so the above expression has a unique maximizer $t^{*}$ that can be derived numerically. The desired conclusion then holds for $\sigma^{2}_{\gamma}:=t_{\gamma}^{*}(1-t_{\gamma}^{*})+A_{\gamma}(t^{*}_{\gamma})$; the fact that $\sigma_{\gamma}^{2}\leq 1$ is immediate by inspection.