正在加载
请稍等

菜单

Home 码农菜园 数据分析 R学习笔记-11 广义线型模型
Home 码农菜园 数据分析 R学习笔记-11 广义线型模型

R学习笔记-11 广义线型模型

数据分析 by   阅读量 10,718

线性模型可以通过一系列连续型和/或类别型预测变量来预测正态分布的响应情况,但在许多情况下,假设因变量为正态分布(甚至连续型变量)并不合理,例如:

  • 结果变量可能是类别型的。二值变量(是/否、通过/失败、存活/死亡)和多类别变量(优/良/可/差)都显然不是正态分布;
  • 结果变量可能是计数型的(一周交通事故的数目、每日酒水消耗的数量),这类变量都是非负的有限值,而它们的均值和方差通常都是相关的(正态分布变量间则是相互独立)。

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

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

许多广泛应用流行的数据分析方法其实都归属于广义线性模型框架。

glm()函数

R中可以使用glm()函数(当然还有其他专门的函数)拟合广义线性模型,形式如下:

其中family为拟合所属的函数族,function为对应的连接函数。glm()可以拟合许多流行的模型,例如Logistic回归、泊松回归和生存分析,以下详细讨论前两个模型。假设有一个相应变量(Y)、三个预测变量(X1、X2、X3)和一个包含数据的数据框(mydata)。

Logistic回归适用于二值响应变量,连接函数为logit函数,概率分布为二项分布,故拟合代码为:

泊松回归适用于在给定时间内响应变量为事件发生数目的情形,其假设Y服从泊松分布,连接函数为log(λ),概率分布为泊松分布,故拟合代码为:

值得注意的是,标准线性模型也是广义线性模型的一个特例,连接函数为恒等函数,概率分布为正态高斯分布,则拟合代码为:

生成的结果与以下代码相同:

总之,广义线性模型通过拟合响应变量的条件均值的一个函数(不是响应变量的条件均值),并假设响应变量服从指数分布族中的某个分布(不限于正态分布),从而极大地扩展了标准线性模型。模型参数估计的推导依据是极大似然估计,而非最小二乘法。

 类似的函数

与分析标准线性模型时lm()类似的许多函数在glm()中都有对应的形式,常见的函数如下所示:

  • summary():展示拟合模型的细节;
  • coefficients()、coef():列出拟合模型的参数(截距项和斜率);
  • confint():给出模型参数的置信区间(默认为95%);
  • residuals():列出拟合模型的残差值;
  • anova():生成两个拟合模型的方差分析表;
  • plot():生成评价拟合模型的诊断图;
  • predict():用拟合模型对新数据集进行预测。

模型拟合和回归诊断

与标准(OLS)线性模型一样,模型适用性的评价对于广义线性模型也非常重要。当评价模型的适用性时,可以绘制初始响应变量的预测值与残差的图形。使用以下代码绘制一个常见的诊断图,其中model为glm()函数返回的对象。

R将列出帽子值(hat value)、学生化残差值和Cook距离统计量的近似值。使用以下代码可创建三幅诊断图:

或者使用以下代码,可以创建一个综合性的诊断图,横轴代表杠杆值,纵轴代表学生化残差值,绘制的符号大小与Cook距离大小成正比。

11.2 Logistic回归

当通过一系列连续型和/或类别型预测变量来预测二值型结果变量时,Logistic回归是一个非常有用的工具。以ARE包中的数据框Affairs为例,该数据从601个参与者身上收集了9个变量,包括一年来出轨的频率以及参与者的性别、年龄、婚龄、是否有小孩、宗教信仰程度(5分制,1分表示非常反对,5分表示非常信仰)、学历、职业(逆向编号的戈登7种分类)、对婚姻的自我评分(5分制,1分表示非常不幸福,5分表示非常幸福)。

首先来看一些描述性的统计信息:

可以得到一些统计信息,如52%的调查对象是女性,72%的人有孩子。我们不关注出轨的次数,而只关注是否出轨,即一个二值型结果(出轨过/从未出轨),使用以下代码将affaris转化为二值型因子ynaffair,该二值型因子即可以作为Logistic回归的结果变量。

从回归系数的p值(最后一栏)可以看到,性别、是否有孩子、学历和职业对方程的贡献都不显著(p值较大)。去除这些变量重新拟合模型,检验新模型是否拟合得更好。

新模型的每个回归系数都非常显著(p<0.05),由于两模型嵌套,故可以使用anova()函数对它们进行比较。对于广义线性回归,可用卡方检验。检验结果的卡方值不显著(p=0.21较大),表明四个预测变量的新模型和九个预测变量的旧模型模拟程度其实差不多。

解释模型参数

首先查看回归系数:

在Logistic回归中,响应变量是Y=1的对数优势比。回归系数的含义是,当其他预测变量不变时,一单位预测变量的变化将引起的响应变量对数优势比的变化。可以将回归系数指数化,以便更易于分析。

可以看到婚龄每增加一年,出轨的优势比将乘以1.106;年龄每增加一岁,出轨的优势比将乘以0.965。因此,随着婚龄的增加和年龄、宗教信仰与婚姻评分的降低,出轨的优势比将上升。因为预测变量不能为0,故截距项在此处没有特定含义。

如果有需要,还可以使用confint()函数获取系数的置信区间,使用exp(confint(fit.reduced))可在优势比尺度上得到系数95%的置信区间。

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

既然已经有了模型,相对于优势比,更有吸引力的结论也许是,给定一个用户行为,预测其出轨的概率。可以使用predict()函数进行预测,以下代码生成了一些虚拟数据并进行了预测,主要评测婚姻评分对出轨概率的影响。

从结果可以看到,当婚姻评分从1分(很不幸福)增加到5分(非常幸福)时,出轨概率从0.53降低至0.15,接下来再看看年龄的影响。

可以看出,当其他变量不变、年龄从17增加到57时,出轨的概率从0.34降低至0.11。

过度离势

过度离势指的是观测到的响应变量方差大于期望的二项分布方差,过度离势会导致奇异的标准误检验和不精确的显著性检验。当出现过度离势时,仍可使用glm()函数拟合Logistic回归,但此时需要将二项分布改为类二项分布。

检测过度离势的一种方法是比较二项分布模型的残差偏差与残差自由度(在四变量模型中分别是615.36和596),如果两者比值比1大很多,则可以认为存在过度离势。

扩展

R中扩展的Logistic回归和变种如下所示:

  • 稳健Logistic回归:当拟合Logistic回归模型数据出现离群点和强影响点时,robust包中的glmRob()函数可用来拟合稳健的广义线性模型;
  • 多项分布回归:当响应变量包含两个以上的无序类别(例如已婚、寡居、离婚)时,可使用mlogit包中的mlogit()函数拟合多项Logistic回归;
  • 序数Logistic回归:当响应变量是一组有序的类别(例如信用为好/良/差)时,可使用rms包中的lrm()函数拟合序数Logistic回归。

11.3 泊松回归

当通过一系列连续型和/或类别型预测变量来预测计数型结果变量时,泊松回归是一个非常有用的工具。以下讨论使用robust包中的Breslow癫痫数据,数据记录了在治疗初期的八周内,抗癫痫药物对癫痫发病数的影响。

响应变量为sumY(随机化后八周内癫痫发病数),预测变量为治疗条件(Trt)、年龄(Age)和前八周内的基础癫痫发病数(Base),我们感兴趣的是药物治疗能否减少癫痫发病数。

首先来看数据的统计汇总信息,虽然一共有12个变量,但是我们只关注之前描述的四个变量。

使用以下代码详细地考察响应变量,基础和随机化后的癫痫发病数都有很高的偏度,从图中可以清楚地看到因变量的偏倚特性以及可能的离群点。

癫痫发病数分布

接下来拟合泊松回归,输出结果列出了偏差、回归参数、标准误和参数为0的检验,需要注意的是,此处预测变量在p<0.05的水平下都非常显著。

解释模型参数

使用coef()函数获取模型参数,或者调用summary()函数输出结果中的Coefficients表格。在泊松回归中,因变量以条件均值的对数形式来建模。年龄的回归参数为0.0227,表明保持其他预测变量不变,年龄增加一岁,癫痫发病数的对数均值将相应增加0.0227。我们同样可以将参数对数化,以更加易于分析。

结果显示,年龄增加一岁,期望的癫痫发病数将乘以1.023,即年龄的增加将导致更高的发病数;一单位Trt的变化(即从安慰剂到治疗组),期望的发病数将乘以0.86,也即保持基础发病数和年龄不变,服药组相对于安慰剂组的发病数降低了14%。

过度离势

在泊松分布中同样可能存在过度离势,泊松分布的方差和均值相等,当响应变量观测的方差比依据泊松分布预测的方差大时,可能存在过度离势。

与Logistic回归类似,但此例中残差偏差(559.44)与残差自由度(55)的比例为10.17远大于1,故存在过度离势。

还可以使用qcc包来检验泊松模型是否存在过度离势。

结果显示,显著性检验的p值为0(小于0.05),进一步表明确实存在过度离势。在这种情况下,可以考虑使用类泊松回归替代泊松回归。

扩展

R中提供了基本泊松回归模型的一些有用扩展:

  • 时间段变化的泊松回归:如果每个观测的时间段长短不同,则因变量需要从计数更改为比率,即发生次数除以观测时间,使用glm()的offset选项即可,offset=log(time),其中time是每名病人的观测时间;
  • 零膨胀的泊松回归:在一个数据集中,0计数的数目时常比用泊松模型预测得多。在Affairs数据集中,很有可能有一群从未出轨的对象,他们称为数据集的结构零值。可以用零膨胀的泊松回归来分析这种情况,它将同时拟合两个模型,一个用来预测哪些人又会出轨,另一个用来预测排除了结构零值之后的调查对象会出轨多少次。可以将其看作是Logistic回归(预测结构零值)和泊松回归(预测无结构零值观测的计数)的组合,pscl包的zeroinfl()函数可做零膨胀泊松回归;
  • 稳健泊松回归:robust包中的glmRob()函数可以拟合稳健广义线性模型,包括稳健泊松回归。

28 2015-04

2条评论

  1. 匿名说道:

    加个reference啊,R in action

发表评论