回归分析第15讲——连续正值型回归模型

在这个系列中，我们一开始学习的是一般的线性回归模型，它处理的是连续的实数值；到后来我们又学习了逻辑回归模型和泊松回归模型，使得我们可以处理离散的正值。最后还剩一种情况我们没有讨论，就是连续的正值，也就是说，响应变量只能取大于0的连续实数。要处理这一类响应变量，我们需要引入新的分布和模型。

注意，本讲我们只介绍相应的分布以及回归的拟合方法，对于模型的检验和诊断方法与之前的内容完全一致，这里不再反复叙述。

一、Beta回归

1.1 Beta分布

要讲清楚Beta回归，首先要了解Beta分布。

Beta分布要从二项分布开始说起。我们已经知道，如果随机变量$X$服从二项分布，即$X\sim B(n,\pi)$，那么其概率质量函数可以写作：

$$\Pr(X=x)=C_n^x\pi^x(1-\pi)^{n-x}$$

事实上，二项分布假定了$\pi$为一个定值，也就是说，二项分布是一个条件分布$\Pr(X=x|\pi)$。那么自然想到，对于参数$\pi$来说，在给定观测下也会存在波动，也满足一个分布。根据贝叶斯公式可以发现：

$$\Pr(X=x|\pi)=\frac{\Pr(\Pi=\pi|X=x)\Pr(X=x)}{\Pr(\pi)}$$

由于$\Pr(X=x)$和$\Pr(\pi)$都是常数，所以：

$$f(\Pi)=\Pr(\Pi=\pi|X=x)\varpropto\pi^a(1-\pi)^b$$

其中$a,b$均为常数，且$\pi\in(0,1)$。

既然存在正比关系，那么我们待定系数$k$，使得$f(\Pi)=k\pi^a(1-\pi)^b$。要使得$f(\Pi$)是个分布，那我们只需要令函数在$(0,1)$之间的积分为1即可：

$$\int_0^1k\pi^a(1-\pi)^bd\pi=1$$

不难求得$k=1/\int_0^1\pi^a(1-\pi)^bd\pi$，带入原式即可。

为了保持和教材的表达一致性，我们令$\alpha=a+1$，$\beta=b+1$，$x=\pi$，积分变量$dt=d\pi$，并且令Beta函数为：

$$B(\alpha,\beta)=\int_0^1t^{\alpha-1}(1-t)^{\beta-1}dt$$

就有分布函数：

$$f(X=x;\alpha,\beta)=\frac{1}{B(\alpha,\beta)}x^{\alpha-1}(1-x)^{\beta-1}$$

这就是Beta分布。为了帮助你理解Beta分布的性质，我们绘制不同$\alpha,\beta$取值下Beta分布的性状：

从图中你可以窥见Beta分布的性质：

• 自变量取值在$(0,1)$之间
• 连续的取值
• 形状多样（U型、倒U型、钟形、均匀分布、逼近正态），可以有偏态

因此，Beta分布可以用于描述0~1之间的连续值，例如转化率、点击率等。

最后，Beta分布的期望和方差分别为：

$$\begin{aligned} EX&=\frac{\alpha}{\alpha+\beta}\\ Var(X)&=\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)} \end{aligned}$$

1.2 Beta回归的拟合

如果响应变量$Y$服从Beta分布，那么就可以搭建Beta回归了。为了描述方便，我们令：

$$\mu=\frac{\alpha}{\alpha+\beta},\phi=\alpha+\beta$$

此时有：

$$Y\sim Beta(\mu\phi,(1-\mu)\phi)$$

我们称$\mu$为均值参数，$\phi$为精度参数。

接下来我们需要选择联系函数。Beta回归选择的联系函数与二项分布相关回归一样，最常见是取逻辑函数，即：

$$g(\mu)=\log\frac{\mu}{1-\mu}$$

当然，你也可以选择probit连接、log-log连接等。

接下来我们以逻辑连接为例，展示参数估计过程。为了方便叙述，我们对Beta函数进行等价转换，即：

$$B(\alpha,\beta)=\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}$$

其中$\Gamma()$是Gamma函数。因此Beta分布可以重写为：

$$f(y;\alpha,\beta)=\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}y^{\alpha-1}(1-y)^{\beta-1}=\frac{\Gamma(\phi)}{\Gamma(\mu\phi)\Gamma[(1-\mu)\phi]}y^{\mu\phi-1}(1-y)^{(1-\mu)\phi-1}$$

你或许对Gamma函数较为陌生，这里我们做一点说明。

Gamma函数是为了把阶乘的数域从整数扩展为实数（甚至复数），如果$n$是正整数，那么有：

$$\Gamma(n)=(n-1)!$$

Gamma函数的定义为：

$$\Gamma(\alpha)=\int_0^\infty x^{\alpha-1}e^{-x}dx,\alpha>0$$

这样定义以后，你就能对实数进行阶乘了。

Gamma函数有如下性质：

• $\Gamma(1)=1$
• $\Gamma(1/2)=\sqrt\pi$
• $\Gamma(\alpha+1)=\alpha\Gamma(\alpha)$

由于Beta分布不属于典型的指数分布族，因此我们在拟合Beta回归时不能完全按照之前的思路进行，需要做一点改动。

极大似然法是一定能用的。首先不难求得对数似然函数为：

$$l_i=\log\Gamma(\phi)-\log\Gamma(\mu_i\phi)-\log\Gamma[(1-\mu_i)\phi]+(\mu_i\phi-1)\log y_i+[(1-\mu_i)\phi-1]\log(1-y_i)$$

由于我们要估计两个参数，即线性预测子的$\beta$和精度参数$\phi$，所以要分别求导，然后联立方程：

$$\begin{aligned} U_{\beta_j}&=\phi\left[\psi[(1-\mu_i)\phi]-\psi(\mu_i\phi)+\log\frac{y_i}{1-y_i}\right]\mu_i(1-\mu_i)x_{ij}\\ U_\phi&=\sum_{i=1}^n\left[\psi(\phi)-\mu_i\psi(\mu_i\phi)-(1-\mu_i)\psi[(1-\mu_i)\phi]+\mu_i\log y_i+(1-\mu_i)\log(1-y_i)\right] \end{aligned}$$

其中$\psi()$称为双Gamma函数，它表示对数Gamma函数的一阶导数，即$(\log\Gamma())^{'}$。

其次是Fisher信息阵。由于有两个参数，所以信息阵需要写成块的形式：

$$\mathfrak I(\beta,\phi)= \begin{pmatrix} \mathfrak I_{\beta\beta} & \mathfrak I_{\beta\phi}\\ \mathfrak I_{\phi\beta} & \mathfrak I_{\phi\phi} \end{pmatrix}$$

可以求得：

$$\begin{aligned} &\mathfrak I_{\beta\beta}=\sum_{i=1}^n\phi^2[\psi^{'}(\mu_i\phi)+\psi^{'}[(1-\mu_i)\phi]]\mu_i^2(1-\mu_i)^2x_ix_i^T\\ &\mathfrak I_{\phi\phi}=\sum_{i=1}^n\left[-\psi^{'}(\phi)+\mu_i^2\psi^{'}(\mu_i\phi)+(1-\mu_i)^2\psi^{'}[(1-\mu_i)\phi]\right]\\ &\mathfrak I_{\beta\phi}=-\sum_{i=1}^n\phi[(1-\mu_i)\psi^{'}[(1-\mu_i)\phi]-\mu_i\psi^{'}(\mu_i\phi)]\mu_i(1-\mu_i)x_i\\ &\mathfrak I_{\phi\beta}=\mathfrak I_{\beta\phi}^T \end{aligned}$$

其中$\psi^{'}()$是双Gamma函数的导数，又可以称为三Gamma函数。

如果你还想用IRWLS算法来拟合Beta回归，也是可以的，但有一定的差别。为了叙述方便，我们令

$$\begin{aligned} &r_i=\psi[(1-\mu_i)\phi]-\psi(\mu_i\phi)+\log\frac{y_i}{1-y_i}\\ &w_i=\psi^{'}(\mu_i\phi)+\psi^{'}[(1-\mu_i)\phi] \end{aligned}$$

工作变量为：

$$z_i^{(t)}=x_i^T\beta^{(t)}+\frac{r_i}{w_i\mu_i(1-\mu_i)}$$

权重为：

$$W_{ii}^{(t)}=\phi^2w_i\mu_i^2(1-\mu_i)^2$$

于是回归系数的最小二乘解为：

$$\beta^{(t+1)}=(X^TW^{(t)}X)^-X^TW^{(t)}z^{(t)}$$

得到回归系数后，再带入更新$\phi$：

$$\phi^{(t+1)}=\phi^{(t)}+\frac{U_\phi(\beta^{(t+1)},\phi^{(t)})}{\mathfrak I_{\phi\phi}(\beta^{(t+1)},\phi^{(t)})}$$

二、指数回归

2.1 指数分布

指数分布要从泊松分布开始讲起。我们已经知道，泊松分布描述的是单位时间或空间内事件发生的次数，期望是$\lambda$。那么容易想到，$1/\lambda$就应该表示事件发生一次需要的时间，即事件发生的时间间隔。时间间隔的长短是否也服从概率分布？

我们尝试从泊松分布出发来导出分布。我们可以把事件发生的时间间隔划分为$t$个单位时间，那么事件经过$t$个单位时间才发生意味着事件的发生时间不小于t，并且每个时间单位里事件的发生次数为0，于是我们求出概率：

$$\Pr(T>t)=\left[Poisson(X=0)\right]^t=e^{-\lambda t}$$

有了这个概率，很自然就能求出事件在$t$个时间单位内发生的概率，即：

$$\Pr(T\le t)=1-\Pr(T>t)=1-e^{-\lambda t}$$

这是一个累计概率密度函数，要求分布只需对$t$求导：

$$f(t)=\frac{d\Pr(T\le t)}{dt}=\lambda e^{-\lambda t}$$

这就是指数分布，指数分布描述了事件下一次发生的时间，在日常生活中非常常见，例如灯泡的寿命、等待下一辆公交车的时间。

一般地，我们认为指数分布的形式是：

$$f(y,\lambda)= \begin{cases} \lambda e^{-\lambda y},&y\ge0\\ 0,&y<0 \end{cases}$$

指数分布的期望和方差分别是：

$$\begin{aligned} &EY=\frac{1}{\lambda}=\mu\\ &Var(Y)=\frac{1}{\lambda^2}=\mu^2 \end{aligned}$$

指数分布肯定属于指数分布族，所以可以写成指数分布族的形式：

$$exp\left\{-\frac{y}{\mu}+\log\frac{1}{\mu}\right\}$$

2.2 指数回归的拟合

如果响应变量服从指数分布，那么可以考虑指数回归。指数回归的联系函数是：

$$g(\mu)=-\frac{1}{\mu}$$

所以不难求出得分统计量和Fisher信息阵：

$$\begin{aligned} &U_j=\sum_{i=1}^n(y_i-\mu_i)x_{ij}\\ &\mathfrak I_{jk}=\sum_{i=1}^n\mu_i^2x_{ij}x_{ik} \end{aligned}$$

也可以直接使用IRWLS算法，工作变量为：

$$z_i=x_i^T\beta^{(t)}+\frac{y_i-\mu_i}{\mu_i^2},\mu_i=-\frac{1}{x_i^T\beta^{(t)}}$$

权重为：

$$w_i=\mu_i^2$$

更新：

$$\beta^{(t+1)}=(X^TWX)^-X^TWz$$

三、Gamma回归

3.1 Gamma分布

上一节我们介绍了指数分布，指数分布描述的是事件下一次发生的等待时间，是通过泊松分布导出的。而进一步的一个问题是，我们能不能描述事件第$k$次发生所等待的时间？

要做到这一点，我们还是从泊松分布来导出。泊松分布描述的是单位时间内事件发生的次数，即：

$$Poisson(y=k,t=1)=\frac{\lambda^k}{k!}e^{-\lambda}$$

如果要描述$t$个单位时间内事件发生的次数，只需要简单地将$\lambda$换成$\lambda t$即可：

$$Poisson(y=k,t)=\frac{(\lambda t)^k}{k!}e^{-\lambda t}$$

此时期望为$E=\lambda t$。

类似于指数分布的推导，我们也可以求出事件发生$k$次的等待时间不少于$t$的概率，即$\Pr(T>t,k)$。要求出这个概率，我们需要把问题等价转换成在$t$个单位时间内事件最多发生$k-1$次，这样就能保证事件发生$k$次时等待时间大于$t$了。既然是至多发生$k-1$次，那么只需要把事件发生0次到$k-1$次的概率相加即可：

$$\Pr(T>t,k)=\sum_{i=0}^{k-1}\frac{(\lambda t)^i}{i!}e^{-\lambda t}$$

于是，$t$个单位时间内事件发生$k$次的累计概率为：

$$\Pr(T\le t,k)=1-\sum_{i=0}^{k-1}\frac{(\lambda t)^i}{i!}e^{-\lambda t}$$

对$t$求导就得到了概率密度：

$$f(t)=\frac{d\Pr(T\le t,k)}{dt}=\frac{\lambda^ke^{-\lambda t}t^{k-1}}{\Gamma(k)}$$

这就是Gamma分布，它描述了事件发生$k$次的等待时间。

在Gamma分布中，$k$是事件发生的次数，又被称为形状参数；$\lambda$是单位时间内事件发生的平均次数，又被称为速率参数。下面展示两张图，你可以从图中直观感受形状参数和速率参数对分布的影响：

Gamma分布的期望和方差分别为：

$$\begin{aligned} &ET=\frac{k}{\lambda}=\mu\\ &Var(T)=\frac{k}{\lambda^2}=\mu^2\phi,\ \ \ \phi=\frac{1}{k} \end{aligned}$$

Gamma分布的期望等于事件发生的次数除以单位时间内事件发生的次数，得到的就是等待时间，符合实际情况。并且，Gamma分布在$k=1$时就退化为指数分布。

3.2 Gamma回归的拟合

如果响应变量服从Gamma分布，那么可以考虑使用Gamma回归。为了获得联系函数，我们把Gamma分布写成指数分布族的形式：

$$\begin{aligned} f(y,\mu,\phi)&=\frac{(\mu\phi)^{-1/\phi}e^{-y/\mu\phi}y^{1/\phi-1}}{\Gamma(1/\phi)}\\ &=exp\left\{\frac{-y}{\mu\phi}-\frac{\log(\mu\phi)}{\phi}+(\frac{1}{\phi}-1)\log y-\log\Gamma(\frac{1}{\phi})\right\}\\ &=exp\left\{\frac{y(1/\mu)-(1-\log\mu)}{-\phi}+\frac{1-\phi}{\phi}\log y-\frac{\log\phi}{\phi}-\log\Gamma(\frac{1}{\phi})\right\} \end{aligned}$$

不难看出，$\theta=1/\mu$，$b(\theta)=-\log\mu$，分散参数$a(\phi)=\phi$。

Gamma回归的联系函数可以取：

$$g(\mu)=\frac{1}{\mu}$$

Gamma回归还可以取对数连接（类比泊松回归）、恒等连接（类比一般线性回归）等。

于是我们可以计算得分统计量：

$$U_j=\sum_{i=1}^n\frac{y_i-\mu_i}{-\phi}x_{ij}$$

Fisher信息阵为：

$$\mathfrak I_{jk}=\sum_{i=1}^n\frac{\mu_i^2x_{ij}x_{ik}}{\phi}$$

如果要使用IRWLS算法，那么先求出工作变量：

$$z_i=x_i^T\beta^{(t)}-\frac{y_i-\mu_i}{\mu_i^2},\mu_i^2=\frac{1}{x_i^T\beta^{(t)}}$$

再求权重：

$$w_i=-\mu_i^2/\phi$$

最后更新：

$$\beta^{(t+1)}=(X^TWX)^-X^TWz$$

你或许注意到了这里有一个分散参数$\phi$，我们似乎还没有估计它。有什么办法可以估计？

如果你仔细比较指数分布和Gamma分布后你就会发现，两个分布的期望都是$\mu$，方差都含有$\mu^2$，只不过Gamma分布多出来一个因子$\phi$。你一定想到了，这就是我们在讨论过度分散问题时提到的放松方差结构的方法。类似于泊松分布和quasi泊松，Gamma分布也可以看成quasi版的指数分布。当$\phi=1$，即$k=1$时，Gamma分布退化为指数分布，就没有过度分散了。

因此，我们可以先假定$k=1$，也就是使用指数回归来拟合模型，之后再通过分散参数的估计方法：

$$\phi=\frac{Pearson \ \chi^2}{n-p}$$

这样就能得到分散参数了。分散参数不为1的指数回归就是Gamma回归，指数回归是Gamma回归的特例。

四、逆高斯回归

4.1 逆高斯分布

在连续正值分布的最后，我们再学习一个逆高斯分布。这里我们不推导逆高斯分布，直接给出表达式：

$$f(y,\mu,\lambda)=\left(\frac{\lambda}{2\pi y^3}\right)^{\frac{1}{2}}exp\left\{-\frac{\lambda(y-\mu)^2}{2\mu^2 y}\right\}$$

其中$\mu$为均值参数，$\lambda>0$为形状参数。当$\lambda\rightarrow+\infty$时，逆高斯分布就接近正态分布（高斯分布）。

逆高斯分布的特点是：极右偏非负分布；具有异方差性；小值数据偏多、大值数据极少而造成严重拖尾。这一点在生存分析中非常重要，已经成为传统生存分析模型的一个优化升级。之后我们在学习生存分析时还会分析逆高斯分布的作用。

为了展示逆高斯分布的特点，我们绘制不同参数组合下的逆高斯分布：

逆高斯分布的均值和方差分别为：

$$\begin{aligned} &EY=\mu\\ &Var(Y)=\frac{\mu^3}{\lambda} \end{aligned}$$

4.2 逆高斯回归的拟合

如果响应变量服从逆高斯分布，那么我们可以考虑使用逆高斯回归。在GLM框架下，我们一般取$\sigma^2=1/\lambda$，所以逆高斯分布可以重写为：

$$f(y,\mu,\sigma^2)=\frac{1}{\sqrt{2\pi y^3\sigma^2}}exr\left\{-\frac{(y-\mu)^2}{2(\mu\sigma)^2y}\right\}$$

这样方差就可以写作$\mu^3\sigma^2$。

逆高斯分布属于指数分布族，可以写成指数分布族的形式：

$$exp\left\{\frac{y/(2\mu^2)-1/\mu}{-\sigma^2}-\frac{1}{2y\sigma^2}-\frac{1}{2}\log(2\pi y^3\sigma^2)\right\}$$

不难看出，$\theta=1/(2\mu^2)$，$b(\theta)=1/\mu$，分散参数$a(\theta)=-\sigma^2$。

逆高斯回归的联系函数可取：

$$g(\mu)=\frac{1}{2\mu^2}$$

得分统计量为：

$$U_j=-\sum_{i=1}^n\frac{y_i-\mu_i}{\sigma^2}x_{ij}$$

如果你尝试对一般线性回归（高斯分布）计算得分统计量，你会发现二者仅仅相差前面的负号。

Fisher信息阵为：

$$\mathfrak I_{jk}=\sum_{i=1}^n\frac{\mu_i^3x_{ij}x_{ik}}{\sigma^2}$$

同理，如果要使用IRWLS算法，需要先求出工作变量：

$$z_i=x_i^T\beta^{(t)}-\frac{y_i-\mu_i}{\mu_i^3},\mu_i=\frac{1}{\sqrt{2x_i^T\beta^{(t)}}}$$

再求出权重：

$$w_i=\frac{\mu_i^3}{\sigma^2}$$

最后更新：

$$\beta^{(t+1)}=(X^TWX)^-X^TWz$$

逆高斯回归还有一个形状参数$\sigma^2$需要估计。形状参数的估计策略是，先给定初始值，更新$\beta^{(t+1)}$，然后更新$\sigma^2$，计算方法为：

$$\sigma^2(\mu_i)=\frac{\sum_{i=1}^n\frac{(y_i-\mu_i)^2}{\mu_i^2 y_i}}{n},\mu_i=\frac{1}{\sqrt{2x_i^T\beta^{(t+1)}}}$$

重复直至收敛即可。