R语言实战第十三章

R语言实战第13章

广义线性模型

在普通的回归分析中，因许多因变量的假设都符合.假设因变量不符合正太分布或不符合连续型变量，如以下几种情况：

（1）：结果变量可能是类别型的，例如二值变量.

（2）：结果变量是计数型的，此类变量是肺腑的有限制，且均值与方差通常是相关的.

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

采用了glm()函数，采用的模型为logistic模型和泊松回归（因变量为计数型模型）.采用的例子为将个人信息作为变量去预测个人是否出轨，药物治疗对于病的发作次数是否有着影响.

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

广义线性模型中的响应变量是条件均值函数（成为连接函数），且相应变量的要求变松即要求为服从指数分布族中的一种分布即可.设好连接函数和概率分布后，便可通过最大似然估计来多次迭代推导各种参数值.

1.1glm()函数

    其函数的基本形式为：

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

    下表为概率分布(family)和相应默认的连接函数(function).

表1 glm()的参数

分布族

默认的连接函数

binomial

link=”logit”

gaussian

link=”identity”

gamma

link=”inverse”

inverse.gaussian

link=”1/mu^2”

poisson

link=”log”

quasi

link=”identity”,variance=”constant”

quasibinomial

link=”logit”

quasipoisson

link=”log”

    逻辑回归适用于二值相应变量，模型假设Y服从二项分布，线性模型的拟合式为：

=

    其中是Y的条件均值，（即给定一系列X的值时Y的概率），是Y=1的优势比，为对数优势比，此处为连接函数，概率分布为二项式，逻辑回归的代码范例为：

    glm(y~x1+x2+x3,family=binomial(link=“logit”)，data=mydata)

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

    其中是Y的均值（也就是方差），此时，连接函数为log()，概率分布为泊松分布.泊松分布的范例为：

       glm(Y~X1+x2+x3,family=poisson(link=”log”_,data=mydata)

       总之，广义线性模型通过拟合相应变量的条件均值的一个函数（非相应变量的条件均值），假设相应变量服从指数分布族的某个分布，极大的拓展了标准线性模型.

13.1.2 连用的函数

       常见的函数间下表.

与glm()连用的函数

函数

描述

summary()

展示拟合的细节

coefficients(),coef()

列出拟合的参数

confint()

给出模型参数的置信区间

residuals()

列出拟合模型的残差值

anova()

生成拟合模型的方差分析表

plot()

生成诊断图

predict()

用拟合模型对数据进行预测

deviance()

拟合模型的偏差

df.residuals()

拟合模型的残差自由度

13.1.3模型拟合和回归诊断

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

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

其中Model为glm()函数返回的对象.

采用帽子值，学生化残差，cook距离统计量的近似值，画出三幅诊断图

plot(hatvalues(model));

plot(rstudent(model))

plot(cook.distant(model))

或者用

Library(car)

influencePlot(model)

可以创造一个综合性的诊断图

13.2逻辑回归

library(AER)

data(Affairs,package="AER")

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"))

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.5

Number of Fisher Scoring iterations: 4

可见性别，是否有孩子，学历及职业对是否出轨的贡献不显著，去除掉该因素再对回归模型进行拟合

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

summary(fit.full)

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

采用anova()函数进行两个回归方程的比对，对于广义线性模型，采用卡方检验来比对，看结果的P值，若p值显著小于0.。05，则表明，新的方程有提高作用.

anova(fit.full,fit.reduct,test="Chisq")

Analysis of Deviance Table

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

    education + occupation + rating

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

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

1       592     609.51                     

2       596     615.36 -4  -5.8474   0.2108

13.2.1参数解释

coef(fit.reduct)

 (Intercept)           age  yearsmarried religiousness        rating 

  1.93083017   -0.03527112    0.10062274   -0.32902386   -0.46136144 

在逻辑回归中，显影变量是Y=1的对数优势比（log），回归系数的含义是当其他预测变量不变时，一单位预测变量的变化可引起相应变量对数优势比的变化.

exp(coef(fit.reduct))

  (Intercept)           age  yearsmarried religiousness        rating 

    6.8952321     0.9653437     1.1058594     0.7196258     0.6304248 

即表明，魂灵增加一年，婚外情的优势比将乘以1.106，相反，年龄增加一岁，婚外情的优势比乘以0.965，即随着婚龄的增加和年龄，总价信仰及婚姻评分的降低，婚外情优势比上升，由预测变量不为0，故截距项无特定含义，也可用confint()函数去获取系数的置信区间，后某预测变量N单元的变化引起较高值优势比的变化为我们所关注的，例如，婚龄增加一年，婚外情优势乘以1.106，若婚龄增加10年，优势比将乘以10.

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

        使用Predict（）函数，可以观察某个越策变量在个水平对结果概率的影响，首先创造某个感兴趣的预测变量值的虚拟数据集，然后对该数据使用predict()函数，预测该值的结果概率.

        采用该方法评价婚姻评分对婚外情概率的影响，先创建虚拟数据集，设定年龄，婚龄，宗教信仰为均值，婚姻评分为1~5.

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

testdata$prob<-predict(fit.reduct,newdata=ttest.data,type="respon")

testdata

  rating      age yearsmarried religiousness      prob

1      1 32.48752     8.177696      3.116473 0.5302296

2      2 32.48752     8.177696      3.116473 0.4157377

3      3 32.48752     8.177696      3.116473 0.3096712

4      4 32.48752     8.177696      3.116473 0.2204547

5      5 32.48752     8.177696      3.116473 0.1513079

    可见，在其他变量不变的情况下，婚姻幸福评分为1时，婚外情概率为0.53，婚姻幸福评分为5时，婚外情概率降为0.15

        现评价下年龄：

testdata<-data.frame(rating=mean(Affairs$rating),age=seq(17,57,10),yearsmarried=mean(Affairs$yearsmarried),religiousness=mean(Affairs$religiousness))

testdata

   rating age yearsmarried religiousness

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.reduct,newdata=testdata,type="respon")

testdata

   rating age yearsmarried 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过渡离势

     过度离势即是因变量的方差大于期望二项分布的方差，过度离势会导致奇异的标准误检验和不精确的显著性检验.若离势过高，仍可用glm()函数进行逻辑回归，但此时需要将二项分布改为类二项分布.

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

     比1大很多，便可认为存在过度离势，在婚外情的例子中：

deviance(fit.reduct)/df.residual(fit.reduct)

[1] 1.03248

        可见其比值接近于1，故可认为不存在过度离势.对过度离势进行检验：

pchisq(summary(fit.od)$dispersion*fit$df.residual,fit$df.residual,lower=F)

[1] 0.4292978

13.3泊松回归

        使用癫痫病数据，讨论治疗初期的霸州内，抗癫痫药物对癫痫病发病次数的影响.

        因变量为sumY,自变量为治疗条件（trt）,年龄(age)，前八周发病次数（Base），探索药物治疗是否能减少癫痫发病次数.

data(bresslow.dat,package="robust")

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

opar<-par(no.readonly=T)

par(mfrow=c(1,2))

attach(breslow.dat)

hist(sumY,breaks=20,xlab="seizure count",main="distribution")

boxplot(sumY~Trt,xlab="treatment",main="group comparisons")

        可见因变量的偏倚性及离群点，药物治疗下癫痫病的发病数变小，且方差也变小（泊松分布中，较小的方差意味着较小的均值），但与普通回归不同，泊松回归不关注其方差性质.

        泊松分布的拟合：      

fito<-glm(sumY~Trt+Age+Base,data=breslow.dat,family=poisson())

summary(fito)

Call:

glm(formula = sumY ~ Trt + Age + Base, family = poisson(), data = breslow.dat)

Deviance Residuals: 

    Min       1Q   Median       3Q      Max  

-6.0569  -2.0433  -0.9397   0.7929  11.0061  

Coefficients:

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

(Intercept)   1.9488259  0.1356191  14.370  < 2e-16 ***

Trtprogabide -0.1527009  0.0478051  -3.194   0.0014 ** 

Age           0.0227401  0.0040240   5.651 1.59e-08 ***

Base          0.0226517  0.0005093  44.476  < 2e-16 ***

---

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

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 2122.73  on 58  degrees of freedom

Residual deviance:  559.44  on 55  degrees of freedom

AIC: 850.71

Number of Fisher Scoring iterations: 5

     模型出现了偏差，回归参数，标准误和参数为0的检验.

13.3.1解释模型参数

exp(coef(fito))

 (Intercept) Trtprogabide          Age         Base 

7.0204403    0.8583864    1.0230007    1.0229102 

一单位的trt的变化，期望的癫痫发病数将乘以0.86，即表示，在发病数和年龄不变的情况下，付耀祖相对于安慰组的发病数降低了20%.13.3.2过度离势

泊松分布的方差和均值相等，当相应变量观测的方差比一句泊松分布预测的方差大时，泊松回归可能发生过度离势，可能产生过度离势的原因有如下:

(1):遗漏某个重要的预测变量

(2):可能因为事件相关，在泊松回归中，会认为每次事件都是独立发生的，以癫痫病为例，发病较多次的病人发生多一次病的概率较大，其概率会大于其他病人.

(3):纵向数据分析，重复测量的数据由于内在群聚特性可能导致过度离势.

若产生过度离势，可能会得到很小的标准误和置信区间，并且显著性检验也会相对教育宽松，及不糊发现并不真实存在的效应.其计算与逻辑回归分析一样，若残差偏差与残差自由度的比例大于1，则表明存在过度离势.

deviance(fito)/df.residual(fito)

[1] 10.1717

即可得比例大一1，可采用qcc包检验泊松分布的过度离势

qcc.overdispersion.test(breslow.dat$sumY,type="poisson")                  

Overdispersion test Obs.Var/Theor.Var Statistic p-value

poisson data          62.87013  3646.468       0

 即表示存在过度离势，处理泊松回归的过度离势的方法是讲family=”poisson”改为family=”quasipoisson”.

fit.od<-glm(sumY~Age+Base+Trt,data=breslow.dat,family=quasipoisson())

Call:

glm(formula = sumY ~ Age + Base + Trt, family = quasipoisson(), 

data = breslow.dat)

Deviance Residuals: 

 Min       1Q   Median       3Q      Max  

-6.0569  -2.0433  -0.9397   0.7929  11.0061  

Coefficients:

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

(Intercept)   1.948826   0.465091   4.190 0.000102 ***

Age           0.022740   0.013800   1.648 0.105085    

Base          0.022652   0.001747  12.969  < 2e-16 ***

Trtprogabide -0.152701   0.163943  -0.931 0.355702    

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

(Dispersion parameter for quasipoisson family taken to be 11.76075)

    Null deviance: 2122.73  on 58  degrees of freedom

Residual deviance:  559.44  on 55  degrees of freedom

AIC: NA

Number of Fisher Scoring iterations: 5

类泊松分布方法所得的参数估计与泊松方法相同，但标准误变大，当考虑过度离势时，会导致没有充分的证据证明，药物对癫痫病的疗效有作用.但应结合具体实际，不可随意删除变量.

13.3.3拓展

13.3.3.1有关泊松回归

泊松回归的公式：

Poisson(=ln(=

其中为Y的均值及方差，指数化后公式为：

即变量每增加1，即为乘以系数，如

13.3.3.2时间段变化的泊松回归

     即将固定的时间变为可变化的时间段，此处假设结果变量是一个比率. 在普通的泊松回归时，左侧为ln(λ)，时间段变化时，当为时间变化的泊松回归，左侧变成了=ln,将多余部分移至等式右边，即为参与时间段变化的泊松回归，在R语言中,采用glm()中的offset选项，在癫痫病案例中，将癫痫病的发病率做因变量，（假设记录了病人的发病时），拟合模型

fit<-glm(sumY~Base+Age+Trt,data=breslow.dat,offset=log(time),family = poisson)

13.3.3.3零膨胀的泊松回归

        零膨胀，即是原数据集中，因变量为0值占得比重很大,采用零膨胀泊松回归即是会同时拟合逻辑回归和泊松回归，一个是预测哪些人会发生婚外情（非零），另一个是预测排除了结构零值之后，那些发生了婚外情的人会有多少次.可采用pscl包中的zeroinf1()函数进行零膨胀的泊松回归.

       关于0膨胀泊松回归分析的例子，数据框DF中有ABCDEY六个变量，其中Y为零膨胀变量，即Y中含有大量0值，需要用其它五个变量来拟合Y变量，于是函数这么写：

zeroinfl(Y~A+B+C|D+E,data=DF) 
zeroinfl(Y~A+B+C+D|A+B+C+E,data=DF) 
zeroinfl(Y~A+B+C|1,data=DF)

　　第一条命令的意思为由A、B、C决定Y的次数，D和E决定Y是否非0，用专业术语来说，就是ABC三个是点模型内变量，DE是零膨胀模型内变量，一个变量既可以属于点模型变量，也属于零膨胀模型变量！如第二条命令。第三条命令的意思是，当Y取得零值是等可能的，此时可以使用1来代替。

13.3.3.4稳健泊松回归

        Robust包中的glmRob()函数可以拟合稳健广义线性回归，用于存在离群点和强影响点中.