13章广义线性模型

线性模型只适应于一系列连续性和/或类别型变量来预测正态分布的相应变量。
但是，许多情况下，假设因变量为正态分布（甚至连续型变量）并不合理：
1：结果变量可能使类别型。二值变量（比如：是/否，活着/死亡）和多分类变量（差/良好/优秀）都显然不是正态分布
2：结果变量可能是计数型的。（比如，一周交通事故的数目，每日酒水消耗的数量）

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

本章重点：logistic回归和泊松回归（计数型）
13.1.1 glm()函数

R语言中用glm函数来拟合广义线性模型

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

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

#Logistic回归适应二值响应变量（0,1）
glm(Y~x1+x2+x3,binomial(link = "logit"),data = mydata)

#泊松分布适应于在给定时间内响应变量为时间发生数目的情形
glm(Y~x1+x2+x3,poisson(link="log"),data = mydata)

##与glm利用的函数
summary() #展示拟合模型的细节
coefficients(),coef() #列出拟合模型的参数（截距项和斜率）
confint() #给出模型参数的置信区间
residuals() #列出拟合模型的残差值
anova() #生成两个拟合模型的方差分析表
plot() #生成评价拟合模型的诊断图
predict() #用拟合模型对新数据进行预测

13.1.3 模型的拟合和回归诊断

#当评价模型的适用性，可以绘制初始相应变量的预测值与残差的图形，例如，如下代码可绘制一个常见的诊断图
plot(predict(model,type="response"),residuals(model,type="deviance"))
#特别注意：当相应变量有很多值时，诊断图非常有用，而当相应变量只有有限的个值时（比如Logistic回归），诊断图的功效就会降低很多


13.2 Logistic实例分析

data(Affairs,package = "AER")
summary(Affairs)
table(Affairs$affairs)
#我们感兴趣的是二值型结果（有过一次/没有货婚外情）
Affairs$yaffairs[Affairs$affairs>0]=1
Affairs$yaffairs[Affairs$affairs==0]=0
Affairs$yaffairs=factor(Affairs$yaffairs,levels = c(0,1),labels = c("No","Yes"))

table(Affairs$yaffairs)

#调用logistic回归模型
attach(Affairs)
fit.full=glm(yaffairs~gender+age+yearsmarried+children+religiousness+education+occupation+rating,data=Affairs,family = binomial())

summary(fit.full)
***********************************************
  结果显示
Call:
  glm(formula = yaffairs ~ 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值可以看出，性别 是否有孩子 学历和职业对方程的贡献都不显著 去除这些变量重新拟合模型 检验新模型是否拟合很好



fit.reduced=glm(yaffairs~age+yearsmarried+religiousness+rating,data=Affairs,family = binomial())
summary(fit.reduced)

***********
  结果分析：新模型的每个回归系数都非常显著。由于两个模型嵌套，可以使用anova() 函数对他们进行比较，对于广义线性回归，可以用卡方检验
anova(fit.reduced,fit.full,test="Chisq")

****************************结果
Analysis of Deviance Table

Model 1: yaffairs ~ age + yearsmarried + religiousness + rating
Model 2: yaffairs ~ 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
*******************************************
卡方值不显著，表明四个预测变量的新模型与九个完整预测变量的模型拟合程度一样好



13.2.1 解释模型参数

先看看回归系数
coef(fit.reduced)
(Intercept)           age  yearsmarried religiousness        rating 
1.93083017   -0.03527112    0.10062274   -0.32902386   -0.46136144 

由于对数优势比解释性差，可对结果进行指数化
exp(coef(fit.reduced))

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

使用predict方法对婚姻平分对婚外情概率的影响
testdata=data.frame(rating=c(1,2,3,4,5),age=mean(Affairs$age),yearsmarried=mean(Affairs$yearsmarried),religiousness=mean(Affairs$religiousness))
testdata

##接下来，使用测试集预测相应的概率
testdata$prob=predict(fit.reduced,newdata = testdata,type = "response")
testdata


13.2.4 扩展
'''
R中扩展的Logistic回归和变种如下所示：
1.稳健的Logistic回归，robust包中的glmRob()函数可用来拟合稳健的广义线性模型，包括稳健的Logistic回归，当拟合Logistic回归模型数据出现利群点和强影响点时，稳健Logistic回归便可排上永昌
2.多项分布回归  若相应变量包含两个或者两个以上的无序类别（比如，已婚/寡居/离婚），便可使用mlogit包中的mlogit()函数拟合多项回归
3.序数Logistic回归  若相应变量时一组有序的类别（比如，差，良，好） 便可使用rms包中的lrm() 函数来拟合
'''

13.3 泊松分布
当通过一系列连续性或者类别型预测变量来预测计数型结果变量时，泊松分布式一个非常好的工具

opar=par(no.readonly = TRUE)
par(mfrow=c(1,2))
attach(breslow.cat)
hist(sumY,breaks=20)
boxplot(sumY~Trt)
par(opar)

#####拟合模型
fit=glm(sumY~.,dta=breslow,family = poisson())
summary(fit)

13.3 扩展

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