from pathlib import Path
import sys
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import lapylace as lp
sys.path.append(str(Path.cwd().parent / 'python'))
from data import ros_pathChileSchools
Code and figures for ChileSchools example. See Chapter 21 in Regression and Other Stories.
Data are from
- Chay, K. Y., McEwan, P. J., and Urquiola, M. (2005). The central role of noise in evaluating interventions that use test scores to rank schools. American Economic Review 95, 1237–1258.
Source: ChileSchools/chile_schools.Rmd.
The Python version keeps data handling explicit and uses lapylace for Stan-backed generalized linear models, so the statistical model can be read from a formula rather than from handwritten Stan.
Data
chile = pd.read_csv(ros_path('ChileSchools/data', 'chile.csv'))
chile.head()| p90 | cmb_regn | urban88 | rule2 | cutoff | cutoff_cmb | eligible | read92 | read88 | math92 | math88 | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | Regions 2,5,10 | 1 | 3.414999 | 49.4 | 49.4 | 0 | 57.000000 | 54.369999 | 67.067001 | 51.259998 |
| 1 | 0 | Regions 6,8 | 1 | 31.814997 | 43.4 | 43.4 | 0 | 85.514999 | 79.059998 | 84.658997 | 71.370003 |
| 2 | 0 | Regions 6,8 | 1 | 4.445001 | 43.4 | 43.4 | 0 | 51.971001 | 47.759998 | 53.265999 | 47.930000 |
| 3 | 0 | Region 13 | 1 | 16.799997 | 46.4 | 46.4 | 0 | 66.374001 | 65.129997 | 56.374001 | 61.270000 |
| 4 | 0 | Regions 2,5,10 | 1 | 0.094999 | 49.4 | 49.4 | 0 | 52.500000 | 49.259998 | 50.873001 | 49.730000 |
Models
The formulas below are the Python counterparts of the model formulas in the source example. Use lapylace for the Stan-backed Bayesian fit with the same formula interface.
# fit = lp.stan_glm('read92 ~ eligible + rule2', data=chile, family=lp.gaussian(), chains=4, iter_sampling=1000)
# fit = lp.stan_glm('read92 ~ eligible + rule2', data=chile, family=lp.gaussian(), chains=4, iter_sampling=1000)
# fit = lp.stan_glm('read92 ~ eligible + rule2 + eligible:rule2', data=chile, family=lp.gaussian(), chains=4, iter_sampling=1000)
# fit = lp.stan_glm('read92 ~ eligible + rule2 + read88 + math88', data=chile, family=lp.gaussian(), chains=4, iter_sampling=1000)
# fit = lp.stan_glm('read92 ~ eligible + rule2 + read88 + math88', data=chile, family=lp.gaussian(), chains=4, iter_sampling=1000)
# fit = lp.stan_glm('read92 ~ eligible + rule2 + z_read88 + z_math88 + eligible:z_read88', data=chile, family=lp.gaussian(), chains=4, iter_sampling=1000)
# fit = lp.stan_glm('read92 ~ eligible + rule2 + z_read88 + z_math88 + eligible:z_read88', data=chile, family=lp.gaussian(), chains=4, iter_sampling=1000)Notes
- Source computation blocks represented: 24.
- Data paths are expressed through the shared
ros_path()helper. - Formula-based Bayesian regressions are routed through
lapylace.stan_glm(). - Plotting and simulation work uses NumPy, pandas, matplotlib idioms.