# Earnings: two-part model
Source: `Earnings/earnings_compound.Rmd`
Annual earnings combine a probability of positive earnings with the level of positive earnings. The Python page uses a Bernoulli `lapylace` model for the first part and a Gaussian `lapylace` model for log earnings among earners.
## Setup
```{python}
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import lapylace as lp
root = Path("../../ROS-Examples")
def coef_medians(fit):
return pd.Series(
[np.median(fit.alpha_draws()), *np.median(fit.beta_draws(), axis=0)],
index=["Intercept", *fit.columns],
)
def fit_glm(formula, data, family=None, seed=1, prior_scale=2.5, intercept_scale=5, aux_scale=10, **kwargs):
return lp.stan_glm(
formula,
data=data,
family=family or lp.gaussian(),
prior=lp.normal(0, prior_scale),
prior_intercept=lp.normal(0, intercept_scale),
prior_aux=lp.exponential(aux_scale),
chains=2,
parallel_chains=2,
iter_warmup=300,
iter_sampling=500,
seed=seed,
refresh=100,
**kwargs,
)
```
```{python}
earnings = pd.read_csv(root / "Earnings/data/earnings.csv").dropna(subset=["earn", "height", "male"])
earnings["positive_earn"] = (earnings["earn"] > 0).astype(int)
pos = earnings.loc[earnings.positive_earn == 1].copy()
pos["log_earn"] = np.log(pos["earn"])
earnings.head()
```
```{python}
fit_positive = fit_glm("positive_earn ~ height + male", earnings, family=lp.bernoulli(), seed=17601, prior_scale=2.5, intercept_scale=10)
fit_logearn = fit_glm("log_earn ~ height + male", pos, seed=17602, prior_scale=2.5, intercept_scale=10, aux_scale=5)
pd.concat({"positive": coef_medians(fit_positive), "log earnings": coef_medians(fit_logearn)}, axis=1).round(3)
```
```{python}
grid = pd.DataFrame({"height": [60, 66, 72], "male": [0, 0, 1]})
prob_positive = fit_positive.posterior_epred(grid).mean(axis=0)
log_level = fit_logearn.posterior_epred(grid).mean(axis=0)
pd.DataFrame({"height": grid.height, "male": grid.male, "Pr(earn>0)": prob_positive, "E(log earn | earn>0)": log_level})
```