Discussion¶
In CmdStanPy, fitting a model to data is straightforward.
Instantiate a CmdStanModel object
Assemble a data dictionary or JSON file which contains definitions for all data variables declared in the model's data block.
Run one of the available inference methods:
sample to do exact Bayesian estimationvariational to do approximate Bayesian estimationoptimize to do penalized maximum likelihood estimation
Extract the estimates for parameters and quantities of interest. CmdStanPy's accessor functions make it easy to get the outputs in whatever data format is appropriate for downstream analysis.
Fitting the model to the data is just one component of the analysis.
The bulk of the code in this notebook is devoted to visualizing both the raw data and the model estimates;
data visualization drives both the model specification process and the model testing process.
With the plotnine package we can create multi-layer plots which overly the raw data with the model estimates
and predictions or multi-layer plots which demonstrate the behaviors of different models.
These activities play a central role in developing trustworthy data analysis pipelines:
Plotting the raw or simulated data before modeling informs the model design process
Plotting prior and posterior predictive checks drives testing
Plotting model estimates and predictions drives documentation and dissemination of results.
In this notebook we work through the multilevel regression model from Gelman and Hill, chapter 12,
restricting our analysis to just the data from the state of Minnesota.
Even if we weren't writing this up as part of an introduction to multi-level modeling,
that is, if we were tasked with this analysis qua analysis, we would proceed similarly -
starting with a very simple model.
There are many possible next steps:
- Expanding the analysis to other states or the entire US.
- Expanding the model to use the county-level soil uranium level measurements
- as a varying-intercept / varying-slope model (Chapter 12)
- as a spatial-effects model
References and Resources¶
This notebook is based on Chris Fonnesbeck's excellent A Primer on Bayesian Multilevel Modeling using PyStan, which was developed as part of a Stan workshop for biomedical statisticians at Vanderbilt University.
Stan Tutorials YouTube Playlist Maggie Lieu - a series of introductory videos on Bayesian modeling with Stan
Stan User's Guide
Visualization in Bayesian workflow
Jonah Gabry, Daniel Simpson, Aki Vehtari, Michael Betancourt, and Andrew Gelman, 2019
Designing for interactive exploratory data analysis requires theories of graphical inference Jessica Hullman and Andrew Gelman, 2021
Data analysis using regression and multilevel/hierarchical models Andrew Gelman and Jennifer Hill, 2007
Statistical rethinking by Richard McElreath -
an intro-stats/linear models course taught from a Bayesian perspective.
Making Plots With plotnine - plotnine tutorial notebook.
Papers by Price and Gelman on the radon data and multilevel models:
Centralized analysis of local data, with dollars and lives on the line: Lessons from the home radon experience. Phillip N. Price and Andrew Gelman, 2015
Analysis of local decisions using hierarchical modeling, applied to home radon measurement and remediation (with discussion) Phillip N. Price and Andrew Gelman, 1999
Acknowledgement and thanks!¶
- Bob Carpenter
- Chris Fonnesbeck
- Jonah Gabry
- Andrew Gelman
- Reshama Shaikh
- Brian Ward
Appendix A: Linear regression review (chapters 2 and 3 of Gelman and Hill)¶
Linear regression models the relationship between a scalar response and one or more prediexplanatory variables.
A simple linear regression where the observed data (red) are assumed to be
the result of random deviations (green) from an underlying relationship (blue)
between the dependent variable (y) and an independent variable (x).
The model fits the data by finding the linear function (a non-vertical straight line)
which minimizes the distance between "$y_i$ at offset $x_i$ and the line at offset $x_i$.
$y_i = \alpha \, + \beta\,x_i + {\epsilon}_i$
- $\alpha$ is the intercept, the offset from zero on the x-axis
- $\beta$ is the slope, the multiplier applied to x.
- ${\epsilon}_i$ is the error term, assuming independent errors drawn from a normal distribution with mean $0$, standard deviation $\sigma$.
In a simple linear regression, there is only a single predictor. In an ordinary linear regression,
there are many predictors, i.e., the model fits a vector of coefficients $\beta$ to the vector
of independent variables $X$.
Interpretation
In Chapter 3, Gelman and Hill write:
Linear regression is a method that summarizes how the average values of a numerical outcome variable vary over subpopulations defined by linear functions of predictors. ... Regression can be used to predict an outcome given a linear function of these predictors, and regression coefficients can be thought of as comparisons across predicted values or as comparisons among averages in the data.
Linear regression: two ways to write the model
The goal of inference is to learn from incomplete or imperfect data.
In the simple linear regression model, the error term $\epsilon$ accounts for imperfect measurements of the data.
$
y_i = \alpha \, + {\beta}\,x_i + {\epsilon}_i
$
where the errors ${\epsilon}_i$ have independent normal distributions with mean $0$ and standard deviation $\sigma$.
An equivalent representation is
$
y_i ∼ \mathrm{N}(\alpha + \beta\,\mathrm{X}_i,\, \sigma^2), \ \ \ \mathrm{for}\ i=1, ..., n
$
The corresponding Stan statement is
y ~ normal(alpha + beta * x, sigma);
Stan provides vectorized versions of all univariate probability distrions.
This statement is far more efficient than using a for loop over all $x_i$ and $y_i$ pairs.
Simple linear regression model in Stan
The simple linear regression with a single predictor and a slope and intercept coefficient and normally distributed noise is the first model discussed in the Stan User's Guide Regression Models chapter.
data {
int<lower=0> N;
vector[N] x;
vector[N] y;
}
parameters {
real alpha;
real beta;
real<lower=0> sigma;
}
model {
y ~ normal(alpha + beta * x, sigma);
}
This model is the minimal possible model.
It consists of three named program blocks:
the data and parameters are declared in the respectively named data and parameters block;
the model block specifies the likelihood, i.e., the probability of the data given the model.
Because no priors are specified on the model parameters, they are given the default prior
distribution, which is uniform from $-\infty$ to $+\infty$.
Appendix B: Data preparation using pandas.¶
The steps required to convert the CSV files from the Gelman and Hill ARM website
into the data structures used in this analysis:
- Merge the county-level soil uranium level measurement into the home radon data.
- Put both outcome and predictors on the log scale, following Gelman and Hill, chapter 4, section 12.
- Restrict the dataset to Minnesota.
- Aggregate county-level information
Pandas objects contain structured arrays which are labeled by
Index objects.
A Series object
manages 1-D arrays and a DataFrame
is a 2-D size-mutable, table, allowing for heterogenous columns.
These index labels are used to perform database-like select, join, and group-by operations.
However, sometimes we need to access just the data, not the index labels.
Most data is backed by an NumPy ndarray.
The array property
returns the underlying numpy.ndarray of the Index or Series object.
Extract relevant columns from CSV as pandas DataFrame
We leverage the pandas.read_csv function
to extract the information we need from the raw CSV files with the following non-default arguments
- parameter
usecols allows us to extract just the relevant columns for this analysis.
- parameter
skipinitialspace strips out initial whitespace from the data.
Once instantiated, we call the convert_dtypes method on the newly instantiated DataFrame so that we can do merge and join operations on all columns.
Image from https://www.health.state.mn.us/communities/environment/air/radon/index.html
