➄ 广义线性模型

其中是一个函数，我们将调用反向链接函数。有许多反向链接函数可供选择；可能最简单的是恒等函数。这是一个返回与其参数相同的值的函数。第3章“线性回归建模”中的所有模型都使用了单位函数，为简单起见，我们只是省略了它。身份功能本身可能不是很有用，但它允许我们以更统一的方式考虑几种不同的模型。

广义线性模型

在上一章中，我们使用输入变量的线性组合来预测输出变量的平均值。我们假设后者为高斯分布。在许多情况下都可以使用高斯分布，但对于其他许多情况，选择不同的分布可能更明智；当我们用 $t$ 分布替换高斯分布时，我们已经看到了一个这样的例子。在本章中，我们将看到更多使用高斯分布以外分布的明智例子。正如我们将了解到的，存在一个通用的主题或模式，可将线性模型推广到许多问题。在本章中，我们将探讨：

• 广义线性模型
• Logistic回归和逆链接函数
• 简单Logistic回归
• 多元Logistic回归
• Softmax函数和多项Logistic回归
• Poisson回归
• 零膨胀Poisson回归

4.1 广义线性模型

本章的核心思想之一相当简单：为了预测输出变量的平均值，我们可以对输入变量的线性组合应用任意函数。

$$\mu = f ( \alpha + X \beta )$$

其中，$f$ 称作逆连结函数。为什么这里把 $f$ 称作逆连结函数而不是连结函数？原因是传统上人们认为此类函数是用来连结输出变量和线性模型的，不过在构建贝叶斯模型的时候，反过来思考可能更容易一些。因此，为了避免疑惑，我们这里统称逆连结函数。前一章中的所有线性模型其实都包含一个逆连结函数，不过我们书写的时候将其省略了，因为它其实是一个恒等函数（函数的返回值和输入值相同）。恒等函数在这里也许没什么用，不过它可以让我们用一种更一般的形式思考不同的模型。

希望使用逆链接函数的一种情况是在处理定类变量时，例如：颜色名称、性别、生物物种或政党/从属关系。所有这些变量都没有被高斯学派很好地建模。想想看，高斯函数原则上适用于在实数集上取任何值的连续变量，而定类变量是离散的，只取几个值(例如红色、绿色或蓝色)。如果改变用来对数据建模的分布，通常也需要改变对这些分布平均值的似是而非的建模方式。例如，如果我们使用二项式分布，就像在第1章和第2章中一样，我们将需要一个返回 [0，1] 区间平均值的线性模型；实现这一点的一种方法是保留线性模型，但使用逆链接函数将输出限制在所需的区间内。这一技巧并不局限于离散变量；我们可能希望对只能取正值的数据进行建模，因此可能希望将线性模型限制为返回分布平均值的正值，例如Gamma或指数。

在继续之前，请注意一些变量可以被编码为定量或定性，这是您必须根据您的问题的上下文做出的决定；例如，如果我们谈论颜色名称，我们可以讨论红色和绿色分类变量，或者如果我们使用波长，我们可以讨论650 nm和510 nm连续变量。

4.2 Logistic回归

回归问题是关于在给定一个或多个输入变量的值的情况下，预测输出变量的连续值。于此不同，分类是在给定一些输入变量时，将离散值（代表离散类）赋给输出变量。在两种情况下，任务都是获得一个模型，该模型能够正确地模拟输出和输入变量之间的映射；为了做到这一点，需要一个具有正确的输出-输入变量对的样本。从机器学习角度看，回归和分类都是有监督学习算法的实例。

我妈妈准备了一道名为 sopa seca 的美味菜肴，基本上是以意大利面为基础的食谱，字面意思是干汤。虽然这听起来可能是用词不当，甚至是自相矛盾，但当我们了解它是如何烹调的时候，这道菜的名字就完全有意义了。Logistic回归也发生了类似的情况，这是一种尽管名为Logistic回归的模型，但通常被认为是解决分类问题的一种方法。

我母亲有一道拿手好菜叫 sopa seca ，其做法源于意大利面，直译过来叫做干汤面，这个名字听起来似乎有点矛盾，不过一旦了解它的做法之后，你会觉得这个叫法其实挺合理的。同样，逻辑斯蒂回归也是如此，它主要用来解决分类问题。

逻辑斯蒂回归模型是线性回归模型的一种扩展，这也是其名称来源。其模型是将 4. 1中的逆连接函数设定为逻辑斯蒂函数而形成，逻辑斯蒂函数为以下形式：

$$\text{logistic}(z)=\frac{1} { 1 + e ^ { - z } }$$

从分类的角度来看，逻辑斯蒂函数最重要的特点之一是不论参数 $z$ 的值为多少，其输出值总是介于0到1之间。因此，该函数将整个实轴压缩到了区间 [0,1] 内。逻辑斯蒂函数也称作 $S$ 型函数（sigmoid function），因为它的形状看起来像 $S$ ，可以运行以下几行代码来看一下。

z = np.linspace(-8, 8)
plt.plot(z, 1 / (1 + np.exp(-z)))
plt.xlabel('z')
plt.ylabel('logistic(z)')

（1）逻辑斯蒂模型

我们几乎具备了将简单线性回归转化为简单Logistic回归的所有要素，接下来将继续学习如何用来对结果进行分类。先从简单问题开始，假设类别只有两类，比如正常邮件/垃圾邮件、安全/不安全、阴天/晴天、健康/生病等。首先对类别进行编码，假设变量 $y$ 只能有两个值 0 或 1 。这样描述之后，问题就有点像抛硬币问题了，在那个例子中我们用到了伯努利分布作为似然。这里的区别是：现在 $q$ 不再是从 $beta$ 分布中生成，而是由一个线性模型定义的。该线性模型可以返回实轴上的任意值，但是伯努利分布的值限定在 [0,1] 区间内。因此我们使用了一个逆连结函数将线性模型的返回值映射到一个适合伯努利分布的区间内，从而将一个线性回归模型转换成了分类模型：

$$\theta=\operatorname{logistic}(\alpha+x \beta) \\ y=\operatorname{Bern}(\theta)$$

请注意，与第3章中的简单线性回归的主要区别在于：使用了伯努利分布而不是高斯分布，并且使用了Logistic函数而不是恒等函数。

下面的Kruschke图显示了逻辑斯蒂回归模型（应该）包含的先验。注意与单参数线性回归模型的区别是：这里用到了伯努利分布而不是高斯分布（或者t分布），并使用了逻辑函数生成区间[0,1]范围内的参数θ，从而适于输入到伯努利分布。

（2）鸢尾花数据集

我们将逻辑斯蒂回归应用到鸢尾花数据集上。在构建模型之前，我们先了解下该数据集。鸢尾花数据集是一个很经典的数据集，包含有Setosa、Versicolour和Virginica 3个种类，这3个类别标签就是我们想要预测的量，即因变量。其中每个种类包含有50个数据，每个数据包含4种变量（或者称为特征，机器学习中通常这么称呼），这4种变量就是我们要分析的自变量，分别是：花瓣长度、花瓣宽度、花萼长度和花萼宽度。花萼有点类似小叶，在花还是芽时包围着花，有保护作用。seaborn中包含鸢尾花数据集，我们可以用如下代码将其导入成Pandas里的Dataframe格式：

iris = pd.read_csv('../data/iris.csv')
iris.head()

现在，可以使用 seaborn 中的 striplot 函数绘制这三个物种与 sepal_length的关系图：

sns.stripplot(x="species", y="sepal_length", data=iris, jitter=True)

上图中y轴是连续的，而x轴是离散的类别变量；图中的点在x轴上是分散开来的，这并没有什么实际意义，只是一个画图的小技巧而已，通过将jitter参数设为True，能够避免所有点都重叠在一条直线上，你可以试着将jitter参数设为False。这里唯一重要的是x轴的含义，即分别代表setosa、versicolour和virginica 3个类别。你还可以用seaborn中的其他作图函数来画这些点（比如violin plot），只需要一行代码就能完成。

另外一种观察数据的方法是用pairplot画出散点图矩阵，用该函数可以得到一个 4×4的网格（因为我们有4种特征）。网格是对称的，上三角和下三角表示的是同样的信息。由于对角线上的散点图其实是变量本身，因此这里用一个特征的KDE图代替了散点图。可以看到，每个子图中，分别用3种颜色表示3种不同的类别标签，这一点与前面图中的表示一致。

sns.pairplot(iris, hue='species', diag_kind='kde')

深入学习之前，花点时间分析下前面这幅图，进一步熟悉这个数据集并了解变量与标签之间的关系。

（3）将逻辑斯蒂回归模型应用到鸢尾花数据集

我们先从一个简单的问题开始：用花萼长度这一特征（自变量）来区分setosa 和 versicolor这两个种类。和前面一样，这里用0和1对类别变量进行编码，利用Pandas可以这么做：

df = iris.query(species == ('setosa', 'versicolor'))
y_0 = pd.Categorical(df['species']).codes
x_n = 'sepal_length'
x_0 = df[x_n].values

现在数据已经表示成了合适的格式，终于可以用 PyMC3 来建模了。留意下面代码中的第一部分与线性回归模型的相似之处。此外还留意两个确定变量：$\theta$ 和 bd 。 $\theta$ 是对变量 $\mu$ 应用逻辑函数之后的值，bd 是一个有边界的值，用于确定分类结果，稍后会详细讨论。还有一点值得提的是除了像下面这样明确写出逻辑函数的完整形式之外，还可以简单地使用 PyMC3 中的 pm.math.sigmoid 函数，该函数是 Theano 中 sigmoid 函数的一个别名。

with pm.Model() as model_0:
    α = pm.Normal('α', mu=0, sd=10)
    β = pm.Normal('β', mu=0, sd=10)
    μ = α + pm.math.dot(x_c, β)
    θ = pm.Deterministic('θ', pm.math.sigmoid(μ))
    bd = pm.Deterministic('bd', -α/β)
    yl = pm.Bernoulli('yl', p=θ, observed=y_0)
    trace_0 = pm.sample(1000)

为节省页数，同时避免您对同一类型图件反复出现感到厌烦，我将省略迹图和其他类似的总结图，但鼓励您制作自己的图和总结，以进一步探索本书中的例子。我们将直接跳到如何生成下图，这是一个数据曲线图，以及拟合的 sigmoid 曲线和决策边界：

theta = trace_0['θ'].mean(axis=0)
idx = np.argsort(x_c)
plt.plot(x_c[idx], theta[idx], color='C2', lw=3)
plt.vlines(trace_0['bd'].mean(), 0, 1, color='k')
bd_hpd = az.hpd(trace_0['bd'])
plt.fill_betweenx([0, 1], bd_hpd[0], bd_hpd[1], color='k', alpha=0.5)
plt.scatter(x_c, np.random.normal(y_0, 0.02),
            marker='.', color=[f'C{x}' for x in y_0])
az.plot_hpd(x_c, trace_0['θ'], color='C2')
plt.xlabel(x_n)
plt.ylabel('θ', rotation=0)
# use original scale for xticks
locs, _ = plt.xticks()
plt.xticks(locs, np.round(locs + x_0.mean(), 1))

前面这张图表示了花萼长度与花的种类（setosa = 0, versicolor=1）之间的关系。蓝色的S型曲线表示theta的均值，这条线可以解释为：在知道花萼长度的情况下花的种类是versicolor的概率，即半透明的蓝色区间是 95% HPD 区间。注意：在这种情况下，逻辑回归其实就是回归，因为我们做的就是在给定特征（或者特征的线性组合）的条件下，回归出一个数据点的类别为 1 的概率。不过需要时刻注意观察到的变量是二值变量而推断的是连续的概率，因此我们需要引入一个规则将连续的概率转换成一个 0-1 结果。决策边界在图中用红色来表示，此外还有一个 95% HPD 区间。根据决策边界， x 值位于其左侧的（这里对应花萼长度）属于类别 0（setosa），而位于其右侧的属于类别 1（versicolor）。决策边界在 x 轴上的取值为对应 y=0.5 时的点，可以证明其结果是 $-\frac{\alpha}{\beta}$ ，推导过程如下。

根据模型的定义，我们有如下关系：

$$\theta=\operatorname{logistic}(\alpha+x \beta)$$

根据逻辑函数的定义，当 $\theta=0.5$ 时，对应的输入为 0，则有：

$$0.5=\operatorname{logistic}\left(\alpha+x_{i} \beta\right) \Leftrightarrow 0=\alpha+x_{i} \beta$$

移项后可以得出，当 $ \theta=0.5 $ 时，对应有：

$$x_{i}=-\frac{\alpha}{\beta}$$

值得一提的是：

• 一般来说，$\theta$ 的价值是 $p(y=1|x)$ 。从这点上来说，Logistic回归是一种真正的回归；关键细节是，在给定特征的线性组合的情况下，我们回归的是数据点属于类别 1 的概率。
• 我们正在模拟一个二分变量的均值，即 [0-1] 区间内的一个数字。然后，我们引入一条规则，将这种概率转化为两类赋值。在这种情况下， 如果 $p(y=1)>=0.5$ ，我们分配类 1 ，否则分配类 0 。
• 这里选取的0.5并不是什么特殊值，你完全可以选其他0到1之间的值。只有当我们认为将标签0（setosa）错误地标为标签1（versicolor）时的代价，与反过来将标签1（versicolor）错误地标为标签0（setosa）时的代价相同时，选取0.5作为决策边界才是可行的。不过大多数情况下，分类出错的代价并不是对称的。

4.3 多元逻辑斯蒂回归

与多变量线性回归类似，多变量逻辑回归使用了多个自变量。这里将花萼长度与花萼宽度结合在一起，注意这里需要对数据做一些预处理。

df = iris.query("species == ('setosa', 'versicolor')")
y_1 = pd.Categorical(df['species']).codes
x_n = ['sepal_length', 'sepal_width']
x_1 = df[x_n].values

（1）决策边界

如果你对如何推导决策边界不感兴趣的话，可以略过这个部分直接跳到模型实现部分。 根据模型，我们有：

$$\theta=\operatorname{logistic}\left(\alpha+\beta_{1} x_{1}+\beta_{2} x_{2}\right)$$

从Logistic函数的定义出发，当Logistic回归的自变量为零时，有 $\theta=0.5 $ 。对应：

$$0.5=\operatorname{logistic}\left(\alpha+\beta_{1} x_{1}+\beta_{2} x_{2}\right) \Leftrightarrow 0=\alpha+\beta_{1} x_{1}+\beta_{2} x_{2}$$

通过移项，我们找到 $\theta=0.5$ 时 $x_2$ 的值：

$$x_{2}=-\frac{\alpha}{\beta_{2}}+\left(-\frac{\beta_{1}}{\beta_{2}} x_{1}\right)$$

这个决策边界的表达式与直线的表达式在数学形式上是一样的，其中第1项表示截距，第2项表示斜率，这里的括号只是为了表达上更清晰，如果你愿意的话完全可以去掉。为什么决策边界是直线呢？想想看，如果我们有一个特征，还有一维的数据，可以用一个点将数据分成两组；如果有两个特征，也就有一个2维的数据空间，从而我们可以用一条直线来对其分割；对于3维的情况，边界是一个平面，对于更高的维度，我们对应有一个超平面。事实上，从概念上讲，超平面可以大致定义为n维空间中n-1维的子空间，因此我们总是可以将决策边界称为超平面。

（2）模型实现

如果要用PyMC3写出多元逻辑回归模型，可以借助其向量化表示的优势，只需要对前面的单参数逻辑回归模型做一些简单的修改即可。

with pm.Model() as model_1:
    α = pm.Normal('α', mu=0, sd=10)
    β = pm.Normal('β', mu=0, sd=2, shape=len(x_n))
    μ = α + pm.math.dot(x_1, β)
    θ = pm.Deterministic('θ', 1 / (1 + pm.math.exp(-μ)))
    bd = pm.Deterministic('bd', -α/β[1] - β[0]/β[1] * x_1[:,0])
    yl = pm.Bernoulli('yl', p=θ, observed=y_1)
    trace_1 = pm.sample(2000)

正如对单个预测变量所做的那样，我们将绘制数据和决策边界：

idx = np.argsort(x_1[:,0])
bd = trace_1['bd'].mean(0)[idx]
plt.scatter(x_1[:,0], x_1[:,1], c=[f'C{x}' for x in y_0])
plt.plot(x_1[:,0][idx], bd, color='k');
az.plot_hpd(x_1[:,0], trace_1['bd'], color='k')
plt.xlabel(x_n[0])
plt.ylabel(x_n[1])

我们已经看到过的，决策边界现在是一条直线，不要被95%HPD区间的曲线给误导了。图中半透明的曲线是由于在中间部分多条直线绕中心区域旋转的结果（大致围绕 $x$ 的平均值 和 $y$ 的平均值）。

（3）关于Logistic回归系数的解释

在解释Logistic回归的 $\beta$ 系数时，我们必须小心。解释并不像我们在第3章“线性回归建模”中讨论的线性模型那样简单。使用Logistic逆连接函数引入了我们必须考虑的非线性。如果为 $\beta$ 为正，则增加 $x$ 会增加一些 $p(y=1)$ 的量，但该量值不是的 $x$ 线性函数；相反，它非线性地依赖于的 $x$ 值。我们可以在图中直观地看到这一事实；我们看到的不是具有恒定斜率的直线，而是具有随 $x$ 变化而不断调整斜率的 $S$ 形直线。一点代数知识可以让我们更深入地了解 $p(y=1)$ 是如何随 $\beta$ 变化的。基础模型是：

$$\theta=\operatorname{logistic}(\alpha+X \beta)$$

逻辑斯蒂的逆函数是 logit函数，它是：

$$\operatorname{logit}(z)=\log \left(\frac{z}{1-z}\right)$$

因此，如果我们取本节中的第一个方程式，并将logit函数应用于这两个项，我们会得到这个方程式：

$$\text{logit}(\theta)=\alpha+X \beta$$

或等价的：

$$\log \left(\frac{\theta}{1-\theta}\right)=\alpha+X \beta$$

记住模型中的 $ \theta $ 是 $ p(y=1) $ ：

$$\log \left(\frac{p(y=1)}{1-p(y=1)}\right)=\alpha+X \beta$$

$ \frac{p(y=1)}{1-p(y=1)} $ 的量值即为赔率。

成功的几率被定义为成功概率与不成功概率之比。虽然掷一个公平骰子得到 2 的概率是 1/6 ，但同一事件的几率是 $\frac{1/6}{5/6}=0.2$ 或1个有利事件对5个不利事件。赌徒经常使用赔率，主要是因为在考虑正确的下注方式时，赔率提供了一种比原始概率更直观的工具。

在Logistic回归中，系数 $\beta$ 编码了对数赔率单位随变量$x$ 单位增加而增加的情况。

从概率到赔率的转换是一个单调的转换，这意味着赔率随着概率的增加而增加，反之亦然。概率被限制在 [0，1] 区间，而赔率则在 [0，∞) 区间内。对数是另一个单调变换，对数赔率在 (-∞，∞) 区间内。下图显示了概率与赔率和对数赔率的关系：

probability = np.linspace(0.01, 1, 100)
odds = probability / (1 - probability)
_, ax1 = plt.subplots()
ax2 = ax1.twinx()
ax1.plot(probability, odds, 'C0')
ax2.plot(probability, np.log(odds), 'C1')
ax1.set_xlabel('probability')
ax1.set_ylabel('odds', color='C0')
ax2.set_ylabel('log-odds', color='C1')
ax1.grid(False)
ax2.grid(False)

因此，汇总数据提供的系数的值是对数赔率标度的：

df = az.summary(trace_1, var_names=varnames)

理解模型的一种非常实用的方法是更改参数，看看会发生什么。在下面的代码块中，我们正在计算支持杂色AS的对数赔率 \text{log_odds_versicolor_i}=\alpha+\beta_1x_1+\beta_2x_2 ，然后用Logistic函数计算杂色的概率。然后，我们通过固定 $x_2$并 $x_1$ 增加1来重复计算：

x_1 = 4.5  # sepal_length
x_2 = 3   # sepal_width
log_odds_versicolor_i = (df['mean'] * [1, x_1, x_2]).sum()
probability_versicolor_i = logistic(log_odds_versicolor_i)
log_odds_versicolor_f = (df['mean'] * [1, x_1 + 1, x_2]).sum()
probability_versicolor_f = logistic(log_odds_versicolor_f)
log_odds_versicolor_f - log_odds_versicolor_i, probability_versicolor_f -
probability_versicolor_i

如果运行代码，您会发现log-odds的增加值约为 4.66，这正是 $\beta_0$ (请查看 $trace\_1$ 的摘要)的值。这与我们之前的发现是一致的，即 $\beta$ 编码了对数赔率单位随变量 $x$ 单位增加而增加的情况。概率的增加约为 0.70 。

（4）处理相关变量

前一章中我们了解了在处理高度相关的变量时可能会遇到一些奇怪的问题。例如，如果用花瓣长度和花瓣宽度作为特征重跑前面的模型，会得到什么样的结果呢？

如果完成了上面的练习，你可能会注意到beta系数要比以前分布得更广，而且95%HPD区间（红色的区间）也更宽。下面的热力图显示了4个变量之间的相关性，可以看到（前一个例子中用到的）花萼长度与花萼宽度之间的相关性要远低于（后一个例子中用到的）花瓣长度与花瓣宽度之间的相关性。我们知道，相关的变量会得到更广的系数组合，因而能更好地解释数据，或者从另外一个角度来说，相关的变量限制模型的能力更弱。

当不同类别的数据能被完美区分时（即用线性模型分隔之后不同类的数据没有重叠），类似的问题就会出现。一个解决办法是不使用相关的变量，不过这个办法可能有时候并不适用。另外一个办法是给先验增加更多的信息，如果我们掌握了一些有用的信息，可以使用一些携带信息的先验，或者更一般地，可以使用一些弱信息的先验。Andrew Gelman和Stan的开发团队建议在进行逻辑回归时使用如下先验：

$ \beta \sim \operatorname{Student} T(0, \nu, s d) $

这里s的取值可以根据期望的尺度引入弱信息，正态参数v的值为3到7附近。该先验的含义是：我们期望系数比较小，同时引入了重尾，从而得到一个比高斯分布更鲁棒的模型。如果你忘记了的话，可以回忆一下前两章中的内容。