R in action读书笔记（18）第十三章

本章内容

建立广义线性模型

预测类别型变量

计数型数据建模

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

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

其中g(μY)是条件均值的函数（称为连接函数）。另外，可放松Y为正态分布的假设，改为Y

服从指数分布族中的一种分布即可。设定连接函数和概率分布后，便可以通过最大似然估计的

多次迭代推导出各参数值。

13.1.1 glm()函数

R中可通过glm函数拟合广义线性模型。

glm(formula,family=family(link=function),data=)

glm()函数可以拟合许多流行的模型，比如Logistic回归、泊松回归和生存分析

13.1.2 连用的函数

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

13.1.3 模型拟合和回归诊断

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

plot(predict(model,type=”response”),residuals(model,type=”deviance”))

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

创建三幅诊断图:

> library(car)

> influncePlot(model)

13.2 Logistic 回归

Logistic回归适用于二值响应变量（0,1）。模型假设Y服从二项分布，线性模型的拟合形式为：

其中π= μY 是Y的条件均值（即给定一系列X的值时Y = 1的概率），(π/1-π)为Y = 1 时的优势比，log(π/1-π)为对数优势比，或logit。本例中，log(π/1-π)为连接函数，概率分布为二项分布。

> data(Affairs,package="AER")

> summary(Affairs)

    affairs          gender         age         yearsmarried    children 

 Min.   : 0.000   female:315   Min.   :17.50   Min.   : 0.125   no :171  

 1st Qu.: 0.000   male  :286   1st Qu.:27.00   1st Qu.: 4.000   yes:430  

 Median : 0.000                Median :32.00   Median : 7.000            

 Mean   : 1.456                Mean   :32.49   Mean   : 8.178            

 3rd Qu.: 0.000                3rd Qu.:37.00   3rd Qu.:15.000            

 Max.   :12.000                Max.   :57.00   Max.   :15.000            

 religiousness     education       occupation        rating     

 Min.   :1.000   Min.   : 9.00   Min.   :1.000   Min.   :1.000  

 1st Qu.:2.000   1st Qu.:14.00   1st Qu.:3.000   1st Qu.:3.000  

 Median :3.000   Median :16.00   Median :5.000   Median :4.000  

 Mean   :3.116   Mean   :16.17   Mean   :4.195   Mean   :3.932  

 3rd Qu.:4.000   3rd Qu.:18.00   3rd Qu.:6.000   3rd Qu.:5.000  

 Max.   :5.000   Max.   :20.00   Max.   :7.000   Max.   :5.000  

> table(Affairs$affairs)

  0   1   2   3   7  12 

451  34  17  19  42  38 

将affairs转化为二值型因子ynaffair

> Affairs$ynaffair[Affairs$affairs>0]<-1

> Affairs$ynaffair[Affairs$affairs==0]<-0

> Affairs$ynaffair<-factor(Affairs$ynaffair,levels=c(0,1),labels=c("No","Yes"))

> table(Affairs$ynaffair)

 No Yes 

451  150 

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

> fit.full<-glm(ynaffair~gender+age+yearsmarried+children+religiousness+education+occupation+rating,data=Affairs,family=binomial())

> summary(fit.full)

Call:

glm(formula = ynaffair ~ gender + age + yearsmarried +children +

religiousness +education + occupation + rating, family = binomial(),

data = Affairs)

Deviance Residuals:

Min 1Q Median 3Q Max

-1.5713 -0.7499 -0.5690 -0.2539 2.5191

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) 1.37726 0.88776 1.551 0.120807

gendermale 0.28029 0.23909 1.172 0.241083

age -0.04426 0.01825 -2.425 0.015301 *

yearsmarried 0.09477 0.03221 2.942 0.003262 **

childrenyes 0.39767 0.29151 1.364 0.172508

religiousness -0.32472 0.08975 -3.618 0.000297 ***

education 0.02105 0.05051 0.417 0.676851

occupation 0.03092 0.07178 0.431 0.666630

rating -0.46845 0.09091 -5.153 2.56e-07 ***

---

Signif. codes: 0‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance:675.38 on 600 degrees of freedom

Residual deviance: 609.51 on 592 degrees of freedom

AIC: 627.51

Number of Fisher Scoring iterations: 4

从回归系数的p值可以看到，性别、是否有孩子、学历和职业对方程的贡献都不显著（无法拒绝参数为0的假设）。去除这些变量重新拟合模型，检验新模型是否拟合得

> fit.reduced<-glm(ynaffair~age+yearsmarried+religiousness+rating,data=Affairs,family=binomial())

> summary(fit.reduced)

Call:

glm(formula = ynaffair ~ age + yearsmarried +religiousness +

rating, family= binomial(), data = Affairs)

Deviance Residuals:

Min 1Q Median 3Q Max

-1.6278 -0.7550 -0.5701 -0.2624 2.3998

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) 1.93083 0.61032 3.164 0.001558 **

age -0.03527 0.01736 -2.032 0.042127 *

yearsmarried 0.10062 0.02921 3.445 0.000571 ***

religiousness -0.32902 0.08945 -3.678 0.000235 ***

rating -0.46136 0.08884 -5.193 2.06e-07 ***

---

Signif. codes: 0‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance:675.38 on 600 degrees of freedom

Residual deviance: 615.36 on 596 degrees of freedom

AIC: 625.36

Number of Fisher Scoring iterations: 4

新模型的每个回归系数都非常显著（p<0.05）。由于两模型嵌套（fit.reduced是fit.full的一个子集），你可以使用anova()函数对它们进行比较，对于广义线性回归，可用卡方检验。

> anova(fit.reduced,fit.full,test="Chisq")

Analysis of Deviance Table

Model 1: ynaffair ~ age + yearsmarried + religiousness +rating

Model 2: ynaffair ~ gender + age + yearsmarried +children + religiousness +

education +occupation + rating

Resid. Df Resid.Dev Df Deviance Pr(>Chi)

1 596 615.36

2 592 609.51 4 5.8474 0.2108

> coef(fit.reduced)

(Intercept) age yearsmarried religiousness rating

1.93083017 -0.03527112 0.10062274 -0.32902386 -0.46136144

13.2.1 解释模型参数

看回归系数：

> coef(fit.reduced)

(Intercept) age yearsmarriedreligiousness rating

1.93083017 -0.03527112 0.10062274 -0.32902386 -0.46136144

在Logistic回归中，响应变量是Y=1的对数优势比（log）。回归系数含义是当其他预测变量不变时，一单位预测变量的变化可引起的响应变量对数优势比的变化。由于对数优势比解释性差，可对结果进行指数化：

> exp(coef(fit.reduced))

(Intercept) age yearsmarriedreligiousness rating

6.8952321 0.9653437 1.1058594 0.7196258 0.6304248

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

使用predict()函数，可观察某个预测变量在各个水平时对结果概率的影响。首先创建一个包含你感兴趣预测变量值的虚拟数据集，然后对该数据集使用predict()函数，以预测这些值的结果概率。

> testdata<-data.frame(rating=c(1,2,3,4,5),age=mean(Affairs$age),yearsmarried=mean(Affairs$yearsmarried),religiousness=mean(Affairs$religiousness))

> testdata

rating age yearsmarried religiousness

1 132.48752 8.177696 3.116473

2 232.48752 8.177696 3.116473

3 332.48752 8.177696 3.116473

4 432.48752 8.177696 3.116473

5 532.48752 8.177696 3.116473

使用测试数据集预测相应的概率：

rating age yearsmarried religiousness prob

1 132.48752 8.177696 3.116473 0.5302296

2 232.48752 8.177696 3.116473 0.4157377

3 332.48752 8.177696 3.116473 0.3096712

4 4 32.48752 8.177696 3.116473 0.2204547

5 532.48752 8.177696 3.116473 0.1513079

再看看年龄的影响：

> testdata<-data.frame(rating=mean(Affairs$rating),

+ age=seq(17,57,10),

+ yearsmarried=mean(Affairs$yearsmarried),

+ religiouseness=mean(Affairs$religiousness))

> testdata

rating ageyearsmarried religiouseness

1 3.93178 17 8.177696 3.116473

2 3.93178 27 8.177696 3.116473

3 3.93178 37 8.177696 3.116473

4 3.93178 47 8.177696 3.116473

5 3.93178 57 8.177696 3.116473

>testdata$prob<-predict(fit.reduced,newdata=testdata,type="response")

> testdata

rating ageyearsmarried religiousness prob

1 3.93178 17 8.177696 3.116473 0.3350834

2 3.93178 27 8.177696 3.116473 0.2615373

3 3.93178 37 8.177696 3.116473 0.1992953

4 3.93178 47 8.177696 3.116473 0.1488796

5 3.93178 57 8.177696 3.116473 0.1094738

13.2.3 过度离势

抽样于二项分布的数据的期望方差是σ2 = nπ(1-π)，n为观测数，π为属于Y = 1组的概率。所谓过度离势，即观测到的响应变量的方差大于期望的二项分布的方差。过度离势会导致奇异的标准误检验和不精确的显著性检验。

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

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

13.2.4 扩展

稳健Logistic回归robust包中的glmRob()函数可用来拟合稳健的广义线性模型，包括稳健Logistic回归。当拟合Logistic回归模型数据出现离群点和强影响点时，稳健Logistic回归便可派上用场。

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

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

13.3 泊松回归

泊松回归适用于在给定时间内响应变量为事件发生数目的情形。它假设Y服从泊松分布，线性模型的拟合形式为：

其中λ是Y的均值（也等于方差）。此时，连接函数为log(λ)，概率分布为泊松分布。

数据集的统计汇总信息：

> library(robust)

> data(breslow.dat,package="robust")

> names(breslow.dat)

[1]"ID" "Y1" "Y2" "Y3" "Y4" "Base" "Age" "Trt"

[9]"Ysum" "sumY" "Age10" "Base4"

> summary(breslow.dat[c(6,7,8,10)])

Base Age Trt sumY

Min. : 6.00 Min. :18.00 placebo :28 Min. : 0.00

1st Qu.:12.00 1st Qu.:23.00 progabide:31 1st Qu.: 11.50

Median :22.00 Median :28.00 Median : 16.00

Mean : 31.22 Mean :28.34 Mean : 33.05

3rd Qu.:41.00 3rd Qu.:32.00 3rd Qu.: 36.00

Max. :151.00 Max. :42.00 Max. :302.00

生成图形：

> par(opar)

> opar<-par(no.readonly=TRUE)

> par(mfrow=c(1,2))

> attach(breslow.dat)

> hist(sumY,breaks=20,xlab="Seizure Count",

+ main="Distribution of Seizures")

>boxplot(sumY~Trt,xlab="Treatment",main="Group Comparisons")

> par(opar)

拟合泊松回归：

输出结果列出了偏差、回归参数、标准误和参数为0的检验。注意，此处预测变量在p<0.05的水平下都非常显著。

13.3.1 解释模型参数

在泊松回归中，因变量以条件均值的对数形式ln(λ)来建模。与Logistic回归中的指数化参数相似，泊松模型中的指数化参数对响应变量的影响都是成倍增加的，而不是线性相加。同样，还需要评价泊松模型的过度离势。可能发生过度离势的原因有如下几个：

遗漏了某个重要的预测变量。

可能因为事件相关。

在纵向数据分析中，重复测量的数据由于内在群聚特性可导致过度离势

13.3.3 扩展

1. 时间段变化的泊松回归

2. 零膨胀的泊松回归

3.稳健泊松回归