Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
77,665 characters · 22 sections · 59 citation commands
Compound Estimation for Binomials
\onehalfspacing
\onehalfspacing
Estimating mean of multiple binomial outcomes is a common task in applied economics, public policy, and experimental analysis. In settings ranging from labor market discrimination and education interventions to development programs and health policy, researchers often observe binomial outcomes—such as employment rate, school attendance, or treatment uptake—and aim to estimate the mean parameters (one-sample setting) or treatment effects (two-sample setting). For instance, bell2019becomes link the tax records with the patent records and report the innovation rates across hundreds of American colleges. These rates can be viewed as the estimated mean parameters from multiple binomials, each corresponding to the number of inventors in a college. As an example of the two-sample setting, kline2024discrimination analyze a large-scale correspondence experiment, sending up to 1,000 fictitious job applications with randomly assigned race and gender indicators to 108 firms. Differences in callback rates across race–gender groups yield firm-level binomial outcomes indicative of discriminatory practices, analyzed via Empirical Bayes estimation to produce adjusted discrimination measures.
Existing methods for estimating multiple binomial means can be broadly classified into three categories. The most common approach is to treat multiple binomials independently. This includes the simple estimator via the empirical average and prediction-based estimator when covariates are available. However, these methods do not allow for information sharing across different binomial units, limiting their efficiency when estimating a large number of binomial outcomes. The second category treats multiple binomial outcomes within a Bayesian framework, using either fully Bayesian or Empirical Bayes (EB) methods. When sample sizes $n_i$ are equal, both Bayes and EB can be applied with the beta priors on the mean parameters griffin1971empirical. However, when $n_i$ varies across observations fienberg1973simultaneous, inference requires specifying a joint prior over $(n, p)$ (where $p$ is the binomial parameter), which complicates computation due to discreteness and the lack of a conjugate prior for $n$, making the effectiveness of Bayes or EB methods less clear in this setting. The third category treat the binomials as approximate Gaussian variables with variances imputed as the large-sample approximation Brown2008inseason,xie2012sure,chen2022empirical,chen2025compound, and develop theoretical guarantees assuming the measurements are exactly Gaussian. Nevertheless, the Gaussian approximation tends to be inaccurate when the sample sizes or the mean parameters are small. In addition, under the Bayesian or EB framework where the size and mean parameters are sampled from a distribution, the prior is only partially identified due to discreteness of binomial measurements kline2021reasonable. By contrast, the prior is point identified with Gaussian measurements under fairly mild regularity conditions efron2012large,chen2022empirical. This fundamental difference poses another threat to the Gaussian-approximation-based EB methods.
Our paper directly addresses the limitation of the Gaussian approximation by directly accommodating binomial outcomes through their exact sampling distribution. Unlike most previous work that develop EB methods, we formulate the problem as a compound estimation problem robbins1951asymptotic, where we treat all size and mean parameters as fixed, rather than random variables generated from a prior distribution. In particular, we allow arbitrarily heterogeneous and bounded sample sizes across observations. We show that, in terms of average mean squared error, our proposed estimator outperforms the maximum likelihood estimator and any single machine-learning estimator under mild regularity conditions. When covariate information is available, the estimator naturally incorporates arbitrary machine learning predictors. In addition, the proposed estimator satisfies a reporting consistency property: its weighted average coincides exactly with the naïve weighted average, a property not shared by existing methods. By contrast, many alternative shrinkage estimators fail to satisfy this property, leading to accounting inconsistency that could be misleading in high-stake settings.
Our methodology proceeds as follows. First, we derive a Stein’s unbiased risk estimator (SURE) for squared-error loss specifically tailored to binomial parameters. Using this, we propose a family of shrinkage estimators that combine the maximum likelihood estimator (MLE), the grand mean, and predictions from machine learning models. Our approach is closely related to the shrinkage estimator introduced by xie2012sure, which shrinks estimates toward either the grand mean or a data-driven location under a Gaussian model assumption with known variances. In contrast, our shrinkage estimator minimizes the SURE associated with average mean squared error under a binomial model assumption. We derive the regret bound, establish the asymptotic normality of our proposed estimator, develop a corresponding statistical inference procedure, and provide practical guidance for constructing confidence regions.
\paragraph{Stein's Identity for Binomial Distributions} Our estimator is inspired by Stein's identity for binomial distributions. Stein's identity was originally introduced by stein1972bound for approximating the distribution of sums of dependent random variables by a normal distribution. Later, barbour2001compound extended Stein's identity to binomial and related discrete distributions. Other related developments include work by ehm1991binomial and soon1996binomial, etc.
\paragraph{Stein's Unbiased Risk Estimator} Additionally, our estimator is based on the Stein's unbiased risk estimator (SURE), which was first proposed by stein1981estimation, as an unbiased estimator of the mean-squared error for estimating the mean of a multivariate normal distribution. Since then, there is a considerable amount of literature that studies the minimization of a SURE-type risk estimate via relatively simple estimators (e.g. linear smoothers) li1985stein,li1986asymptotic,li1987asymptotic,johnstone1988inadmissibility,kneip1994ordered,donoho1995adapting. More recently, there has been a line of work applying SURE for tuning regularization parameters for high-dimensional methods such as the Lasso, reduced-rank regression, and singular value thresholding tibshirani2012degrees,candes2013unbiased,mukherjee2015degrees. Further, xie2012sure proposed a class of SURE-based shrinkage estimators and showed a uniform consistency property for SURE in a hierarchical model.
Notably, both the earlier work and recent studies xie2012sure,ghosh2025stein,nobel2023tractable,karamikabir2021wavelet,kim2021noise2score, chen2025compound rely on the Gaussian distribution assumption. While eldar2008generalized derived the SURE for mean squared error within general exponential families, their primary goal was to select regularization parameters for rank-deficient Gaussian models and linear Gaussian models.
\paragraph{Bayes and Empirical Bayes Methods} A substantial body of literature also addresses the estimation of binomial outcomes through Bayes and empirical Bayes methods. griffin1971empirical derived a Bayes estimator expressed in terms of the marginal probabilities rather than the prior, thereby obtaining an empirical Bayes estimator. berry1979empirical applied the theory of Dirichlet processes to the empirical Bayes estimation for the binomial outcomes. albert1984empirical proposed a empirical Bayes method by defining a class of prior distributions for a set of binomial probabilities to reflect the user's prior belief about the similarity of the probabilities. sobel1993bayes constructed ranking procedures for comparing multiple binomial parameters via both Bayes and empirical Bayes approaches. sivaganesan1993robust proposed an empirical Bayes approach that partially identifies the posterior means of binomial outcomes by imposing moment conditions on the unknown prior. consonni1995bayesian adapted a Bayesian approach to estimate the binomial parameters by imposing prior information about the partitions of the binomial experiments. weiss2010bayesian used a Bayesian hierarchical approach to simultaneously estimate the parameters of multiple binomial distributions. kline2021reasonable employed an empirical Bayes approach to estimate binomial parameters while treating sample sizes as fixed, and develop partial-identification methods for moments in the two-sample setting, with an application to job-level discrimination detection. gu2025reasonable further advanced this line of work by incorporating both partial-identification and sampling uncertainty to construct valid confidence intervals for empirical Bayes estimators.
For any positive integer $n$, we use $[n]$ to denote $\{1,2,\ldots,n\}$. Given a random variale $A$ and a distribution $\mathcal{F}$, we write $A\sim\mathcal{F}$ to imply that $A$ follows distribution $\mathcal{F}$. For any two random variables $A$ and $B$, $A\ \raisebox{0.05em}{\rotatebox[origin=c]{90}{$\models$}} \ B$ means $A$ is independent of $B$. We use $\mathbf{1}(\cdot)$ to denote the indicator function. For any matrix or vector $\boldsymbol{\nu}$, we use $\boldsymbol{\nu}^T$ to denote the transpose of $\boldsymbol{\nu}$. Given any matrix $\mathbf{M}$, we use $\|\mathbf{M}\|_{\infty,\infty}$ to denote the maximum absolute value for the entries of $\mathbf{M}$. For any vector $\boldsymbol{\nu}=(\nu_1,\ldots,\nu_k)\in\mathbb{R}^d$, $\|\boldsymbol{\nu}\|_{\infty}=\max_{k\in[d]}|\nu_k|$ and $\|\boldsymbol{\nu}\|_{2}=\sqrt{\sum_{k=1}^d\nu_k^2}$. Given any two matrix $\mathbf{M}_1$ and $\mathbf{M}_2$, we write $\mathbf{M}_1\preceq\mathbf{M}_2$ to mean that $\mathbf{M}_2-\mathbf{M}_1$ is positive semi-definite. For any sequence of random variables $X_n$, we say that $X_n=\mathrm{o}_p(1)$ if $X_n$ converges in probability to $0$, and we say that $X_n=\mathrm{O}_p(R_n)$ for some real sequence $R_n$, if $X_n=Y_nR_n$, and $\{Y_n\}$ is uniformly tight. We use $\rightsquigarrow$ to imply “converges in distribution to” and $\overset{p}{\rightarrow}$ to imply “converges in probability to”. For any function $f(\cdot)$, we use $\nabla f(\cdot)$ to denote the gradient of $f$. For any $d\in\mathbb{Z}_{+}$, we use $\mathbf{I}_d$ to denote the $d$-by-$d$ identity matrix.
We aim to estimate the unknown binomial parameters $\{\theta_i^{\mathrm{o}}\}_{i\in[N]}$ in the one-sample setting and $\{\theta_{i1}^{\mathrm{t}},\theta_{i2}^{\mathrm{t}}\}_{i\in[N]}$ in the two-sample setting.
In the one-sample setting, we observe $\{n_i,\mathbf{X}_i,Y_i\}_{i\in[N]}$, where $n_i\in\mathbb{Z}_{+}$, $\mathbf{X}_i\in\mathbb{R}^d$ are fixed for any $i\in[N]$, $\{Y_i\}_{i\in[N]}$ are independently but not necessarily identically distributed (i.n.i.d.) random variables such that $Y_{i}\sim \mathrm{Bin}(n_i, \theta_i^{\textrm{o}})$, and
Let $\hat{\theta}_i^{\textrm{o}}$ denote the estimator for $\theta_i^{\textrm{o}}$, and $\hat{\theta}_i^{\textrm{o}}$ could depend on all observations.
In the two-sample setting, we observe $\{n_{i1},\mathbf{X}_{i1},Y_{i1}\}_{i\in[N]}$ and $\{n_{i2},\mathbf{X}_{i2},Y_{i2}\}_{i\in[N]}$ from group one and two, where $n_{i\ell}\in\mathbb{Z}_{+}$ and $\mathbf{X}_{i\ell}\in\mathbb{R}^d$ are fixed, $\forall i\in[N], \ell\in\{1,2\}$. The group one binomial outcomes $\{Y_{i1}\}_{i\in[N]}$ are independent from $\{Y_{i2}\}_{i\in[N]}$ of group two. Here $Y_{i1}\sim \mathrm{Bin}(n_{i1}, \theta_{i1}^{\textrm{t}})$ and $Y_{i2}\sim \mathrm{Bin}(n_{i2},\theta_{i2}^{\textrm{t}})$, where
Suppose $\{Y_{i1}\}_{i\in[N]}$ and $\{Y_{i2}\}_{i\in[N]}$ are both i.n.i.d. across $i\in[N]$. We estimate $\theta_{i1}^{\textrm{t}}-\theta_{i2}^{\textrm{t}}$ using the two-sample estimator $\hat{\theta}_i^{\textrm{t}}$ which could also depend on all observations.
Denote the vectors of the unknown estimands as $\theta^{\mathrm{o}} := (\theta_1^{\mathrm{o}}, \ldots, \theta_N^{\mathrm{o}})$, $\theta_{\ell}^{\mathrm{t}}:=(\theta_{1\ell}^{\mathrm{t}},\ldots,\theta_{N\ell}^{\mathrm{t}})$, $\ell\in\{1,2\}$. Denote the vectors of the one-sample and two-sample estimators as $\hat{\theta}^{\mathrm{o}} := (\hat{\theta}_1^{\mathrm{o}}, \ldots, \hat{\theta}_N^{\mathrm{o}})$ and $\hat{\theta}^{\mathrm{t}}:=(\hat{\theta}_1^{\mathrm{t}},\ldots,\hat{\theta}_N^{\mathrm{t}})$. The objective is to assess the performance of the estimators in both one-sample and two-sample settings through their $L_2$ risks, specified in (ref) and (ref), respectively:
Since $\sum_{i=1}^N \left(\theta_i^{\mathrm{o}}\right)^2$ and $\sum_{i=1}^N \left(\theta_{i1}^{\mathrm{t}} - \theta_{i2}^{\mathrm{t}}\right)^2$ are constants, minimizing the one-sample $L_2$ risk is equivalent to minimizing $L_{\mathrm{o}}$ as follows:
and minimizing the two-sample $L_2$ risk is equivalent to minimizing $L_{\mathrm{t}}$ defined as follows:
The first terms in (ref) and (ref) can be unbiasedly estimated by the plug-in estimator $(1/N) \sum_{i=1}^N \left(\hat{\theta}_i^{\mathrm{o}}\right)^2$ and $(1/N) \sum_{i=1}^N \left(\hat{\theta}_i^{\mathrm{t}}\right)^2$, respectively. The main challenge stems from the second terms because of the unknown parameters $\theta_i^\mathrm{o}$ and $(\theta_{i1}^\mathrm{t} - \theta_{i2}^\mathrm{t})$. In the following subsections, we derive a Stein identity for the binomial distribution and use it to construct Stein’s unbiased risk estimators (SUREs) for the second terms in the one-sample and two-sample objectives, respectively.
We start with a simpler problem of estimating a single Binomial parameter $\theta$ for $Y \sim \mathrm{Bin}(n,\theta)$, without incorporating covariates, by applying Stein’s identity for Binomial distributions.
In what follows, we show by Theorem (ref) that a SURE does exist for the binomial case when the estimator is constructed via polynomial functions. For any function $h$ on $\{0, \ldots, n\}$, define
and
For any $a\in \{0, 1, \ldots, n\}$, define
Theorem (ref) implies that if $h(\cdot)$ is a polynomial of degree less than $n$, then $\mathcal{T}h(Y;n,a)$ is an unbiased estimator of $\theta\mathbb{E}[h(Y)]$. Throughout the remainder of the paper, we set $a=\lfloor n/2\rfloor$, which is a robust choice according to statement (i) of Theorem (ref).
We propose a class of estimators for the one-sample estimands (ref), where for any $i\in[N]$,
and $\hat{g}(\mathbf{X}_i)$ is the machine learning (ML) estimator for $g(\mathbf{X}_i)=\mathbb{E}[\theta_i^{\mathrm{o}}|\mathbf{X}_i]$ defined as (ref). Similarly, for two-sample case, for any $i\in[N]$, $\ell\in\{1,2\}$ and $\boldsymbol{\lambda}=(\lambda_1,\lambda_2)\in[0,1]\times\mathbb{R}$, define
where for any $\ell\in\{1,2\}$, $\hat{g}_{\ell}(\mathbf{X}_{i1})$ is the ML estimator for $g_{\ell}(\mathbf{X}_{i1})=\mathbb{E}[\theta_{i\ell}^{\mathrm{t}}|\mathbf{X}_{i1}]$ defined as (ref). We propose a class of estimators for the two-sample estimands (ref) as $\hat{\theta}_i^{\mathrm{t}}(\boldsymbol\lambda) := \hat{\theta}_{i1}^{\mathrm{o}}(\boldsymbol{\lambda}) - \hat{\theta}_{i2}^{\mathrm{o}}(\boldsymbol{\lambda})$.
It is straightforward to see that both the one-sample and two-sample estimators satisfy a reporting consistency property: the aggregate of the unit-level estimates exactly equals the overall empirical proportion.
This feature of reporting consistency is typically not satisfied by standard estimators (e.g., machine learning models, empirical Bayes), whose fitted unit-level values need not aggregate back to the raw overall proportion. Such inconsistency could be problematic in litigation and regulatory settings, where differences between the reported overall rate and the model-based unit-level estimates may call the analysis into question. In applications such as discrimination or pay‐equity studies, where a single firm-level figure is reported to regulators or courts, our estimator avoids this problem: the job-level estimates always average exactly to the overall firm-level rate used for reporting.
We use the $K$-fold cross-fitting method to estimate $\hat{g}(\cdot)$ for the one-sample estimator or $\hat{g}_{\ell}(\cdot)$, $\forall\ell\in\{1,2\}$ for the two-sample estimator. In particular, we split the sample index $[N]$ into $K$ disjoint folds $\mathcal{I}_1,\ldots,\mathcal{I}_K$. For any $i\in[N]$, let $k(i)$ denote the fold that the $i$-th sample belongs to. Notice that it is not required that the folds have exactly equal size, for example each $k(i)$ could be drawn uniformly from $[K]$. We proceed with having equal-size folds for simplicity, without loss of generality. We let $\hat{g}^{-k}(\cdot)$ (resp. $\hat{g}_{\ell}^{-k}(\cdot)$, $\forall\ell\in\{1,2\}$) denote $\hat{g}(\cdot)$ (resp. $\hat{g}_{\ell}(\cdot), \forall\ell\in\{1,2\}$) computed without using observations from fold $k$. Further, we clip $\hat{g}(\cdot)$, $\hat{g}_{\ell}(\cdot), \forall\ell\in\{1,2\}$ between $0$ and $1$ as the final ML outputs. Specifically, given any $k\in[K]$, for any $i\in\mathcal{I}_k$, we set $\hat{g}(\mathbf{X}_i)=\hat{g}^{-k}(\mathbf{X}_i)$ in (ref), and we set $\hat{g}_{\ell}(\mathbf{X}_{i\ell})=\hat{g}_{\ell}^{-k}(\mathbf{X}_{i\ell})$ in (ref).
Then the one-sample binomial estimator parametrized by $\boldsymbol{\lambda}=(\lambda_1,\lambda_2)\in[0,1]\times\mathbb{R}$ can be rewritten as
Similarly, both components (ref) of the two-sample binomial estimator parametrized by $\boldsymbol{\lambda}=(\lambda_1,\lambda_2)\in[0,1]\times\mathbb{R}$ for each $\ell\in\{1,2\}$ can be rewritten as
and
We refer to the estimators (ref) and (ref) as the one-sample binomial-shrinkage estimator and two-sample binomial-shrinkage estimator respectively. Both estimators interpolate the maximum likelihood estimator (MLE), the grand mean and the machine learning (ML) model estimator. They can also be viewed as shrinking MLE towards the grand mean and the ML estimates. In the absence of covariates, we simply set $\lambda_2 = 0$, so the estimators for this case reduce to special instances of (ref) and (ref).
For any $\boldsymbol{\lambda}\in[0,1]\times\mathbb{R}$, define
Thus, (ref) is equivalent to the one-sample objective function (ref) with respect to the binomial shrinkage estimator $\hat{\theta}_i^{\mathrm{o}}(\boldsymbol{\lambda})$, and (ref) is equivalent to the two-sample objective function (ref) with respect to the binomial shrinkage estimator $\hat{\theta}_i^{\mathrm{t}}(\boldsymbol{\lambda})$.
Recall that the functional $\mathcal{T}$ is defined as (ref), where we omit $a$ in the notation since throughout we set $a=\lfloor n/2\rfloor$ in (ref). In order to align the form with the definition of the operator $\mathcal{T}$ in (ref), we write the one-sample binomial shrinkage estimator in (ref) as $\hat{\theta}_i^{\mathrm{o}}(Y_i;n_i| \boldsymbol{\lambda})$ to emphasize its dependence on $Y_i$ and $n_i$. Specifically, we define $\mathcal{T}\hat{\theta}_i^{\mathrm{o}}(\boldsymbol{\lambda}):=\mathcal{T}\hat{\theta}_i^{\mathrm{o}}(Y_i;n_i|\boldsymbol{\lambda})$, where $\mathcal{T}\hat{\theta}_i^{\mathrm{o}}(Y_i;n_i|\boldsymbol{\lambda})$ denotes the term where we apply functional $\mathcal{T}$ to $\hat{\theta}_i^{\mathrm{o}}(\boldsymbol{\lambda})$ by fixing the grand mean $(\sum_{i=1}^NY_i)/(\sum_{i=1}^Nn_i)$ and the ML model outputs $\hat{g}^{-k(j)}(\mathbf{X}_j)$, $\forall j\in[N]$.
For the two-sample binomial-shrinkage estimator (ref), we write it as $\hat{\theta}_i^{\mathrm{t}}(Y_{i1};n_{i1}|\boldsymbol{\lambda})$ (resp. $\hat{\theta}_i^{\mathrm{t}}(Y_{i2};n_{i2}|\boldsymbol{\lambda})$) to emphasize it dependence on the parameter $\boldsymbol{\lambda}$, $Y_{i1}$ and $n_{i1}$ (resp. $Y_{i2}$ and $n_{i2}$). Specifically, for any $\ell\in\{1,2\}$, $\mathcal{T}\hat{\theta}_i^{\mathrm{t}}(Y_{i\ell};n_{i\ell}|\boldsymbol{\lambda})$ denotes the term where we apply the functional $\mathcal{T}$ to $\hat{\theta}_i^{\mathrm{t}}(\boldsymbol{\lambda})$ by fixing $\{Y_{im}$, $n_{im}\}$ for $m=3-\ell$, the grand means $(\sum_{i=1}^NY_{i\ell})/(\sum_{i=1}^Nn_{i\ell}),\ \forall\ell\in\{1,2\}$ and the ML model outputs $\hat{g}^{-k(j)}(\mathbf{X}_{j\ell})$ for any $\ell\in\{1,2\}, j\in[N]$. Define
The explicit expressions obtained by expanding (ref) and (ref) are given in (ref) and (ref) in the Appendix (ref).
To derive the bias bounds in the one- and two-sample settings, we impose the following assumptions.
Proposition (ref) below implies that $\hat{L}_{\mathrm{o}}(\boldsymbol{\lambda})$ is an approximate SURE for $L_{\mathrm{o}}(\boldsymbol{\lambda})$ defined in (ref), and that $\hat{L}_{\mathrm{t}}(\boldsymbol{\lambda})$ is an approximate SURE for $L_{\mathrm{t}}(\boldsymbol{\lambda})$ defined in (ref).
Proposition (ref) shows that, as $N\rightarrow \infty$, $\mathbb{E}[\hat{L}_{\mathrm{o}}(\boldsymbol{\lambda})]$ and $\mathbb{E}[\hat{L}_{\mathrm{t}}(\boldsymbol{\lambda})]$ closely approximate the true one-sample and two-sample objectives $L_{\mathrm{o}}(\boldsymbol{\lambda})$ and $L_{\mathrm{t}}(\boldsymbol{\lambda})$ up to the constant terms, respectively. Hence, $\hat{L}_{\mathrm{o}}(\boldsymbol{\lambda})$ and $\hat{L}_{\mathrm{t}}(\boldsymbol{\lambda})$ are approximate SUREs for their corresponding objectives. A natural approach, therefore, is to use the minimizers of $\hat{L}_{\mathrm{o}}(\boldsymbol{\lambda})$ and $\hat{L}_{\mathrm{t}}(\boldsymbol{\lambda})$ as approximations to the minimizers of $L_{\mathrm{o}}(\boldsymbol{\lambda})$ and $L_{\mathrm{t}}(\boldsymbol{\lambda})$, respectively.
We note that using the class of estimators as proposed in (ref) and (ref) the one-sample and two-sample objectives (ref), (ref), together with their approximate SUREs, are quadratic functions of $\boldsymbol{\lambda}$, as summarized in the following proposition:
By first-order condition, Proposition (ref) implies that the unconstrained one-sample parameter estimator $\hat{\boldsymbol{\lambda}}_{\mathrm{o}}$\footnote{In the following we use subscript $\mathrm{o}$ (resp. $\mathrm{t}$) with $\hat{\boldsymbol{\lambda}}_\mathrm{o}$ (resp. $\hat{\boldsymbol{\lambda}}_\mathrm{t}$) or $\boldsymbol{\lambda}_\mathrm{o}^*$ (resp. $\boldsymbol{\lambda}_\mathrm{t}^*$) to imply the one-sample (resp. two-sample) estimated $\hat{\boldsymbol{\lambda}}$ or the $\boldsymbol{\lambda}^*$ as the true minimizer related to the one-sample (resp. two-sample) $L_2$ loss function.} that minimizes $\hat{L}_{\mathrm{o}}(\boldsymbol{\lambda})$ over $\boldsymbol{\lambda}\in\mathbb{R}\times\mathbb{R}$ is
The unconstrained two-sample parameter estimator $\hat{\boldsymbol{\lambda}}_{\mathrm{t}}$ that minimizes $\hat{L}_{\mathrm{t}}(\boldsymbol{\lambda})$ over $\boldsymbol{\lambda}\in\mathbb{R}\times\mathbb{R}$ is
Let $\boldsymbol{\lambda}_{\mathrm{o}}^*$ and $\boldsymbol{\lambda}_{\mathrm{t}}^*$ denote the minimizers of $L_{\mathrm{o}}(\boldsymbol{\lambda})$ (ref) and $L_{\mathrm{t}}(\boldsymbol{\lambda})$ (ref) over $[0,1]\times\mathbb{R}$, respectively. We refer to $\hat{\boldsymbol{\lambda}}_\mathrm{o}$ and $\hat{\boldsymbol{\lambda}}_\mathrm{t}$ as the one-sample approximate SURE and the two-sample approximate SURE, respectively.
When $\boldsymbol{\lambda}_{\mathrm{o}}^*$ and $\boldsymbol{\lambda}_{\mathrm{t}}^*$ are unconstrained, i.e. $\lambda_{\mathrm{o}1}^*\in(0,1)$ and $\lambda_{\mathrm{t}1}^*\in(0,1)$, the asymptotic distributions of $\hat{\boldsymbol{\lambda}}_{\mathrm{o}}$ and $\hat{\boldsymbol{\lambda}}_{\mathrm{t}}$ are as follows:
Assumption (ref) and Assumption (ref) include Lindeberg–Feller type conditions, which are used to establish a multivariate central limit theorem for i.n.i.d. samples. Recall from (ref) that the feasible parameter $\boldsymbol{\lambda}=(\lambda_1,\lambda_2)\in[0,1]\times\mathbb{R}$. The asymptotic normality result in Theorem (ref) holds only when $\boldsymbol\lambda_{\mathrm{o}}^*$ and $\boldsymbol\lambda_{\mathrm{t}}^*$ lie in the interior of the feasible region. So the validity of statistical inference based on Theorem (ref) fails when $\boldsymbol\lambda_{\mathrm{o}}^*$ or $\boldsymbol\lambda_{\mathrm{t}}^*$ is on the boundary (i.e., $\boldsymbol{\lambda}_{\mathrm{o}1}^*\in\{0,1\}$ or $\boldsymbol{\lambda}_{\mathrm{t}1}^*\in\{0,1\}$). In such cases, we adopt the inference method for constrained extremum estimators proposed by li2024inference, with details provided in (ref).
The regret bounds for both one-sample and two-sample objectives follow immediately from Theorem (ref) as follows:
In some applications, it may be costly to replace a simple status-quo estimator, such as the empirical average, with the more sophisticated estimator given by our SURE method. As a result, the policy maker would adopt the new estimator only if it has a sufficient efficiency gain over the incumbent. Suppose the class of estimators nests the current estimator with $\boldsymbol{\lambda} = \boldsymbol{\lambda}_0$, we can formulate the problem of whether to use $\hat{\boldsymbol{\lambda}}$ in place of $\boldsymbol{\lambda}_0$ through a hypothesis test with the null hypothesis $H_0: \boldsymbol{\lambda}^* = \boldsymbol{\lambda}_0$.
In this section, we consider the more general problem of on constructing confidence regions for $\boldsymbol{\lambda}^*$ based on the approximate SURE estimators $\hat{\boldsymbol{\lambda}}_{\mathrm{o}}$ and $\hat{\boldsymbol{\lambda}}_{\mathrm{t}}$. Specifically, we consider two scenarios. First, if the true value of $\boldsymbol{\lambda}_{\mathrm{o}}^*$ or $\boldsymbol{\lambda}_{\mathrm{t}}^*$ is believed to lie in the interior of the parameter space $[0,1]\times\mathbb{R}$, we construct the confidence region using standard inference techniques for unconstrained estimators, based upon the asymptotic normality results of Theorem (ref) and Theorem (ref). Alternatively, if the true value of $\boldsymbol{\lambda}_{\mathrm{o}}^*$ or $\boldsymbol{\lambda}_{\mathrm{t}}^*$ is suspected to lie on the boundary of $[0,1]\times\mathbb{R}$, we apply inference methods tailored to constrained estimators following li2024inference.
To perform statistical inference on the unconstrained $\boldsymbol{\lambda}^*$ (where $\boldsymbol{\lambda}^*=\boldsymbol{\lambda}_{\mathrm{o}}^*$ for one-sample case and $\boldsymbol{\lambda}^*=\boldsymbol{\lambda}_{\mathrm{t}}^*$ for two-sample case), we utilize the asymptotic normality results derived for unconstrained $\boldsymbol{\lambda}$ in Theorem (ref) and Theorem (ref), and estimate the variance of $\hat{\boldsymbol{\lambda}}$ using bootstrap methods (where $\hat{\boldsymbol{\lambda}}=\hat{\boldsymbol{\lambda}}_{\mathrm{o}}$ for one-sample case and $\hat{\boldsymbol{\lambda}}=\hat{\boldsymbol{\lambda}}_{\mathrm{t}}$ for two-sample case). Specifically, we generate $B$ bootstrap samples, and compute $\{\hat{\boldsymbol{\lambda}}^{(b)}\}_{b\in[B]}$ for each bootstrap sample, and then compute the covariance matrix $\hat{\boldsymbol{V}}$ of the resulting bootstrap estimates $\hat{\boldsymbol{\lambda}}^{(b)}$. Then we construct the level $1-\alpha$ confidence set for $\boldsymbol{\lambda}^*$ as
where $\chi_{2,1-\alpha}^2$ is the chi-square critical value with $2$ degrees of freedom.
For statistical inference on the constrained $\boldsymbol{\lambda}^*$, we follow the procedure outlined in Section 2 of li2024inference to construct the confidence set. Specifically, we detail the steps for computing the confidence set for the constrained estimator in the two-sample case, as this method is only applied to the two-sample discrimination report application in the empirical application in Section (ref), whose $\hat{\boldsymbol{\lambda}}=\hat{\boldsymbol{\lambda}}_{\mathrm{t}}$ (ref) falls on the boundary of the feasible region (thus it is believed that the true $\boldsymbol{\lambda}_{\mathrm{t}1}^*=0$ for this application). The procedure for the one-sample case is very similar, and thus is omitted here for brevity.
Suppose we have already computed $\mathbf{D}_{N,1}, \mathbf{D}_{N,2}$ and $\hat{\boldsymbol{\lambda}}=\hat{\boldsymbol{\lambda}}_{\mathrm{t}}$ according to Lemma (ref) and (ref), the steps to construct the confidence set are as follows:
Under general regularity conditions, li2024inference shows that the confidence set constructed as described above for the constrained estimator is uniformly valid. In our applications, we choose $\gamma=0.5$, $\delta_n=1/\sqrt{n}$.
We close this section by introducing a procedure that splits a single binomial observation into two independent binomial observations that share the same parameter. This is useful in empirical settings where only aggregated binomial outcomes are available. For example, the total number of students in a school and how many pass a test, or how many are classified as innovators. We use the data thinning/fission method neufeld2024data,leiner2025data to construct separate training and holdout samples from such aggregated outcomes, which allows us to compare our binomial-shrinkage estimators with alternative methods, as we do extensively in Section (ref). This is in similar spirit to the Coupled Bootstrap method that validates the EB methods with Gaussian measurements chen2022empirical.
For the one-sample case, suppose each observation indexed by $i\in[N]$ corresponds to $n_i$ binary samples, among which $Y_i$ of them have outcome $1$, so $Y_i\sim\mathrm{Binomial}(n_i,\theta_i^{\mathrm{o}})$. We select $m_i$ samples from the original $n_i$ samples, where $m_i < n_i$. We then generate $Y_i^{(1)}$ according to a hyper-geometric distribution with parameters $(m_i, n_i - m_i, Y_i)$, and define $Y_i^{(2)} = Y_i - Y_i^{(1)}$. We take $\{Y_i^{(2)},n_i - m_i\}_{i\in[N]}$ as the training set and $\{Y_i^{(1)}, m_i\}_{i\in[N]}$ as the holdout set. Consequently, according to neufeld2024data,leiner2025data, we have $$Y_i^{(1)}\sim\mathrm{Binomial}(m_i,\theta_i^{\mathrm{o}}) \mbox{ and } Y_i^{(2)}\sim\mathrm{Binomial}(n_i - m_i,\theta_i^{\mathrm{o}}).$$ We might then compute the binomial estimators on the training set $\mathcal{F}_T^{\mathrm{o}}:=\{Y_i^{(2)},n_i - m_i\}_{i\in[N]}$ and honestly evaluate them with the holdout set $\mathcal{F}_H^{\mathrm{o}}:=\{Y_i^{(1)}, m_i\}_{i\in[N]}$.
Similarly, in the two-sample application, suppose each observation indexed by $i\in[N]$ corresponds to two independent populations: population one with $n_{i1}$ samples (of which $Y_{i1}$ have outcome 1), and population two with $n_{i2}$ samples (of which $Y_{i2}$ have outcome 1). Thus, $Y_{i1}\sim\textrm{Binomial}(n_{i1},\theta_{i1})$, $Y_{i2}\sim\textrm{Binomial}(n_{i2},\theta_{i2})$, and the two-sample binomial parameter $\theta_i^{\mathrm{t}}=\theta_{i1}-\theta_{i2}$. To create the holdout set, we choose $m_{i1}<n_{i1}$ samples from population one and $m_{i2}<n_{i2}$ samples from population two. We then generate $Y_{i1}^{(1)}$ and $Y_{i2}^{(1)}$ independently from hypergeometric distributions with parameters $(m_{i1},n_{i1}-m_{i1},Y_{i1})$ and $(m_{i2},n_{i2}-m_{i2},Y_{i2})$, respectively, and define $Y_{i1}^{(2)}=Y_{i1}-Y_{i1}^{(1)}$ and $Y_{i2}^{(2)}=Y_{i2}-Y_{i2}^{(1)}$. As a result, we have independence: $Y_{i1}^{(1)}\ \raisebox{0.05em}{\rotatebox[origin=c]{90}{$\models$}} \ Y_{i1}^{(2)}$ and $Y_{i2}^{(1)}\ \raisebox{0.05em}{\rotatebox[origin=c]{90}{$\models$}} \ Y_{i2}^{(2)}$. According to neufeld2024data,leiner2025data, we have $$Y_{i1}^{(1)}\sim\mathrm{Binomial}(m_{i1},\theta_{i1}),\ Y_{i1}^{(2)}\sim\mathrm{Binomial}(n_{i1}-m_{i1},\theta_{i1}),$$ $$Y_{i2}^{(1)}\sim\mathrm{Binomial}(m_{i2},\theta_{i2}),\ Y_{i2}^{(2)}\sim\mathrm{Binomial}(n_{i2}-m_{i2},\theta_{i2}).$$ Hence, we might take $\mathcal{F}_T^{\mathrm{t}}:=\{n_{i1}-m_{i1},Y_{i1}^{(2)},n_{i2}-m_{i2},Y_{i2}^{(2)}\}_{i\in[N]}$ as the training set, while $\mathcal{F}_H^{\mathrm{t}}:=\{m_{i1},Y_{i1}^{(1)},m_{i2},Y_{i2}^{(1)}\}_{i\in[N]}$ is taken as the holdout set.
The following proposition shows that the data thinning/fission procedures yield unbiased estimators of the $L_2$ risk function for a generic estimator of the binomial parameters.
Proposition (ref) enables out-of-sample evaluation of binomial estimators even when separate training and holdout datasets are not directly available. We later apply this method to validate the performance of our estimators in the innovation and education empirical applications.
To examine the practical performance of the one-sample and two-sample binomial-shrinkage estimators, we now present three empirical applications related to bell2019becomes, kline2024discrimination and gang2023ranking.
Using data on inventors from patent records linked to tax records, bell2019becomes investigated the determinants of becoming a successful inventor and published an “Opportunity Atlas” detailing patent rates across various population groups, segmented by neighborhood, college attendance, parental income level, and racial background.
Leveraging data from a large-scale resume correspondence experiment, which signaled race and gender to employers through randomly assigned distinctive names, kline2024discrimination measured disparities in contact rates across race and gender categories, yielding noisy estimates of discriminatory behavior at the firm level. They subsequently constructed a "discrimination report card" summarizing experimental evidence on biases exhibited by a broad range of Fortune 500 companies using an Empirical Bayes approach.
gang2023ranking proposed a ranking and selection framework as an alternative to conventional false-discovery-rate analyses. They validated their approach using an empirical study of K-12 school test performance data, identifying significant differences in passing rates between students from socioeconomically advantaged (SEA) and disadvantaged (SED) backgrounds within each school. Specifically, they computed $p$-values based on a normal approximation, despite over $30\%$ of the schools having fewer than $100$ data points in at least one of the two groups.
We demonstrate our methods on three applications, one in the one-sample setting and two in the two-sample settings. The one-sample problem we consider is the Opportunity Atlas innovation dataset\footnote{The dataset description is available at \url{https://opportunityinsights.org/wp-content/uploads/2018/04/Inventors-Codebook-Table-3.pdf}.}, which documents innovation rates across colleges by linking tax records to inventor information from patent records bell2019becomes. The dataset contains the fraction of inventors, defined as individuals who were listed on a patent application between 2001 and 2012 or granted a patent between 1996 and 2014, from $423$ colleges. For college $i\in [423]$, we will treat this share measurement as $Y_i$ and the total number of students as $n_i$. In this application, $Y_i$s are generally close to zero.
The first two-sample application is to detect employment discrimation using a large-scale resume correspondence experiment dataset\footnote{The dataset is available at \url{https://www.aeaweb.org/articles?id=10.1257/aer.20230700}} from kline2022systemic. Each observation corresponds to a job applicant and contains his/her demographic information such as race and gender, the job identifier for which the candidate applied, an indicator of whether the candidate received a callback, and additional characteristics. In this dataset, every job receives up to four resumes from each of the two racial groups. We analyze $9821$ jobs, each with exactly four white applicants and four black applicants. Following kline2021reasonable, kline2022systemic, kline2024discrimination, we measure the extent of discrimination for a job as the callback rate gap white and black job applicants. For any job $i\in [9821]$, $n_{i1} = n_{i2} = 4$, $Y_{i1}$ and $Y_{i2}$ are the fractions of callbacks among the four applications in each racial group.
The second two-sample application utilizes the dataset\footnote{The dataset we use is `AYP_05.csv' from \url{https://github.com/bgang92/rankingselection}} on K-12 school test performance data drawn from the 2005 Annual Yearly Performance (AYP) study, which is analyzed in gang2023ranking. Schools in this dataset are categorized into three types: 'H' for high schools, 'M' for middle schools, and 'E' for elementary schools. For each school $i\in[6398]$, $Y_{i1}$ and $Y_{i2}$ measure number of socially-economically advantaged (SEA) and socially-economically disadvantaged (SED) students who took and passed the test, and $n_{i1}$ and $n_{i2}$ measure the number of SEA and SED students. In this application, $(n_{i1}, n_{i2})$ are vastly heterogeneous.
\paragraph{Reporting Innovation Rates.} We perform a linear regression to construct $\hat{g}(x)$ with 10-fold cross-fitting. In particular, we regress $Y_i$ on two covariates: `patent' (the total number of patents granted to students) and `total_cites' (the total number of patent citations received by students). The regression model is specified as follows: $$\textrm{inventor}=\beta_1\cdot\textrm{patent}+\beta_2\cdot\textrm{total\_cites}+\epsilon.$$ We use the one-sample binomial-shrinkage estimator (ref) to estimate the inventor fraction for each institution. The constrained estimate $\hat{\boldsymbol\lambda}$ coincides with the unconstrained one, which is
The resulting estimated inventor fractions for each college are illustrated in Figure (ref) and the $95\%$ confidence region is plotted in Figure (ref). While $\lambda_1\approx 1$ and $\lambda_2\approx 0$, we can reject that the MLE or global mean is MSE-optimal.
\paragraph{Estimating Employment Discrimination.} We apply a gradient boosted tree learner to classify the callback status (`cb') with 10-fold cross-fitting. The feature set includes job and market characteristics --- such as the state and census region of the employer, the experimental wave, and the industry of the job posting --- as well as application-level information, including the applicant’s implied age at submission, education credentials, and the order in which applications were sent. We additionally incorporate demographic and identity signals embedded in the resume, including gender indicators, age-over-40 status, involvement in LGBTQ, political, or academic organizations, and the use of gender-related pronouns. Finally, the model controls for experiment-design variables, including paired-application identifiers, balance indicators, occupational skill categories, and the number of experimental waves associated with each posting. We define $\hat{g}_{i1}$ and $\hat{g}_{i2}$ as the average predicted callback rates from the trained classifier for white and black applicants for each job $i$, respectively. We then interpolate between the maximum likelihood estimate (MLE), the machine learning predictions, and the grand mean within the estimator framework defined in (ref). The constrained estimate, $\hat{\boldsymbol\lambda}$, is computed accordingly as
The resulting estimates of the callback rate differences for each job are depicted in Figure (ref) and the $95\%$ confidence region is plotted in Figure (ref). In this case, we cannot reject the null that the global mean is MSE-optimal, though we can confidently claim that the MLE is so.
\paragraph{Estimating Test Passing Rate Gaps.} We use the school type as the only covariate and define $\hat{g}_{i1}$ and $\hat{g}_{i2}$ as the average passing rate of SEA and SED students within the same school type as school $i$. We then apply the same procedure as in the second application. The resulting constrained estimate $\hat{\boldsymbol\lambda}$ matches the unconstrained solution:
Figure (ref) plots the estimated differences in test passing rates between SEA and SED students across individual schools. The $95\%$ confidence region is plotted in Figure (ref). Unlike the previous two cases, we find strong evidence that the MSE-optimal estimator within the class must utilize both the shrinkage to global mean and the assistance by the prediction model.
To empirically validate the effectiveness of the binomial-shrinkage estimator, we plot both the one-sample and two-sample SURE surface, defined by (ref) and (ref), respectively, alongside the corresponding holdout mean squared error (MSE). Specifically, we randomly split each dataset into training and holdout subsets. On the training set, we compute the SURE values across a range of $\boldsymbol{\lambda}$ values; on the holdout set, we evaluate the holdout MSE for our binomial-shrinkage estimators at these same $\boldsymbol{\lambda}$ values.
Note that the one-sample and two-sample SUREs are unbiased estimators for (ref) and (ref), respectively, up to constants from (ref) and (ref). Consequently, if our method is valid, the SURE surface and holdout MSE surface should be parallel up to a constant.
To approximate the true risk, we apply the standard data splitting to the employment discrimination application where the individual-level measurements are available. Specifically, for each job position $i$, we randomly select one white applicant and one black applicant as the holdout observation and leave the other three in each group as the training observations. For the other two applications, we only have access to the aggregate data. Thus, we apply the data-thinning procedure described in (ref) with $m_i=\lfloor 0.2n_i\rfloor$ for the application on reporting the innovation rates and $m_{i1}=\lfloor 0.2n_{i1}\rfloor, m_{i2}=\lfloor 0.2n_{i2}\rfloor$ for the application on estimating the test passing rate gaps.
The SURE surfaces based on the training data and the holdout MSE surface are plotted in Figure (ref) - (ref). Clearly, the surfaces are closed to parallel, justifying the validity of our SURE approach.
We also compare our binomial-shrinkage estimators with existing methods proposed by xie2012sure and chen2022empirical for all three empirical applications. Both papers work with Gaussian measurements $Z_i|\theta_i\sim\mathcal{N}(\theta_i,\sigma_i^2)$, where $Z_i$ are independent observations with known variances $\sigma_i^2$. Under this assumption of normality, xie2012sure derive a Stein's unbiased risk estimator (SURE) for the squared-error loss and propose two classes of shrinkage estimators: one is a linear interpolation between the maximum likelihood estimator (MLE) and the grand mean; the other is a linear interpolation between the MLE and a data-driven location. Another recent work by chen2022empirical highlights that Gaussian empirical Bayes methods commonly rely on a precision-independence assumption, which posits that the parameters of interest are independent of their known standard errors—an assumption that is often theoretically problematic and empirically rejected. Consequently, chen2022empirical propose the CLOSE framework to estimate Gaussian means with known variances that accommodates precision dependence.
For the method proposed by xie2012sure, we construct three distinct estimators. The first, labeled SURE (grand mean), linearly interpolates between the MLE and the grand mean. The second estimator, labeled SURE (weighted ML mean), is a linear interpolation between the MLE and the weighted average of cross-fitted machine learning predictions, where each prediction is weighted by the total sample size $n_i$ of observation $i$. The third estimator, labeled SURE (data-driven location), shrinks toward a data-driven location computed following the procedure proposed by xie2012sure. For the approach proposed by chen2022empirical, we implement two estimators: CLOSE-NPMLE and CLOSE-Gauss. The CLOSE-NPMLE estimator employs a nonparametric maximum likelihood (NPMLE) approach to estimate the prior distribution of $\theta_i|\sigma_i$, while CLOSE-Gauss assumes a standard Gaussian prior for the same distribution. We implemented the “close” R package to compute the posterior means using both of the CLOSE estimators. We compute posterior means under both specifications using the \texttt{close} R package closeRpackage.
We first compare the holdout mean squared error (MSE) of all the aforementioned estimators with that of our binomial-shrinkage estimator across the three empirical applications. The methods proposed by xie2012sure and chen2022empirical assume known variances; hence, following common practice, we plug in the estimated variance $Y_i/n_i(1 - Y_i/n_i)/n_i$ for the one-sample case, where $n_i$ is the total sample size and $Y_i$ is the number of positive outcomes (value $1$) for observation $i$. Similarly, for the two-sample case, we use the variance estimator $Y_{i1}/n_{i1}(1-Y_{i1}/n_{i1})/n_{i1}+Y_{i2}/n_{i2}(1-Y_{i2}/n_{i2})/n_{i2}$, where $n_{i1}$ and $Y_{i1}$ represent the sample size and count of positive outcomes from population one for observation $i$, and similarly, $n_{i2}$ and $Y_{i2}$ correspond to population two. One important note is that the CLOSE method is not applicable to the employment discrimination application. Specifically, in the training set for this application, we have $n_{i1}=n_{i2}=3$, implying that $Y_{i1}, Y_{i2}\in\{0,1,2,3\}$. Consequently, there are only a few distinct values for the variance, making it unsuitable for the nonparametric estimation required by the CLOSE framework.
Additionally, we also compute the weighted average of the estimates for all three applications. According to (ref), the weighted average of the one-sample binomial-shrinkage estimator coincides with the grand mean, i.e.
Further, note that the two-sample binomial-shrinkage estimator defined in equation (ref) can be equivalently expressed as the difference between two one-sample estimators applied separately to populations one and two, denoted as $\hat{\theta}_{i1}(\mathbf{Y};\boldsymbol{\lambda})$ and $\hat{\theta}_{i2}(\mathbf{Y};\boldsymbol{\lambda})$, respectively. Similarly, the difference between the weighted averages of these estimators, $\hat{\theta}_{i1}(\mathbf{Y};\boldsymbol{\lambda})$ and $\hat{\theta}_{i2}(\mathbf{Y};\boldsymbol{\lambda})$, corresponds exactly to the grand-mean difference between the two populations, i.e.
\paragraph{Results. } Table (ref) presents the MSE values for all estimators. Our binomial-shrinkage estimator achieves the lowest MSE across all three applications. Table (ref) reports the weighted-average estimates from the one-sample innovation application, as well as the differences between the weighted-average estimates for two populations in the two-sample applications (discrimination report and education applications), computed over the entire dataset. Notably, the weighted-average estimates from all the considered estimators differ from the grand mean calculated from the full dataset. In contrast, our binomial-shrinkage estimator satisfies this property exactly, as shown in equations (ref) and (ref).
This paper develops a compound decision approach for estimating many binomial parameters by working directly with the exact binomial distribution. We construct an approximate Stein’s unbiased risk estimator for the average mean squared error that remains valid under heterogeneous sample sizes and small counts, without relying on Gaussian approximations. For a class of machine-learning–assisted linear shrinkage estimators, we establish asymptotic optimality, regret bounds relative to the oracle, and valid inference. The estimators satisfy a reporting-consistency property and are guaranteed to improve upon the maximum likelihood estimator and any single predictive model under mild conditions. More broadly, the framework offers a computationally tractable method for combining information pooling with covariate adjustment across large collections of binomial problems. Empirical applications to discrimination detection, education outcomes, and innovation rates show that the proposed estimators deliver stable and interpretable gains in accuracy in settings with heterogeneous sample sizes, small samples, or small binomial means.
For future work, we plan to incorporate Bayes estimators into the shrinkage framework and to extend the analysis to settings with selection decisions chen2025compound, which will require a more sophisticated bias calculation for our risk estimator. In addition, if most $n_i$ are much larger than $2$, we can also add higher-order polynomials of $Y_i$ as covariates.