{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# GLM: Poisson Regression"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Running on PyMC3 v3.6\n"
]
}
],
"source": [
"## Interactive magics\n",
"%matplotlib inline\n",
"import sys\n",
"import re\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"plt.style.use('seaborn-darkgrid')\n",
"import seaborn as sns\n",
"import patsy as pt\n",
"import pymc3 as pm\n",
"\n",
"plt.rcParams['figure.figsize'] = 14, 6\n",
"np.random.seed(0)\n",
"print('Running on PyMC3 v{}'.format(pm.__version__))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is a minimal reproducible example of Poisson regression to predict counts using dummy data.\n",
"\n",
"This Notebook is basically an excuse to demo Poisson regression using PyMC3, both manually and using the `glm` library to demo interactions using the `patsy` library. We will create some dummy data, Poisson distributed according to a linear model, and try to recover the coefficients of that linear model through inference.\n",
"\n",
"For more statistical detail see:\n",
"\n",
"+ Basic info on [Wikipedia](https://en.wikipedia.org/wiki/Poisson_regression)\n",
"+ GLMs: Poisson regression, exposure, and overdispersion in Chapter 6.2 of [ARM, Gelmann & Hill 2006](http://www.stat.columbia.edu/%7Egelman/arm/)\n",
"+ This worked example from ARM 6.2 by [Clay Ford](http://www.clayford.net/statistics/poisson-regression-ch-6-of-gelman-and-hill/)\n",
"\n",
"This very basic model is inspired by [a project by Ian Osvald](http://ianozsvald.com/2016/05/07/statistically-solving-sneezes-and-sniffles-a-work-in-progress-report-at-pydatalondon-2016/), which is concerned with understanding the various effects of external environmental factors upon the allergic sneezing of a test subject.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Local Functions"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def strip_derived_rvs(rvs):\n",
" '''Convenience fn: remove PyMC3-generated RVs from a list'''\n",
" ret_rvs = []\n",
" for rv in rvs:\n",
" if not (re.search('_log',rv.name) or re.search('_interval',rv.name)):\n",
" ret_rvs.append(rv) \n",
" return ret_rvs\n",
"\n",
"\n",
"def plot_traces_pymc(trcs, varnames=None):\n",
" ''' Convenience fn: plot traces with overlaid means and values '''\n",
"\n",
" nrows = len(trcs.varnames)\n",
" if varnames is not None:\n",
" nrows = len(varnames)\n",
" \n",
" ax = pm.traceplot(trcs, var_names=varnames, figsize=(12,nrows*1.4),\n",
" lines=tuple([(k, {}, v['mean']) \n",
" for k, v in pm.summary(trcs, varnames=varnames).iterrows()]))\n",
"\n",
" for i, mn in enumerate(pm.summary(trcs, varnames=varnames)['mean']):\n",
" ax[i,0].annotate('{:.2f}'.format(mn), xy=(mn,0), xycoords='data',\n",
" xytext=(5,10), textcoords='offset points', rotation=90,\n",
" va='bottom', fontsize='large', color='#AA0022') "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Generate Data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This dummy dataset is created to emulate some data created as part of a study into quantified self, and the real data is more complicated than this. Ask Ian Osvald if you'd like to know more https://twitter.com/ianozsvald\n",
"\n",
"\n",
"### Assumptions:\n",
"\n",
"+ The subject sneezes N times per day, recorded as `nsneeze (int)`\n",
"+ The subject may or may not drink alcohol during that day, recorded as `alcohol (boolean)`\n",
"+ The subject may or may not take an antihistamine medication during that day, recorded as the negative action `nomeds (boolean)`\n",
"+ I postulate (probably incorrectly) that sneezing occurs at some baseline rate, which increases if an antihistamine is not taken, and further increased after alcohol is consumed.\n",
"+ The data is aggregated per day, to yield a total count of sneezes on that day, with a boolean flag for alcohol and antihistamine usage, with the big assumption that nsneezes have a direct causal relationship.\n",
"\n",
"\n",
"Create 4000 days of data: daily counts of sneezes which are Poisson distributed w.r.t alcohol consumption and antihistamine usage"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# decide poisson theta values\n",
"theta_noalcohol_meds = 1 # no alcohol, took an antihist\n",
"theta_alcohol_meds = 3 # alcohol, took an antihist\n",
"theta_noalcohol_nomeds = 6 # no alcohol, no antihist\n",
"theta_alcohol_nomeds = 36 # alcohol, no antihist\n",
"\n",
"# create samples\n",
"q = 1000\n",
"df = pd.DataFrame({\n",
" 'nsneeze': np.concatenate((np.random.poisson(theta_noalcohol_meds, q),\n",
" np.random.poisson(theta_alcohol_meds, q),\n",
" np.random.poisson(theta_noalcohol_nomeds, q), \n",
" np.random.poisson(theta_alcohol_nomeds, q))),\n",
" 'alcohol': np.concatenate((np.repeat(False, q),\n",
" np.repeat(True, q),\n",
" np.repeat(False, q),\n",
" np.repeat(True, q))),\n",
" 'nomeds': np.concatenate((np.repeat(False, q),\n",
" np.repeat(False, q),\n",
" np.repeat(True, q),\n",
" np.repeat(True, q)))})"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" nsneeze \n",
" alcohol \n",
" nomeds \n",
" \n",
" \n",
" \n",
" \n",
" 3995 \n",
" 38 \n",
" True \n",
" True \n",
" \n",
" \n",
" 3996 \n",
" 31 \n",
" True \n",
" True \n",
" \n",
" \n",
" 3997 \n",
" 30 \n",
" True \n",
" True \n",
" \n",
" \n",
" 3998 \n",
" 34 \n",
" True \n",
" True \n",
" \n",
" \n",
" 3999 \n",
" 36 \n",
" True \n",
" True \n",
" \n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" nsneeze \n",
" \n",
" \n",
" nomeds \n",
" False \n",
" True \n",
" \n",
" \n",
" alcohol \n",
" \n",
" \n",
" \n",
" \n",
" False \n",
" 1.018 \n",
" 5.866 \n",
" \n",
" \n",
" True \n",
" 2.938 \n",
" 35.889 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"g = sns.catplot(x='nsneeze', row='nomeds', col='alcohol', data=df,\n",
" kind='count', size=4, aspect=1.5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Observe:**\n",
"\n",
"+ This looks a lot like poisson-distributed count data (because it is)\n",
"+ With `nomeds == False` and `alcohol == False` (top-left, akak antihistamines WERE used, alcohol was NOT drunk) the mean of the poisson distribution of sneeze counts is low.\n",
"+ Changing `alcohol == True` (top-right) increases the sneeze count `nsneeze` slightly\n",
"+ Changing `nomeds == True` (lower-left) increases the sneeze count `nsneeze` further\n",
"+ Changing both `alcohol == True and nomeds == True` (lower-right) increases the sneeze count `nsneeze` a lot, increasing both the mean and variance."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Poisson Regression"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our model here is a very simple Poisson regression, allowing for interaction of terms:\n",
"\n",
"$$ \\theta = exp(\\beta X)$$\n",
"\n",
"$$ Y_{sneeze\\_count} ~ Poisson(\\theta)$$\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Create linear model for interaction of terms**"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"fml = 'nsneeze ~ alcohol + antihist + alcohol:antihist' # full patsy formulation"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"fml = 'nsneeze ~ alcohol * nomeds' # lazy, alternative patsy formulation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 1. Manual method, create design matrices and manually specify model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Create Design Matrices**"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"(mx_en, mx_ex) = pt.dmatrices(fml, df, return_type='dataframe', NA_action='raise')"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" Intercept \n",
" alcohol[T.True] \n",
" nomeds[T.True] \n",
" alcohol[T.True]:nomeds[T.True] \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 1.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" \n",
" \n",
" 1 \n",
" 1.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" \n",
" \n",
" 2 \n",
" 1.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" \n",
" \n",
" 3997 \n",
" 1.0 \n",
" 1.0 \n",
" 1.0 \n",
" 1.0 \n",
" \n",
" \n",
" 3998 \n",
" 1.0 \n",
" 1.0 \n",
" 1.0 \n",
" 1.0 \n",
" \n",
" \n",
" 3999 \n",
" 1.0 \n",
" 1.0 \n",
" 1.0 \n",
" 1.0 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rvs_fish = [rv.name for rv in strip_derived_rvs(mdl_fish.unobserved_RVs)]\n",
"plot_traces_pymc(trc_fish, varnames=rvs_fish)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Observe:**\n",
"\n",
"+ The model converges quickly and traceplots looks pretty well mixed"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Transform coeffs and recover theta values"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" hpd_2.5 \n",
" hpd_97.5 \n",
" \n",
" \n",
" \n",
" \n",
" b0_intercept \n",
" 1.017360 \n",
" 0.956431 \n",
" 1.075815 \n",
" \n",
" \n",
" b1_alcohol[T.True] \n",
" 2.886177 \n",
" 2.696525 \n",
" 3.100204 \n",
" \n",
" \n",
" b2_nomeds[T.True] \n",
" 5.765253 \n",
" 5.404684 \n",
" 6.128869 \n",
" \n",
" \n",
" b3_alcohol[T.True]:nomeds[T.True] \n",
" 2.120270 \n",
" 1.973526 \n",
" 2.280686 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rvs_fish_alt = [rv.name for rv in strip_derived_rvs(mdl_fish_alt.unobserved_RVs)]\n",
"plot_traces_pymc(trc_fish_alt, varnames=rvs_fish_alt)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Transform coeffs"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" hpd_2.5 \n",
" hpd_97.5 \n",
" \n",
" \n",
" \n",
" \n",
" Intercept \n",
" 1.018445e+00 \n",
" 0.956916 \n",
" 1.082025e+00 \n",
" \n",
" \n",
" alcohol[T.True] \n",
" 2.884207e+00 \n",
" 2.693380 \n",
" 3.100136e+00 \n",
" \n",
" \n",
" nomeds[T.True] \n",
" 5.757942e+00 \n",
" 5.368011 \n",
" 6.129526e+00 \n",
" \n",
" \n",
" alcohol[T.True]:nomeds[T.True] \n",
" 2.121873e+00 \n",
" 1.972392 \n",
" 2.295317e+00 \n",
" \n",
" \n",
" mu \n",
" 3.340292e+18 \n",
" 1.005653 \n",
" 3.927977e+54 \n",
" \n",
" \n",
"
\n",
"