这里是佳奥！让我们开始广义线性模型篇的学习吧！

广义线性模型扩展了线性模型的框架，它包含了非正态因变量的分析。

在本章中，我们将首先简要概述广义线性模型，并介绍如何使用glm()函数来进行估计。

然后我们将重点关注该框架中两种流行的模型：Logistic回归（因变量为类别型）和泊松回归（因变量为计数型）。

1 广义线性模型和glm( )函数

现假设想要对响应变量Y和p个预测变量X1...Xp间的关系进行建模。

在标准线性模型中，可以假设Y呈正态分布，关系的形式为：

给定一系列X变量的值，赋予X变量合适的权重，然后将它们加起来，便可预测Y观测值分布的均值。

而广义线性模型拟合的形式为：

其中g(μY)是条件均值的函数（称为连接函数）。设定好连接函数和概率分布后，便可以通过最大似然估计的多次迭代推导出各参数值。

1.1 glm( )函数

函数的基本形式为：

glm( )的参数：

glm( )函数可以拟合许多流行的模型，比如Logistic回归、泊松回归。后续内容针对这两个展开：

假设有一个响应变量（Y）、三个预测变量（X1、X2、X3）和一个包含数据的数据框（mydata）。

可用如下代码拟合Logistic回归模型：

拟合泊松回归模型：

1.2 连用的函数

与分析标准线性模型时lm( )连用的许多函数在glm( )中都有对应的形式，与glm( )连用的函数见下表：

1.3 模型拟合和回归诊断

当评价模型的适用性时，可以绘制初始响应变量的预测值与残差的图形：

其中，model为glm( )函数返回的对象。

R将列出帽子值（hat value）、学生化残差值和Cook距离统计量的近似值。

对于识别异常点的阈值，现在并没统一答案，其中一个方法就是绘制各统计量的参考图，然后找出异常大的值。例如，如下代码可创建三幅诊断图：

它可以创建一个综合性的诊断图。在后面的图形中，横轴代表杠杆值，纵轴代表学生化残差值，而绘制的符号大小与Cook距离大小成正比。

2 Logistic 回归

以AER包中的数据框Affairs为例，我们将通过探究婚外情的数据来阐述Logistic回归的过程。

该数据从601个参与者身上收集了9个变量，包括一年来婚外私通的频率以及参与者性别、年龄、婚龄、是否有小孩、宗教信仰程度（5分制，1分表示反对，5分表示非常信仰）、学历、职业（逆向编号的戈登7种分类），还有对婚姻的自我评分（5分制，1表示非常不幸福，5表示非常幸福）。

我们先来看一下数据：

从这些统计信息可以看到，52%的调查对象是女性，72%的人有孩子，样本年龄的中位数为32岁。

对于响应变量，72%的调查对象表示过去一年中没有婚外情（451/601），而婚外偷腥的最多次数为12（占了6%）。

但此处我们感兴趣的是二值型结果（有过一次婚外情/没有过婚外情），将affairs转化为二值型因子ynaffair：

该二值型因子现可作为Logistic回归的结果变量：

从回归系数的p值（最后一栏）可以看到，性别、是否有孩子、学历和职业对方程的贡献都不显著（无法拒绝参数为0的假设）。

去除这些变量重新拟合模型，检验新模型是否拟合得好：

新模型的每个回归系数都非常显著（p<0.05）。

由于两模型嵌套（fit.reduced是fit.full的一个子集），可以使用anova( )函数对它们进行比较，对于广义线性回归，可用卡方检验：

结果的卡方值不显著（p=0.21），表明四个预测变量的新模型与九个完整预测变量的模型拟合程度一样好。

结论：添加性别、孩子、学历和职业变量不会显著提高方程的预测精度，因此可以依据更简单的模型进行解释。

2.1 解释模型参数

先看看回归系数：

在Logistic回归中，响应变量是Y=1的对数优势比（log）。

回归系数含义是当其他预测变量不变时，一单位预测变量的变化可引起的响应变量对数优势比的变化。

由于对数优势比解释性差，可对结果进行指数化：

可以看到婚龄增加一年，婚外情的优势比将乘以1.106（保持年龄、宗教信仰和婚姻评定不变）。而如果婚龄增加10年，优势比将乘以1.106^10，即2.7。

2.2 评价预测变量对结果概率的影响

使用predict( )函数，可观察某个预测变量在各个水平时对结果概率的影响。

首先创建一个包含你感兴趣预测变量值的虚拟数据集，然后对该数据集使用predict( )函数，以预测这些值的结果概率：

从这些结果可以看到，当婚姻评分从1（很不幸福）变为5（非常幸福）时，婚外情概率从0.53降低到了0.15（假定年龄、婚龄和宗教信仰不变）。

下面我们再看看年龄的影响：

此处可以看到，当其他变量不变，年龄从17增加到57时，婚外情的概率将从0.34降低到0.11。

2.3 过度离势

过度离势会导致奇异的标准误检验和不精确的显著性检验。

检测过度离势的一种方法是比较二项分布模型的残差偏差与残差自由度，如果比值：

比1大很多，便可认为存在过度离势。

回到婚外情的例子，可得：

它非常接近于1，表明没有过度离势。

还可以对过度离势进行检验。为此，需要拟合模型两次，第一次使用family ="binomial"，第二次使用family = "quasibinomial"。

假设第一次glm( )返回对象记为fit，第二次返回对象记为fit.od，那么，

提供的p值即可对零假设H0：Φ = 1 与备择假设H1：Φ ≠ 1进行检验。若p很小（小于0.05），便可拒绝零假设。

又报错了，，，

一行一行排查，原来是拼写错误，，

此处p值（0.34）显然不显著（p > 0.05），我们认为不存在过度离势。

2.4 扩展

R中扩展的Logistic回归和变种：

1、稳健Logistic回归 robust包中的glmRob( )函数可用来拟合稳健的广义线性模型。

2、多项分布回归 若响应变量包含两个以上的无序类别（比如，已婚/寡居/离婚），便可使用mlogit包中的mlogit( )函数拟合多项Logistic回归。

3、序数Logistic回归 若响应变量是一组有序的类别（比如，信用风险为差/良/好），便可使用rms包中的lrm( )函数拟合序数Logistic回归。

3 泊松回归

当通过一系列连续型和/或类别型预测变量来预测计数型结果变量时，泊松回归是一个非常有用的工具。

我们将使用robust包中的Breslow癫痫数据。

我们就遭受轻微或严重间歇性癫痫的病人的年龄和癫痫发病数收集了数据，包含病人被随机分配到药物组或者安慰剂组前八周和随机分配后八周两种情况。响应变量为sumY（随机化后八周内癫痫发病数），预测变量为治疗条件（Trt）、年龄（Age）和前八周内的基础癫痫发病数（Base）。

在解释这些协变量后，我们感兴趣的是药物治疗是否能减少癫痫发病数。

首先看看数据集的统计汇总信息：

我们只关注之前描述的四个变量。基础和随机化后的癫痫发病数都有很高的偏度。

现在，我们更详细地考察响应变量：

可以清楚地看到因量的偏倚特性以及可能的离群点。

初看图形时，药物治疗下癫痫发病数似乎变小了，且方差也变小了（泊松分布中，较小的方差伴随着较小的均值）。

与标准最小二乘回归不同，泊松回归并不关注方差异质性。

接下来拟合泊松回归：

输出结果列出了偏差、回归参数、标准误和参数为0的检验。

此处预测变量在p<0.05的水平下都非常显著。

3.1 解释模型参数

使用coef( )函数可获取模型系数，或者调用summary( )函数的输出结果中的Coefficients表格：

年龄的回归参数为0.0227，表明保持其他预测变量不变，年龄增加一岁，癫痫发病数的对数均值将相应增加0.03。

指数化系数：

保持其他变量不变，年龄增加一岁，期望的癫痫发病数将乘以1.023。

3.2 过度离势

泊松分布的方差和均值相等。当响应变量观测的方差比依据泊松分布预测的方差大时，泊松回归可能发生过度离势。

与Logistic回归类似，此处如果残差偏差与残差自由度的比例远远大于1，那么表明存在过度离势。对于癫痫数据，它的比例为：

很显然，比例远远大于1。

qcc包提供了一个对泊松模型过度离势的检验方法：

后续步骤同Logistic，这里不再赘述。

扩展

R提供了基本泊松回归模型的一些有用扩展，包括允许时间段变化、存在过多0时会自动修正的模型，以及当数据存在离群点和强影响点时有用的稳健模型。包括：

1时间段变化的泊松回归

2零膨胀的泊松回归

3稳健泊松回归

4 小结

本篇在简短介绍了通用方法后，探究了Logistic回归和泊松回归。

随后，我们讨论了过度离势问题，包括如何检测以及依据它进行调整等方法。

在下一篇的学习中，我们将学习如何使用因子分析方法检测和检验这些无法被观测到的变量的假设。

我们下一篇再见！