GLM: 负二项回归#
注意
本笔记本使用了一些非 PyMC 依赖库，因此需要专门安装才能运行。请打开下面的下拉菜单以获取更多指导。
额外的依赖安装说明
为了运行此笔记本（在本地或 Binder 上），您不仅需要安装可用的 PyMC 以及所有可选依赖项，还需要安装一些额外的依赖项。有关 PyMC 本身安装的建议，请参阅安装
您可以使用您偏好的包管理器安装这些依赖项，我们在此处提供 pip 和 conda 命令作为示例。
$ pip install seaborn
请注意，如果您想（或需要）从笔记本内部而不是命令行安装软件包，您可以通过运行 pip 命令的变体来安装软件包
import sys
!{sys.executable} -m pip install seaborn
您不应运行 !pip install，因为它可能会将软件包安装在不同的环境中，即使安装了也无法从 Jupyter notebook 中使用。
另一种替代方法是使用 conda
$ conda install seaborn
当使用 conda 安装科学 Python 软件包时，我们建议使用 conda forge
RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)
%config InlineBackend.figure_format = "retina"
az.style.use("arviz-darkgrid")
本笔记本紧密跟随 Jonathan Sedar 的 GLM 泊松回归示例（该示例又受到 Ian Osvald 的一个项目 的启发），不同之处在于，此处的数据是负二项分布而不是泊松分布。
负二项回归用于建模方差高于均值的计数数据。负二项分布 可以被认为是泊松分布，其速率参数是伽马分布的，因此可以调整速率参数以解释增加的方差。
生成数据#
与泊松回归示例中一样，我们假设打喷嚏以某个基线速率发生，并且饮酒、不服用抗组胺药或两者兼有都会增加其频率。
泊松数据#
首先，让我们看一下泊松回归示例中的一些泊松分布数据。
# Mean Poisson 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_pois = 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_pois.groupby(["nomeds", "alcohol"])["nsneeze"].agg(["mean", "var"])
| 均值 | 方差 | ||
|---|---|---|---|
| nomeds | 酒精 | ||
| 否 | 否 | 1.047 | 1.127919 |
| 是 | 2.986 | 2.960765 | |
| 是 | 否 | 5.981 | 6.218858 |
| 是 | 35.929 | 36.064023 |
由于泊松分布随机变量的均值和方差相等，因此样本均值和方差非常接近。
负二项数据#
现在，假设数据集中的每个受试者都患有流感，从而增加了他们打喷嚏的方差（并导致少数不幸的人每天打喷嚏超过 70 次）。如果平均打喷嚏次数保持不变但方差增加，则数据可能服从负二项分布。
# Gamma shape parameter
alpha = 10
def get_nb_vals(mu, alpha, size):
"""Generate negative binomially distributed samples by
drawing a sample from a gamma distribution with mean `mu` and
shape parameter `alpha', then drawing from a Poisson
distribution whose rate parameter is given by the sampled
gamma variable.
"""
g = stats.gamma.rvs(alpha, scale=mu / alpha, size=size)
return stats.poisson.rvs(g)
# Create samples
n = 1000
df = pd.DataFrame(
{
"nsneeze": np.concatenate(
(
get_nb_vals(theta_noalcohol_meds, alpha, n),
get_nb_vals(theta_alcohol_meds, alpha, n),
get_nb_vals(theta_noalcohol_nomeds, alpha, n),
get_nb_vals(theta_alcohol_nomeds, alpha, n),
)
),
"alcohol": np.concatenate(
(
np.repeat(False, n),
np.repeat(True, n),
np.repeat(False, n),
np.repeat(True, n),
)
),
"nomeds": np.concatenate(
(
np.repeat(False, n),
np.repeat(False, n),
np.repeat(True, n),
np.repeat(True, n),
)
),
}
)
df
| nsneeze | 酒精 | nomeds | |
|---|---|---|---|
| 0 | 2 | 否 | 否 |
| 1 | 1 | 否 | 否 |
| 2 | 2 | 否 | 否 |
| 3 | 2 | 否 | 否 |
| 4 | 3 | 否 | 否 |
| ... | ... | ... | ... |
| 3995 | 28 | 是 | 是 |
| 3996 | 63 | 是 | 是 |
| 3997 | 19 | 是 | 是 |
| 3998 | 34 | 是 | 是 |
| 3999 | 40 | 是 | 是 |
4000 rows × 3 columns
df.groupby(["nomeds", "alcohol"])["nsneeze"].agg(["mean", "var"])
| 均值 | 方差 | ||
|---|---|---|---|
| nomeds | 酒精 | ||
| 否 | 否 | 0.990 | 1.137037 |
| 是 | 2.983 | 4.030742 | |
| 是 | 否 | 6.004 | 8.658643 |
| 是 | 35.918 | 169.018294 |
与泊松回归示例中一样，我们看到饮酒和/或不服用抗组胺药会在不同程度上增加打喷嚏的频率。与该示例不同，对于 alcohol 和 nomeds 的每种组合，nsneeze 的方差都高于均值。这表明泊松分布不适合该数据，因为泊松分布的均值和方差相等。
可视化数据#
g = sns.catplot(x="nsneeze", row="nomeds", col="alcohol", data=df, kind="count", aspect=1.5)
# Make x-axis ticklabels less crowded
ax = g.axes[1, 0]
labels = range(len(ax.get_xticklabels(which="both")))
ax.set_xticks(labels[::5])
ax.set_xticklabels(labels[::5]);
/opt/homebrew/Caskroom/miniconda/base/envs/pymc-dev/lib/python3.11/site-packages/seaborn/axisgrid.py:123: UserWarning: The figure layout has changed to tight
self._figure.tight_layout(*args, **kwargs)
/opt/homebrew/Caskroom/miniconda/base/envs/pymc-dev/lib/python3.11/site-packages/seaborn/axisgrid.py:123: UserWarning: The figure layout has changed to tight
self._figure.tight_layout(*args, **kwargs)
负二项回归#
创建 GLM 模型#
COORDS = {"regressor": ["nomeds", "alcohol", "nomeds:alcohol"], "obs_idx": df.index}
with pm.Model(coords=COORDS) as m_sneeze_inter:
a = pm.Normal("intercept", mu=0, sigma=5)
b = pm.Normal("slopes", mu=0, sigma=1, dims="regressor")
alpha = pm.Exponential("alpha", 0.5)
M = pm.ConstantData("nomeds", df.nomeds.to_numpy(), dims="obs_idx")
A = pm.ConstantData("alcohol", df.alcohol.to_numpy(), dims="obs_idx")
S = pm.ConstantData("nsneeze", df.nsneeze.to_numpy(), dims="obs_idx")
λ = pm.math.exp(a + b[0] * M + b[1] * A + b[2] * M * A)
y = pm.NegativeBinomial("y", mu=λ, alpha=alpha, observed=S, dims="obs_idx")
idata = pm.sample()
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [intercept, slopes, alpha]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 11 seconds.
查看结果#
az.plot_trace(idata, compact=False);
# Transform coefficients to recover parameter values
az.summary(np.exp(idata.posterior), kind="stats", var_names=["intercept", "slopes"])
| 均值 | 标准差 | hdi_3% | hdi_97% | |
|---|---|---|---|---|
| 截距 | 0.993 | 0.033 | 0.935 | 1.056 |
| 斜率[nomeds] | 6.048 | 0.221 | 5.635 | 6.438 |
| 斜率[酒精] | 3.006 | 0.115 | 2.789 | 3.215 |
| 斜率[nomeds:酒精] | 1.995 | 0.086 | 1.842 | 2.163 |
az.summary(idata.posterior, kind="stats", var_names="alpha")
| 均值 | 标准差 | hdi_3% | hdi_97% | |
|---|---|---|---|---|
| alpha | 9.909 | 0.487 | 9.05 | 10.887 |
平均值接近我们在生成数据时指定的数值
基准率是一个常数 1。
饮酒使基准率增加三倍。
不服用抗组胺药使基准率增加 6 倍。
饮酒和不服用抗组胺药使预期速率（如果它们的速率是独立的）翻倍。如果它们是独立的，那么两者都做会将基准率增加 3*6=18 倍，但实际上基准率增加了 3*6*2=36 倍。
最后，nsneeze_alpha 的均值也非常接近其 10 的实际值。
另请参阅 bambi's 负二项示例 以供进一步参考。
许可声明#
本示例库中的所有笔记本均根据 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"
}
渲染后可能如下所示