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.
64,506 characters · 19 sections · 54 citation commands
Counterfactual Density Effects and the German East--West Income Gap
\doublespacing
\noindentKeywords: Density regression, decomposition methods, causal inference.
The most prevalent approaches in causal inference are based on the study of mean-based quantities such as the average treatment effect (ATE) and the average treatment effect on the treated (ATT). While it is true that ATEs are easy to interpret and can be estimated reliably in many situations, they are nonetheless restricted to identifying location-shifts in the data while ignoring effects that go beyond the first moment of the distribution. Recently, there has been growing interest in more nuanced approaches based on quantiles or other distributional characteristics chernozhukov2005iv, firpo2007efficient, chernozhukov2013inference.
In this paper, we advocate for a density-focused approach for causal inference. Similar to quantiles, densities reflect the entire distribution of the variable of interest. However, densities offer some important advantages. The computational burden is reduced, as estimating conditional densities involves a single estimate, whereas distributional and quantile regressions must be run separately for different values of the threshold index and quantile level, respectively. There is also no monotonicity issue in the estimation unlike in quantile-based approaches. Consequently, rearrangement methods to avoid quantile crossing, such as those proposed in chernozhukov2010quantile, are not required. Additionally, it can be argued that densities are better suited to display and intuitively understand the shape of the data compared to distribution and quantile functions, e.g., in the presence of bimodalities or shifts in the probability mass. Finally, methods based on quantile regression usually require continuous distributions, while our density-based approach is also applicable in the case of discrete and mixed-type outcome variables. These problems are highly relevant in the case of income distributions considered in our empirical application, which exhibit both bimodalities and probability masses at zero.
Our definition of counterfactual densities and the resulting density effects build on the decomposition literature, originating with blinder1973wage and oaxaca1973male. Specifically, we define counterfactual densities based on the conditional density of group A, evaluated as if it were exposed to the covariate distribution of another group, B. The original Oaxaca--Blinder decomposition was introduced for the purpose of explaining differences in means of two groups. A special case of the decomposition focusing on proportions and categorical covariates had already been proposed by kitagawa1955components. For approaches beyond the mean, similar decompositions based on counterfactual quantities have been studied for quantiles and other distributional characteristics dinardo1996labor, firpo2018decomposing, chernozhukov2013inference. We refer to fortin2011decomposition for a comprehensive overview of decomposition methods in economics. These approaches are based on an additive Oaxaca--Blinder type decomposition. As a crucial difference, we instead focus on a multiplicative decomposition. We consider two types of counterfactual density effects that admit a causal interpretation under the standard unconfoundedness and overlap assumptions. First, the distribution effect captures the counterfactual impact of changing the conditional density from that of the control group to that of the treatment group, while holding the covariate distribution fixed at the treatment group level. Second, the covariate effect represents the effect of a hypothetical change in the covariate distribution, with the conditional density held fixed at the control group.
The estimation of counterfactual densities depends critically on accurate estimates of conditional densities. However, many existing approaches have important limitations. On the one hand, fully nonparametric estimators suffer from the curse of dimensionality, making them impractical even in moderately dimensional settings. Examples include the log-spline approach of stone1991asymptotics, stone1994use and local polynomial-based estimators fan1996estimation, cattaneo2024boundary. On the other hand, fully parametric models rely on strong distributional assumptions, which can lead to model misspecification when these assumptions are violated. In this paper, we instead rely on the Bayes Hilbert space approach for conditional density estimation of maier2024conditional. The authors propose a flexible structured additive regression model that obeys the logic of densities, i.e., the estimated densities are non-negative and integrate to one. Estimation relies on the use of basis functions (e.g., splines) and is based on a Poisson approximation to the Bayes Hilbert space likelihood.
Our counterfactual density methodology is motivated by the analysis of the East--West income gap in Germany. This empirical application illustrates the advantages of density-based approaches for several reasons. First, some of the estimated densities exhibit heavy skewness and bimodality, features that are easily detectable in density plots but may be difficult to identify using estimated quantile or distribution functions. Second, the income variable is zero-inflated, with the fraction of zeros varying with the covariates, which poses challenges for quantile-based methods. This is particularly relevant since a focus on the positive part of the income distribution would avoid this problem only at the expense of losing information about the unemployed. Empirically, we find that the East--West gap has narrowed over the past 30 years after reunification; however, notable differences still persist. Our results suggest that these differences are largely driven by the conditional distribution rather than by differences in the composition of covariates. I.e., differences in the covariate distributions explain only a small part of the observed differences in the income distributions. Finally, we find that these differences are much more pronounced when focusing on the male subpopulation. We therefore conclude that the East--West income gap in Germany is, to a large extent, a male-specific issue.
Methods for counterfactual distributions based on quantile and distributional regression instead of density regression have been discussed before. chernozhukov2013inference is one of the most closely related papers. They focus on distributional (i.e., for the cumulative distribution function) and quantile effects, but consider a similar Oaxaca--Blinder decomposition of effects. Instead of estimating the conditional density, their approach relies on quantile and distribution function regression. For both approaches, they rely on known basis functions of covariates. Similarly, machado2005counterfactual consider an Oaxaca--Blinder decomposition based on estimating conditional quantiles. In particular, they estimate counterfactual densities based on linear quantile regression fits for several quantile levels and a simulation-based procedure. I.e., they generate counterfactual outcomes by using the estimated quantile regression coefficients and by drawing from the covariate distribution (to integrate-out the effect of the covariates) and by sampling quantile ranks from a uniform distribution. The resulting counterfactual densities need to be estimated by kernel density estimation. We consider our density regression-focused approach to be complementary to these existing approaches. This is particularly the case in settings with mixed-type dependent variables, as in our empirical application. Quantile-based methods are not applicable in such scenarios.
The issue of counterfactual density estimation has been studied in the existing literature. dinardo1996labor considers a similar decomposition of effects and a similar plug-in estimator for counterfactual densities. However, (i) they do not explicitly discuss causal implications of their estimated densities, (ii) they use kernel density estimates, which restricts the applicability in many settings and (iii) they consider differences between densities instead of ratios, which might be hard to interpret in low-density regions of the support. martinez2024counterfactual model counterfactual densities using Kernel--Stein discrepancies, but rely on parametric assumptions. melnychuk2023normalizing propose a deep learning method called `Interventional Normalizing Flows' for estimating fully parametric counterfactual densities. Similarly, kennedy2023semiparametric approximate counterfactual densities by parametric models; for instance, relying on the exponential family, Gaussian mixture models, or on truncated series regression. Theoretically, they provide asymptotic results for the estimator, such as root-$n$ consistency and semiparametric efficiency bounds. They also consider `density effects', which they define based on distances or other measures of discrepancy. However, the above-mentioned approximations suffer from several limitations, namely potential model misspecification and the risk that the estimated densities can take negative values and do not integrate to one. Further, the specific setting of our empirical application with bimodalities and point masses cannot be adequately addressed by fixed parametric distributions. The focus of the present paper is therefore different from that of the above literature for two reasons. First, we explicitly embed our framework for causal inference on counterfactual densities within the decomposition methods literature. Second, while the above papers try to detect discrepancies between two densities in the form of scalar quantities, our method allows the localization of the regions in the support of the dependent variable where the discrepancies are substantial, by relying on ratios instead of distances.
Our contributions are fourfold. First, we propose a causal inference framework based on a multiplicative Oaxaca--Blinder decomposition of counterfactual densities. Compared with conventional additive decompositions, we thus focus on relative differences instead of absolute differences between densities, which has clear advantages for the analysis of low-density regions. A multiplicative approach further obeys the logic of Bayes Hilbert spaces, a suitable functional space for density functions. Second, we propose an estimation procedure that relies on a flexible additive model specification for the conditional densities that avoids both the restrictiveness of parametric models and the curse of dimensionality of fully nonparametric models. Further, the approach does not require continuous outcome variables but also allows for discrete and mixed-type distributions, as is the case in our empirical application. In practice, estimation can be carried out as an approximate Poisson regression problem. Third, the additive specification of the conditional densities allows us to further isolate the effect of different covariates on the covariate effect. Finally, on the empirical side, we use our counterfactual density methodology to gain more detailed insights into the East--West income gap in Germany that would otherwise have been lost when relying on mean-based procedures.
The remainder of this paper is structured as follows. Section (ref) introduces the model setup and notation and provides a definition of counterfactual densities as well as the corresponding counterfactual density effects. The Bayes Hilbert space approach for the estimation of conditional densities is discussed in Section (ref). We provide a short simulation study to analyze the finite-sample properties of our estimation method in Section (ref). In Section (ref) we employ our counterfactual density methodology to study the East--West income gap in Germany. Section (ref) concludes.
In this section, we introduce our framework for counterfactual densities. For this purpose, we will use the following notation. Let $Y\in\mathbb{R}$ be the dependent variable, $X\in\mathbb{R}^d$ the vector of covariates, and let $D\in\{0,1\}$ be the treatment variable, which takes the value one for the treated and zero otherwise. For the sake of exposition, we focus on binary treatment variables. However, extensions to more than two treatment groups will be discussed later. We consider the following potential outcome setting rubin1974estimating,
where $Y_1$ is the potential outcome for the treated and $Y_0$ the potential outcome for the control group. Since only one potential outcome is observed, $Y_0$ is a latent variable for the treatment group and $Y_1$ for the control group.
Let $f_{Y_{\langle1,1\rangle}}$ and $f_{Y_{\langle0,0\rangle}}$ denote the unconditional densities of $Y_1$ and $Y_0$, respectively. Further define the conditional densities for treatment and control group as $f_{Y_1|X}(y|x)$ and $f_{Y_0|X}(y|x)$. In many applications, it may be of interest to define densities for the counterfactual quantities. For this purpose, $f_{Y_{\langle1,0\rangle}}$ denotes the (unconditional) counterfactual density of the treated if they had faced the covariate distribution of the control group,
where $F_{X_0}$ is the marginal cumulative distribution function (cdf) and $\mathcal{X}_0$ is the support of $X$ for the control group. Vice versa we can define $f_{Y_{\langle0,1\rangle}}$, the counterfactual density of the untreated if they had faced the covariate distribution of the treated group. These quantities are the object of study in dinardo1996labor. chernozhukov2013inference introduce a similar relationship for counterfactual distribution functions. Equation ((ref)) reveals that the counterfactual density is completely determined by the conditional density of the treated group and the marginal covariate distribution function of the control group. The crucial part in the estimation of counterfactual densities thus is to devise a suitable estimator of the conditional density. We will discuss this in detail in Section (ref).
To attribute a causal interpretation to the counterfactual densities introduced in the previous subsection, we have to impose the following assumptions.
Both assumptions are commonly used in the causal inference literature (e.g., rosenbaum1983central), and are not particular to our density-focused framework. However, they need to be carefully checked and discussed in any application of our framework.
A central issue is to define causal quantities which (i) are relevant for practitioners and (ii) are based on the counterfactual densities. One possibility is to consider the decomposition by oaxaca1973male originally introduced for the mean, but studied by dinardo1996labor in the context of densities and by chernozhukov2013inference for distribution and quantile effects. The decomposition for densities is given by
This gives three different effects, (i) the effect of changing the conditional density, (ii) the effect of changing the covariate distribution, and (iii) a combination of both.
Recently, kennedy2023semiparametric studied another quantity, which they labeled as `density effects'. Consider two densities, $g_1(y)$ and $g_0(y)$, and some discrepancy function $h:\mathbb{R}^2\to\mathbb{R}^{+}$. Then, their density effects are scalar quantities defined as
Examples for $D_h$ include the total variation distance and the KL divergence. While kennedy2023semiparametric did not explicitly consider a decomposition approach, their measure can easily be incorporated under an Oaxaca-Blinder type decomposition. Potential downsides of this approach are limited interpretability and relevance for practitioners, as well as the loss of information entailed by aggregating discrepancies between densities into a single number. Even though it might be a suitable measure to detect the existence of possible discrepancies, it does not provide any information on the direction of the effect. In general, focusing on a scalar quantity leads to a loss of a large proportion of the information. We argue that when studying counterfactual densities, the focus should lie on the analysis of the heterogeneity of possible effects, i.e., how the treatment affects certain regions of the distribution. E.g., the treatment can have an effect on the lower or upper tail of the distribution or it can affect the entire distribution with a shift in the mean. For example, the introduction of a minimum wage will primarily affect the lower tail of the income distribution.
Going into a similar direction as the decomposition of oaxaca1973male and blinder1973wage, we can consider an alternative decomposition of effects. Instead of looking at differences between densities, we consider ratios,
The motivation for the use of a multiplicative decomposition is twofold. First, for the estimation of conditional densities we rely on the use of Bayes Hilbert spaces, which are suitable spaces for density functions. These are vector spaces in which addition corresponds to multiplication and subtraction corresponds to taking ratios (for details see Section (ref)). Related to this point, the use of ratios offers clear advantages over differences for the analysis of densities. Our proposal can detect discrepancies along the whole domain of the density, whereas differences can only find minor differences in regions with low density values, e.g., tail regions can never show important differences in contrast to high density regions. Similar to chernozhukov2013inference, we can define three kinds of effects.
An advantage of the classical mean-based Oaxaca-Blinder approach is the possibility to further decompose the distribution and covariate effects into individual contributions of covariates. Unfortunately, this advantage does not directly translate to the study of other nonlinear distributional quantities, such as quantile or distributional effects. rothe2015decomposing demonstrates this problem for quantile treatment effects, even in the case of a linear model for the conditional quantiles. This is due to potential dependence among different covariates. The same issue complicates any attempt to decompose our counterfactual density effects additively. As for quantile effects, the problem even persists in the case of a purely multiplicative density regression model as we will present in Section (ref). A possible but not entirely satisfactory solution for this issue is the sequential conditioning approach of chernozhukov2013inference, which suffers from the issue of path-dependence: in most applications the order of covariates is chosen arbitrarily. Further, rothe2012partial argues that such an approach is unable to accurately reflect the impact of group differences in the marginal distribution of a single covariate.
Instead, we propose to isolate the contribution of individual covariates by considering marginal covariate effects, holding all remaining variables fixed at their control group distribution. This approach avoids path-dependence and has a transparent interpretation: each quantity measures the effect on the outcome density of shifting a single covariate's distribution from the control to the treatment group, while leaving the joint distribution of the remaining covariates unchanged. While the resulting quantities do not multiply to the total covariate effect, they nonetheless provide interpretable summaries of which covariates drive the overall covariate effect.
Recall the definition of the covariate effect,
Consider the marginal effect of the $j$-th variable on the density of $Y_k$, after integrating out the effect of all other variables using the marginal covariate distribution of group $l$,
for $k,l=0,1$, where $F_{X_{l,-j}}$ denotes group $l$'s cdf of all covariates except variable $j$. Related quantities have been studied in the context of nonparametric estimation of the conditional mean, see e.g., linton1995kernel and hardle2004nonparametric. Then we can consider the following quantity to measure variable $X_j$'s contribution to the covariate effect, by changing its distribution (instead of that of all $X$) from control to treatment group,
Under the additive model specification for the conditional density that we will introduce in Section (ref), the partial effect $\widetilde{h}_{j,0|0}$ in ((ref)) corresponds to the additive effect of variable $X_j$. The expression is therefore independent of the marginal distribution of $X_{-j}$ in this special case. The quantity $\operatorname{CE}_j$ answers the following question: how would the outcome density change if only the marginal distribution of covariate $X_j$ were shifted from that of the control group to that of the treatment group, while all other covariates remain distributed as in the control group? It provides a direct answer to this question at each point $y$ of the support of the dependent variable, rather than compressing it into a scalar summary as in the mean case.
Complementary to (ref), one can further analyze the distribution effect by looking at the contribution of a given covariate $X_j$ to the distribution effect. We define the contribution to the distribution effect as follows,
The quantity $\operatorname{DE}_j$ compares the partial effect of covariate $X_j$ on the conditional density across the two groups. Specifically, the numerator evaluates the treatment group's structural relationship between $X_j$ and $Y$, while the denominator evaluates the control group's, both integrated over the common (control group) distribution of $X_j$. Thus, $\operatorname{DE}_j(y)$ measures how much the density at $y$ would change if only the way $X_j$ shapes the conditional density were switched from the control to the treatment group specification, while the distribution of $X_j$ itself is held fixed. This isolates the role of $X_j$ in driving the distribution effect, and can be a useful tool for identifying the set of variables that are most responsible for the distribution effect.
In this section, we present our estimation procedure for the counterfactual densities. The critical step in the estimation is to obtain a suitable estimate of the conditional densities. Let $\{(Y_{ki},X_{ki})\}_{i=1}^{n_k}$ denote a sample of $n_k$ i.i.d. copies of $(Y_k,X_k)$, for $k=0,1$. We consider a plug-in estimator of the form
where $\widehat{f}_{Y_1|X}(y|x)$ is an estimator of the conditional density and
is the empirical distribution function of $X_0$, with the inequality $X_{0i}\leq x$ interpreted entry-wise for each $x_j$. $\widehat{f}_{Y_{\langle0,1\rangle}}(y)$ and $\widehat{F}_{X_1}(x)$ can be defined analogously. dinardo1996labor propose a kernel density estimation procedure for the conditional densities. However, such fully nonparametric methods suffer from the curse of dimensionality and do not work well in settings with moderate to large dimensions of covariates. Other approaches suffer from restrictive parametric assumptions or do not provide suitable estimates for densities, i.e., the non-negativity or integrating-to-one constraints might be violated. Instead, for the estimation of the conditional densities we follow the Bayes Hilbert space approach of maier2024conditional, who developed a flexible additive model framework for modeling conditional densities.
Before setting up the density regression model, we first provide a concise introduction to the used Bayes Hilbert spaces boogart2010bayes,boogart2014bayes. For a more detailed introduction we refer the reader to maier2025additive and maier2024conditional. The Bayes Hilbert space on the measurable space $(\mathcal{T},\mathcal{A})$ with reference measure $\mu$ is defined by $B^2(\mu)=B^2(\mathcal{T},\mathcal{A},\mu):=\{f\in B(\mu)|\int_{\mathcal{T}}(\log f)^2d\mu<\infty\}$, where $B(\mu)$ is a Bayes space (i.e., a set of equivalence classes of $\mu$-densities that are $\mu$-a.e. positive and unique) with reference measure $\mu$. This is a vector space with addition, $f_1\oplus f_2=_{\mathcal{B}}f_1f_2$, and scalar multiplication, $\alpha\odot f_1=_{\mathcal{B}}(f_1)^{\alpha}$, for $f_1,f_2\in B^2(\mu)$ and $\alpha\in\mathbb{R}$, where $=_{\mathcal{B}}$ denotes equality up to scale. A crucial concept in this framework is the centered log-ratio (clr) transformation, $\operatorname{clr}(f):=\log(f)-\frac{1}{\mu(\mathcal{T})}\int_{\mathcal{T}}\log(f)d\mu$, which maps an element in $B^2(\mu)$ to a function in $L^2_0(\mu)=L_0^{2}(\mathcal{T},\mathcal{A},\mu):=\{\tilde{f}\in L_0^2(\mu)|\int_{\mathcal{T}}\tilde{f}d\mu=0\}$, a closed subspace of $L^2(\mu)$. Transforming the data is beneficial for practical implementation and computational reasons as it enables the use of tools established for $L^2$-spaces. The clr transform is an isometric isomorphism, and it is bijective with inverse transformation $\operatorname{clr}^{-1}(\tilde{f})=_{\mathcal{B}}\exp(\tilde{f})$. $B^2(\mu)$ is a Hilbert space, with inner product $\langle f_1,f_2\rangle_{B^2(\mu)}=\int_{\mathcal{T}}\operatorname{clr}(f_1)\cdot\operatorname{clr}(f_2)d\mu$. Although Bayes Hilbert spaces are defined for arbitrary measurable spaces, we will focus on $\mathcal{T}\subset\mathbb{R}$. Following maier2024conditional, we distinguish between three cases. For the continuous case, we have $\mathcal{T}=[a,b]$ and $\mu$ is the Lebesgue measure $\lambda$ on $\mathcal{T}$. For the discrete case, $\mathcal{T}=\{t_1,\ldots,t_D\}$ and $\mu$ is a weighted sum $\delta$ of Dirac measures. Finally, we consider the mixed case with $\mathcal{T}=[a,b]\cup\{t_1,\ldots,t_D\}$ and $\mu=\lambda+\delta$. To demonstrate the importance of the latter case, note that in our empirical application we will look at income distributions in Germany, which have mixed-type densities with additional point masses at zero and for incomes above a certain threshold.
We now present the flexible additive density regression setup of maier2024conditional. In the following, we suppress the group-specific index of the data for the sake of exposition. Consider an i.i.d. sample of observations, $(y_i,x_i)\in\mathcal{T}\times\mathcal{X}$, $\mathcal{X}\subseteq\mathbb{R}^d$, $i=1,\ldots,n$. The conditional density of $Y$ given $X=x_i$, denoted by $f_i:=f(Y|X=x_i)$, is assumed to be an element of the Bayes Hilbert space $B^2(\mu)$. As described in the previous subsection, the framework is flexible enough to handle continuous, discrete, as well as mixed distributions and data. We assume the following additive structure,
where the partial effects $h_j$ are also elements of the Bayes Hilbert space, $h_j(x_i)\in B^2(\mu)$. Each partial effect can depend on one, several (for interactions) or no covariates (the intercept) and can be linear or nonlinear in $x_i$. Each effect is assumed to be represented by the following tensor product basis,
where $b_{j,l}:\mathbb{R}^{d}\to\mathbb{R}$, $b_{\mathcal{T},m}\in B^2(\mu)$ are basis functions over the covariates and over $\mathcal{T}$, respectively, and $\theta_{j,l,m}\in\mathbb{R}$ are the corresponding coefficients. For identification, we center smooth main effects around the intercept $\beta_0$ and interactions around corresponding main effects, see maier2024conditional.
Applying the centered log-ratio (clr) transformation to ((ref)) and using the basis function representation in ((ref)) yields
where $\tilde{b}_{\mathcal{T},m}=\operatorname{clr}(b_{\mathcal{T},m})$, $b(x)=(b_{1,1}(x),\ldots,b_{J,d_J}(x))^{\top}\in\mathbb{R}^{\sum_{j=1}^{J}d_j}$, and $\tilde{b}_{\mathcal{T}}=(\tilde{b}_1,\ldots,\tilde{b}_{d_{\mathcal{T}}})^{\top}\in(L_0^2(\mu))^{d_{\mathcal{T}}}$ are vectors of basis functions, and the corresponding parameter vector $\boldsymbol{\theta}=(\boldsymbol{\theta}_1^{\top},\ldots,\boldsymbol{\theta}_J^{\top})^{\top}\in\mathbb{R}^R$ with $\boldsymbol{\theta}_j=(\theta_{j,1,1},\ldots,\theta_{j,d_j,d_{\mathcal{T}}})^{\top}$ and the dimension of $\boldsymbol{\theta}$ is $R=\sum_{j=1}^J d_jd_{\mathcal{T}}$. The choice of the covariate-specific basis functions $b_j(x_i)$ depends on the type of the considered effect. For instance, a smooth non-linear effect can be modeled via B-splines. Similarly, the choice of the basis functions $\tilde{b}_{\mathcal{T}}$ depends on the reference measure $\mu$. maier2025additive describe constructions based on transformations of B-splines and of indicators for the continuous and the discrete part, respectively.
In principle, $\boldsymbol{\theta}$ can be estimated via maximum likelihood estimation, with likelihood and log-likelihood functions given by
To increase the smoothness of the estimated densities, it is also possible to include additional penalty terms for the coefficients, and thus to consider a penalized log-likelihood. The estimation can be computationally challenging due to the presence of the integral term in the log-likelihood. Following maier2024conditional, we instead estimate $\boldsymbol{\theta}$ by approximating the problem via additive Poisson regression.
For reducing the computational burden, the estimation procedure can be approximated by using a (shifted) multinomial log-likelihood. In this section we focus on the continuous case for notational simplicity, although maier2024conditional show that the discrete and mixed cases can also be covered. For this purpose, we need to partition the support of $Y$ into discrete histogram bins. Let $a=a_0<a_1<\ldots<a_G=b$. Then $U_g=[a_{g-1},a_g)$ for $g=1,\ldots,G-1$ and $U_G=[a_{G-1},a_G]$ partition the interval $[a,b]$. The values of the histogram are $n_g^i=\mathbf{1}_{\{y_i\in U_g\}}$ and the corresponding histogram widths are $\Delta_g=a_g-a_{g-1}$ for $g=1,\ldots,G$. Further, denote the bin center of histogram bin $U_g$ as $u_g$. The vector $(n_1^i,\ldots,n_G^i)$ can be viewed as a realization of a multinomial variable with sample size $1$ and with class probabilities
The multinomial log-likelihood up to constants is
maier2024conditional show that the multinomial log-likelihood converges to the Bayes Hilbert space log-likelihood as the maximal bin size approaches zero, as well as the convergence of the corresponding maximum likelihood estimator and the inverse Fisher information used for inference. Further, they show the equivalence of the multinomial and a certain Poisson likelihood. In particular, for computational reasons it is beneficial to pool observations that share the same combination of covariates, and fit the histogram counts using a Poisson model with an additional intercept parameter for each unique covariate combination. We can thus rely on additive Poisson regression for estimation of the parameter vector $\boldsymbol{\theta}$.
Having introduced the additive density regression framework and the associated estimation procedure, a natural next question concerns the issue of uncertainty quantification. In particular, it is of major interest to empirical researchers whether the distribution and covariate effects are significant, i.e., different from one. For this purpose, we rely on the asymptotic results of maier2024conditional. We propose drawing values $\boldsymbol{\theta}_k^{(b)}$, $b=1,\ldots,B$, of the regression parameters from the $(1-\alpha)$ Wald confidence regions defined in Lemma A.15 of the above paper, for a given significance level $\alpha$. For each simulation iteration, we obtain the corresponding conditional density based on the simulated parameter values and the fixed basis functions, $\widehat{f}_{Y_k|X}^{(b)}(y|x)=(b(x)\otimes \tilde{b}_{\mathcal{T}})^{\top}\boldsymbol{\theta^{(b)}_k}$, for $k=0,1$. Further, we calculate the respective counterfactual densities by integrating with respect to the covariate distributions of the treatment and control groups, $\widehat{f}^{(b)}_{Y_{\langle k,l\rangle}}(y)=\int_{\mathcal{X}_l} \widehat{f}^{(b)}_{Y_k|X}(y|x)d\widehat{F}_{X_l}(x)$, for $k,l=0,1$. Since the conditional densities of both groups are estimated independently, we can follow the above procedure separately for each group. We thus obtain estimates of the respective density effects based on simulations from the asymptotic confidence regions of the model parameters. These estimates can be plotted alongside the original estimate of the respective density effect and thus serve to quantify estimation uncertainty.
In this section, we study the finite sample performance of our proposed estimator for the counterfactual densities using a simulation study. To enable the comparison with an alternative estimator for the conditional densities, we restrict our attention on a data-generating process with categorical covariates and continuous dependent variables. We consider a setting with $d=3$ covariates, each taking two possible values. We assume the additive model specification for the conditional densities introduced in (ref). The partial effects $h_j$ are generated as beta density functions with different parameter values, $h_j(x)=x_j\odot \beta_j$, $j=1,2,3$. As a consequence, the conditional densities are also beta densities. We consider sample sizes ranging from $n=500$ to $n=100{,}000$. Further details on the data-generating process and other aspects of the simulation study are provided in Section (ref) of the supplementary material. It should be noted that our model is not correctly specified, as the spline basis functions only approximate the true conditional densities.
The performance is evaluated by the total variation (TV) distance between the true and estimated counterfactual densities. We further evaluate the estimation accuracy of the conditional densities using the same metric. To provide a benchmark, we compare the estimation accuracy with an alternative approach based on kernel density estimation with a Gaussian kernel and Silverman's rule of thumb bandwidth, carried out separately for every covariate combination, which is possible in this case with only binary covariates. All simulation results are based on $1{,}000$ Monte Carlo iterations.
Table (ref) reports the estimation accuracy of the four counterfactual densities in terms of the TV distance between estimated and true densities. As expected, the estimation becomes more accurate with increasing sample size. The table further reports the estimation accuracy of the benchmark method based on kernel density estimates of the conditional densities. The performance of both methods is quite similar in small samples, whereas the Bayes Hilbert space approach has a slight advantage in settings with larger samples ($n=10{,000},20{,}000$). However, for the largest sample size ($n=100{,}000$), the performance of the two methods appears to converge.
Since this is the key step of our estimation methodology, we separately analyze the estimation accuracy of the conditional densities using the TV distance between the true and estimated densities, comparing the Bayes–Hilbert approach with kernel density estimation. To make the estimation results easier to interpret, we aggregate the results by averaging the TV distances over all covariate combinations. Interestingly, the results in Table (ref) tell a different story from the estimation results of the unconditional counterfactual densities. In fact, the estimation accuracy of the Bayes Hilbert space approach is higher in all the settings we consider. The reason for this is that the method makes explicit use of the multiplicative structure of the conditional densities. To illustrate this point, for the estimation of one particular conditional density, $f_{Y_j|X}(y|X=x)$, the Bayes Hilbert space approach can borrow strength from data that do not belong to this conditional density. I.e., an observation $i$ with $X_i\neq x$ can still impact the estimation of the conditional density. In contrast, the kernel density estimator can only make use of observations for which $X_i=x$. However, since the estimation of the counterfactual densities involves taking averages over the estimated conditional densities, the kernel density estimator has the advantage that the estimates of the different conditional densities are statistically independent. Therefore, the results for the counterfactual densities in Table (ref) are less clear than the results in Table (ref).
We want to point out that for comparison purposes the simulation study is restricted to discrete regressors and a continuous outcome variable. An additional advantage of the Bayes Hilbert space approach is that it can be easily applied to settings with continuous regressors as well as discrete and mixed-type outcome variables. To account for continuous regressors, the kernel density benchmark would also need to rely on smoothing in the covariate dimension, which would make the estimation problem of the conditional densities much more difficult. In contrast, the Bayes Hilbert space approach can be readily applied in these settings as well.
We apply our counterfactual density framework to analyze the East--West gap in gross incomes in Germany. Most existing studies on this subject consider a location-based definition of Easterners and Westerners, irrespective of the place of birth and socialization burda1997getting,kluge2018decomposing. However, this approach is likely to suffer from endogeneity bias. Place of residence and place of work are to a large extent choice variables that might depend on latent factors, thereby distorting the effects of interest. In contrast to these above studies, dickey2021persistent are interested both in a location-based and an origin-based East--West gap. These existing papers use Oaxaca--Blinder type decomposition methods, which are based either on mean regression, quantile regression, or the unconditional quantile regression approach of fortin2011decomposition. We use our decomposition approach for densities, which allows to look at effects on the whole distribution in an interpretable way while also being able to cover zero as well as positive incomes.
Due to the aforementioned endogeneity issues, we restrict our analysis to an origin-based definition of the East--West gap. The data are obtained from the Socio-Economic Panel (SOEP), which provides person-specific information on demographic and socio-economic aspects (see goebel2019german). We identify Easterners using the variable `loc1989', which provides information on the place of residence immediately before the fall of the Berlin Wall and German reunification. If this information is not available, we further classify a person as an Easterner if their birth region lies geographically in the East. This secondary identification is essential for categorizing individuals born after 1989. Formally, this is done using the variable `birthregion_ew'. We restrict our analysis to individuals aged 18–67. We further exclude retirees, students (both at school and at university), as well as persons with disabilities. To account for over- and under-representation of certain demographic groups, we use the corresponding cross-sectional individual weights provided in the dataset.
The dependent variable is monthly gross labor income, which includes both primary and secondary income sources. To ensure comparability across time, all income values are adjusted to 2021 price levels. We note that the density of this variable is of a mixed type, with a point mass at zero. Additionally, we use an upper bound for the support of the income variable at EUR $10{,}000$, setting values above to the category $10{,}000+$ represented by a second point mass. This introduces another point mass at the threshold value. As covariates, we consider sex ($x_{sex}$); a categorical variable for level of education ($x_{edu}$); an indicator variable for whether the place of residence is urban or rural ($x_{rural}$); and a categorical variable for employer size ($x_{size}$). Apart from these categorical explanatory variables, we also consider the age of the individual as a smooth effect. By the logic of the Oaxaca--Blinder approach, we run two separate regressions for Easterners and Westerners, with the following additive density regression model
where $k=1$ for the East and $k=0$ for the West. The partial effects $\beta_{j,k}$ denote group-specific intercepts for the $j$-th categorical variables, where $\beta_{j,k}=0_{B^2}$, the additive neutral element of the Bayes Hilbert space for the respective reference category, and $g_k(\cdot)$ represents a smooth effect for age. Since East Germans are our treatment group and West Germans the control group, the distribution and covariate effects are defined as $\operatorname{DE}(y)=f_{Y_{\langle\text{East},\text{East}\rangle}}(y)/f_{Y_{\langle\text{West},\text{East}\rangle}}(y)$ and $\operatorname{CE}(y)=f_{Y_{\langle\text{West},\text{East}\rangle}}(y)/f_{Y_{\langle\text{West},\text{West}\rangle}}(y)$, respectively. To account for changes over time, we conduct separate analyses for the years 1991 ($n_0=3{,}132$, $n_1=6{,}764$), 2001 ($n_0=3{,}936$, $n_1=9{,}929$), 2011 ($n_0=4{,}893$, $n_1=12{,781}$) and 2021 ($n_0=2{,}447$, $n_1=6{,}433$). For the years 1991 and 2001, we choose smaller upper bounds for the continuous part of the income distribution: EUR $5{,}000$ and EUR $7{,}000$, respectively.
This section presents the estimation results for the counterfactual densities, the corresponding distribution and covariate effects, and their development over time. The results for the East--West gap for 2001 and 2021 are visualized in Figures (ref) and (ref). We refer to Figures (ref) and (ref) in Section (ref) of the supplementary material for additional results for 1991 and 2011. We add 100 draws from the 95% confidence regions for the distribution and covariate effects to quantify the estimation uncertainty, following the procedure outlined in Subsection (ref). First, we observe that the discrepancies between income densities for East and West have decreased substantially over recent decades. This confirms both previous empirical findings and theoretical predictions about East--West wage convergence. Still, the distribution effect for 2021 shows that differences persist in the lower as well as upper regions of the income distribution. For instance, even if West Germans were to have the same age structure and demographic characteristics as East Germans, they would be far more likely to be in the upper tail of the income distribution. As a second major empirical result, we find that the density-based East--West gap can be overwhelmingly attributed to the distribution effect, with the covariate effect only playing a minor role. Thus, the observed differences can not be well explained by factors such as education or urbanization.
To illustrate the benefit of our density-focused analysis, we include a comparison with a classical (additive) Oaxaca--Blinder decomposition for the mean incomes of the treatment (East) and control (West) groups. When looking at the left side of Table (ref), it becomes clear that a mean-based analysis provides a simple, scalar summary that is straightforward to interpret. However, this interpretability comes at the expense of not being able to capture the nuances of the differences. For example, the higher average income of Westerners can stem either from a large share of high incomes or from a very low share of low incomes. Empirically, the results reaffirm our finding that the covariate effect is dominated by the distribution effect, and that the latter effect is decreasing over time. Interestingly, the covariate effect is positive in the years following the reunification, but the sign switches between the years 2001 and 2011.
We also compare our results with the `density effects' proposed by kennedy2023semiparametric, which are again a scalar measure for the discrepancy between two (counterfactual) densities. As a metric we choose the total variation distance, and we also decompose the effect additively into a distribution and covariate effect. For simplicity, we use the same estimated counterfactual densities as in our main analysis. The results on the right hand side of Table (ref) show the fundamental limitation of this approach. Even if the approach can reliably detect discrepancies between the estimated densities beyond the case of simple mean shifts, it is still incapable of identifying the regions of the income distributions that are responsible for the discrepancy. An additional disadvantage is the lack of information about the direction of the effect across the distribution. This information is of course of utmost importance for policymakers. We therefore argue that scalar measures for discrepancies between densities can be beneficial as an additional tool for analysis, but important information about the direction and location of the effect is lost in the process.
In the previous analysis, we assumed that the sex variable enters additively in the two income regressions. However, the decomposition of the East--West income gap may be fundamentally different for men and women. We therefore present additional estimation results for counterfactual densities and density effects based on separate estimations for men and women for 2001. Figures (ref) and (ref) indeed show that the East--West gap is a more pronounced issue for the male population. This can partly be explained by the larger share of part-time work for women in the West compared to the East. Indeed, one can see that the distribution effect is below one in the lower regions of the income distribution. I.e., relatively more women have a low income in West Germany than in East Germany. In contrast, we see similar but less strong effects for women in the upper tail regions of the income distribution.
In the following, we want to further analyze the covariate effect by isolating the impact of single variables using equation ((ref)). We refer to Figure (ref) in the supplementary material for the contributions to the covariate effect in 2001 of the variables education, $\operatorname{CE}_{\text{edu}}(y)$, and age, $\operatorname{CE}_{\text{age}}(y)$. For education, the contribution to the covariate effect is significantly below one in the lower income regions, which implies that the share of East Germans in these regions is higher despite and not because of differences in education. For age, we do not observe any significant effects in any direction.
Finally, in Section (ref) of the supplement, we include a robustness check by including an additional categorical variable controlling for the industrial sector of the main job. The inclusion leads to an endogeneity issue, since having an industry code presupposes the employment of the person, which is why we do not include it in our main empirical analysis. Due to this issue, we need to restrict the analysis to the continuous part of the income distribution. We show that the results are indeed robust towards controlling for industry, i.e., the inclusion of the variable cannot explain the remaining discrepancies in the income distributions between East and West.
In this paper, we presented a new framework for conducting causal inference based on counterfactual densities. The approach is based on a multiplicative Oaxaca--Blinder type decomposition of the densities of the treatment and control groups into a distribution and covariate effect. To estimate the conditional densities, we rely on the Bayes Hilbert space additive density regression model of maier2024conditional. As an application of our approach, we analyze the German East--West income gap. We find that differences between income densities decreased over time and that the decline can mainly be attributed to changes in the conditional distribution. In contrast, differences in the covariate distribution between East and West only play a minor role. Additionally, we find that the East--West gap is much more pronounced for the male sub-population.
A major advantage of our density-based approach is that it is able to capture high degrees of heterogeneity in the causal effects. Further, visualization of the counterfactual densities and corresponding effects allows for an intuitive interpretation of the modeled effects. Compared to alternative approaches based on quantiles, our approach allows for mixed types of the dependent variable. There are some limitations to our approach, which can be the subject of future research. First, the computational complexity of the Poisson estimation problem can be substantial in the case of large sample sizes and continuous explanatory variables. Second, it is currently assumed that the underlying additive density regression model is correctly specified. Future research could analyze the effect of misspecification on the estimation results.
The authors declare no conflict of interests.
Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Project number 513634041.
{ \onehalfspacing }