广义线性模型

回归适⽤于⼆值响应变量，连接函数为logit函数，概率分布为⼆项分布：glm(Y ~ X1 + X2 + X3, family=binomial(link="logit"), data=data)
Logistics回归模型的表达形式为：
模型解读：为暴露于某种状态下的结局概率，是⼀种变量变换⽅式，表⽰对进⾏变换，为偏回归系数，表⽰在其他⾃变量不变的条件下，每变化⼀个单位的估计值。

Logistics回归是通过最⼤似然估计求解常数项和偏回归系数，基本思想时当从总体中随机抽取n个样本后，最合理的参数估计量应该使得这n个样本观测值的概率最⼤。

最⼤似然法的基本思想是先建⽴似然函数与对数似然函数，再通过使对数似然函数最⼤求解相应的参数值，所得到的估计值称为参数的最⼤似然估计值。

1. 数据准备（类别型变量进⾏0/1量化）
⾸先，我们选⽤AER包中的Affairs数据集来构建Logistics回归模型，这个数据集记录了⼀组婚外情数据，其中包括参与者性别、年龄、婚龄、是否有⼩孩、宗教信仰程度（5分制，1表⽰反对，5表⽰⾮常信仰）、学历、职业和婚姻的⾃我评分（5分制，1表⽰⾮常不幸福，5表⽰⾮常幸福）。

在使⽤数据集之前，载⼊AER包
> library(AER)
> data(Affairs,package="AER")
对于这个数据集，我们关注是否出轨，即这是⼀个⼆值型结果（出轨过/从未出轨）。

因此，我们接下来将'affaris'特征转化为⼆值型因⼦'ynaffair'，该⼆值型因⼦即可以作为Logistic回归的结果变量。

> 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"))
2. Logistics模型构建：
> myfit <- glm(ynaffair ~ gender + age + yearsmarried + children + religiousness + education + occupation + rating, data=Affairs, family=binomial())
> summary(myfit)
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
从回归系数的p值（最后⼀栏）可以看到，性别、是否有孩⼦、学历和职业对⽅程的贡献都不显著（p值较⼤）。

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

> myfit_reduced <- glm(ynaffair ~ age + yearsmarried + religiousness + rating, data=Affairs, family=binomial())
> summary(myfit_reduced)
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 ***
从最后⼀列可以看到，新模型的每个回归系数都⾮常显著（p<0.05）
3. 多模型⽐较
由于这个两模型是嵌套关系，所以我们可以使⽤anova()函数对它们进⾏⽐较。

⽽对于⼴义线性回归，可⽤卡⽅检验进⾏判断。

> anova(myfit,myfit_reduced,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 59
2 609.51
2 596 615.36 -4 -5.8474 0.2108
检验结果的卡⽅值不显著（p=0.21），这表明，四个预测变量的新模型和九个预测变量的旧模型模拟程度其实差不多。

⼀般情况下我们更倾向于选择四个预测变量的新模型。

4. 评价预测变量对结果概率的影响（predict函数）
既然已经有了模型，那么给定⼀个⽤户⾏为，预测其出轨的概率。

可以使⽤predict()函数进⾏预测，以下代码⽣成了⼀些虚拟数据并进⾏了预
测，主要评测婚姻评分对出轨概率的影响。

# 构建dataframe
> 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(myfit_reduced, newdata=testdata, type="response")
> testdata
rating age yearsmarried religiousness prob
1 1 32.4875
2 8.177696 3.11647
3 0.5302296
2 2 32.48752 8.177696 3.11647
3 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.17769
6 3.116473 0.1513079
从结果可以看到，当婚姻评分从1分（很不幸福）增加到5分（⾮常幸福）时，出轨概率从0.53降低⾄0.15。

> testdata <- data.frame(rating=mean(Affairs$rating), age=seq(17, 57, 10), yearsmarried=mean(Affairs$yearsmarried), religiousness=mean(Affairs$religiousness))
> testdata$prob <- predict(myfit_reduced, newdata=testdata, type="response")
> testdata
rating age yearsmarried religiousness prob
1 3.93178 17 8.177696 3.116473 0.3350834
2 3.93178 27 8.177696 3.11647
3 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.17769
6 3.116473 0.1094738
在年龄的影响⽅⾯，可以看出，当其他变量不变、年龄从17增加到57时，出轨的概率从0.34降低⾄0.11。

5. 过度离势
过度离势指的是观测到的响应变量⽅差⼤于期望的⼆项分布⽅差，过度离势会导致奇异的标准误检验和不精确的显著性检验。

当出现过度离势时，仍可使⽤glm()函数拟合Logistic回归，但此时需要将⼆项分布改为类⼆项分布。

检测过度离势的⼀种⽅法是⽐较⼆项分布模型的残差偏差与残差⾃由度（在四变量模型中分别是615.36和596），如果两者的⽐值⽐1⼤很多，则可以认为存在过度离势。

6. Logistic回归模型的进阶学习
扩展的Logistic回归和变种如下所⽰：
稳健Logistic回归：当拟合Logistic回归模型数据出现离群点和强影响点时，robust包中的glmRob()函数可⽤来拟合稳健的⼴义线性模型
多项分布回归：当响应变量包含两个以上的⽆序类别（例如已婚、寡居、离婚）时，可使⽤mlogit包中的mlogit()函数拟合多项Logistic回归
多项分布回归：当响应变量包含两个以上的⽆序类别（例如已婚、寡居、离婚）时，可使⽤mlogit包中的mlogit()函数拟合多项Logistic回归
序数Logistic回归：当响应变量是⼀组有序的类别（例如信⽤为好/良/差）时，可使⽤rms包中的lrm()函数拟合序数Logistic回归。

三、泊松回归
泊松回归的⾃变量是连续性或类别型变量，因变量是计数型的变量。

泊松回归因变量通常局限在⼀个固定长度时间段内进⾏测量（如过去⼀年
交通事故数），整个观测集中时间长度都是不变的。

Poisson回归主要有两个假设，⾸先，具有相同特征和同时的不同对象的⼈时风险是同质
的，其次，当样本量越来越⼤时，频数的均数趋近于⽅差。

调⽤模型的公式如下：
> myfit <- glm(y~x1+x2+...+xn,data=,family = poisson)
1. 数据准备
robust包中Breslow癫痫数据记录了治疗初期⼋周内，抗癫痫药物对癫痫发病数的影响。

响应变量为sumY（随机化后⼋周内癫痫发病数），预测变量为治疗条件（Trt）、年龄（Age）和前⼋周内的基础癫痫发病数（Base），在这个数据集中，我们感兴趣的是药物治疗能否减少癫痫发病数。

> data(breslow.dat, package="robust")
接下来我们考察响应变量，基础和随机化后的癫痫发病数都有很⾼的偏度，从图中可以清楚地看到因变量的偏倚特性以及可能的离群点。

> library(ggplot2)
> p1 <- ggplot(breslow.dat,aes(sumY)) + geom_histogram(color='green')
> p2 <- ggplot(breslow.dat,aes(Trt,sumY)) + geom_bar()
# 组合图形
> library(gridExtra)
> grid.arrange(p1,p2,nrow=1)
从图中可以清楚的看到因变量的偏移特性及可能的离群点。

药物治疗下癫痫的发病数似乎变⼩，且⽅差也变⼩了。

2. 构建模型
> myfit <- glm(sumY ~ Base + Age + Trt, data=breslow.dat, family=poisson())
> summary(myfit)
3. 过度离势
在泊松分布中同样可能存在过度离势，泊松分布的⽅差和均值相等，当响应变量观测的⽅差⽐依据泊松分布预测的⽅差⼤时，可能存在过度离势。

> summary(myfit)
(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
与Logistic回归类似，但此例中残差偏差（559.44）与残差⾃由度（55）的⽐例为10.17，这⼀值远⼤于1，故存在过度离势。

在这种情况下，可以考虑使⽤类泊松回归替代泊松回归（family=quasipoisson()）。

除此之外，我们还可以使⽤qcc包来检验泊松模型是否存在过度离势。

> library(qcc)
> qcc.overdispersion.test(breslow.dat$sumY, type="poisson")
4. 泊松回归模型的进阶学习：
时间段变化的泊松回归：如果每个观测的时间段长短不同，则因变量需要从计数更改为⽐率，即发⽣次数除以观测时间，使⽤glm()的offset选项即
可，offset=log(time)，其中time是每名病⼈的观测时间
零膨胀的泊松回归：在⼀个数据集中，0计数的数⽬时常⽐⽤泊松模型预测得多。

在Affairs数据集中，很有可能有⼀群从未出轨的对象，他们称为数据集的结构零值。

可以⽤零膨胀的泊松回归来分析这种情况，它将同时拟合两个模型，⼀个⽤来预测哪些⼈⼜会出轨，另⼀个⽤来预测排除了结构零值之后的调查对象会出轨多少次。

可以将其看作是Logistic回归（预测结构零值）和泊松回归（预测⽆结构零值观测的计数）的组合，pscl包的zeroinfl()函数可做零膨胀泊松回归
稳健泊松回归：robust包中的glmRob()函数可以拟合稳健⼴义线性模型，包括稳健泊松回归。