NO.29.tip： 广义线性模型

背景

广义线性模型是指数族分布的条件版本。更准确地说，该模型具有以下形式：

$$p(y_n\mid x_n,w,{\sigma }^2)=\exp \left[\frac{y_n{\eta }_n-A({\eta }_n)}{{\sigma }^2}+\log h(y_n,{\sigma }^2)\right]\text{\hspace{1em}(1)}$$

其中，${\eta }_n=w^{\top }{x}_n$ 是分布的自然参数，$A({\eta }_n)$ 是对数归一化器，$T(y)=y$ 是充分统计量，${\sigma }^2$ 是离散项。我们可以得出响应变量的均值和方差如下：

$${\mu }_n\triangleq \mathbb{E}[y_n\mid x_n,w,{\sigma }^2]=A^{\prime }({\eta }_n)\triangleq {\ell }^{-1}({\eta }_n)\text{\hspace{1em}(2)}$$

$$\mathbb{V}[y_n\mid x_n,w,{\sigma }^2]=A^{\prime \prime }({\eta }_n){\sigma }^2\text{\hspace{1em}(3)}$$

我们将使用 ${\mu }_n={\ell }^{-1}({\eta }_n)$ 表示从线性输入到输出均值的映射，其中函数 $\ell$ 称为链接函数(link function)，${\ell }^{-1}$ 称为均值函数(mean function)。这种关系通常记作：

$$\ell ({\mu }_n)={\eta }_n=w^{\top }{x}_n\text{\hspace{1em}(4)}$$

由于线性的假设，广义线性模型的预测能力非常有限(尽管我们总是可以在 $x_n$ 上使用基函数展开来提高灵活性)。然而，广义线性模型在统计学文献中的主要用途并不是预测，而是假设检验。这取决于计算后验 $p(w\mid \mathcal{D})$ 的能力。可以利用这一点得出结论，即判断任何输入(例如，代表不同的群组)是否对输出有显著影响。

源与流

一些流行的广义线性模型

我们将给出一些广泛使用的广义线性模型示例。

线性回归

回想一下，线性回归具有以下形式：

$$p(y_n\mid x_n,w,{\sigma }^2)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{1}{2{\sigma }^2}(y_n-w^{\top }{x}_n{)}^2\right)\text{\hspace{1em}(5)}$$

因此，

$$\log p(y_n\mid x_n,w,{\sigma }^2)=-\frac{1}{2{\sigma }^2}(y_n-{\eta }_n{)}^2-\frac{1}{2}\log (2\pi {\sigma }^2)\text{\hspace{1em}(6)}$$

其中，${\eta }_n=w^{\top }{x}_n$。我们可以将其记作广义线性模型的形式，如下所示：

$$\log p(y_n\mid x_n,w,{\sigma }^2)=\frac{y_n{\eta }_n-\frac{{\eta }_n^2}{2}}{{\sigma }^2}-\frac{1}{2}\left(\frac{y_n^2}{{\sigma }^2}+\log (2\pi {\sigma }^2)\right)\text{\hspace{1em}(7)}$$

我们发现 $A({\eta }_n)={\eta }_n^2/2$，因此，

$$\mathbb{E}[y_n]={\eta }_n=w^{\top }{x}_n\text{\hspace{1em}(8)}$$

$$\mathbb{V}[y_n]={\sigma }^2\text{\hspace{1em}(9)}$$

二项式回归

如果响应变量是 $N_n$ 次试验中的成功次数，即 $y_n\in \{0,\cdots ,N_n\}$，那么我们可以使用二项式回归(binomial regression)，其定义为：

$$p(y_n\mid x_n,N_n,w)=\text{Bin}(y_n\mid \sigma (w^{\top }{x}_n),N_n)\text{\hspace{1em}(10)}$$

我们发现二元逻辑回归是 $N_n=1$ 时的一种特殊情况。

其对数概率密度函数由下式给出：

$$\log p(y_n\mid x_n,N_n,w)=y_n\log {\mu }_n+(N_n-y_n)\log (1-{\mu }_n)+\log \left(\genfrac{}{}{0ex}{}{N_n}{y_n}\right)\text{\hspace{1em}(11)}$$

$$=y_n\log \left(\frac{{\mu }_n}{1-{\mu }_n}\right)+N_n\log (1-{\mu }_n)+\log \left(\genfrac{}{}{0ex}{}{N_n}{y_n}\right)\text{\hspace{1em}(12)}$$

其中，${\mu }_n=\sigma ({\eta }_n)$。为了将其记作广义线性模型形式，让我们定义：

$${\eta }_n\triangleq \log \left[\frac{{\mu }_n}{(1-{\mu }_n)}\right]=\log \left[\frac{1}{1+e^{-w^{\top }{x}_n}}\frac{1+e^{-w^{\top }{x}_n}}{e^{-w^{\top }{x}_n}}\right]=\log \frac{1}{e^{-w^{\top }{x}_n}}=w^{\top }{x}_n\text{\hspace{1em}(13)}$$

因此，我们可以将二项式回归写成广义线性模型形式，如下所示：

$$\log p(y_n\mid x_n,N_n,w)=y_n{\eta }_n-A({\eta }_n)+h(y_n)\text{\hspace{1em}(14)}$$

其中，$h(y_n)=\log \left(\genfrac{}{}{0ex}{}{N_n}{y_n}\right)$，并且有

$$A({\eta }_n)=-N_n\log (1-{\mu }_n)=N_n\log (1+e^{{\eta }_n})\text{\hspace{1em}(15)}$$

因此，

$$\mathbb{E}[y_n]=\frac{\mathrm{d}A}{\mathrm{d}{\eta }_n}=\frac{N_n{e}^{{\eta }_n}}{1+e^{{\eta }_n}}=\frac{N_n}{1+e^{-{\eta }_n}}=N_n{\mu }_n\text{\hspace{1em}(16)}$$

并且有

$$\mathbb{V}[y_n]=\frac{{\mathrm{d}}^{2}A}{\mathrm{d}{\eta }_n^2}=N_n{\mu }_n(1-{\mu }_n)\text{\hspace{1em}(17)}$$

泊松回归

如果响应变量是一个整数计数，即 $y_n\in \{0,1,\cdots \}$，那么我们可以使用泊松回归(Poisson regression)，其定义为：

$$p(y_n\mid x_n,w)=\text{Poi}(y_n\mid \exp (w^{\top }{x}_n))\text{\hspace{1em}(18)}$$

其中，

$$\text{Poi}(y\mid \mu )=e^{-\mu }\frac{{\mu }^y}{y!}\text{\hspace{1em}(19)}$$

是泊松分布。泊松回归广泛用于生物统计应用，其中 $y_n$ 可以表示给定个人或地点的疾病数量，或高通量测序环境中基因组位置的读取次数。

其对数概率密度函数由下式给出：

$$\log p(y_n\mid x_n,w)=y_n\log {\mu }_n-{\mu }_n-\log (y_n!)\text{\hspace{1em}(20)}$$

其中，${\mu }_n=\exp (w^{\top }{x}_n)$。因此，在广义线性模型形式中，我们有

$$\log p(y_n\mid x_n,w)=y_n{\eta }_n-A({\eta }_n)+h(y_n)\text{\hspace{1em}(21)}$$

其中，${\eta }_n=\log ({\mu }_n)=w^{\top }{x}_n$，$A({\eta }_n)={\mu }_n=e^{{\eta }_n}$，$h(y_n)=-\log (y_n!)$。因此

$$\mathbb{E}[y_n]=\frac{\mathrm{d}A}{\mathrm{d}{\eta }_n}=e^{{\eta }_n}={\mu }_n\text{\hspace{1em}(22)}$$

并且有

$$\mathbb{V}[y_n]=\frac{{\mathrm{d}}^{2}A}{\mathrm{d}{\eta }_n^2}=e^{{\eta }_n}={\mu }_n\text{\hspace{1em}(23)}$$

零膨胀泊松回归

在许多形式的计数数据中，即使考虑了预测因素，观察到的0的数量也比模型预测的要多。直观上来说，这是因为可能有很多方法不会产生输出结果。例如，考虑预测产品的销售数据。如果销售额是0，这意味着该产品不受欢迎(因此需求也非常低)，还是只是因为该产品销售一空(意味着需求很高，但供应量为零)？类似的问题也出现在基因组学、流行病学等领域。

为了处理这种情况，通常使用零膨胀泊松(Zero-Inflated Poisson, ZIP)模型。该模型的似然是两种分布的混合：0处的尖峰和标准泊松分布。形式上，我们定义

$$\text{ZIP}(y\mid \rho ,\lambda )=\{\begin{array}{llc}\rho +(1-\rho )\exp (-\lambda )& \text{如果 }y=0& \\ (1-\rho )\frac{{\lambda }^y\exp (-\lambda )}{y!}& \text{如果 }y& gt;0\end{array}\text{\hspace{1em}(24)}$$

其中，$\rho$ 是拾取尖峰的先验概率，$\lambda$ 是泊松速率。我们发现，生成0有两种"机制"：一种方法(以概率 $\rho$)是选择尖峰，另一种方法(以概率 $1-\rho$)则仅仅因为泊松速率很低而生成零计数。(后一个事件的概率为 ${\lambda }^0{e}^{-\lambda }/0!=e^{-\lambda }$。)

具有非规范链接函数的广义线性模型

我们已经看到输出分布的均值参数是如何由 $\mu ={\ell }^{-1}(\eta )$ 给出的，其中函数 $\ell$ 是链接函数。链接函数可以有若干种选择，接下来我们将展开讨论。

规范链接函数(canonical link function) $\ell$ 满足 $\theta =\ell (\mu )$ 的性质，其中 $\theta$ 是规范(自然)参数。因此

$$\theta =\ell (\mu )=\ell ({\ell }^{-1}(\eta ))=\eta \text{\hspace{1em}(25)}$$

这就是我们迄今为止的假设。例如，对于伯努利分布，其规范参数是对数几率 $\eta =\log (\mu /(1-\mu ))$，由logit变换给出：

$$\eta =\ell (\mu )=\text{logit}(\mu )=\log \left(\frac{\mu }{1-\mu }\right)\text{\hspace{1em}(26)}$$

其逆函数为sigmoid函数或逻辑函数：

$$\mu ={\ell }^{-1}(\eta )=\sigma (\eta )=1/(1+e^{-\eta })\text{\hspace{1em}(27)}$$

然而，我们可以自由使用其他类型的链接函数。例如，我们使用

$$\eta =\ell (\mu )={\Phi }^{-1}(\mu )\text{\hspace{1em}(28)}$$

$$\mu ={\ell }^{-1}(\eta )=\Phi (\eta )\text{\hspace{1em}(29)}$$

这被称为probit链接函数。

另一个链接函数是互补对数-对数函数(complementary log-log)，它有时用于二元响应的情况：

$$\eta =\ell (\mu )=\log (-\log (1-\mu ))\text{\hspace{1em}(30)}$$

这可以用于观察到0个事件(用 $y=0$ 表示)或一个或多个事件(由 $y=1$ 表示)的应用中，其中假设事件由速率为 $\lambda$ 的泊松分布控制。设 $E$ 为事件数。泊松假设意味着 $p(E=0)=\exp (-\lambda )$，因此

$$p(y=0)=(1-\mu )=p(E=0)=\exp (-\lambda )\text{\hspace{1em}(31)}$$

因此，$\lambda =-\log (1-\mu )$。当 $\lambda$ 是协变量函数时，需要确保函数的值是正数，所以我们使用 $\lambda =e^{\eta }$，因此有

$$\eta =\log (\lambda )=\log (-\log (1-\mu ))\text{\hspace{1em}(32)}$$

最大似然估计

广义线性模型可以使用与我们拟合逻辑回归相似的方法来拟合。特别是，负对数似然具有以下形式(忽略常数项)：

$$\text{NLL}(w)=-\log p(\mathcal{D}\mid w)=-\frac{1}{{\sigma }^2}\sum _{n=1}^N{\ell }_n\text{\hspace{1em}(33)}$$

其中，

$${\ell }_n\triangleq {\eta }_n{y}_n-A({\eta }_n)\text{\hspace{1em}(34)}$$

其中，${\mu }_n=w^{\top }{x}_n$。为了简单起见，我们假设 ${\sigma }^2=1$。

我们可以按如下方式计算单个项的梯度：

$$g_n\triangleq \frac{\partial {\ell }_n}{\partial w}=\frac{\partial {\ell }_n}{\partial {\eta }_n}\frac{\partial {\eta }_n}{\partial w}=(y_n-A^{\prime }({\eta }_n))x_n=(y_n-{\mu }_n)x_n\text{\hspace{1em}(35)}$$

其中，${\eta }_n=f(w^{\top }{x}_n)$，并且 $f$ 是从规范参数映射到均值参数的逆链接函数。(例如，在逻辑回归的情况下，我们有 ${\mu }_n=\sigma (w^{\top }x)$。)

海森矩阵由下式给出：

$$H=\frac{{\partial }^2}{\partial w\partial w^{\top }}\text{NLL}(w)=-\sum _{n=1}^N\frac{\partial g_n}{\partial w^{\top }}\text{\hspace{1em}(36)}$$

其中，

$$\frac{\partial g_n}{\partial w^{\top }}=\frac{\partial g_n}{\partial {\mu }_n}\frac{\partial {\mu }_n}{\partial w^{\top }}=-x_n{f}^{\prime }(w^{\top }{x}_n)x_n^{\top }\text{\hspace{1em}(37)}$$

因此，

$$H=\sum _{n=1}^N{f}^{\prime }({\eta }_n)x_n{x}_n^{\top }\text{\hspace{1em}(38)}$$

例如，在逻辑回归的情况下，$f({\eta }_n)=\sigma ({\eta }_n)={\mu }_n$，$f^{\prime }({\eta }_n)={\mu }_n(1-{\mu }_n)$。一般而言，我们发现海森矩阵是正定的，因为 $f^{\prime }({\eta }_n)>0$；因此负对数似然是凸函数，广义线性模型的最大似然估计是唯一的(假设对所有的 $n$，$f({\eta }_n)>0$)。

对于较小的数据集，我们可以使用迭代重新加权最小二乘法(Iteratively Reweighted Least Squares, IRLS)算法(这是牛顿法的一种形式)来计算最大似然估计。对于较大的数据集，可以使用随机梯度下降。(在实践中，将随机梯度下降与自动调整步长的方法相结合通常非常有效。)

贝叶斯推理

最大似然估计提供了参数的点估计，但不传达任何不确定性的概念，这对于假设检验以及避免过拟合都很重要。为了计算不确定性，我们将对参数进行贝叶斯推理。为了做到这一点，我们首先需要指定一个先验。选择合适

的先验取决于链接函数的形式。例如，如果将 $\alpha$ 传递给sigmoid函数，偏移项 $\alpha \in \mathbb{R}$ 上的"平坦"或"无信息性"先验在概率尺度上就不会转化为无信息先验。

一旦我们选择了先验，就可以使用各种近似推理方法来计算后验。对于小型数据集，最简单的方法是哈密尔顿蒙特卡罗，因为只需要记下对数似然和对数先验，然后可以使用自动梯度来计算可以传递给哈密尔顿蒙特卡罗引擎的导数。

对于大型数据集，哈密尔顿蒙特卡罗方法可能很慢，因为该方法是一种全批次处理算法。在这样的设置中，变分贝叶斯、期望传播或者更专业的算法是最佳选择。