## Motivation

In regression analyses, we often wonder about "how much covariates matter?" for explaining the relationship between a target variable $D$ and an outcome variable $Y$.

For example, we might start analysing the gender wage gap with a simple regression model as `log(wage) on gender`. But arguably, men and women differ in many socio-economic characteristics: they might have different (average) levels of education or career experience, and they might work in different industries and select into different higher- or lower-paying industries. So which fraction of the gender wage gap can be explained by these observable characteristics?

In this notebook, we will compute and decompose the gender wage gab based on a subset of the PSID data set using a method commonly known as the
"Gelbach Decomposition" ([Gelbach, JoLE 2016](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=1425737)).

We start with loading a subset of the PSID data provided by the AER R package.

```{python}
import re

import pandas as pd

import pyfixest as pf

psid = pd.read_csv(
    "https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/refs/heads/master/csv/AER/PSID7682.csv"
)
psid["experience"] = pd.Categorical(psid["experience"])
psid["year"] = pd.Categorical(psid["year"])
psid.head()
```

Computing a first correlation between gender and wage, we find that males earn on average 0.474 log points
more than women.

```{python}
fit_base = pf.feols("log(wage) ~ gender", data=psid, vcov="hetero")
fit_base.summary()
```

To examine the impact of observable on the relationship between wage and gender, a common strategy in applied research is to incrementally add a set of covariates to the baseline regression model of `log(wage) on gender`. Here, we will incrementally add the following covariates:

- education,
- experience
- occupation,
- industry
- year
- ethnicity

We can do so by using **multiple estimation syntax**:

```{python}
fit_stepwise1 = pf.feols(
    "log(wage) ~ gender + csw0(education, experience, occupation, industry, year, ethnicity)",
    data=psid,
)
pf.etable(fit_stepwise1)
```

Because the table is so long that it's hard to see anything, we restrict it to display only a few variables:

```{python}
pf.etable(fit_stepwise1, keep=["gender", "ethnicity", "education"])
```

We see that the coefficient on gender is roughly the same in all models. Tentatively, we might already conclude that the observable characteristics in the data do not explain a large part of the gender wage gap.

But how much do differences in education matter? We have computed 6 additional models that contain education as a covariate. The obtained point estimates
vary between $0.059$ and $0.075$. Which of these numbers should we report?

Additionally, note that while we have only computed 6 additional models with covariates, the number of possible models is much larger.
If I did the math correctly, simply by additively and incrementally adding covariates, we could have computed $57$ different models (not all of which would have included `education` as a control).

As it turns out, different models **will lead to different point estimates**. The order of incrementally adding covariates **might** impact our conclusion. To illustrate this, we keep the same ordering as before, but start with `ethnicity` as our first variable:

```{python}
fit_stepwise2 = pf.feols(
    "log(wage) ~ gender + csw0(ethnicity, education, experience, occupation, industry, year)",
    data=psid,
)
pf.etable(fit_stepwise2, keep=["gender", "ethnicity", "education"])
```

We obtain 5 new coefficients on `education` that vary between 0.074 and 0.059.

So, which share of the "raw" gender wage gap can be attributed to differences in education between men and women? Should we report a statisics based on the 0.075 estimate? Or
on the 0.059 estimate? Which value should we pick?

To help us with this problem, Gelbach (2016, JoLE) develops a decomposition procedure building on the omitted variable bias formula that produces a single value for the contribution of a given covariate, say education, to the gender wage gap.

## Notation and Gelbach's Algorithm

Before we dive into a code example, let us first introduce the notation and Gelbach's algorithm. We are interested in "decomposing" the effect of
a variable $X_{1} \in \mathbb{R}$ on an outcome $Y \in \mathbb{R}$ into a part explained by covariates $X_{2} \in \mathbb{R}^{k_{2}}$ and an unexplained part.

Thus we can specify two regression models:

- The **short** model
    $$
        Y = X_{1} \beta_{1} + u_{1}
    $$

- the **long** (or full) model

    $$
        Y = X_{1} \beta_{1} + X_{2} \beta_{2} + e
    $$

By fitting the **short** regression, we obtain an estimate $\hat{\beta}_{1}$, which we will denote as the **direct effect**, and by estimating the **long** regression, we obtain an estimate of the regression coefficients $\hat{\beta}_{2} \in \mathbb{R}^{k_2}$. We will denote the estimate on $X_1$ in the long regression as the **full** effect.

We can then compute the contribution of an individual covariate $\hat{\delta}_{k}$ via the following algorithm:

- Step 1: we compute coefficients from $k_{2}$ auxiliary regression models $\hat{\Gamma}$ as
    $$
        \hat{\Gamma} = (X_{1}'X_{1})^{-1} X_{1}'X_{2}
    $$

    In words, we regress the target variable $X_{1}$ on each covariate in $X_{2}$. In practice, we can easily do this in one line of code via `scipy.linalg.lstsq()`.

- Step 2: We can compute the total effect **explained** by the covariates, which we denote by $\delta$, as

    $$
        \hat{\delta} = \sum_{k=1}^{k_2} \hat{\Gamma}_{k} \hat{\beta}_{2,k}
    $$

    where $\hat{\Gamma}_{k}$ are the coefficients from an auxiliary regression $X_1$ on covariate $X_{2,k}$ and $\hat{\beta}_{2,k}$ is the associated estimate on $X_{2,k}$ from the **full** model.

    The individual **contribution of covariate $k$** is then defined as

    $$
        \hat{\delta}_{k} = \hat{\Gamma}_{k} \hat{\beta}_{2,k}.
    $$


After having obtained $\delta_{k}$ for each auxiliary variable $k$, we can easily aggregate multiple variables into a single groups of interest. For example, if $X_{2}$ contains a set of dummies from industry fixed effects, we could compute the explained part of "industry" by summing over all the dummies:

$$
        \hat{\delta}_{\textit{industry}} = \sum_{k \in \textit{industry dummies}} \hat{\Gamma}_{k} \hat{\beta}_{2,k}
$$

## `PyFixest` Example

To employ Gelbach's decomposition in `pyfixest`, we start with the **full** regression model that contains **all variables of interest**:

```{python}
fit_full = pf.feols(
    "log(wage) ~ gender + ethnicity + education + experience + occupation + industry +year",
    data=psid,
)
```

After fitting the **full model**, we can run the decomposition procedure by calling the `decompose()` method. The only required argument is to specify the focal variable `decomp_var`, which in this case is "gender". Inference is conducted via a non-parametric bootstrap and can optionally be turned off.

```{python}
gb = fit_full.decompose(decomp_var="gender[T.male]", digits=5)
```

As before, this produces a pretty big output table that reports
- the **direct effect** of the regression of `log(wage) ~ gender`
- the **full effect** of gender on log wage using the **full regression** with all control variables
- the **explained effect** as the difference between the full and direct effect
- a **single scalar value** for the individual contributions of a covariate to overall **explained effect**

For our example at hand, the additional covariates only explain a tiny fraction of the differences in log wages between men and women - 0.064 points.
Of these, around one third can be attributed to ethnicity, 0.00064 to years of eduaction, etc.

Note that for now, we ask `etable()` to only report effects in levels. By switching to `panels = "all"`, we would also report normalized coefficient; but we decided not to do so here as otherwise the table would have turned out even longer than it already has.

```{python}
gb.etable(
    panels="levels",
)
```

Note: we can also return the decomposition results as a pd.DataFrame via the `tidy` method.

```{python}
gb.tidy().head()
```

Because experience is a categorical variable, the table gets pretty unhandy: we produce one estimate for "each" level. Luckily, Gelbach's decomposition
allows us to group individual contributions into a single number. In the `decompose()` method, we can combine variables via the `combine_covariates` argument:

```{python}
gb2 = fit_full.decompose(
    decomp_var="gender[T.male]",
    combine_covariates={
        "experience": re.compile("experience"),
        "occupation": re.compile("occupation"),
        "industry": re.compile("industry"),
        "year": re.compile("year"),
        "ethnicity": re.compile("ethnicity"),
    },
    only_coef=True,  # suppress bootstrap for inference
)
```

We now report a single value for "experience", which explains a good chunk - around half - of the explained part of the gender wage gap.

```{python}
gb2.etable(
    panels="levels",
)
```

We can aggregate even more to "individual level" and "job" level variables:

```{python}
gb3 = fit_full.decompose(
    decomp_var="gender[T.male]",
    combine_covariates={
        "job": re.compile(r".*(occupation|industry).*"),
        "personal": re.compile(r".*(education|experience|ethnicity).*"),
        "year": re.compile("year"),
    },
    only_coef=True,  # suppress inference
)
```

```{python}
gb3.etable(panels="levels")
```

We can easily change multiple aspects of the GT table that `etable` returns.
If we set `panels = "all"`, `etable()` will also report normalized coefficients.

```{python}
gb3.etable(
    column_heads=["Column A", "Column B", "Column C"],
    panel_heads=["Panel A", "Panel B", "Panel C"],
    panels="all",
    add_notes="We add more notes.",
    caption="Gelbach Decomposition",
)
```

We can visualise the Gelbach decomposition using the `coefplot` method.

```{python}
gb3.coefplot()
```

## Literature

- ["When do Covariates Matter? And Which Ones, and How Much?" by Gelbach, Jonah B. (2016), Journal of Labor Economics](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=1425737)