一般线性模型与广义线性模型
由于统计模型的多样性和各种模型的适应性，针对因变量（也就是我们常用于回归方程中的y）和解释变量（就是自变量，也就是我们常用于回归的x）的取值性质，可以将统计模型分为多种类型。通常自变量为定性变量的线性模型称为一般线性模型（general linear model，GLM），如实验设计模型、方差分析模型。因变量为非正态分布的线性模型为广义线性模型（Generalized linear model,GENMOD），如Logistic回归模型、对数线性模型和Cox比例风险模型。
对于一般线性线性模型，其基本假定是y服从正态分布，或至少y的方差$\sigma^2$为有限常数。但是在实际研究中，情况并非总是如此，例如当y是发病率(y=k/n)时，y服从二项分布，期望值和方差分别为E(y)=π，Var(y)=1/n x π(1-π)，其方差与倒数呈反比，且是π的函数。又如，当y是单位时间内的放射性计数时，y服从Poisson分布，期望值和方差分别为E(y)=μ,Var(y)=μ，方差是μ的函数。实际数据中有很多资料均不符合一般线性模型的基本假定。尽管也可以将频率或频数作为y代入一般线性模型，但拟合结果往往不能令人满意，例如会出现频率的拟合值y>1、频数的拟合值y<0这种奇怪现象。
后来有人在一般线性模型的基础上，对$\sigma^2$为有限常数的假定作了进一步的推广，提出了广义线性模型的概念和拟似然函数(quasi-likelihood function)的方法，用于求解满足下列条件的线性模型：
$E(y) = \mu $
$m(\mu) = X\beta$
$cov(y) = \sigma^2 V(\mu) \quad （一）$
其中m为连接函数组成的向量，将$\mu$转化为$\beta$的线性表达式，$V(\mu)$为$n \times n$矩阵，其中每个元素均为$\mu$的函数，当处$y_{i}$值相互独立时，$V(\mu)$为对角矩阵。当$m(\mu)=\mu$，$V(\mu)=I$时，(一)中的三个公式为一般线性模型，也就是说，（一）中包括了一般线性模型。在广义线性模型中，均假定观察值y具有指数族概率密度函数，如下所示：
$f(y|\theta,\phi) = exp{[y\theta -b(\theta)]/a(\phi)=c(y,\phi)} \quad (二)$
其中，a、b和c是三种函数形式，θ为典则参数，如果给定φ（散布参数，有时写作$\sigma^2$），（二）式就是具有参数θ的指数族密度函数，以正式分布为例：
与（二）式对照，如下所示：
$\theta=\mu$
$b(\theta)=\mu^2/2$
$\phi=\sigma^2$
$a(\phi)=\sigma^2$
$c(y,\phi)=-\frac{1}{2}[y^2/\sigma^2 + ln(2\pi\sigma^2) ]$
根据样本和y的函数可建立对数似然函数，并可导出y的期望值和方差。在广义线性模型中，（二）式中的典则参数不仅仅是μ的函数，还是参数β0，β1,..βp的线性表达式，因此对μ作变换，则可得到下面三种分布连接函数的形式：
正态分布：$m(\mu)=\mu=\sum\beta_{j}x_{j}$
二项分布：$m(\mu)=ln(\frac{\mu}{1-\mu})=\sum\beta_{j}x_{j}$
Poisson分布：$mu(\mu)=ln(\mu)=\sum\beta_{j}x_{j}$
Logistic属于广义线性模式的一种，它是通常的正态线性模型的推广，它要求应变量只能通过线性形式依赖于自变量。上述推广体现在以下两个方面：
通过一个连接函数，将应变量与自变量建立线性关系
$m(E(y))=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+\cdots+\beta_{p}x_{p}$
通过一个误差函数，说明广义线性模型的最后一部分随机项。
因此，Logistic是关于应变量为0到1定性变量的广义线性回归问题，且广义线性模型的分布族为二项分布。
广义线性模型中常用的分布族如下所示：
| 分布 | 函数 | 模型 |
|---|---|---|
| 正态(Gaussian) | $E(y)=X’\beta$ | 普通线性模型 |
| 二项(Binomial) | $E(y)=\frac{exp(X’\beta)}{1+exp(X’\beta)}$ | Logistic模型和概率模型单位(probit)模型 |
| 泊松(Poisson) | $E(y)=exp(X’\beta)$ | 对数线性模型 |
在R语言中，正态分布（也叫高斯分布）的广义线性模型事实上同线性模型是相同的，即
|
|
与线性模型fm <- lm(formula, data)得到的结论是一致的，当然，其效率会差很多。下面是广义线性模型函数glm()的用法：
|
|
这样，在广义线性意义下，我们不仅知道了一般线性模型是广义线性模型的一个特例，而且导出了处理频率资料的Logistic模型和处理频数资料的对数线性模型。这个重要结果还说明，虽然Logistic模型和对数线性模型都是非线性模型，即μ和β呈非线性关系，但通过连接函数使m(μ)和β呈线性关系，从而使我们可以用线性拟合的方法求解这类非线性模型。更有意义的是，实际研究中的主要数据形式无非是计量资料、频率资料和频数资料（半计量资料实际上可以看做有序的频数资料），因此，掌握了广义线性模型的思想和方法，就可以使用统一的方法处理各种类型的统计数据。
Logistic回归简介
在回归中，如果因变量不是数量型的，而是取两个可能值的分类变量，而自变量可以为定性变量及定量变量，那么此时普通的线性回归模型就不适用了，至少回归方程的左边不能是分类变量。我们回忆一下二项分布Bin(n,p)，取两个值的二分布变量和二项分布有关，满足二项分布的变量取两种可能的值，处理过程可以看成为n次独立的随机试验（Bernoulli试验），每次试验只可能取两个可能的值之一，如果称一个试验结果为“成功”，而另一个为“不成功”，并记每次试验成功的概率为p，那么不成功的概率则为1-p。在这里，在所有n次试验中，成功概率p保持不变，然而在实际情况中，概率在所有情况下都不变的假定往往不成立，概率p可能是一些自变量的函数，此时，回归方程的右边与普通的回归方程的类似，而左边则是概率p的函数，而不是p（因变量成功概率）本身，这种回归就是Logistic回归。
Logistic回归（logistic regression）属于概率型线性回归，它是研究二分类（可以扩展到多类）观察结果与一些影响因素之间关系的一种多变量分析方法。在流行病学研究中，经常需要分析疾病与各危险因素之间的定量关系，如食管癌的发生与吸烟、饮酒、不良包含习惯等危险因素的关系，为了正确说明这种关系，需要排除一些混杂因素的影响，传统上常使用Mantel-Haenszel分层分析方法，但这一方法适用于样本含量大、分析因素较少的情况。如果用线性回归，由于应变量Y是一个二值变量（通常取值为1或0），不满足应用条件，尤其当各因素都处于低水平或高水平时，预测值$\hat{Y}$可能超出0~1的范围，出现不合理的现象，用Logistic回归分析就是为了解决上述问题的。
基本概念
logistic回归模型
设应变量是一个二值变量，取值为0或1，如下所示：
$Y=
\begin{cases}
1 &出现阳性结果（发病、有效、死亡等） \\
0 &出现阴性结果（未发病、无效、存活等）
\end{cases}$
另有影响Y取值的m个自变变量，分别为$X_{1},X_{2},\dots,X_{m}$，记作$P=P(Y=1|X_{1},X_{2},\dots,X_{m})$，表示在m个自变量作用下阳性结果发生的概率，Logistic回归模型就可以表示为下面的公式，如下所示：
$P=\frac{1}{1+exp[-(\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\cdots+\beta_{m}X_{m})]} \quad (16-1)$
其中$\beta_{0}$是常数项，$\beta_{1},\beta_{2},\cdots,\beta_{m}$是回归系数，如果用Z表示m个自变量的线性组合，就是下面的样子：
$Z=\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\cdots+\beta_{m}X_{m}$
则Z与P之间的Logistic曲线就是下面的图形，如下所示：
从图中可以看出：当Z趋于$+\infty$时，P值逐渐趋于1；当Z趋于$-\infty$时，P值逐渐趋于0；P值的变化在0~1这个范围内，并且随着Z的增加或减少，以点$(0,0.5)$为中心呈对称S形变化，logistic模型的这些特点能够较好地配合生物学反应资料，对(16-1)的公式作对数变量，logistic回归模型就成了下面的线性形式，如下所示：
$ln(\frac{P}{1-P})=\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\cdots+\beta_{m}X_{m} \quad (16-2)$
上述公式的左端为阳性与阴性结果发生概率之比的自然对数，称为P的Logit变换，称为logitP，这就可以看出，虽然概率P的取值范围是0~1，但是logitP却没有数值界限，这里我们使用了ln(p/1-p)将p/1-p变成了线性方程，这里的ln(x)就是我们前面提到的连接函数。
模型参数意义
以流行病学来说明一下模型参数的意义，从公式(16-2)可以看出，常数项$\beta_{0}$表示暴露剂量为0时，个体发病与不发病概率比之的自然对数，回归系数$\beta_{j}(j=1,2,\cdots,m)$表示自变量$X_{j}$改变一个单位时，logitP的改变量，它与衡量危险因素作用大小的比数变称优势比（odds ratio，OR，也有教材译为比数比）有一个对应的关系，对比某一危险因素两个不同暴露水平$X_{j}=c_{1}$与$X_{j}=c_{0}$的发病情况（假定其他因素的水平相同），其优势比的自然对数如下所示：
$ \begin{align}lnOR_{j}&=ln[\frac{P_{1}/(1-P_{1})}{P_{0}/(1-P_{0})}]\\
&=logitP_{1}-logitP_{0}\\&=(\beta_{0}+\beta_{j}c_{1}+\sum \limits_{t\neq j}^m \beta_{t}X_{t}-(\beta_{0}+\beta_{j}c_{0}+\sum \limits_{t\neq j}^m \beta_{t}X_{t}) \\&=\beta_{j}(c_{1}-c_{0}) \end{align}\quad(16-3)$
转换成真数，则是如下的样子：
在上面的这些公式中，$P_{1}$和$P_{2}$分别表示在$X_{j}$取值为$c_{1}$及$c_{0}$时的发病概率，$OR_{j}$称作多变量调整后的优势比（adjusted odds ratio），表示扣除了其他自变量影响后危险因素的作用，如果$X_{j}$赋值上为1，如下所示：
那么，暴露组与非暴露组发病的优势比则为
$OR_{j}=exp(\beta_{j}) \qquad $(16-5)
当$\beta_{j}=0$时，$OR_{j}=1$，说明因素$X_{j}$对产生的发生不起作用；当$\beta_{j}>0$时，$OR_{j}>1$，说明$X_{j}$是一个危险因子；当$\beta_{j}<1$时，$OR_{j}<1$，说明$X_{j}$是一个保护因子。
由于$OR_{j}$值与模型中的常数项无关，$\beta_{0}$在危险因素分析中视为无效参数，对于发病率较低的慢性疾病如心脑血管疾病、恶性中肿瘤等，由于P很小，优势比可以作用相对危险度（relative risk，RR）的近似估计，如：
这是logistic回归用于流行病学调查资料的优势之一，即得到某一因素的回归系数估计值后，便可以得到该因素不同水平下优势比的估计值和相对危险度的近似估计值，关于比数比，可以参考这篇笔记《比数与比数比》。
Logistic回归模型的参数估计
参数估计
根据一组实际观察资料估计Logistic回归模型的参数时，通常采用最大似然法估计（maximum likelihood estimate, MLE），关于最大似然法，可以参考这篇笔记《最大似然法详解》。
最大似然法即建立一个样本似然函数，如下所示：
式中$P_{i}$表示第$i$例观察对象在暴露条件下阳性结果发生的概率，如果实际出现的是阳性结果，取$Y_{i}=1$，否则取$Y_{i}=0$，根据最大似然原理，在一次抽样中，获得现有样本的概率应该最大，即似然函数L应该达到最大值，为了简化计算，通常取似然函数的对数形式，如下所示：
形成要计算的目标函数lnL，然后采用Newton-Raphson迭代方法使对数似然函数达到极大值，此时参数的取值$b_{0},b_{1},b_{2},\cdots,b_{m}$即为$\beta_{0},\beta_{1},\beta_{2},\cdots,\beta_{m}$的最大似然估计值，同时得到参数估计值的方差，即协方差矩阵（对角线元素开平方为标准误$S_{b_{0}},S_{b_{1}},S_{b_{2}},\cdots,S_{b_{m}}$。
如果我们有下述的分层资料，如下所示：
高有g层，每层有$n_{k}(k=1,2,\cdots,g)$例观察以对象，其中有$d_{k}$个病例，则似然函数为：
其中，$P_{k}$表示第k层阳性结果发生的概率，上述求解过程现在都是通过软件来完成的。
优势比估计
根据公式(16-4)，其一因素两个不同水平$(c_{1},c_{0})$优势比的估计值如下所示：
$OR_{j}$的可信区间可以利用$b_{j}$的抽样分布来估计，在样本含量比较大的情况下，它近似服从正态分布，特殊情况下，若自变量$X_{j}$只有暴露和非暴露两个水平，则优势比$OR_{j}$的$1-\alpha$可信区间计算公式如下所示：
Logistic案例计算：
对前面所述的表格（16-1）做Logistic回归分析，也就是下面的这个表格
| 分层 | 吸烟$X_{1}$ | 饮酒$X_{2}$ | 观察例数$n_{k}$ | 阳性数$d_{k}$ | 阴性数$n_{k}-d_{k}$ |
|---|---|---|---|---|---|
| 1 | 0 | 0 | 199 | 63 | 136 |
| 2 | 0 | 1 | 170 | 63 | 107 |
| 3 | 1 | 0 | 101 | 44 | 57 |
| 4 | 1 | 1 | 416 | 265 | 151 |
确定各变量的赋值或编码，如下所示：
进行logistic回归计算，如下所示：
|
|
结果如下所示：
|
|
从上面结果可以看出：
$b_{0}=-0.9099,S_{b_{0}}=0.1358,b_{1}=0.8856,S_{b_{1}}=0.1500,b_{2}=0.5261,S_{b_{2}}=0.1572$
吸烟与不吸烟的优势比为：
$ \widehat{OR_{1}}=exp(b_{1})=exp(0.8856)=2.42 $
$OR_{1}$的95%可信区间如下所示：
饮酒与不饮酒的优势比为：
$ \widehat{OR_{2}}=exp(b_{2})=exp(0.5261)=1.69 $
$OR_{2}$的95%可信区间如下所示：
Logistic回归模型的假设检验
得到logistic回归方程后，还需要对其回归系数进行假设检验，用于说明所研究的自变量对于应变量Y的是影响是否有统计学意义，为此，需要对模型中的回归系数是否不全为0作出检验，这个跟多重线性回归的检验是一样的，检验假设条件如下所示：
$H_{0}:\beta_{1}=\beta_{2}=\cdots=\beta_{m}=0$
$H_{1}:$各$\beta_{j}(j=1,2,\cdots,m)$不全为0
更典型的问题是对一个回归系数的检验，检验假设为$H_{0}:\beta_{j}=0,H_{1}:\beta \neq 0$。常用的检验验方法有似然比检验（likelihood ratio test），Wald检验和计分检验（score test）。
似然比检验
似然比检验的基本思想是比较在两种不同假设条件下的对数似然函数值，看其差别大小。具体做法是，先拟合一个不包含准备检测的变量在内的logistic回归模型，求出它的对数似然函数值$lnL_{0}$，然后把需要检验的变量加入模型中去，再检验配合，得到一个新的对数似然函数值$lnL_{1}$，假设前后两个模型分别包含$l$个自变量和$p$个自变量，似然比统计量G的计算公式为：
当样本含量较大时，在零假设下得到的G统计量近似服从自由度为$d(d=p-l)$的$\chi^2$分布，若$G\geq{\chi^2_{\alpha,d}}$时，则新加入的$d$个自变量对回归方程有统计学意义，如果只对一个回归系数检验，则$d=1$，由上面的案例可以计算出来：
$lnL(X_{1})=-585.326,lnL(X_{2})=-597.436,lnL(X_{1},X_{2})=-579.711$
符号$L(X_{1})$和$L(X_{2})$分别表示模型中只含有$X_{1}$和$X_{2}$的最大似然函数值，而$L(X_{1},X_{2})$则表示模型中同时含有$X_{1}$和$X_{2}$的最大似然函数。
对于$H_{0}:\beta_{1}=0,H_{1}:\beta \neq 0$，则有：
$G=2[lnL(X_{1},X_{2})-lnL(X_{2})]=2[-579.711-(-597.436)]=35.45$
查$\chi^2$界值表可以得，$\chi^2_{0.05,1}=3.84$，$G>3.84$，因此在$\alpha=0.05$检验水平上拒绝$H_{0}$，接受$H_{1}$，说明控制了饮酒因素是呼后，食管癌与吸烟有显著性关系。
同理，对于$H_{0}:\beta_{2}=0,H_{1}:\beta_{2} \neq 0$，则有
$G=2[lnL(X_{1},X_{2})-lnL(X_{2})]=2[-579.711-(-585.326)]=11.23$
$G>3.84$，拒绝$H_{0}$，接受$H_{1}$，说明控制了吸烟因素的影响后，食管癌与饮酒有显著性关系。
看一下R对于似然比卡方的检验，因变量服从二项分布，因此对模型的方差检验需要使用似然比卡方检验，anova()函数提供了检验Logistic模型的馀在卡方比的功能，它也能逐一地检验模型的自变量是否显著。
|
|
从结果可能看出来，wine能够产生11.23的偏差，其p值小于0.01，对模型的影响是十分显著的。模型偏差度量了自变量对模型的影响程度，该值越小，就说明自变量造的偏差越小，对因变量的分类越没有帮助；该值越大，就说明自变量造成的偏差越大，对因变量的分类帮助越多，根据模型偏差的显著与否，能合理地判断自变量是否适合引物Logistic回归模型。
Wald检验
Wald检验的原理可以再看一下这篇笔记《比数与比数比》。
Wald只需要将各参数$\beta_{j}$的估计值$b_{j}$与0比较，并且它的标准误$S_{b_{j}}$作为参照；为检验$H_{0}:\beta_{j}=0,H_{1}:\beta_{j} \neq0$，计算如下统计量
$\mu=\frac{b_{j}}{S_{b_{j}}} \qquad (16-13)$
或
$\chi^2=(\frac{b_{j}}{S_{b_{j}}})^2 \qquad (16-14)$
对于大样本资料，在零假设下$\mu$近似服从标准正态分布，而$\chi^2$则近似服从自由度$\nu=1$的$\chi^2$分布，在本例中，对于$H_{0}:\beta_{1}=0,H_{1}:\beta_{2} \neq 0$，则有：
$\chi^2=(\frac{0.8856}{0.1500})^2=34.86$
对于$H_{0}:\beta_{2}=0,H_{1}:\beta_{2} \neq 0$，则有：
$\chi^2=(\frac{0.5261}{0.1572})^2=11.20$
$\chi^2$值均大于3.84，说明食管癌与吸烟、吸烟有显著性关系，结论同前。
Logistic回归变量筛选
与多元线性回归分析类似，当自变量的个数较多时，为了使建立的logistic回归模型比较稳定和便于解释，应尽可能将回归效果显著的自变量选入模型中，将作用不显著的自变量排除在外，具体的算法有前进法、后退法和逐步法。logistic逐步回归与线性逐步回归过程非常相似，使用其中所用的检验统计量不再是F统计量，而是似然比统计量、Wald统计量和计分统计量之一，例如使用似然比统计量，即利用如下的公式：
$G=2(lnL_{1}^{(l)}-lnL_{0}^{(l)})$
在进行到第$l$步时，通过比较含有某一自变量$X_{j}$和不含有$X_{j}$的模型，决定$X_{j}$是否引入模型或从模型中剔除。
看下面的案例（《医学统计学》，第四版，孙振球）：
案例B
为了探讨冠心病的有关危险因素，对26例冠心病患者和28便对照者进行病例对照研究，各因素的说明及资料见下表1和表2，试用logistic逐步回归分析法筛选危险因素（$\alpha_{入}=0.10$，$\alpha_{出}=0.15$），如下所示：
| 表1 | ||
|---|---|---|
| 因素 | 变量名 | 赋值说明 |
| 年龄（岁） | X1 | <45=1,45~=2,55~=3,65~=4 |
| 高血压史 | X2 | 无=0,有=1 |
| 高血压家族史 | X3 | 无=0,有=1 |
| 吸烟 | X4 | 不吸=0,吸=1 |
| 高血脂史 | X5 | 无=0,有=1 |
| 动物脂肪摄入 | X6 | 低=0,高=1 |
| 体重指数（BMI） | X7 | <24~=1,24~=2,26~=3 |
| A型性格 | X8 | 否=0,是=1 |
| 冠心病 | Y | 对照=0,病例=1 |
| 表2 | |||||||
|---|---|---|---|---|---|---|---|
| 序号 | X1 | X2 | X3 | … | X7 | X8 | Y |
| 1 | 3 | 1 | 0 | … | 1 | 1 | 0 |
| 2 | 2 | 0 | 1 | … | 1 | 0 | 0 |
| 3 | 2 | 1 | 0 | … | 1 | 0 | 0 |
| … | … | … | … | … | … | … | … |
| 54 | 3 | 1 | 1 | … | 3 | 1 | 1 |
（原始数据已放github上）
计算过程如下所示：
|
|
计算结果如下所示：
|
|
从结果来看，X6，X8这2项的P值小于0.05，现在再对这2项进行回归拟合，如下所示：
|
参考资料
- 孙振球. 医学统计学.第4版[M]. 人民卫生出版社, 2014.
- 王斌会. 多元统计分析及R语言建模[M]. 暨南大学出版社, 2010.
- 逻辑回归Logistics regression公式推导