GLM: 泊松回归#
这是一个使用虚拟数据预测计数的泊松回归的最小可重现示例。
本笔记本基本上是演示使用 PyMC 进行泊松回归的借口，包括手动方式和使用 bambi 演示使用 formulae 库的交互作用。我们将创建一些虚拟数据，这些数据根据线性模型进行泊松分布，并尝试通过推断恢复该线性模型的系数。
有关更多统计详细信息，请参阅
维基百科上关于泊松回归的基本信息
ARM, Gelmann & Hill 2006 第 6.2 章中的 GLM：泊松回归、暴露和过度离散
Clay Ford 提供的 ARM 6.2 中的工作示例
这个非常基本的模型灵感来自 Ian Osvald 的一个项目，该项目关注于理解外部环境因素对测试对象过敏性喷嚏的各种影响。
#!pip install seaborn
import arviz as az
import bambi as bmb
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import seaborn as sns
from formulae import design_matrices
RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)
%config InlineBackend.figure_format = 'retina'
az.style.use("arviz-darkgrid")
本地函数#
生成数据#
此虚拟数据集的创建是为了模拟作为量化自我研究一部分创建的一些数据，而真实数据比这更复杂。如果您想了解更多信息，请咨询 Ian Osvald @ianozsvald。
假设：#
受试者每天打喷嚏 N 次，记录为
nsneeze (int)受试者当天可能饮酒也可能不饮酒，记录为
alcohol (boolean)受试者当天可能服用也可能不服用抗组胺药物，记录为负面行为
nomeds (boolean)我们假设（可能不正确）喷嚏以某个基线速率发生，如果不服用抗组胺药，则会增加，饮酒后会进一步增加。
数据按天聚合，以产生当天的喷嚏总数，以及饮酒和抗组胺药物使用的布尔标志，并假设喷嚏数具有直接的因果关系。
创建 4000 天的数据：每日喷嚏计数，这些计数相对于饮酒和抗组胺药物使用情况呈泊松分布
# decide poisson theta values
theta_noalcohol_meds = 1 # no alcohol, took an antihist
theta_alcohol_meds = 3 # alcohol, took an antihist
theta_noalcohol_nomeds = 6 # no alcohol, no antihist
theta_alcohol_nomeds = 36 # alcohol, no antihist
# create samples
q = 1000
df = pd.DataFrame(
{
"nsneeze": np.concatenate(
(
rng.poisson(theta_noalcohol_meds, q),
rng.poisson(theta_alcohol_meds, q),
rng.poisson(theta_noalcohol_nomeds, q),
rng.poisson(theta_alcohol_nomeds, q),
)
),
"alcohol": np.concatenate(
(
np.repeat(False, q),
np.repeat(True, q),
np.repeat(False, q),
np.repeat(True, q),
)
),
"nomeds": np.concatenate(
(
np.repeat(False, q),
np.repeat(False, q),
np.repeat(True, q),
np.repeat(True, q),
)
),
}
)
df.tail()
| nsneeze | alcohol | nomeds | |
|---|---|---|---|
| 3995 | 26 | True | True |
| 3996 | 35 | True | True |
| 3997 | 36 | True | True |
| 3998 | 32 | True | True |
| 3999 | 35 | True | True |
查看各种组合的均值（泊松均值）#
df.groupby(["alcohol", "nomeds"]).mean().unstack()
| nsneeze | ||
|---|---|---|
| nomeds | False | True |
| alcohol | ||
| False | 1.047 | 5.981 |
| True | 2.986 | 35.929 |
简要描述数据集#
g = sns.catplot(
x="nsneeze",
row="nomeds",
col="alcohol",
data=df,
kind="count",
height=4,
aspect=1.5,
)
for ax in (g.axes[1, 0], g.axes[1, 1]):
for n, label in enumerate(ax.xaxis.get_ticklabels()):
label.set_visible(n % 5 == 0)
/Users/reshamashaikh/miniforge3/envs/pymc-ex/lib/python3.10/site-packages/seaborn/axisgrid.py:118: UserWarning: This figure was using constrained_layout, but that is incompatible with subplots_adjust and/or tight_layout; disabling constrained_layout.
self._figure.tight_layout(*args, **kwargs)
观察
这看起来很像泊松分布的计数数据（因为它就是）
当
nomeds == False且alcohol == False时（左上角，即使用了抗组胺药，没有饮酒），喷嚏计数的泊松分布的均值较低。更改
alcohol == True（右上角）会稍微增加喷嚏计数nsneeze更改
nomeds == True（左下角）会进一步增加喷嚏计数nsneeze同时更改
alcohol == True and nomeds == True（右下角）会大大增加喷嚏计数nsneeze，同时增加均值和方差。
泊松回归#
我们这里的模型是一个非常简单的泊松回归，允许项的交互
为项的交互创建线性模型
fml = "nsneeze ~ alcohol + nomeds + alcohol:nomeds" # full formulae formulation
fml = "nsneeze ~ alcohol * nomeds" # lazy, alternative formulae formulation
1. 手动方法，创建设计矩阵并手动指定模型#
创建设计矩阵
dm = design_matrices(fml, df, na_action="error")
mx_ex = dm.common.as_dataframe()
mx_en = dm.response.as_dataframe()
mx_ex
| 截距 | alcohol | nomeds | alcohol:nomeds | |
|---|---|---|---|---|
| 0 | 1 | 0 | 0 | 0 |
| 1 | 1 | 0 | 0 | 0 |
| 2 | 1 | 0 | 0 | 0 |
| 3 | 1 | 0 | 0 | 0 |
| 4 | 1 | 0 | 0 | 0 |
| ... | ... | ... | ... | ... |
| 3995 | 1 | 1 | 1 | 1 |
| 3996 | 1 | 1 | 1 | 1 |
| 3997 | 1 | 1 | 1 | 1 |
| 3998 | 1 | 1 | 1 | 1 |
| 3999 | 1 | 1 | 1 | 1 |
4000 行 × 4 列
创建模型
with pm.Model() as mdl_fish:
# define priors, weakly informative Normal
b0 = pm.Normal("Intercept", mu=0, sigma=10)
b1 = pm.Normal("alcohol", mu=0, sigma=10)
b2 = pm.Normal("nomeds", mu=0, sigma=10)
b3 = pm.Normal("alcohol:nomeds", mu=0, sigma=10)
# define linear model and exp link function
theta = (
b0
+ b1 * mx_ex["alcohol"].values
+ b2 * mx_ex["nomeds"].values
+ b3 * mx_ex["alcohol:nomeds"].values
)
## Define Poisson likelihood
y = pm.Poisson("y", mu=pm.math.exp(theta), observed=mx_en["nsneeze"].values)
采样模型
with mdl_fish:
inf_fish = pm.sample()
# inf_fish.extend(pm.sample_posterior_predictive(inf_fish))
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [Intercept, alcohol, nomeds, alcohol:nomeds]
Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 36 seconds.
查看诊断
az.plot_trace(inf_fish);
观察
模型快速收敛，迹图看起来混合良好
转换系数并恢复 theta 值#
az.summary(np.exp(inf_fish.posterior), kind="stats")
| 均值 | 标准差 | hdi_3% | hdi_97% | |
|---|---|---|---|---|
| 截距 | 1.045 | 0.034 | 0.981 | 1.106 |
| alcohol | 2.862 | 0.109 | 2.668 | 3.065 |
| nomeds | 5.734 | 0.205 | 5.341 | 6.115 |
| alcohol:nomeds | 2.102 | 0.086 | 1.948 | 2.266 |
观察
每个特征作为基线喷嚏计数的乘数的贡献似乎与数据生成一致
exp(Intercept): 均值=1.05 cr=[0.98, 1.10]
当没有饮酒和服药时，大致线性基线计数，与生成的数据一致
theta_noalcohol_meds = 1 (如上设置) theta_noalcohol_meds = exp(Intercept) = 1
exp(alcohol): 均值=2.86 cr=[2.67, 3.07]
添加酒精的非零正效应，基线喷嚏计数的约 3 倍乘数，与生成的数据一致
theta_alcohol_meds = 3 (如上设置) theta_alcohol_meds = exp(Intercept + alcohol) = exp(Intercept) * exp(alcohol) = 1 * 3 = 3
exp(nomeds): 均值=5.73 cr=[5.34, 6.08]
添加 nomeds 的更大、非零正效应，基线喷嚏计数的约 6 倍乘数，与生成的数据一致
theta_noalcohol_nomeds = 6 (如上设置) theta_noalcohol_nomeds = exp(Intercept + nomeds) = exp(Intercept) * exp(nomeds) = 1 * 6 = 6
exp(alcohol:nomeds): 均值=2.10 cr=[1.96, 2.28]
酒精和药物的小而正的交互作用效应，基线喷嚏计数的约 2 倍乘数，与生成的数据一致
theta_alcohol_nomeds = 36 (如上设置) theta_alcohol_nomeds = exp(Intercept + alcohol + nomeds + alcohol:nomeds) = exp(Intercept) * exp(alcohol) * exp(nomeds * alcohol:nomeds) = 1 * 3 * 6 * 2 = 36
2. 替代方法，使用 bambi#
创建模型
使用 bambi 的替代自动公式
model = bmb.Model(fml, df, family="poisson")
拟合模型
inf_fish_alt = model.fit()
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [Intercept, alcohol, nomeds, alcohol:nomeds]
Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 36 seconds.
查看迹图
az.plot_trace(inf_fish_alt);
转换系数#
az.summary(np.exp(inf_fish_alt.posterior), kind="stats")
| 均值 | 标准差 | hdi_3% | hdi_97% | |
|---|---|---|---|---|
| 截距 | 1.049 | 0.033 | 0.987 | 1.107 |
| alcohol | 2.849 | 0.105 | 2.659 | 3.047 |
| nomeds | 5.708 | 0.190 | 5.362 | 6.058 |
| alcohol:nomeds | 2.111 | 0.083 | 1.960 | 2.258 |
观察
迹图看起来混合良好
转换后的模型系数看起来或多或少与手动模型生成的系数相同
请注意，后验预测样本具有极端的偏斜
posterior_predictive = model.predict(inf_fish_alt, kind="pps")
我们可以使用 az.plot_ppc() 来检查后验预测样本是否与观察到的数据相似。
有关后验预测检查的更多信息，我们可以参考 先验和后验预测检查。
az.plot_ppc(inf_fish_alt);
水印#
%load_ext watermark
%watermark -n -u -v -iv -w -p pytensor,aeppl
Last updated: Wed Nov 30 2022
Python implementation: CPython
Python version : 3.10.5
IPython version : 8.4.0
pytensor: 2.7.5
aeppl : 0.0.32
matplotlib: 3.5.2
seaborn : 0.12.1
pymc : 4.1.2
numpy : 1.21.6
pandas : 1.4.3
arviz : 0.12.1
bambi : 0.9.1
Watermark: 2.3.1
许可声明#
本示例库中的所有笔记本均根据 MIT 许可证提供，该许可证允许修改和再分发用于任何用途，前提是保留版权和许可声明。
引用 PyMC 示例#
要引用此笔记本，请使用 Zenodo 为 pymc-examples 存储库提供的 DOI。
重要提示
许多笔记本改编自其他来源：博客、书籍... 在这种情况下，您也应该引用原始来源。
还请记住引用您的代码使用的相关库。
这是一个 bibtex 中的引用模板
@incollection{citekey,
author = "<notebook authors, see above>",
title = "<notebook title>",
editor = "PyMC Team",
booktitle = "PyMC examples",
doi = "10.5281/zenodo.5654871"
}
渲染后可能看起来像