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.
87,414 characters · 13 sections · 7 citation commands
The Sorted Effects Method: Discovering Heterogeneous Effects Beyond Their Averages
In nonlinear and interactive linear models the partial (ceteris paribus) effects of interest often vary with respect to the underlying covariates. For example, consider a binary response model with conditional choice probability ${\mathrm{P}}(Y = 1 \mid X) = F(X^\mathsf{T} \beta)$, where $Y$ is a binary response variable, $X$ is a vector of covariates, $F$ is a distribution function such as the standard normal or logistic, and $\beta$ is a vector of coefficients. The partial or predictive effect (PE) of a marginal change in a continuous covariate $X_j$ with coefficient $\beta_j$ on the conditional choice probability is $$\Delta(X) = f(X^\mathsf{T} \beta) \beta_j, \quad f(v) = \partial F(v)/\partial v,$$ which generally varies in the population of interest with the covariate vector $X$, as $X$ varies according to some distribution, say $\mu$. A common empirical practice is to report the average partial effect (APE), $${\mathrm{E}} [\Delta (X)] = \int \Delta(x) d\mu(x),$$ as a single summary measure of the PE (e.g., \citeasnoun[Chap. 2]{wooldridge:text}), or to report effects for some groups (e.g., \citeasnoun{angrist2008mostly}). However, the APE completely disregards the heterogeneity of the PE and may give a very incomplete picture of the impact of the covariates.
In this paper we propose complementing the APE by reporting the entire set of PEs sorted in increasing order and indexed by a ranking with respect to the distribution of the covariates in the population of interest. These sorted effects correspond to percentiles of the PE,
and provide a more complete representation of the heterogeneity of $\Delta(X)$. We shall call these effects as sorted predictive or partial effects (SPE) by default, as most models are predictive.\footnote{When the underlying model has a structural or causal interpretation, we may use the name sorted structural effects or sorted treatment effects.} We also show how to use the SPEs to carry out classifications analysis (CA). This analysis consists of classifying the observational units into most or least affected depending on whether their PEs are above or below some tail SPE, and then comparing the moments or distribution of the covariates of the most and least affected groups.
Heterogeneous effects also arise in the most basic linear models with interactions oaxaca73,cox84. Consider a conditional mean model for the Mincer earnings function: $$ Y = P(T,W)^\mathsf{T} \beta + \epsilon, \quad {\mathrm{E}}[\epsilon \mid T, W] = 0, \quad X = (T, W), $$ where $Y$ is log wage, $T$ is an indicator of gender (or race, treatment, or program participation), and $W$ is a vector of labor market characteristics. The vector $P(T,W)$ is a collection of transformations of $T$ and $W$, involving some interaction between $T$ and $W$. For example, \citeasnoun{oaxaca73} used the specification $P(T,W) = (TW, (1-T)W)$. Then, the PE of changing $T=0$ to $T=1$ is $$ \Delta(X) = P(1,W)^\mathsf{T} \beta - P(0,W)^\mathsf{T} \beta,$$ which is a measure of the gender wage gap conditional on worker characteristics. The function $u \mapsto \Delta^*_\mu(u)$ provides again a complete summary of the entire range of PEs. The left panel of Figure (ref) illustrates the SPE of the conditional gender wage gap for women.
The SPE varies sharply from around $-40$ to $6.5\%$, and does not coincide with the average PE of $-20\%$. The PE is especially (negatively) large for women who have any of the following characteristics: married, low educated, high experience, and working on sales occupations -- this follows from the classification analysis, where we compare the average characteristics of the subpopulations of women with covariate values $X$ such that $\Delta(X)$ is above the 90% percentile and below the 10% percentile. We refer the reader to Section (ref) for a detailed discussion of this example.
The general settings that we deal with in this paper as well as the specific results we obtain are as follows: Let $X$ denote a covariate vector, $\Delta(X)$ denote a generic PE of interest, $\mu$ denote the distribution of $X$ in the population of interest, and $\mathcal{X}$ denote the interior of the support of $X$ in this population. The SPE is obtained by sorting the multivariate function $x \mapsto \Delta(x)$ in increasing order with respect to $\mu$. Using tools from differential geometry, we prove that this multivariate sorting operator is Hadamard differentiable with respect to the PE function $\Delta$ and the distribution $\mu$ at the regular values of $x \mapsto \Delta(x)$ on $\mathcal{X}$. This key and new mathematical result allows us to derive the large sample properties of the empirical SPE, which replace $\Delta$ and $\mu$ by sample analogs, obtained from parametric or semi-parametric estimators, using the functional delta method. In particular, we derive a functional central limit theorem and a bootstrap functional central limit theorem for the empirical SPE. The main requirement of these theorems is that the empirical $\Delta$ and $\mu$ also satisfy functional central limit theorems, which hold for many estimators used in empirical economics under general sampling conditions. We use the properties of the empirical SPE to construct confidence sets for the SPE that hold uniformly over quantile indices. We also show under the same conditions that the empirical version of the objects in the classification analysis follow functional central limit theorems and bootstrap functional central limit theorems. We derive these result by establishing the Hadamard differentiability of a classification operator related to the multivariate sorting operator.
\paragraph{Related technical literature: } Previously, \citeasnoun{CFG-10} derived the properties of the rearrangement (sorting) operator in the univariate case with known $\mu$ (standard uniform distribution). Those results were motivated by a completely different problem -- namely, restoration of monotonicity in conditional quantile estimation -- rather than the problem of summarizing heterogeneous effects by the SPEs. These prior technical results are not applicable to our case as soon as the dimension of $X$ is greater than one, which is the case in all modern applications where effects are of interest. Moreover, the previous results are not applicable even in the univariate case since the measure $\mu$ is not known in all envisioned applications. The properties of the sorting operator are different in the multivariate case and require tools from differential geometry: computation of functional (Hadamard) derivatives of the sorting operator with respect to perturbations of $\Delta$ require us to work with integration on $(d_x-1)$-dimensional manifolds of the type $\{\Delta(x) = \delta\}$, where $d_x = \dim X$. Moreover, we also need to compute functional derivatives with respect to suitable perturbations of the measure $\mu$. In econometrics or statistics, \citeasnoun{sasaki-15} also used differential geometry to characterize the structural properties of derivatives of conditional quantile functions in nonseparable models; and \citeasnoun{kim:pollard} used tools from differential geometry to derive the large sample properties of the maximum score and other cube root consistent estimators. Relative to these papers, we share the use of differential geometry tools as a general proof strategy, but we apply these tools to establish the analytical properties of different functionals -- namely, the SPEs. Moreover, our results on the functional differentiability of the sorting and classification operators in the multivariate case constitute new mathematical results, which are of interest in their own right.
\paragraph{Organization of the paper: } In Section (ref) we discuss the quantities of interest in nonlinear and interactive linear models with examples; introduce the SPE and related CA, along with their empirical counterparts; and outline the main inferential results. In Section (ref) we provide an empirical application to the gender wage gap in the U.S. in 2015. We derive the properties of the empirical SPE and CA in large samples and show how to use these properties to make inference in Section (ref). Appendix (ref) provides some key mathematical results on the differentiability of the multivariate sorting and classification operators and Appendix (ref) contains the proof of the main results. All other proofs are given in the online appendix with supplementary material (SM), which also contains additional technical material, and results from Monte Carlo simulations and an empirical application to mortgage denials using binary response models cfy17sup.
We start by discussing the objects of interest in nonlinear and interactive linear models.
We consider a general model characterized by a predictive function $g(X)$, where $X$ is a $d_x$-vector of covariates that may contain unobserved components, as in quantile regression models. The function $g$ usually arises from a model for a response variable $Y$, which can be discrete or continuous. We call the function $g$ predictive because the underlying model can be either predictive or causal under additional assumptions, but we do not insist on estimands having a causal interpretation. For example, in a mean regression model, $g(X) = {\mathrm{E}}[Y \mid X]$ corresponds to the expectation function of $Y$ conditional on $X$; in a binary response model, $g(X) = {\mathrm{P}}[Y = 1 \mid X]$ corresponds to the choice probability of $Y=1$ conditional on $X$; in a quantile regression model, $g(X) = {\mathrm{Q}}_{Y}[\epsilon \mid Z]$, where the covariate $X = (\epsilon,Z)$ consists of the unobservable rank variable $\epsilon$ with a uniform distribution, $\epsilon \mid Z \sim U(0,1)$, and the observed covariate vector $Z$, and where ${\mathrm{Q}}_{Y}[ \tau \mid Z]$ is the conditional $\tau^{th}$-quantile of $Y$ given $Z$.
Let $X = (T,W)$, where $T$ is the key covariate or treatment of interest, and $W$ is a vector of control variables. We are interested in the effects of changes in $T$ on the function $g$ holding $W$ constant. These effects are usually called partial effects, marginal effects, or treatment effects. We call them predictive effects (PE) throughout the paper, as such a name most accurately describes the meaning of the estimand (especially when a causal interpretation is not available). If $T$ is discrete, the PE is
where $t_1$ and $t_0$ are two values of $T$ that might depend on $t$ (e.g., $t_0 = 0$ and $t_1 = 1$, or $t_0 = t$ and $t_1 = t + 1$). This PE measures the effect of changing $T$ from $t_0$ to $t_1$ holding $W$ constant at $w$. If $T$ is continuous and $t \mapsto g(t,w)$ is differentiable, the PE is
where $\partial_{t}$ denotes $\partial/ \partial t$, the partial derivative with respect to $t$. This PE measures the effect of a marginal change of $T$ from the level $t$ holding $W$ constant at $w$.\footnote{We can also consider high-order and crossed effects. For example, $\Delta(x) = \partial^2_{t^2} g(t,w)$ gives the second-order PE of the continuous treatment $T$ if $t \mapsto g(t,w)$ is twice differentiable; and, letting $X=(T,S,W)$ where $T$ and $S$ are discrete, $\Delta(x) = g(t_1,s_1,w) - g(t_0,s_1,w)-g(t_1,s_0,w) + g(t_0,s_0,w)$ gives the crossed effect or interaction of $T$ and $S$.}
We consider the following examples in the empirical applications of Section (ref) and SM.
The set of examples listed above are the most basic, leading cases, arising mostly in predictive analysis and program evaluation. Our theoretical results are rather general and are not limited to these cases. Thus, they allow for both $\Delta$ and $\mu$ to originate from causal or structural models and to be estimated by structural methods. For example, in treatment effects models with selection on observables rr83, the PE is the conditional average treatment effect $\Delta(x) = {\mathrm{E}}(Y_1 - Y_0 \mid X=x),$ where $Y_1$ and $Y_0$ are potential outcomes in the treated and non-treated statuses and $X$ is a vector of covariates. The standard approach is to aggregate the conditional average treatment effects by integration with respect to the distribution of the covariates in the population of interest $\mu$. This yields the average treatment effect if $\mu$ is the distribution in the entire population or the average treatment effect on the treated if $\mu$ is the distribution in the treated subpopulation. The SPE can be used to complement the analysis by reporting the entire range of conditional average treatment effects, and also to determine the optimal treatment allocation with budget constraints. Thus, \citeasnoun{bd12} showed that under some conditions this optimal allocation has a cutoff determined by a tail percentile of the conditional average treatment effects, i.e. by a SPE. Another example is the welfare analysis described in \citeasnoun{hn17}, where $\Delta(x)$ is the compensating or equivalent variation of a price change conditional on covariates such as income and demographic characteristics, and $\mu$ is the distribution of covariates in the population of interest.
In all the previous examples, the PE $\Delta(x)$ is a function of $x$ and therefore can be different for each observational unit. To summarize this effect in a single measure, a common practice in empirical economics is to average the PEs. Averaging, however, masks most of the heterogeneity in the PE allowed by nonlinear or interactive linear models. We propose reporting the entire set of values of the PE sorted in increasing order and indexed by a ranking $u \in [0,1]$ with respect to the population of interest. These sorted effects provide a more complete representation of the heterogeneity in the PE than the average effects.
The $u$-SPE is the $u^{th}$-quantile of $\Delta(X)$ when $X$ is distributed according to $\mu$. As for the average effect, $\mu$ can be chosen to select a target subpopulation from the entire population. For example, when $T$ is a treatment indicator:
By considering $\Delta_{\mu}^*(u) $ at multiple quantile indices, we obtain a one-dimensional representation of the heterogeneity of the PE. Accordingly, our object of interest is the SPE-function $$ \{u \mapsto \Delta_{\mu}^*(u) : u \in \mathcal{U} \}, \ \ \mathcal{U} \subseteq [0,1], $$ where $\mathcal{U}$ is the set of quantile indices of interest.
We also show how to use the $u$-SPE for classification analysis. Let $u \in \mathcal{U},$ with $u < 1/2$, and $Z$ be a $d_z$-dimensional random vector that includes $X$ and possibly other variables such as $Y$ in Examples 1--3. By abuse of notation, we also denote the distribution of $Z$ over its support $\mathcal{Z}$ as $\mu$.
In practice, we replace the PE $\Delta$ and the distribution $\mu$ by sample analogs to construct plug-in estimators of the SPE. Let $\widehat{\Delta}(x)$ and $\widehat{\mu}(x)$ be estimators of $\Delta(x)$ and $\mu(x)$ obtained from $\{(Y_i,T_i,W_i) : 1 \leqslant i \leqslant n \}$, an independent and identically distributed sample of size $n$ from $(Y,T,W)$.
Example 1 (Binary response model, cont.) \ The estimator of the PE is
where $\widehat \beta$ is the maximum likelihood (ML) estimator of $\beta$, $$ \widehat \beta \in \arg \max_{b \in \mathbb{R}^{d_p}} \sum_{i=1}^n [Y_i \log F(P(T_i,W_i) ^\mathsf{T} b) + (1-Y_i)\log \{1 - F(P(T_i,W_i) ^\mathsf{T} b)\}], \ \ d_p = \dim P(T,W). $$\qed
Example 2 (Interactive linear model with additive error, cont.) \ The estimator of the PE is
where $\widehat \beta$ is the ordinary least squares (OLS) estimator of $\beta$, $$ \widehat \beta \in \arg \min_{b \in \mathbb{R}^{d_p}} \sum_{i=1}^n [Y_i - P(T_i,W_i) ^\mathsf{T} b]^2, \ \ d_p = \dim P(T,W). $$\qed
Example 3 (Linear model with non-additive error, cont.) \ The estimator of the PE is
where $\widehat \beta(\tau)$ is the \citeasnoun{Koenker:1978} quantile regression (QR) estimator of $\beta(\tau)$, $$ \widehat \beta(\tau) \in \arg \min_{b \in \mathbb{R}^{d_p}} \sum_{i=1}^n \rho_{\tau}(Y_i - P(T_i,W_i) ^\mathsf{T} b), \ \ d_p = \dim P(T,W), \ \ \rho_{\tau}(v) = (\tau - 1\{v < 0\})v. $$ \qed
The empirical version of the $u$-CA classifies the observations in the sample using the empirical PEs and $u$-SPE, and computes the moments and distributions in the resulting most and least affected subsamples.
The main inferential result for the SPE can be previewed as follows. Assume that the PE function $x \mapsto \Delta(X)$ is not locally flat in the sense that the norm of its gradient does not vanish anywhere over the support, and other regularity conditions stated in Section 4. Then, the empirical SPE-process is $\sqrt{n}$-consistent and converges in distribution to a centered Gaussian process, namely
the metric space of bounded functions on $\mathcal{U}$, as a stochastic process indexed by $u \in \mathcal{U}$, where $\mathcal{U}$ is a compact subset of $(0,1)$. Moreover, the exchangeable bootstrap algorithm specified in Algorithm (ref) estimates consistently the law of $Z_{\infty}(u)$.
The next corollary to Theorem (ref) in Section (ref) provides uniform bands that cover the SPE-function simultaneously over a region of values of $u$ with prespecified probability in large samples. It does cover pointwise confidence bands for the SPE-function at a specific quantile index $u$ as a special case by simply taking $\mathcal{U}$ to be the singleton set $\{u\}$.
We now describe a practical bootstrap algorithm to estimate the quantiles of $t(\mathcal{U})$. Let $(\omega_{1}, \ldots, \omega_{n})$ denote the bootstrap weights, which are nonnegative random variables independent of the data obeying the conditions stated in \citeasnoun{vdV-W}. For example, $(\omega_{1}, \ldots, \omega_{n})$ is a multinomial vector with dimension $n$ and probabilities $(1/n,\ldots,1/n)$ in the empirical bootstrap. In what follows $B$ is the number of bootstrap draws, such that $B \to \infty$. In our experience, setting $B \geqslant 500$ suffices for good accuracy.
Let $\Lambda_{\Delta,\mu}^{u}(t) := [\Lambda_{\Delta,\mu}^{-u}(t),\Lambda_{\Delta,\mu}^{+u}(t)]$ and $\widehat{\Lambda}_{\Delta,\mu}^{u}(t) := [\widehat{\Lambda}_{\Delta,\mu}^{-u}(t), \widehat{\Lambda}_{\Delta,\mu}^{+u}(t)]$. The main inferential result for CA can be previewed as follows: the empirical CA-process converges in distribution to a centered bivariate Gaussian process, namely
as a stochastic process indexed by $t \in \mathbb{R}^{d_z}$. Moreover, exchangeable bootstrap estimates consistently the law of $Z_{\infty}^{u}(t)$.
The next corollary to Theorem (ref) in Section (ref) provides uniform bands that cover $L$ linear combinations of the $2$-dimensional vector $\Lambda_{\Delta,\mu}^{u}(t)$ with coefficients $c_1, \ldots, c_L$ simultaneously over $t \in \mathcal{T}$ with prespecified probability in large samples. It covers pointwise confidence intervals for the mean of the $k^{th}$ component of $Z$ for least affected as a special case with $L=1$ linear combination, $c_1 = (1,0)'$, and $\mathcal{T} = \{e_k\}$, where $e_k$ is a unit vector with a one in the $k^{th}$ position. Joint confidence intervals for $s$ differences of means of the $k_1^{th}$, $\ldots$, $k_s^{th}$ components of $Z$ between most and least affected are a special case with $L=1$ linear combination, $c_1 = (-1,1)'$, and $\mathcal{T} = \{e_{k_1}, \ldots, e_{k_s}\}$. Joint uniform bands for the distribution of the $k^{th}$ component of $Z$ for most and least affected are also a special case with $L=2$ linear combinations, $c_1 =(1,0)$, $c_2 = (0,1)$, and $\mathcal{T} = \{ t \in \mathbb{R}^{d_z}: t_j = \bar T, j \neq k \}$, where $\bar T$ is an arbitrarily large number. By appropriate choice of the linear combinations and the index set $T$, we can therefore conduct multiple tests while preserving the significance level from simultaneous inference problems rsw10b,rsw10,lsx16. We show examples in the empirical application of Section (ref).
In addition to moments and distributions, we can conduct inference on the subpopulations of most and least affected.\footnote{Here we follow the set inference approach described in \citeasnoun{ckm15}, which builds on \citeasnoun{cht:bounds}. In addition our results justify the use of subsapling-based methods as in \citeasnoun{cht:bounds} and \citeasnoun{RSset}.} Let $$ \mathcal{M}^{-u}:=\{ (x,y) \in \mathcal{Z}: \Delta(x) \leqslant \Delta^*_\mu (u) \}, \quad \mathcal{M}^{+u} := \{ (x,y) \in \mathcal{Z} : \Delta(x) \geqslant \Delta^*_\mu(1-u) \}, $$ be the sets representing the $u$-least and $u$-most affected subpopulation, respectively. Here we assume that $\mathcal{Z}$ is compact or that the support of $(X,Y)$ has been intersected with a compact set to form $\mathcal{Z}$. We can construct an outer $(1-\alpha)$-confidence set for $\mathcal{M}^{-u}$ as\footnote{Note that we can also similarly construct an inner confidence region, which is the complement of the outer confidence region of $\mathcal{X}\setminus \mathcal{M}^{-u}$, see \citeasnoun{ckm15} for relevant discussion.} $$ \mathcal{CM}^{-u}(1-\alpha) = \{ (x,y) \in \mathcal{Z}: \widehat \Sigma^{-1/2}(x,u) \sqrt{n}[ \widehat \Delta(x) - \widehat \Delta^*_\mu (u) ] \leqslant \widehat c(1-\alpha) \}, $$ where $\widehat c(1-\alpha)$ is a consistent estimator of $c(1-\alpha)$, the $(1-\alpha)$-quantile of the random variable $$ V_{\infty} = \sup_{ \{x \in \mathcal{X}: \Delta(x) = \Delta^*_\mu (u)\}} \Sigma^{-1/2}(x,u) [ G_\infty (x) - Z_\infty(u) ], $$ and $x \mapsto \widehat \Sigma(x,u)$ is a uniformly consistent estimator of $x \mapsto \Sigma(x,u)$, the variance function of the process $G_\infty(x) - Z_\infty(u)$ defined in Section 4. The estimator $\widehat c(1-\alpha)$ can be obtained as the $(1-\alpha)$-quantile of the bootstrap version of $V_{\infty}$, $$ \widetilde V^*_{\infty} = \sup_{ \{ x \in \mathcal{X}: \widehat \Delta(x) = \widehat \Delta^*_\mu (u) \}} \widehat \Sigma^{-1/2}(x,u) \sqrt{n} \Big ( [\widetilde \Delta (x) - \widetilde{\Delta_{\mu}^*}(u)] - [ \widehat \Delta (x) - \widehat \Delta^*_\mu(u)] \Big ), $$ where $\widetilde \Delta(x)$ and $\widetilde{\Delta_{\mu}^*}(u)$ are defined as in Algorithm (ref). A similar $(1-\alpha)$-confidence set, $\mathcal{CM}^{+u}(1-\alpha),$ can be constructed for $\mathcal{M}^{+u}$. These sets can be visualized by plotting all 2 or 3 dimensional projections of their elements. We provide an example of such plots in Section (ref). An immediate consequence of the set inference results in \citeasnoun{ckm15} and the results of this paper is the following corollary:
We report the main results of the application to the gender wage gap using data from the U.S. March Supplement of the Current Population Survey (CPS) in 2015. In Appendix (ref) of the SM, we complement the analysis with supporting results from a simulation calibrated to this application. There, we find that our estimation and inference methods perform well in finite samples that closely mimic the characteristics of the CPS data. This exercise serves to indirectly verify the plausibility of the main regularity conditions mentioned in Section (ref) and formally stated in Section (ref).
Our sample consists of white, non-hispanic individuals who are aged 25 to 64 years and work more than 35 hours per week during at least 50 weeks of the year. We exclude self-employed workers; individuals living in group quarters; individuals in the military, agricultural or private household sectors; individuals with inconsistent reports on earnings and employment status; individuals with allocated or missing information in any of the variables used in the analysis; and individuals with hourly wage rate below $\$3$. The resulting sample contains $32,523$ workers including $18,137$ men and $14,382$ of women.
We estimate interactive linear models with additive and non-additive errors, using mean and quantile regressions, respectively. The outcome variable $Y$ is the logarithm of the hourly wage rate constructed as the ratio of the annual earnings to the total number of hours worked, which is constructed in turn as the product of number of weeks worked and the usual number of hours worked per week. The key covariate $T$ is an indicator for female worker, and the control variables $W$ include 5 marital status indicators (widowed, divorced, separated, never married, and married); 5 educational attainment indicators (less than high school graduate, high school graduate, some college, college graduate, and advanced degree); 4 region indicators (midwest, south, west, and northeast); a quartic in potential experience constructed as the maximum of age minus years of schooling minus 7 and zero, i.e., $experience = \max(age-education-7,0)$; 5 occupation indicators (management, professional and related; service; sales and office; natural resources, construction and maintenance; and production, transportation and material moving); 12 industry indicators (mining, quarrying, and oil and gas extraction; construction; manufacturing; wholesale and retail trade; transportation and utilities; information; financial services; professional and business services; education and health services; leisure and hospitality; other services; and public administration); and all the two-way interactions between the education, experience, occupation and industry variables except for the occupation-industry interactions.\footnote{The sample selection criteria and the variable construction follow \citeasnoun{Mulligan-Rubinstein-08}. The occupation and industry categories follow the 2010 Census Occupational Classification and 2012 Census Industry Classification, respectively.} All calculations use the CPS sampling weights to account for nonrandom sampling in the March CPS.
Table (ref) reports sample means of the variables used in the analysis. Working women are more highly educated than working men, have about the same potential experience, and are less likely to be married and more likely to be divorced. They work relatively more often in managerial and sales occupations and in the industries providing education and health services. Working men are relatively more likely to work in construction and production occupations within non-service industries. The unconditional gender wage gap is 23%.
Figure (ref) of Section (ref) plots estimates and 90% confidence bands for the APE and SPE-function on the treated (women) of the conditional gender wage gap using additive and non-additive error models. The PEs are obtained as described in Examples (ref) and (ref) with $P(T,W) = (TW,(1-T)W)$. In this case $\dim P(T,W) = 332$, which makes it very difficult to identify any pattern about the gender wage gap just by looking at the regression coefficients. The distribution $F_{T,W}$ is estimated by the empirical distribution of $(T,W)$ for women, and $F_{\epsilon}$ is approximated by a uniform distribution over the grid $\{.02, .03, \ldots, .98\}$. The confidence bands are constructed using Algorithm (ref) with standard exponential weights (weighted bootstrap) and $B=500$, and are uniform for the SPE-function over the grid $\mathcal{U} = \{.01, .02, \ldots, .98\}$. We monotonize the bands using the rearrangement method described in Remark (ref), and implement the finite sample corrections described in Remark (ref). After controlling for worker characteristics, the gender wage gap for women remains on average around 20%. More importantly, we uncover a striking amount of heterogeneity, with the PE ranging between -6.5 and 40% in the additive error model and between -14 and 54% in the non-additive error model.\footnote{In the 2016 version of the paper we found similar patterns of heterogeneity using CPS 2012 data with a specification that did not include occupation and industry indicators.}
Table (ref) shows the results of a classification analysis, exhibiting characteristics of women that are most and least affected by the gender wage gap together with standard errors obtained by weighted bootstrap. We focus here on the non-additive model, but the results from the additive model are similar. Since the PE are predominantly negative, we define the most affected as $\Delta(X) < \Delta^*_{\mu}(u)$ and the lest affected as $\Delta(X) > \Delta^*_{\mu}(1-u)$ to facilitate the interpretation. According to this model the 10% of the women most affected by the gender wage gap on average earn lower wages, are much more likely to be married, much less likely to be never married, have lower education, live in the South, possess much more potential experience, are more likely to have sales and non managerial occupations, and work more often in manufacture and retail and less often in education industries than the 10% least affected women.
Table (ref) tests if the differences found in table (ref) are statistically significant. It reports p-values for the test of equality of means for most and least affected women. The first p-value accounts for simultaneous inference on all variables within a given category. For example, it accounts that we are conducting five tests corresponding to the five categories of marital status. For the non categorical variables log wage and experience the p-values are for one test. The second p-value accounts for simultaneous inference of all the differences displayed in the table.\footnote{We employ the so called "single-step" methods for controlling the family-wise error rate. To generate a (somewhat) higher power, we recommend to employ the p-values generated via "step-down" methods, such as those reported in \citeasnoun{RWpvals} and \citeasnoun{lsx16}.} These p-values are obtained by Algorithm (ref) with the appropriate choice of vectors of linear combinations and set $\mathcal{T}$, and 500 weighted bootstrap repetitions. The p-values show that most of the differences from table (ref) are statistically significant at conventional significant levels after controlling for simultaneous inference. In particular, the most affected women are significantly more likely to be married, high-school graduates, more experienced, and in sales occupations, and less likely to be never married and in managerial occupations under the most strict simultaneous inference correction. \citeasnoun{bk17} have recently documented the importance of differences in occupation and industry to explain the gender wage gap using data from the Panel Study of Income Dynamics (PSID) 1980-2010 and a different methodology based on wage decompositions. Consistent with our findings, they argue that this importance might be due to compensating differentials. Unlike \citeasnoun{bk17} and previous studies in the literature, our analysis uncovers significant heterogeneity in the extent of the gender wage gap and relates this heterogeneity to human capital, occupation, industry and other characteristics.
We further explore these findings by analyzing the APE and SPE on the treated conditional on marital status and unobserved rank in the non-additive error model. Figures (ref) and (ref) show estimates and 90% confidence bands of the APE and SPE-function of the gender wage gap for 2 subpopulations defined by marital status (married and never married) and 3 subpopulations defined by unobserved rank (first decile, median and ninth decile, where the unobserved rank is .1, .5 and .9, respectively). The confidence bands are constructed as in fig. (ref). We find significant heterogeneity in the gender gap within each subpopulation, and also between subpopulations defined by marital status and unobserved rank. The SPE-function is more negative for married women and at the tails of the conditional distribution. Married women at the top decile suffer from the highest gender wage gaps. This pattern is consistent with “glass-ceiling” effects behind the gender wage gap abv03.
Figure (ref) plots simultaneous 90% confidence bands for the distribution of experience and log wage for the most and least affected women. They are obtained by Algorithm (ref) with 500 weighted bootstrap replications. The estimated distribution of experience for the most affected first-order stochastically dominates the same estimated distribution for the least affected women. Moreover, the uniform bands confirm that this dominance is statistically significant at the 90% confidence level for the underlying distributions. The estimated (marginal) distribution of log wage for the least affected first-order dominates the same estimated distribution for most affected, but we cannot reject that the underlying distributions are equal at the 10% significance level. The results of the classification analysis are consistent with preferences that make never married highly educated young women working on managerial occupations be more career-oriented.\footnote{We find similar results using the additive error model. We do not report these results for the sake of brevity.}
Finally, Figure (ref) plots two dimensional projections of experience-log wage and experience-marital status of the confidence sets for the 10% most and least affected subpopulations. We show the results from the additive error model for the conditional expectation. Here we use a simplified specification that excludes the two-way interactions from $W$ to get more precise estimates of all the PEs. We obtain 90% confidence sets for the most and least affected subpopulations by weighted bootstrap with standard exponential weights and 500 repetitions. The sets $\mathcal{CM}^{-0.1}(0.90)$ and $\mathcal{CM}^{+0.1}(0.90)$ include $23\%$ and $19\%$ of the women in the sample, respectively.\footnote{Recall that in this application the set $\mathcal{CM}^{-0.1}(0.90)$ corresponds to most affected women and $\mathcal{CM}^{+0.1}(0.90)$ to least affected women. We drop one woman that is included in both sets.} The projections show that there are relatively more least affected women with low experience at all wage levels, more high affected women with high wages with between 15 and 25 years of experience, and more least affected women which are not married at all experience levels.
For an open set $\mathcal{K}$, let the class $\mathcal{C}^1$ on $\mathcal{K}$ denote the set of continuously differentiable real valued functions on $\mathcal{K}$. We make the following technical assumptions about the PE function $\Delta: \Bbb{R}^{d_x} \mapsto \Bbb{R}$ and the distribution of the covariates:
${\sf S}.1$. The part of the domain of the PE function $x \mapsto \Delta(x)$ of interest, $\mathcal{X}$, is open and its closure $\overline{\mathcal{X}}$ is compact. The distribution $\mu$ is absolutely continuous with respect to the Lebesgue measure with density $\mu'$. There exists an open set $B(\mathcal{X})$ containing $\overline{\mathcal{X}}$ such that $x \mapsto \Delta(x)$ is $\mathcal{C}^1$ on $B(\mathcal{X})$, and $x \mapsto \mu'(x)$ is continuous on $B(\mathcal{X})$ and is zero outside the domain of interest, i.e. $\mu'(x) = 0$ for any $x \in B(\mathcal{X}) \setminus \mathcal{X}$.
${\sf S}.2$. Let $\mathcal{M}_{\Delta}(\delta):= \{x \in {\mathcal{X}}: \Delta(x)=\delta\}$. For any regular value $\delta$ of $\Delta$ on $\overline{\mathcal{X}}$, we assume that the closure of $\mathcal{M}_{\Delta}(\delta)$ has a finite number of connected branches.
The following property of the set $\mathcal{M}_{\Delta}(\delta)$ is a useful implication of Assumptions ${\sf S}.1$ and ${\sf S}.2$ that we will exploit in the analysis.
Assumption ${\sf S}.1$ imposes mild smoothness conditions on the PE function $x \mapsto \Delta(x)$. It also requires that all the components of the covariate $X$ are continuous random variables. We defer the treatment of the case where $X$ has both continuous and discrete components to the SM. As a matter of generalization, our theoretical analysis allows us to replace that $x \mapsto \mu'(x)$ vanishes on $\partial \mathcal{X}$, by the weaker condition that the intersection of $\mathcal{M}_\Delta(\delta)$ and the boundary of $\mathcal{X}$ have zero volume with respect to $\mu$, namely
where $\partial \mathcal{X}$ denotes the boundary of $\overline \mathcal{X}$, $\partial \Delta(x)$ is the gradient of $x \mapsto \Delta(x)$, and $\int_{\mathcal{M}} f(x) d\mathrm{Vol}$ denotes the integral of the function $f$ on the manifold $\mathcal{M}$ with respect to volume; see Appendix (ref) in the SM for a brief review on Differential Geometry. This relaxation is relevant to cover the case where $X$ includes an uniformly distributed component such as the unobserved rank in Example (ref).\footnote{ In the numerical examples of Section (ref) in the SM, the first two designs only satisfy this relaxed condition.}
Assumption ${\sf S}.2$ imposes shape restrictions on $x \mapsto \Delta(x)$ that rule out cases such as infinite cyclical oscillations or flat areas. A simple sufficient condition for ${\sf S}.2$ is that the map $x \mapsto \Delta(x)$ does not have critical points on $\overline \mathcal{X}$. This means that $x \mapsto \Delta (x)$ is not locally flat anywhere on $\overline \mathcal{X}$, which we define to mean that the norm of the gradient, $\| \partial \Delta(x) \|$, does not vanish on $ x\in \overline \mathcal{X}$. In this case, any $\delta$ in the image of $\overline \mathcal{X}$ under $\Delta$ is regular. This condition is probably the most relevant for practice and can be verified in applications, at least informally.
We make the following assumptions about the estimator of the PE. Let $\ell^{\infty}(\mathcal{T})$ denote the set of bounded and measurable functions $g:\mathcal{T} \to \mathbb{R}$ and $\mathcal{F}$ a fixed subset of continuous functions on $B(\mathcal{X})$. Let $\ell^{\infty} (B(\mathcal{X}))$ be the set of bounded and measurable functions on $B(\mathcal{X})$ and $\rightsquigarrow$ denote weak convergence (convergence in distribution).
${\sf S}.3$. $\widehat{\Delta}$, the estimator of $\Delta$, belongs to $\mathcal{F}$ with probability approaching $1$ and obeys a functional central limit theorem, namely, $$a_n(\widehat{\Delta}-\Delta)\rightsquigarrow G_{\infty} \text{ in } \ell^{\infty} (B(\mathcal{X})),$$ where $a_n$ is a sequence such that $a_n \to \infty$ as $n \to \infty$, and $x \mapsto G_{\infty}(x)$ is a tight process that has almost surely uniformly continuous sample paths on $B(\mathcal{X})$.
In the parametric and semiparametric models of Examples 1--3, ${\sf S}.3$ holds under weak conditions that guarantee asymptotic normality of the ML, OLS and QR estimators. For the QR estimator in Example (ref) where the unobserved rank is one of the covariates, these conditions include that the density of $Y$ conditional on $X$ be bounded away from zero koenker:book, which is facilitated by excluding tail quantile indexes.
Let $\widehat \mu$ be the estimator of the distribution $\mu$. It is convenient to identify $\mu$ and $\widehat \mu$ with the operators: $$g \mapsto \mu (g) =\int g(x) d\mu(x), \quad g \mapsto \widehat \mu (g) =\int g(x) d\widehat \mu(x), $$ mapping from the set $\mathcal{G} := \{ x \mapsto 1( f(x) \leqslant \delta): f \in \mathcal{F}, \delta \in \mathcal{V}\}$ to $\Bbb{R}$, where $\mathcal{F}$ is the fixed subset of continuous functions on $B(\mathcal{X})$ containing $\Delta$, and $\mathcal{V}$ is any compact set of $\Bbb{R}$. We require $\mathcal{G}$ to be totally bounded under the $L^2(\mu)$ norm. Define $\mathbb{H}$ as the set of all bounded linear operators $H$ on $\mathcal{G}$ of the form $$ g \mapsto H(g), $$ which are uniformly continuous on $g \in \mathcal{G}$ under the $L^2(\mu)$ norm. We define the boundedness of these operators with respect to the norm: $$\|H\|_{\mathcal{G}}=\sup_{g\in \mathcal{G}}{|H(g)|},$$ and define the corresponding distance between two operators $H$ and $\widetilde H$ in $\mathbb{H}$ as $\|H - \widetilde H \|_{\mathcal{G}} = \sup_{g\in \mathcal{G}} |H(g)-\widetilde H(g)|$. Clearly, $\mu \in \mathbb{H}$.
We make the following assumption about $\widehat \mu$.
${\sf S}.4$. The function $x \mapsto \widehat \mu(x)$ is a distribution over $B(\mathcal{X})$ obeying in $\mathbb{H}$,
where $g\mapsto H_{\infty}(g)$ is a.s. an element of $\mathbb{H}$ (i.e. it has almost surely uniformly continuous sample paths on $\mathcal{G}$ with respect to the $L^2(\mu)$ metric) and $b_n$ is a sequence such that $b_n \to \infty$ as $n \to \infty$.
When $\widehat \mu$ is the empirical distribution based on a random sample from the population with distribution $\mu$, then $b_n = \sqrt{n}$ and $H_{\infty} = B_{\mu}$, where $B_{\mu}$ is a $\mu$-Brownian Bridge, i.e. a Gaussian process with zero mean and covariance function $(g_1,g_2) \mapsto \mu(g_1 g_2) - \mu(g_1) \mu(g_2)$. In this case condition ${\sf S}.4$ imposes that the function class $$\mathcal{G} = \{x \mapsto 1(f(x) \leqslant \delta) : f \in \mathcal{F}, \delta \in \mathcal{V} \}$$ is $\mu$-Donsker. Note that $\mathcal{F}$ is the parameter space that contains $\Delta(x)$ as well as $\widehat \Delta(x)$ in ${\sf S}.3$. In parametric models for the PE where $ \mathcal{F} = \{ f(x,\theta) : \theta \in \Theta \}$, $f$ is known, $\theta \subseteq {\Bbb{R}}^{d_\theta}$ with $d_\theta < \infty$, and $x \mapsto f(x,\theta)$ is $\mathcal{C}^1$ on $\mathcal{X}$ for all $\theta \in \Theta$, the class $\mathcal{G}$ is $\mu$-Donsker under mild conditions specified for example in \citeasnoun[Chap. 19]{vdV}. Examples 1 and 2 specify the PE parametrically. Lemma (ref) in the SM gives other sufficient conditions for the Donsker property.
The following result is derived as a consequence of the new mathematical results on the Hadamard differentiability of the sorting operator, stated in Lemma (ref) in the Appendix (proof given in SM due to space constraints), in conjunction with the functional delta method. It shows that the empirical SPE-function follows a FCLT over sets of quantiles corresponding to $\Delta_{\mu}^*$ pre-images of compact sets of $\mathbb{R}$.
Define $\mathcal{D}$ as a compact set consisting of regular values of $x \mapsto \Delta(x)$ on $\overline \mathcal{X}$, and $\mathcal{U} := \{\widetilde u \in [0,1] : \Delta^*_{\mu}(\widetilde u) \in \mathcal{D}, f_{\Delta,\mu}(\Delta_{\mu}^*(\widetilde u)) > \varepsilon\},$ for a fixed $\varepsilon>0$, where $ f_{\Delta,\mu}(\Delta_{\mu}^*(\widetilde u))$ is the density of $\Delta(X)$ defined in Lemma (ref)(a). Let $r_n := a_n\wedge b_n$, the slowest of the rates of convergence of $\widehat \Delta$ and $\widehat \mu$. Assume $r_n/a_n \to s_{\Delta} \in [0,1]$ and $r_n/b_n \to s_{\mu} \in [0,1]$, where $s_{\Delta} = 0$ when $b_n = o(a_n)$ and $s_{\mu} = 0$ when $a_n = o(b_n)$. For example, $s_{\mu} = 0$ if $\mu$ is treated as known.
It is convenient to modify the notation for the $u$-CA separating the dependence on $\Delta^*_{\mu}(u)$ from $\Delta$ and $\mu$ and specifying the characteristic of interest as $\varphi_t$. Moreover, when $Z = (X,Y)$ we remove the dependence on $Y$ by taking expectations conditional on $X$. Let $\Lambda_{\widehat \Delta,\widehat \mu,\widehat\Delta^*_{\widehat \mu}(u)}(\varphi_t):= \widehat{\Lambda}_{\Delta,\mu}^{u}(t)$ and $\Lambda_{\Delta,\mu,\Delta^*_{\mu}(u)}(\varphi_t) := \Lambda_{\Delta,\mu}^{u}(t)$, where $\varphi_t \in \mathcal{F}_M \cup \mathcal{F}_I$, $t = (t_1, \ldots, t_{d_z}) \in \mathbb{R}^{d_z}$, $u \in \mathcal{U}$, $\mathcal{F}_M:=\{\int z_{1}^{t_1} \cdots z_{d_z}^{t_{d_z}} d\mu(y \mid x): t_1,...,t_{d_z} \in \{0,1,2, \ldots\}, \int |z_{1}^{t_1} \cdots z_{d_z}^{t_{d_z}}| d\mu(z) < \infty , t_1 + \ldots + t_{d_z} \leqslant M\}$, $M$ is some fixed integer, $\mu(y \mid x)$ is the distribution of $Y$ at $y$ conditional on $X=x$, and $\mathcal{F}_I:=\{\int 1(z_{1}\leqslant t_1,...,z_{d_z}\leqslant t_{d_z}) d\mu(y \mid x): t_1,...,t_{d_z} \in \mathbb{R}\}$. For example, $\varphi_t(x) = x^{t_x} {\mathrm{E}}[Y^{t_y} \mid X = x]$ or $\varphi_t(x) = 1(x \leqslant t_x) \mu(t_y \mid x)$ for $t = (t_x,t_y)$. To derive the properties of $\Lambda_{\widehat \Delta,\widehat \mu,\widehat\Delta^*_{\widehat \mu}(u)}(\varphi_t)$, we use that the class of functions $\widetilde{\mathcal{G}}=\{1(f\leqslant \delta)\varphi : \varphi\in \mathcal{F}_M \cup \mathcal{F}_I,\delta\in \mathcal{V}, f\in \mathcal{F}\}$ is $\mu$-Donsker. When $x \mapsto \mu(y \mid x)$ is continuous, this property holds by assumption ${\sf S}.4$ when $\widehat \mu$ is the empirical distribution.\footnote{Lemma (ref) in the SM gives other sufficient conditions for the Donsker property.}
The following result is derived as a consequence of the new mathematical results on the Hadamard differentiability of the classification operator, stated in Lemma (ref) in the Appendix (proof given in SM due to space constraints), in conjunction with the functional delta method.
Assumption ${\sf AS}.1$ is a technical condition stated in Appendix (ref) of the SM to deal with the discontinuity of the indicator functions when $\varphi_t \in \mathcal{F}_I$. A sufficient condition for ${\sf AS}.1$ is that $$\int_{\mathcal{M}_{\Delta}(\delta) \cap \{x : x_k=t_k\}}d \mathrm{Vol}=0$$ holds uniformly over all $\delta\in \mathcal{V}$, $t_k \in \mathbb{R}$ and $k=1,2,...,d_x$. In other words, the manifold $\mathcal{M}_\Delta(\delta)$ and the set of points $\{x : x_k = t_k\}$ can not have an intersection with positive volume of $(d_x-1)$-dimension.
Corollaries (ref) and (ref) use critical values of statistics related to the limit processes $Z_{\infty}$ and $Z^u_{\infty}$ to construct confidence bands and p-values. These critical values can be hard to obtain in practice. In principle one can use simulation, but it might be difficult to numerically locate and parametrize the manifold $\mathcal{M}_{\Delta}(\delta)$, and to evaluate the integrals on $\mathcal{M}_{\Delta}(\delta)$ needed to compute the realizations of $Z_{\infty}(u)$ and $Z^u_{\infty}(t)$. This creates a real challenge to implement our inference methods. To deal with this challenge we employ (exchangeable) bootstrap to compute critical values praestgaard-wellner-93,vdV-W instead of simulation. We show that the bootstrap law is consistent to approximate the distribution of the limit processes of Theorems (ref) and (ref).
To state the bootstrap validity result formally, we follow the notation and definitions in \citeasnoun{vdV-W}. Let $\mathrm{D}_n$ denote the data vector and let $\mathrm{B}_n = (\omega_1, \dots, \omega_n)$ be the vector of bootstrap weights. Consider a random element $\widetilde{Z}_n = Z_n(\mathrm{D}_n,\mathrm{B}_n)$ in a normed space $\mathbb{D}$. We say that the bootstrap law of $\widetilde{Z}_n$ consistently estimates the law of some tight random element $Z_{\infty}$ and write $\widetilde{Z}_n\rightsquigarrow_{{\mathrm{P}}} Z_{\infty}$ if $$ \sup_{h \in \mathrm{BL}_1(\mathbb{D})} |{\mathrm{E}}_{\mathrm{B}_n} h(\widetilde{Z}_n) - {\mathrm{E}}_{{\mathrm{P}}} h(Z_{\infty})| \to_{{\mathrm{P}}} 0, $$ where $ \mathrm{BL}_1(\mathbb{D})$ denotes the space of functions with Lipschitz norm at most 1; ${\mathrm{E}}_{\mathrm{B}_n}$ denotes the conditional expectation with respect to $\mathrm{B}_n$ given the data $\mathrm{D}_n$; ${\mathrm{E}}_{{\mathrm{P}}}$ denotes the expectation with respect to ${\mathrm{P}}$, the distribution of the data $\mathrm{D}_n$; and $\to_{{\mathrm{P}}}$ denotes convergence in (outer) probability.
The next result is a consequence of the functional delta method for the exchangeable bootstrap. Let $\Lambda_{\widetilde \Delta,\widetilde \mu,\widetilde \Delta^*_{\widetilde \mu}(u)}(\varphi_t):= \widetilde{\Lambda}_{\Delta,\mu}^{u}(t)$, the bootstrap draw of $\widehat{\Lambda}_{\Delta,\mu}^{u}(t)$ defined in Algorithm (ref).
Theorem (ref) employs the high-level condition that the bootstrap can approximate consistently the laws of $\widehat \Delta$ and $\widehat \mu$, after suitable rescaling. In Examples 1-3 when $\widehat \mu$ is the empirical measure based on the random sample of size $n$, the exchangeable bootstrap method entails randomly reweighing the sample using the weights $(\omega_{1}, \ldots, \omega_{n})$, which include empirical boostrap and i.i.d. exponential weights, for example. In this case the high level condition holds if the weights satisfy the conditions stated in equation (3.6.8) of \citeasnoun{vdV-W}. We refer to \citeasnoun{vdV-W} and \citeasnoun{CFM} for bootstrap FCLT for parametric and semi parametric estimators of $\Delta$ including least squares, quantile regression, and distribution regression, as well as nonparametric estimators of $\mu$ including the empirical distribution function.