1 Logistic 回归
当通过一系列连续型和/或类别型预测变量来预测二值型结果变量时，Logistic回归是一个非
常有用的工具。以AER包中的数据框Affairs为例，我们将通过探究婚外情的数据来阐述Logistic回归的过程。
婚外情数据即著名的“Fair’s Affairs”，取自于1969年《今日心理》（Psychology Today）所做的一个非常有代表性的调查，该数据从601个参与者身上收集了9个变量，包括一年来婚外私通的频率以及参与者性别、年龄、婚龄、是否有小孩、宗教信仰程度（5分制，1分表示反对，5分表示非常信仰）、学历、职业（逆向编号的戈登7种分类），还有对婚姻的自我评分（5分制，1表示非常不幸福，5表示非常幸福）。
install.packages("AER")
data(Affairs, package="AER")
summary(Affairs)
table(Affairs$affairs) #统计婚外情次数的人数
从这些统计信息可以看到，315/601=52%的调查对象是女性，430/601=72%的人有孩子，样本年龄的中位数为32岁。对于响应变量，72%的调查对象表示过去一年中没有婚外情（451/601），而婚外偷情的最多次数为12（占了6%）。
虽然这些婚姻的轻率举动次数被记录下来，但此处我们感兴趣的是二值型结果（有过一次婚
外情/没有过婚外情）。按照如下代码，可将affairs转化为二值型因子ynaffair：
Affairs$ynaffair[Affairs$affairs > 0] <- 1 #赋值
Affairs$ynaffair[Affairs$affairs == 0] <- 0
Affairs$ynaffair <- factor(Affairs$ynaffair,
levels=c(0,1), #值的结果是0和1
labels=c("No","Yes")) #给二值型因子贴上标签
table(Affairs$ynaffair)
该二值型因子现可作为Logistic回归的结果变量：
fit.full <- glm(ynaffair ~ gender + age + yearsmarried + children +
religiousness + education + occupation +rating,
data=Affairs, family=binomial()) #family是设置概率分布
summary(fit.full)
结果分析：从回归系数的p值（最后一栏）可以看到，性别、是否有孩子、学历和职业对方程的贡献都不显著（后面没有星星），也就是我们无法拒绝参数为0的假设。去除这些变量重新拟合模型，检验新模型是否拟合得好呢？
fit.reduced <- glm(ynaffair ~ age + yearsmarried + religiousness +
rating, data=Affairs, family=binomial())
summary(fit.reduced)
结果分析：新模型的每个回归系数都非常显著（p<0.05）。由于两模型嵌套（fit.reduced是fit.full的一个子集），你可以使用anova()函数对它们进行比较，对于广义线性回归，可用卡方检验。
anova(fit.reduced, fit.full, test="Chisq")
结果分析：结果的卡方值不显著（p=0.21），表明四个预测变量的新模型与九个完整预测变量的模型拟合程度一样好。这使得你更加坚信添加性别、孩子、学历和职业变量不会显著提高方程的预测精度，因此可以依据更简单的模型进行解释。
解释模型参数：先看看回归系数
coef(fit.reduced)
对每个观测值各分解项的值取指数，将结果转化为原始尺度：
exp(coef(fit.reduced))
结果分析：可以看到婚龄增加一年，婚外情的优势比将乘以1.106（保持年龄、宗教信仰和婚姻评定不变）；相反，年龄增加一岁，婚外情的的优势比则乘以0.965。因此，随着婚龄的增加和年龄、宗教信仰与婚姻评分的降低，婚外情优势比将上升。
1.1 评价预测变量对结果概率的影响
使用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")
#predict：从拟合的稳健广义线性模型（GLM）对象中获取预测值，并可选择估计这些预测值的标准误差。Newdata：可选地，一个数据框架，在其中寻找用于预测的变量，如果省略，则使用拟合的线性预测器。
testdata
结果分析：从这些结果可以看到，当婚姻评分从1（很不幸福）变为5（非常幸福）时，婚外情概率从0.53降低到了0.15（假定年龄、婚龄和宗教信仰不变）。下面我们再看看年龄的影响：
testdata <- data.frame(rating=mean(Affairs$rating),
age=seq(17, 57, 10), #从17岁到57岁，每隔10岁生成一个
yearsmarried=mean(Affairs$yearsmarried),
religiousness=mean(Affairs$religiousness))
testdata
testdata$prob <- predict(fit.reduced, newdata=testdata, type="response")
Testdata
结果分析：此处可以看到，当其他变量不变，年龄从17增加到57时，婚外情的概率将从0.34降低到0.11。