R语言统计分析——Logistic回归

参考资料：R语言实战【第2版】

        当通过一系列连续型和/或类别型变量来预测二值型结果变量时，Logistic回归是一个非常有用的工具。以AER包中的数据框Affairs为例。首次使用该数据前，需安装AER包（install.packages("AER")）。

Affairs数据框详细内容见下图：

# 安装AER包
install.packages("AER")
# 加载AER包
library(AER)
# 加载Affairs数据
data(Affairs,package="AER")
# 查看数据基本信息
summary(Affairs)
# 查看出轨的次数
table(Affairs$affairs)

        虽然对婚姻不忠的次数被记录下来，但我们感兴趣的是二值型结果（是否有过婚外情）。下面将affairs数据转化为二值型因子ynaffair。

# 将响应变量转换为二值型因子变量
Affairs$ynaffair[Affairs$affairs>0]<-1
Affairs$ynaffair[Affairs$affairs==0]<-0
Affairs$ynaffair<-factor(Affairs$ynaffair,
                            level=c(0,1),
                            labels=c("no","yes"))
# 查看转换后的数据
table(Affairs$ynaffair)

        下面对二值型因子作为Logistic回归的结果变量进行广义线性拟合：

# glm拟合
fit.full<-glm(ynaffair~gender+age+yearsmarried+children+religiousness+education+occupation+rating,
              data=Affairs, family=binomial(link="logit"))
# 查看拟合结果
summary(fit.full)

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

# 去掉不显著的自变量重新拟合
fit.reduced<-glm(ynaffair~age+yearsmarried+religiousness+rating,
                 data=Affairs,family=binomial(link="logit"))
# 查看拟合结果
summary(fit.reduced)

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

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

        结果显示，卡方值不显著（p=0.21），说明四个预测变量的新模型和九个完整预测边阿玲的模型拟合程度一样好。

1、解释模型参数

        先看看回归系数：

coef(fit.reduced)

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

        由于对数优势比解释性差，我们可以对结果进行指数化：

exp(coef(fit.reduced))

        由商量结果可看出：婚龄增加一年，婚外情的优势比将乘以1.106（保持年龄、宗教信仰和婚姻评定不变）；相反，年龄增加一岁，婚外情的的优势比则乘以0.965。因此，随着婚龄的增加和年龄、宗教信仰与婚姻评分的降低，婚外情优势比将上升。因为预测变量不能等于0，截距项在此处没有什么特定含义。

        如果有需要，我们还可以使用confint()函数获取系数的置信区间。

        最后，预测变量一单位的变化可能并不是我们最想关注的。对于二值型Logistic回归，某预测变量n单位的变化引起的较高值上优势比的变化为exp(βj)^n，它反映的信息可能更为重要。比如，保持其他预测变量不变，婚龄增加一年，婚外情的优势比将乘以1.106，而如果婚龄增加10年，优势比将乘以1.106^10，即2.7。

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
# 使用本数据集对数据进行预测
testdata$prob<-predict(fit.reduced,newdata=testdata,type="response")
testdata

        从这些结果可以看到，当婚姻评分从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(fit.reduced,newdata=testdata,type="response")
testdata

        此处可以看到，当其他变量不变，年龄从17增加到57时，婚外情的概率将从0.34降低到0.11。利用该方法，我们可探究每一个预测变量对结果概率的影响。

3、过度离势

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

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

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

=残差偏差 / 残差自由度

比1大很多，我们便可认为存在过度离势。回到本例中

        此比值非常接近于1，表明没有过度离势。

        我们还可以对过度离势进行检验。为此，我们需要拟合模型两次，第一次使用family=binomial，第二次使用family=quasibinomial。假设第一次glm()返回对象记为fit，第二次返回对象记为fit.od，那么：

fit<-glm(ynaffair~age+yearsmarried+religiousness+rating,
            data=Affairs,family=binomial(link="logit"))
fit.od<-glm(ynaffair~age+yearsmarried+religiousness+rating,
            data=Affairs,family=quasibinomial(link="logit"))
pchisq(summary(fit.od)$dispersion*fit$df.residual,
       fit$df.residual,lower=F)

提供的p值即可对零假设H0：=1与备择假设H1： ≠1进行检验。若p很小（小于0.05），我们便可拒绝零假设。

本例中，p值为0.34，大于0.05，可以确定不存在过度离势。