# Inference

## An inferential problem: The Gender Wage Gap

In the previous lab, we already analyzed data from the March Supplement of the U.S. Current Population Survey (2015) and answered the question how to use job-relevant characteristics, such as education and experience, to best predict wages. Now, we focus on the following inference question:


**What is the difference in predicted wages between men and women with the same job-relevant characteristics?**

Thus, we analyze if there is a difference in the payment of men and women (gender wage gap). The gender wage gap may partly reflect discrimination against women in the labor market or may partly reflect a selection effect, namely that women are relatively more likely to take on occupations that pay somewhat less (for example, school teaching).

To investigate the gender wage gap, we consider the following log-linear regression model

\begin{align} \log(Y) &= \beta'X + \epsilon\\ &= \beta_1 D + \beta_2' W + \epsilon, \end{align}
where $D$ is the indicator of being female ($1$ if female and $0$ otherwise) and the $W$'s are controls explaining variation in wages. Considering transformed wages by the logarithm, we are analyzing the relative difference in the payment of men and women.

---


```python
import numpy as np
import pandas as pd
import statsmodels.api as sm
import sklearn.linear_model as lm
import statsmodels.formula.api as smf
from sklearn.preprocessing import PolynomialFeatures
```

## Data Analysis

We consider the same subsample of the U.S. Current Population Survey (2015) as in the previous lab. Let us load the data set.


```python
file_dir = '/Users/jaumevives/Desktop/MIT 2/14.388/wage2015_subsample_inference.csv'
df = pd.read_csv(file_dir)
```

To start our (causal) analysis, we compare the sample means given gender:



```python
df.head()
```


```python
table = pd.DataFrame()

cols = ["lwage","sex","shs","hsg","scl","clg","ad","ne","mw","so","we","exp1"]

table['Variable'] = ["Log Wage","Sex","Less then High School","High School Graduate","Some College",
                      "Gollage Graduate","Advanced Degree", "Northeast","Midwest","South","West","Experience"]

table['All'] = df[cols].mean().values

table['Men'] = df[df['sex']==0][cols].mean().values
table['Women'] = df[df['sex']==1][cols].mean().values


# Show results
table
```


```python
df.describe()
```


```python
# print to Latex
table.to_latex()
```

In particular, the table above shows that the difference in average logwage between men and women is equal to  0,038.


```python
# compute difference
table.loc[0, 'Women'] - table.loc[0, 'Men'] 
```

Thus, the unconditional gender wage gap is about  $3,8\%$ for the group of never married workers (women get paid less on average in our sample). We also observe that never married working women are relatively more educated than working men and have lower working experience.

This unconditional (predictive) effect of gender equals the coefficient  $\beta$  in the univariate ols regression of  $Y$  on  $D$ :

$$
log(Y)=\beta D+\epsilon.
$$

We verify this by running an ols regression.


```python
nocontrol_fit = smf.ols("lwage ~ sex", data=df).fit()
nocontrol_est = nocontrol_fit.params['sex']
nocontrol_se = nocontrol_fit.HC0_se['sex']

print("The estimated gender coefficient is "\
      + str(nocontrol_est) + " and the corresponding robust standard error is " + str(nocontrol_se))
```

Note that the standard error is computed to be robust to heteroskedasticity.

Next, we run an ols regression of  $Y$  on  $(D,W)$  to control for the effect of covariates summarized in  $W$ :

$$
log(Y)=\beta_1 D + \beta_2' W + \epsilon.
$$
 
Here, we are considering the flexible model from the previous lab. Hence,  $W$ controls for experience, education, region, and occupation and industry indicators plus transformations and two-way interactions.

Let us run the ols regression with controls.


```python
# Ols regression with controls

flex = "lwage ~ sex + (exp1+exp2+exp3+exp4)*(shs+hsg+scl+clg+C(occ2)+C(ind2)+mw+so+we)"

#   Note that ()*() operation in formula objects in statsmodels creates a formula of the sort:
#  (exp1+exp2+exp3+exp4)+ (shs+hsg+scl+clg+occ2+ind2+mw+so+we) +  (exp1+exp2+exp3+exp4)*(shs+hsg+scl+clg+occ2+ind2+mw+so+we)
#  This is not intuitive at all, but that's what it does.

control_fit = smf.ols(flex, data=df).fit()
control_est = control_fit.params['sex']
control_se = control_fit.HC0_se['sex']

control_fit.summary()
```


```python
control_est
```

The estimated regression coefficient  $\beta_1≈-0.06955$  measures how our linear prediction of wage changes if we set the gender variable  $D$  from 0 to 1, holding the controls  $W$  fixed. We can call this the predictive effect (PE), as it measures the impact of a variable on the prediction we make. Overall, we see that the unconditional wage gap of size  4 \% for women increases to about  7 \% after controlling for worker characteristics.

Next, we are using the Frisch-Waugh-Lovell theorem from the lecture partialling-out the linear effect of the controls via ols.


```python
# Partialling-Out using ols

# models
flex_y = "lwage ~  (exp1+exp2+exp3+exp4)*(shs+hsg+scl+clg+C(occ2)+C(ind2)+mw+so+we)" # model for Y
flex_d = "sex ~ (exp1+exp2+exp3+exp4)*(shs+hsg+scl+clg+C(occ2)+C(ind2)+mw+so+we)"   # model for D

# partialling-out the linear effect of W from Y
t_Y = smf.ols(flex_y, data=df).fit().resid
# partialling-out the linear effect of W from D
t_D = smf.ols(flex_d, data=df).fit().resid

# regression of Y on D after partialling-out the effect of W
partial_fit = sm.OLS(t_Y, t_D).fit()
partial_est = partial_fit.params['x1']

print("Coefficient for D via partialling-out " + str(partial_est))

# standard error
partial_se = partial_fit.HC0_se['x1']

# confidence interval
print("2.5 - 97.5 % CI: " + str(partial_fit.conf_int().values[0]))
```

Again, the estimated coefficient measures the linear predictive effect (PE) of  $D$  on  $Y$  after taking out the linear effect of  $W$  on both of these variables. This coefficient equals the estimated coefficient from the ols regression with controls.

We know that the partialling-out approach works well when the dimension of $W$ is low in relation to the sample size  $n$ . When the dimension of  $W$  is relatively high, we need to use variable selection or penalization for regularization purposes.


In the following, we illustrate the partialling-out approach using lasso instead of ols.


```python
# Partialling-Out using lasso

# models
X_d =  smf.ols(flex_d, data=df).data.exog
y_d = smf.ols(flex_d, data=df).data.endog

X_y =  smf.ols(flex_y, data=df).data.exog
y_y = smf.ols(flex_y, data=df).data.endog

# partialling-out the linear effect of W from Y
t_Y_lasso = y_y - lm.LassoCV(cv=10, fit_intercept=False,
                      alphas=np.logspace(-6,2, 100)).fit(X_y,y_y).predict(X_y)

# partialling-out the linear effect of W from D
t_D_lasso = y_d - lm.LassoCV(cv=10, fit_intercept=False,
                      alphas=np.logspace(-6,2, 100)).fit(X_d,y_d).predict(X_d)

# regression of Y on D after partialling-out the effect of W
partial_lasso_fit = sm.OLS(t_Y_lasso, t_D_lasso).fit()
partial_lasso_est = partial_lasso_fit.params[0]

print("Coefficient for D via partialling-out using lasso " + str(partial_lasso_est))

# standard error
partial_lasso_se = partial_lasso_fit.HC0_se[0]

# confidence interval
print("2.5 - 97.5 % CI: " + str(partial_lasso_fit.conf_int()[0]))
```

Using lasso for partialling-out here provides similar results as using ols.


```python
table2 = pd.DataFrame()

table2['Model'] = ["Without controls", "full reg", "partial reg", "partial reg via lasso"]

table2['Estimate'] = [nocontrol_est, control_est, partial_est, partial_lasso_est]

table2['Std. Error'] = [nocontrol_se, control_se, partial_se, partial_lasso_se]

# Show results
table2
```


```python
# print to Latex
table2.to_latex()
```

It it worth to notice that controlling for worker characteristics increases the gender wage gap from less that 4\% to 7\%. The controls we used in our analysis include 5 educational attainment indicators (less than high school graduates, high school graduates, some college, college graduate, and advanced degree), 4 region indicators (midwest, south, west, and northeast); a quartic term (first, second, third, and fourth power) in experience and 22 occupation and 23 industry indicators.

Keep in mind that the predictive effect (PE) does not only measures discrimination (causal effect of being female), it also may reflect selection effects of unobserved differences in covariates between men and women in our sample.

Next we try "extra" flexible model, where we take interactions of all controls, giving us about 1000 controls.


```python
extraflex = "lwage ~ sex + (exp1+exp2+exp3+exp4+shs+hsg+scl+clg+C(occ2)+C(ind2)+mw+so+we)**2"

covs = ["exp1", "exp2", "exp3", "exp4", "shs", "hsg", "scl", "clg", "occ2", "ind2", "mw", "so", "we"]

# we use sklearn polynomial features
poly = PolynomialFeatures(interaction_only=True, degree=4)
#X_poly = poly.fit_transform(df[covs].values)
#X = np.concatenate([df[['sex']].values, X_poly], axis=1)

X = smf.ols(extraflex, data=df).data.exog
y = smf.ols(extraflex, data=df).data.endog
```


```python
# model
extraflex = "lwage ~ sex + (exp1+exp2+exp3+exp4+shs+hsg+scl+clg+C(occ2)+C(ind2)+mw+so+we)**2"

# covs = ["exp1", "exp2", "exp3", "exp4", "shs", "hsg", "scl", "clg", "occ2", "ind2", "mw", "so", "we"]
# alternatively e use sklearn polynomial features
# poly = PolynomialFeatures(interaction_only=True, degree=4)
# X_poly = poly.fit_transform(df[covs].values)
# X = np.concatenate([df[['sex']].values, X_poly], axis=1)

control_fit_extra = smf.ols(extraflex, data=df).fit()
control_extra_est = control_fit_extra.params['sex']

n = df.shape[0]
p = control_fit_extra.params.shape[0]

control_extra_se = control_fit_extra.HC0_se[0]*np.sqrt(n/(n-p))
# crude adjustment for the effect of dimensionality on OLS standard errors, motivated by Cattaneo, Jannson, and Newey (2018)
# for really correct way of doing this, we need to implement Cattaneo, Jannson, and Newey (2018)'s procedure.

print("Number of Extra-Flex Controls " + str( p- 1))

print("Coefficient for OLS with extra flex controls " +  str(control_extra_est))
```


```python
# Design matrices
X = smf.ols(extraflex, data=df).data.exog
y = smf.ols(extraflex, data=df).data.endog

# partialling-out the linear effect of W from Y
t_Y_lasso_extra = y_y - lm.LassoCV(cv=10, fit_intercept=False,
                      alphas=np.logspace(-6,2, 100)).fit(X,y).predict(X)

# partialling-out the linear effect of W from D
t_D_lasso_extra = y_d - lm.LassoCV(cv=10, fit_intercept=False,
                      alphas=np.logspace(-6,2, 100)).fit(X,y).predict(X)

# regression of Y on D after partialling-out the effect of W
partial_lasso_fit_extra = sm.OLS(t_Y_lasso_extra, t_D_lasso_extra).fit()
partial_lasso_est_extra = partial_lasso_fit.params[0]
partial_lasso_se_extra = partial_lasso_fit_extra.HC0_se[0]


print("Coefficient for D via partialling-out using lasso " + str(partial_lasso_est_extra))

```


```python
table3 = pd.DataFrame()

table3['Model'] = ["full reg", "partial reg via lasso"]

table3['Estimate'] = [control_extra_est, partial_lasso_est_extra]

table3['Std. Error'] = [control_extra_se, partial_lasso_se_extra]

# Show results
table3
```


```python
# print to Latex
table3.to_latex()
```

In this case $p/n = 20\%$, that is  $p/n$  is no longer small and we start seeing the differences between unregularized partialling out and regularized partialling out with lasso (double lasso). The results based on double lasso have rigorous guarantees in this non-small $p/n$ regime under approximate sparsity. The results based on OLS still have guarantees in $p/n< 1$ regime under assumptions laid out in Cattaneo, Newey, and Jansson (2018), without approximate sparsity, although other regularity conditions are needed.

abfde9c2c5fd6d95688a7ba933307656bf4e4b2b
ipynb
Jupyter Notebook
PM1/inference.ipynb
CausalAIBook/MetricsMLNotebooks
1740f88b3e9ce355bd3cbafd9f988a9cf79f3ff4
MIT
2022-02-04T21:38:55.000Z
2022-02-26T16:02:27.000Z
PM1/inference.ipynb
CausalAIBook/MetricsMLNotebooks
1740f88b3e9ce355bd3cbafd9f988a9cf79f3ff4
MIT
2022-02-08T15:22:49.000Z
2022-03-07T04:07:47.000Z
PM1/inference.ipynb
CausalAIBook/MetricsMLNotebooks
1740f88b3e9ce355bd3cbafd9f988a9cf79f3ff4
MIT
2022-02-08T15:05:34.000Z
2022-02-27T04:55:13.000Z
Qwen/Qwen-72B
1. YES
2. YES