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_pathElectric Company
Analysis of “Electric company” data. See Chapters 1, 16, 19 and 20 in Regression and Other Stories.
Source: ElectricCompany/electric.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
electric = pd.read_csv(ros_path('ElectricCompany/data', 'electric.csv'))
electric.head()
electric_wide = pd.read_csv(ros_path('ElectricCompany/data', 'electric_wide.txt'), sep=r'\s+')
electric_wide.head()| city | grade | treated_pretest | treated_posttest | control_pretest | control_posttest | supplement | |
|---|---|---|---|---|---|---|---|
| 0 | Fresno | 1 | 13.8 | 48.9 | 12.3 | 52.3 | Supplement |
| 1 | Fresno | 1 | 16.5 | 70.5 | 14.4 | 55.0 | Replace |
| 2 | Fresno | 1 | 18.5 | 89.7 | 17.7 | 80.4 | Supplement |
| 3 | Fresno | 1 | 8.8 | 44.2 | 11.5 | 47.0 | Replace |
| 4 | Fresno | 1 | 15.3 | 77.5 | 16.4 | 69.7 | Supplement |
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('y ~ x + t', data=electric, family=lp.gaussian(), chains=4, iter_sampling=1000)
# fit = lp.stan_glm('y ~ x + t', data=electric, family=lp.gaussian(), chains=4, iter_sampling=1000)
# fit = lp.stan_glm('y ~ x + t + x:t', data=electric, family=lp.gaussian(), chains=4, iter_sampling=1000)
# fit = lp.stan_glm('post_test ~ treatment + pre_test + treatment:pre_test', data=electric, family=lp.gaussian(), chains=4, iter_sampling=1000)
# fit = lp.stan_glm('post_test ~ treatment + pre_test + treatment * pre_test', data=electric, family=lp.gaussian(), chains=4, iter_sampling=1000)
# fit = lp.stan_glm('post_test ~ treatment', data=electric, family=lp.gaussian(), chains=4, iter_sampling=1000)
# fit = lp.stan_glm('post_test ~ treatment + pre_test', data=electric, family=lp.gaussian(), chains=4, iter_sampling=1000)
# fit = lp.stan_glm('supp ~ pre_test', data=electric, family=lp.bernoulli(), chains=4, iter_sampling=1000)
# fit = lp.stan_glm('post_test ~ supp + pre_test', data=electric, family=lp.gaussian(), chains=4, iter_sampling=1000)
# fit = lp.stan_glm('y ~ x + t + x*t)', data=electric, family=lp.gaussian(), chains=4, iter_sampling=1000)Notes
- Source computation blocks represented: 32.
- 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.