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.
27,566 characters · 17 sections · 3 citation commands
Identification and Inference with Min-over-max Estimators for the Measurement of Labor Market Fairness
Min-over-max style estimators arise naturally in problems studying disparities across groups. Consider the estimator that we call Demographic Parity (DP), which is the ratio of a certain outcome in the lowing-performing group to the same in the highest-performing group. One might compute such an estimator on, say, application response rates across groups where equity is desired. It is easy to compute and interpret: DP is less than or equal to 1 and the closer it is to unity, the closer we are to equity between groups. The DP metric makes equity audits simple: a value below 0.8 is indicative of inequity and calls for further investigation.
Because the min and max functions are not differentiable everywhere, min-over-max estimators run into certain issues with statistical inference and asymptotic normality is not available. In these notes, we take an approximation approach to min-over-max estimators. Using smooth and differentiable approximations to the min and max functions, where the level of approximation is decided by a tuning parameter, we control the skewness of the asymptotic distributions of the ratio estimators and recover normality. Since larger values of the tuning parameter mean better approximations of our constituent functions, but also more asymptotic skewness, we provide upper bounds on it as a function of the sample size of the dataset. In this way, as datasets get larger and larger, the approximation matters less and we can get closer to the true ratio estimator, while preserving statistical normality.
The rest of the notes is organized as follows. Section (ref) motivates the problem with the Demographic Parity estimator and formalizes the mathematical notation. Section (ref) provides the basic asymptotic result we build on. Section (ref) introduces the approximation trick that is core to this work and Section (ref) quantifies the approximation error in their asymptotic distribution. Section (ref) presents the central asymptotic result with the approximation estimators, along with how to do inference in practical settings including how to tune the approximation parameter. Sections (ref) and (ref) go over extensions and limitations of our approximation strategy, and Section (ref) concludes.
Given this setup, our metric of interest, Demographic Parity (DP for short), is defined as
and the theoretical estimand is
Note that, by definition, both $\text{DP}_0$ and $\hat{\text{DP}}$ are limited to the interval $[0, 1]$ since we have $0 \leq \min_g \left\{ a_g \right\} \leq \max_g \left\{ a_g \right\} \leq 1$ for any set of nonnegatives $\left\{ a_g \right\}_{g=1}^G$.
The central limit theorem gives us
where $\sigma_g^2 = \frac{s_g (1 - s_g)}{p_g}$.
This variance-covariance matrix can be easily estimated with regression of binary outcomes on group fixed effects and heteroskedastic standard errors. The command to use for those in R is:
To get the asymptotic distribution of the max and min of these sample means, we would typically resort to a delta method. However, the delta method requires a continuously differentiable function, which the max and min are not. For this reason, we use a “smooth” version of the max function, which is
for a parameter $\alpha > 0$. Here $\text{rsmax}$ stands for “real softmax” zhang_lipton_li_smola. This approximation function is also known as the LogSumExp function. The function has the nice property that
while maintaining differentiability. Figure (ref) provides a visualization of the approximation.
Its gradient is
The interpretation of the gradient is very simple when we look at the limiting case of $\alpha \rightarrow \infty$:
The smooth minimum is the min version of the smooth maximum, and is defined as
for $\alpha > 0$. As $\alpha \rightarrow \infty$, it converges to the true $\min$, while maintaining differentiability throughout.
Its gradient is
and the interpretation of the gradient is that it is 1 at the unique minimum and 0 everywhere else. As usual, in the event of multiple equal minima, it places equal weight on each of them (still with 0 elsewhere).
Computing the functions rsmax and softmax has problems with floating point overflow, even though each $\hat{s}_j \in [0, 1]$ because we scale them by $\alpha$, which is a large number. Since we perform an expoentiating operation, even for moderate values of $\alpha$, the term $e^{\alpha s_j}$ could result in overflow. \\
To deal with this problem, we shift all inputs to the functions by their largest value blanchard2019accurate. That is, we compute them as
for $a_{(G)} = \max_g a_g$. This shift ensures every term $\alpha(a_j - a_{(G)}) \leq 0$. Thus, the exponential terms are always bounded by 1, and therefore the computation of the functions is numerically stable. \\
Given, these numerically stable rsmax and softmax functions, we can compute the corresponding rsmin and softmin functions as usual as
Applying the delta method to equation (ref) with function
we get
where $\Sigma$ is the variance-covariance matrix from equation (ref), and $\nabla h$ is the Jacobian of function $h$ evaluated at the true means, $\left\{ a_j \right\}$. \\
The interpretation of this new variance-covariance matrix is straightforward. Consider the limiting case $\alpha \rightarrow \infty$, for the sake of exposition:
Remember that the gradient of the true max function does not exist when two (or more) elements of the input vector are equal. However, the limit of the gradient of the real softmax function, i.e. the limit of the softmax function as $\alpha \rightarrow \infty$ does exist.
The true max and min functions are differentiable except at inputs where there are multiple argmaxes and argmins respectively. In that sense, we could directly apply the delta method on the true max and min functions and obtain the same asymptotic distribution as the ones derived above, because the limit of the softmax and softmin functions is identical to the gradients of the true max and min functions. The caveat is that this equivalence only works at inputs where there aren't multiple argmaxes (or argmins) in $\left\{ s_j \right\}$. \\
What about in the cases when we do have multiple argmaxes (or argmins) in $\left\{ s_j \right\}$? The practical answer is that if we used a small enough $\alpha$, we have a smooth enough function and have perfect inference even in these edge cases. However, the lower the $\alpha$, the worse is the approximation of the true max function. Thus there is a tradeoff between having asymptotical normality (at low $\alpha$) and inferring the right estimand at the cost of asymptotic bias and non-normality (at high $\alpha$). See Appendix (ref) for more on non-differentiability. \\
We now present a formal analysis of the approximation error when using Taylor expansions in the delta method.
We now quantify the error in our asymptotic distribution from using our approximation functions and provide guidance on how to choose the approximation parameter $\alpha$ that trades off the need for a better approximation with the stability properties that come from a well-behaved estimator.
The CLT of Bernoulli sample means gives us
for vector valued $\hat{s}$ and $s$. Consider a function $g$. In our context, $g$ could be the rsmax function or the rsmin. A first-order Taylor expansion of the estimate around the true mean gives us
where $R_1$ is the remainder from the first-order approximation and is equal to $\frac{\nabla^2 g(s^\prime)}{2!} (\hat{s} - s)^2$ for some $s^\prime$ between $s$ and $\hat{s}$ (a convex combination of the two). Thus, our asymptotic distribution has two terms:
The remainder (approximation error) term, $\sqrt{N}R_1$, is
Thus, if we want the approximation error to die down with $N$, we require $q < 1/2$. But what is $q$? Remember that it is the probabilistic order of the Hessian of our function $g$. Consider the case of the rsmax function (similar arguments apply to the rsmin). Its Hessian is the derivative of the softmax function. Denote the softmax function evaluated on the $i$th element as $\mathcal{S}_i$. Then its derivative is given by
Here $\delta_{ij}$ is the Kronecker delta. Thus $\alpha$ determines the error rate of Hessian and from the above analysis, we determine that
For a discussion on using higher-order approximations, see Appendix (ref).
A couple of points to conclude this discussion:
Let $\Sigma^\prime$ denote the variance-covariance matrix in equation (ref). Using the function $g(a, b) = \frac{a}{b}$ and the delta method on equation (ref), we then get
for $\nabla g (a, b) = [1/b, -a/b^2]^\top$ evaluated at $a = \text{rsmin}_\alpha \left\{ s_j \right\}$ and $b = \text{rsmax}_\alpha \left\{ s_j \right\}$. This gives us the desired asymptotic distribution of DP.
Call this final variance $\sigma^2_{\text{DP}}$. Then an appropriate 95% two-sided confidence interval would be
As discussed above, where max and min functions are differentiable, we have perfect asymptotic normality everywhere. However, there can be significant nonnormality when the true effects are identical (when the estimands are equal for at least some groups, i.e. where the min or max functions are nondifferentiable). We present a simple visualization to see how.\\
We work with two equally sized groups from a sample of size $N=1,000,000$. First, we start with distinct true success rates at 10% and 5%. We then plot the distribution of $\sqrt{N} (\hat{\text{DP}} - \text{DP}_0)$ for soft approximations of the DP function (centered on the soft version of the ground truth) and the ratio of the true min and max functions. We plot these histograms using 100,000 simulations each in Figure (ref). As expected, we find perfect asymptotic normality everywhere. \\
Next, we plot the case where the true means are in fact equal (at 10%) in Figure (ref). Remember that this is the case where we no longer have smoothness in the DP function to claim clean delta method inference. We find that we have asymptotic normality as long as $\alpha$ is small, specifically as $\alpha < O(\sqrt{N})$. This is indeed what our theory predicted (Look at equation (ref)). When $\alpha$ is large, we see significant left-skewness is the asymptotic distribution. This arises from the fact that, when the true means are equal, the min estimator is mechanically always lower than the max estimator and no longer identifies the “true” min group.
Recall that $\text{DP} \in [0, 1]$. Therefore testing if $\text{DP} \in [0.8, 1.25]$ (the Equal Employment Opportunity Commission thresholds for parity) is the same as testing if $\text{DP} < 0.8$. We formulate the hypotheses as
The Z-statistic of interest would be
and we test if $Z < -1.645$, the 5% quantile of the standard Normal distribution. In other words, the p-value for this one-sided test is
for a standard normal CDF, $\Phi(.)$, and we reject the null if this p-value is under 0.05.
Suppose we run an A/B test with two variants, $A$ and $B$. We seek to test if the experiment caused a change in the equity distribution across variants. Then we have two sets of statistics, $(\hat{DP}_A, \frac{\hat{\sigma}_A}{\sqrt{N_A}})$ for variant $A$ and $(\hat{DP}_B, \frac{\hat{\sigma}_B}{\sqrt{N_B}})$ for variant $B$. The two-sided hypotheses are:
The test statistic now is
and we test if $|Z| < 1.96$. If not, we reject the null. The p-value for this test is $2 (1 - \Phi(|Z|))$.
If a one-sided test is desired, say to test of variant $B$ is more equal than variant $A$, like so
the test statistic remains the same, but now we test if $Z < 1.645$. If not, we reject the null. The p-value for this test is $1 - \Phi(Z)$. As a reminder, a one-sided test is more powerful because we work under the assumption that $B$ would only be significantly better than $A$, and do not test the opposite direction.
In our analysis, we worked under a Bernoulli model where we assumed iid data and that we observe successes for each trial. In an online jobs marketplace, this means we observe whether each job application was successful or not. That is, we assume we can see all job offers a member has received. In the real world, this may not be the case. We may only observe all applications a member has submitted and the final job they accepted (and updated on their profile). \\
Consider this new model:
Then, from the binomial distribution, we know that $Y_i$ is 1 with probability $1 - (1 - s_g)^{N_i}$ and 0 otherwise. That is,
Thus, within each group, the log likelihood function is
We can then obtain the MLE $\hat{s}_g$, which is asymptotically normal. Given this asymptotic normality, we can reapply the machinery we developed above to do inference on DP.
For the estimand, $DP_0 = \frac{\min_g s_g}{\max_g s_g}$, the sample ratio estimator, $\frac{\min_g \hat{s}_g}{\max_g \hat{s}_g}$ is natural. Indeed, we showed above that it is consistent and has asymptotic normality. The problem of this estimator, though, is that it assumes there is an unambiguous min and max and the only uncertainty to be modeled is in the precise magnitude of this min and max.
In finite sample, when confidence intervals of various $s_g$ may overlap, there is additional uncertainty in whether we picked the right $g$ as the argmin or argmax. This means $\min_g \hat{s}_g$ is biased downward and $\max_g \hat{s}_g$ is biased upward (To see this, recall that the max function is convex. Now apply Jensen's inequality to show the result).
Practically, if the true $DP_0 = 1$, then the estimator $\hat{\text{DP}}$ is biased downward and will have less than perfect coverage. We may also have poor coverage when some groups are especially small and thus have their sample means estimated with much imprecision.
As an alternative to the sample min and max estimators, we may use cross-fitted estimators van2013estimating. They work as follows:
We conclude by reiterating that the DP inference machinery we developed here is much more generic than the motivating Binomial model suggests. It can be applied for the min and maxes of any parameters, for whose estimators we have asymptotic normality, which is several of the most popular econometric estimators.