The PyMC3 workflow¶
PyMC3 aims for intuitive and readable, yet powerful syntax that reflects how statisticians describe models. The modeling process generally follows these five steps:
- Encode a probability model by defining the following:
- The prior distributions that quantify knowledge and uncertainty about latent variables
- The likelihood function that conditions the parameters on observed data
- Analyze the posterior using one of the options:
- Obtain a point estimate using MAP inference
- Sample from the posterior using MCMC methods
- Approximate the posterior using variational Bayes.
- Check the model using various diagnostic tools.
- Generate predictions.
The resulting model can be used for inference to gain detailed insights into parameter values as well as to predict outcomes for new data points.
Logistic regression with PyMC3¶
Logistic regression estimates a linear relationship between a set of features and a binary outcome, mediated by a sigmoid function to ensure the model produces probabilities. The frequentist approach resulted in point estimates for the parameters that measure the influence of each feature on the probability that a data point belongs to the positive class, with confidence intervals based on assumptions about the parameter distribution.
In contrast, Bayesian logistic regression estimates the posterior distribution over the parameters itself. The posterior allows for more robust estimates of what is called a Bayesian credible interval for each parameter with the benefit of more transparency about the model’s uncertainty.
Imports¶
import warnings
warnings.filterwarnings('ignore')
%matplotlib inline
from pathlib import Path
import pickle
from collections import OrderedDict
import pandas as pd
import numpy as np
from scipy import stats
from sklearn import preprocessing
from sklearn.model_selection import train_test_split
from sklearn.metrics import (roc_curve, roc_auc_score, confusion_matrix, accuracy_score, f1_score,
precision_recall_curve)
from mlxtend.plotting import plot_confusion_matrix
import theano
import pymc3 as pm
from pymc3.variational.callbacks import CheckParametersConvergence
import statsmodels.formula.api as smf
import arviz as az
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import seaborn as sns
from IPython.display import HTML
data_path = Path('data')
fig_path = Path('figures')
model_path = Path('models')
for p in [data_path, fig_path, model_path]:
if not p.exists():
p.mkdir()
The Data¶
The goal of this dataset is to create a binary classification model that predicts whether or not over $11,000$ customers will subscribe to a term deposit after a marketing campaign the bank has performed, based on $17$ predictor variables. The target variable is given as deposit and takes on a value of $1$ if the customer has subscribed and $0$ otherwise.
This is a slightly imbalanced class problem because there are more customers that did not subscribe the term deposit than did.
Load from Disk¶
data = pd.read_csv('bank.csv')
data.info()
data.head()
Data Preprocessing¶
categorical_features = data.select_dtypes(include=[np.object])
numeric_features = data.select_dtypes(include=[np.number])
cols = numeric_features.columns
numeric_features.loc[:, cols] = preprocessing.scale(numeric_features.loc[:, cols])
le = preprocessing.LabelEncoder()
categorical_features = data[categorical_features.columns].apply(lambda x: le.fit_transform(x))
data = pd.concat([categorical_features, numeric_features], axis=1)
data.info()
data.head()
Quick EDA¶
total = len(data)
plt.figure(figsize=(7,5))
g = sns.countplot(x='deposit', data=data)
g.set_ylabel('Count', fontsize=14)
for p in g.patches:
height = p.get_height()
g.text(p.get_x()+p.get_width()/2.,
height + 1.5,
'{:1.2f}%'.format(height/total*100),
ha="center", fontsize=14, fontweight='bold')
plt.margins(y=0.1)
plt.show()
plt.figure(figsize=(13, 13))
corr = data.corr()
mask = np.tri(*corr.shape).T
sns.heatmap(corr.abs(), mask=mask, annot=True)
b, t = plt.ylim()
b += 0.5
t -= 0.5
plt.ylim(b, t)
plt.show()
n_fts = len(data.columns)
colors = cm.rainbow(np.linspace(0, 1, n_fts))
data.drop('deposit',axis=1).corrwith(data.deposit).sort_values(ascending=True).plot(kind='barh',
color=colors, figsize=(12, 6))
plt.title('Correlation to Target (deposit)')
plt.show()
print('\n',data.drop('deposit',axis=1).corrwith(data.deposit).sort_values(ascending=False))
data['housing'] = -data['housing']
data['contact'] = -data['contact']
data['campaign'] = -data['campaign']
data['poutcome'] = -data['poutcome']
data['loan'] = -data['loan']
corr = data.corr()
cor_target = corr["deposit"]
relevant_features = cor_target[cor_target > 0.08]
relevant_features.sort_values(ascending=False)
sns.pairplot(data[relevant_features.index]);
Models¶
simple_model = 'deposit ~ duration + pdays + previous'
full_model = 'deposit ~ duration + pdays + previous + education + balance + loan + poutcome + campaign + housing + contact'
MAP Inference¶
A probabilistic program consists of observed and unobserved random variables (RVs). We define the observed RVs via likelihood distributions and unobserved RVs via prior distributions. PyMC3 includes numerous probability distributions for this purpose.
The PyMC3 library makes it very straightforward to perform approximate Bayesian inference for logistic regression. Logistic regression models the probability that individual $i$ subscribes to a deposit based on $k$ features.
Logistic Regression:
$$ p(y_i = 1 \mid \beta) = \sigma \beta_0 + \beta_1 x_{i1} + \dots + \beta_k x_{ik}) $$Where $\sigma$ is the logistic function: $$ \sigma(t) = \frac{1}{1 + e^{-1}} $$
Manual Model Specification¶
We will use the context manager with to define a manual_logistic_model that we can refer to later as a probabilistic model:
- The random variables for the unobserved parameters for intercept and $3$ features are expressed using uninformative priors that assume normal distributions with mean $0$ and standard deviation of $100$.
- The likelihood combines the parameters with the data according to the specification of the logistic regression.
- The outcome is modeled as a Bernoulli RV with success probability given by the likelihood.
with pm.Model() as manual_logistic_model:
# random variables for coefficients with
# uninformative priors for each parameter
intercept = pm.Normal('intercept', 0, sd=100)
beta_1 = pm.Normal('beta_1', 0, sd=100)
beta_2 = pm.Normal('beta_2', 0, sd=100)
beta_3 = pm.Normal('beta_3', 0, sd=100)
# Transform random variables into vector of probabilities p(y_i=1)
# according to logistic regression model specification.
likelihood = pm.invlogit(intercept + beta_1 * data.duration + beta_2 * data.pdays + beta_3 * data.previous)
# Bernoulli random vector with probability of success
# given by sigmoid function and actual data as observed
pm.Bernoulli(name='logit', p=likelihood, observed=data.deposit)
manual_logistic_model.model
The command pm.model_to_graphviz(manual_logistic_model) produces the plate notation displayed below.
It shows the unobserved parameters as light and the observed elements as dark circles. The rectangle indicates the number of repetitions of the observed model element implied by the data included in the model definition.
pm.model_to_graphviz(manual_logistic_model)
MAP Estimate¶
Maximum a posteriori probability (MAP) estimation leverages that the evidence is a constant factor that scales the posterior to meet the requirements for a probability distribution. Since the evidence does not depend on $θ$, the posterior distribution is proportional to the product of the likelihood and the prior. Hence, MAP estimation chooses the value of $θ$ that maximizes the posterior given the observed data and the prior belief, that is, the mode of the posterior.
$$ \theta_{MAP} = \operatorname*{arg\,max}_\theta P(X \mid \theta) P(\theta) $$We obtain point MAP estimates for the three parameters using the just defined model’s .find_MAP() method:
with manual_logistic_model:
# compute maximum a-posteriori estimate
# for logistic regression weights
manual_map_estimate = pm.find_MAP()
def print_map(result):
return pd.Series({k: np.asscalar(v) for k, v in result.items()})
print_map(manual_map_estimate)
GLM Model¶
PyMC3 includes numerous common models so that we can usually leave the manual specification for custom applications.
The following code defines the same logistic regression as a member of the Generalized Linear Models (GLM) family using the formula format inspired by the statistical language R and ported to python by the patsy library:
with pm.Model() as logistic_model:
pm.glm.GLM.from_formula(simple_model,
data,
family=pm.glm.families.Binomial())
pm.model_to_graphviz(logistic_model)
PyMC3 solves the optimization problem of finding the posterior point with the highest density using the quasi-Newton Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm but offers several alternatives provided by the scipy library. The result is virtually identically to the corresponding statsmodels estimate:
model = smf.logit(formula=simple_model, data=data[['deposit', 'duration', 'pdays', 'previous']])
result = model.fit()
print(result.summary())
print_map(manual_map_estimate)
result.params
Markov Chain Monte Carlo¶
Markov chains are stochastic models that describe sequences of possible events. Each event comes from a set of outcomes, and each outcome determines which outcome occurs next, according to a fixed set of probabilities. An important feature of Markov chains is that they are memoryless: everything that you would possibly need to predict the next event is available in the current state, and no new information comes from knowing historical events.
Monte Carlo methods rely on repeated random sampling to approximate results that may be deterministic, but that does not permit an analytic, exact solution.
Many algorithms apply the Monte Carlo method to a Markov Chain, and generally proceed as follows:
- Start at the current position.
- Draw a new position from a proposal distribution.
- Evaluate the probability of the new position in light of data and prior distributions:
- If sufficiently likely, move to the new position
- Otherwise, remain at the current position
- Repeat from step 1.
- After a given number of iterations, return all accepted positions.
Define the Model¶
We will use a slightly more complicated model to illustrate Markov chain Monte Carlo inference:
full_model = 'deposit ~ duration + pdays + previous + education + balance + loan + poutcome + campaign + housing + contact'
with pm.Model() as logistic_model:
pm.glm.GLM.from_formula(formula=full_model,
data=data,
family=pm.glm.families.Binomial())
logistic_model.basic_RVs
Hamiltonian Monte Carlo – going NUTS¶
By default, PyMC3 automatically selects the most efficient sampler and initializes the sampling process for efficient convergence. For a continuous model, PyMC3 chooses the NUTS sampler. It also runs variational inference via ADVI to find good starting parameters for the sampler. One among several alternatives is to use the MAP estimate.
Hamiltonian Monte Carlo (HMC) is a hybrid method that leverages the first-order derivative information of the gradient of the likelihood to propose new states for exploration and overcome some of the challenges of MCMC. In addition, it incorporates momentum to efficiently jump around the posterior. As a result, it converges faster to a high-dimensional target distribution than simpler random-walk Metropolis or Gibbs sampling.
To see what the convergence looks like, we first draw $1,000$ samples after tuning the sampler for $1,000$ iterations that will be discarded. The sampling process can be parallelized for multiple chains using the cores argument (except when using GPU).
with logistic_model:
trace = pm.sample(tune=1000,
draws=1000,
chains=4,
init = 'adapt_diag',
cores=3)
Inspect Trace¶
pm.plot_trace(trace);
pm.trace_to_dataframe(trace).info()
with open(model_path / 'logistic_model_mh.pkl', 'wb') as buff:
pickle.dump({'model': logistic_model, 'trace': trace}, buff)
PyMC3 produces various summary statistics for a sampler. These are available as individual functions in the stats module, or by providing a trace to the pm.summary() function:
pm.summary(trace)
draws = 100
trace_df = pm.trace_to_dataframe(trace).assign(
chain=lambda x: x.index // draws)
trace_df.info()
Persist Results¶
with open(model_path / 'logistic_model_nuts.pkl', 'wb') as buff:
pickle.dump({'model': logistic_model,
'trace': trace}, buff)
with open(model_path / 'logistic_model_nuts.pkl', 'rb') as buff:
data0 = pickle.load(buff)
logistic_model, trace_NUTS = data0['model'], data0['trace']
Combine Traces¶
draws = 10000
df = pm.trace_to_dataframe(trace_NUTS).iloc[200:].reset_index(
drop=True).assign(chain=lambda x: x.index // draws)
trace_df = pd.concat([trace_df.assign(samples=100),
df.assign(samples=10000)])
trace_df.info()
Visualize both traces¶
trace_df_long = pd.melt(trace_df, id_vars=['samples', 'chain'])
trace_df_long.info()
Convergence¶
We can visualize the samples over time and their distributions to check the quality of the results. The following charts show the posterior distributions after an initial $100$ and an additional $10,000$ samples, respectively, and illustrate how convergence implies that multiple chains identify the same distribution. The pm.trace_plot() function shows the evolution of the samples as well.
g = sns.FacetGrid(trace_df_long, col='variable', row='samples',
hue='chain', sharex='col', sharey=False)
g = g.map(sns.distplot, 'value', hist=False, rug=False)
Estimate Odds Ratio¶
An odds ratio (OR) is a measure of association between a certain property $A$ and a second property $B$ in a population. Specifically, it tells you how the presence or absence of property $A$ has an effect on the presence or absence of property $B$. The OR is also used to figure out if a particular exposure (like eating processed meat) is a risk factor for a particular outcome (such as colon cancer), and to compare the various risk factors for that outcome. As long as you have two properties you think are linked, you can calculate the odds.
You have two choices for the formula: $$(a \div c) \div (b \div d)$$
or, equivalently: $$(a \times d) \div (b \times c)$$
- $OR > 1$ indicates increased occurrence of event
- $OR < 1$ indicates decreased occurrence of event (protective exposure)
b = trace['duration']
lb, ub = np.percentile(b, 2.5), np.percentile(b, 97.5)
lb, ub = np.exp(lb), np.exp(ub)
print(f'P({lb:.3f} < Odds Ratio < {ub:.3f}) = 0.95')
We can interpret something along those lines: "With probability $0.95$ the odds ratio is greater than $5.310$ and less than $6.233$, so the duration effect takes place because a person with a longer duration has at least $5.310$ higher probability to subscribe to a term deposit than a person with a lower duration, while holding all the other variables constant.
Computing Credible Intervals¶
We can compute the credible intervals, the Bayesian counterpart of confidence intervals, as percentiles of the trace. The resulting boundaries reflect our confidence about the range of the parameter value for a given probability threshold, as opposed to the number of times the parameter will be within this range for a large number of trials.
fig, ax = plt.subplots(figsize=(8, 4))
sns.distplot(np.exp(b), axlabel='Odds Ratio', ax=ax)
ax.set_title(f'Credible Interval: P({lb:.3f} < Odds Ratio < {ub:.3f}) = 0.95')
ax.axvspan(lb, ub, alpha=0.5, color='gray');
Variational Inference¶
Variational Inference (VI) is a machine learning method that approximates probability densities through optimization. In the Bayesian context, it approximates the posterior distribution as follows:
- Select a parametrized family of probability distributions
- Find the member of this family closest to the target, as measured by Kullback-Leibler divergence
Compared to MCMC, Variational Bayes tends to converge faster and scales to large data better. While MCMC approximates the posterior with samples from the chain that will eventually converge arbitrarily close to the target, variational algorithms approximate the posterior with the result of the optimization, which is not guaranteed to coincide with the target.
Variational Inference is better suited for large datasets and to quickly explore many models.In contrast, MCMC will deliver more accurate results on smaller datasets or when time and computational resources pose fewer constraints.
Run Automatic Differentation Variational Inference (ADVI)¶
The interface for variational inference is very similar to the MCMC implementation. We just use the fit() instead of the sample() function, with the option to include an early stopping CheckParametersConvergence callback if the distribution-fitting process converged up to a given tolerance:
with logistic_model:
callback = CheckParametersConvergence(diff='absolute')
approx = pm.fit(n=100000, callbacks=[callback])
Persist Result¶
with open(model_path / 'logistic_model_advi.pkl', 'wb') as buff:
pickle.dump({'model': logistic_model,
'approx': approx}, buff)
Sample from approximated distribution¶
We can draw samples from the approximated distribution to obtain a trace object as above for the MCMC sampler:
trace_advi = approx.sample(2000)
pm.summary(trace_advi)
# pm.summary(trace_advi).to_csv(model_path / 'trace_advi.csv')
Let us visualize the covariance structure of the model
az.plot_pair(trace_NUTS, figsize=(10, 10));
Let us visualize the covariance structure of the ADVI model
az.plot_pair(trace_advi, figsize=(10, 10));
Clearly, ADVI does not capture (as expected) the interactions between variables because of the mean field approximation, and so it underestimates the overall variance.
Model Diagnostics¶
Bayesian model diagnostics includes validating that the sampling process has converged and consistently samples from high-probability areas of the posterior, and confirming that the model represents the data well.
For high-dimensional models with many variables, it becomes cumbersome to inspect numerous traces. When using NUTS, the energy plot helps to assess problems of convergence. It summarizes how efficiently the random process explores the posterior. The plot shows the energy and the energy transition matrix that should be well matched as in the below example.
Energy Plot¶
When using NUTS, the energy plot helps to assess problems of convergence. It summarizes how efficiently the random process explores the posterior. The plot shows the energy and the energy transition matrix, which should be well-matched.
pm.energyplot(trace_NUTS);
Forest Plot¶
A forest plot, also known as a blobbogram, is a graphical display of estimated results from a number of scientific studies addressing the same question, along with the overall results. It was developed for use in medical research as a means of graphically representing a meta-analysis of the results of randomized controlled trials.
az.plot_forest([trace_advi, trace_NUTS]);
Posterior Plot¶
pm.plot_posterior(trace_NUTS);
Posterior Predictive Checks¶
PPCs are very useful for examining how well a model fits the data. They do so by generating data from the model using parameters from draws from the posterior. We use the function pm.sample_ppc for this purpose and obtain $n$ samples for each observation (the GLM module automatically names the outcome $y$):
ppc = pm.sample_ppc(trace_NUTS, samples=500, model=logistic_model)
ppc['y'].shape
Check AUC Score¶
y_score = np.mean(ppc['y'], axis=0)
pred_scores = dict(y_true=data.deposit,y_score=y_score)
roc_auc_score(**pred_scores)
Prediction¶
Follows PyMC3 docs
Predictions use theano’s shared variables to replace the training data with test data before running posterior predictive checks. To facilitate visualization, we create the train and test datasets, and convert the former to a shared variable. Note that we need to use numpy arrays and provide a list of column labels:
Train-test split¶
X = data.drop('deposit', axis=1)
y = data.deposit
X_train, X_test, y_train, y_test = train_test_split(X, y,
test_size=0.2,
random_state=42)
labels = X_train.columns
Create shared theano variable¶
X_shared = theano.shared(X_train.values)
Define logistic model¶
with pm.Model() as logistic_model_pred:
pm.glm.GLM(x=X_shared, labels=labels,
y=y_train, family=pm.glm.families.Binomial())
Run NUTS sampler¶
with logistic_model_pred:
pred_trace = pm.sample(draws=1000,
tune=1000,
chains=2,
cores=2,
init='adapt_diag')
Replace shared variable with test set¶
We then run the sampler as before, and apply the pm.sample_ppc function to the resulting trace after replacing the train with test data:
X_shared.set_value(X_test)
ppc = pm.sample_ppc(pred_trace,
model=logistic_model_pred,
samples=100)
AUC Score¶
y_score = np.mean(ppc['y'], axis=0)
roc_auc_score(y_score=np.mean(ppc['y'], axis=0),
y_true=y_test)
pred_scores = dict(y_true=y_test, y_score=y_score)
cols = ['False Positive Rate', 'True Positive Rate', 'threshold']
roc = pd.DataFrame(dict(zip(cols, roc_curve(**pred_scores))))
precision_recall_curve¶
precision, recall, ts = precision_recall_curve(y_true=y_test, probas_pred=y_score)
pr_curve = pd.DataFrame({'Precision': precision, 'Recall': recall})
f1_score¶
f1 = pd.Series({t: f1_score(y_true=y_test, y_pred=y_score>t) for t in ts})
best_threshold = f1.idxmax()
Plot results¶
fig, axes = plt.subplots(ncols=3, figsize=(13, 5))
ax = sns.scatterplot(x='False Positive Rate', y='True Positive Rate', data=roc, s=50, legend=False, ax=axes[0])
axes[0].plot('False Positive Rate', 'True Positive Rate', data=roc, lw=1, color='k')
axes[0].plot(np.linspace(0,1,100), np.linspace(0,1,100), color='k', ls='--', lw=1)
axes[0].fill_between(y1=roc['True Positive Rate'], x=roc['False Positive Rate'], alpha=.3, color='red')
axes[0].set_title('Receiver Operating Characteristic')
sns.scatterplot(x='Recall', y='Precision', data=pr_curve, ax=axes[1])
axes[1].set_ylim(0,1)
axes[1].set_title('Precision-Recall Curve')
f1.plot(ax=axes[2], title='F1 Scores', ylim=(0,1))
axes[2].set_xlabel('Threshold')
axes[2].axvline(best_threshold, lw=1, ls='--', color='k')
axes[2].text(text=f'Max F1 @ {best_threshold:.2f}', x=.60, y=.95, s=5)
fig.suptitle(f'roc_auc_score = {round(roc_auc_score(**pred_scores),2)}', fontsize=24)
fig.tight_layout()
plt.subplots_adjust(top=.8)
plt.show();
