title: "第一章-二分类logistic回归预测模型"
author: "大鹏统计工作室R语言"
date: "2020/4/28"
案例1
二分类logsitic回归
➢以AER包中的数据框Affairs为例，我们将通过探究婚外情的数据来阐述Logistic回归分析的过程。
➢婚外情数据取自1969年（Psychology Today）所做一个非常有代表性的调查。
该数据从601个参与者收集了9个变量，包括一年来婚外私通的频率以及参与者性别、年龄、婚龄、是否有小孩、宗教信仰程度（5分制，1分表示反对，5分表示非常信仰）、学历、职业，还有对婚姻自我评分（5 分制，1表示非常不幸福，5表示非常幸福）。
使用软件版本信息R 3.6.1
R version 3.6.1 (2019-07-05)
统计分析代做，联系微信，1265950003
课程购买，联系微信，1265950003
下载.png
安装本小节需要的R包，如已安装请忽略
#install.packages("AER")
#调用AER包中的Affairs数据,并将Affairs数据传递给mydata
data(Affairs, package="AER")
mydata<-Affairs
#查看数据集mydata中各变量基本统计信息
summary(mydata)
image.png
#查看因变量affairs的分布
table(mydata$affairs)
image.png
其中0为过去几年中没有婚外性行为，1、2、3分别为有1、2和3次婚外性行为；当婚外性行为的次数在4-10之间时，affairs定义为7；当婚外性行为频率较高，为每周，每月或每日时，affairs定义为12
#将婚外情次数大于0的修改成1
mydata$affairs[mydata$affairs> 0] <-1
table(mydata$affairs)
image.png
0表示没有婚外情，1表示存在婚外情
#将affairs处理成分类变量
mydata$affairs<-factor(mydata$affairs, levels=c(0,1), labels=c("No","Yes"))
table(mydata$affairs)
#拟合模型，模型名称model
model<-glm(affairs~ gender + age + yearsmarried+ children +religiousness + education + occupation +rating, data=mydata, family=binomial())
#查看模型结果
summary(model)$coefficients
image.png
存在统计学意义的变量age、yearsmarried、religiousness、rating
年龄、婚龄、宗教信仰程、对婚姻自我评分
#第一列Estimate → coef
coef<-summary(model)$coefficients[,1]
#第二列Std. Error → se
se<-summary(model)$coefficients[,2]
#第四列Pr(>|z|) → pvalue
pvalue<-summary(model)$coefficients[,4]
第一种方法计算OR及其可信区间
Results<-cbind(exp(coef),exp(coef-1.96*se),exp(coef+1.96*se),pvalue)
dimnames(Results)[[2]]<-c("OR","LL","UL","P value")
Results
image.png
存在统计学意义的变量age、yearsmarried、religiousness、rating
年龄、婚龄、宗教信仰程、对婚姻自我评分
第二种方法计算OR及其可信区间
#exp(cbind("OR"=coef(model),confint(model)))
confint(model)
exp(confint(model))
image.png
image.png
#####逐步回归
step(model,direction ="backward")
image.png
AIC赤池信息准则
衡量统计模型拟合优良性(Goodness of fit)的一种标准，通常最小AIC值相对应的模型作为最优对象
#根据逐步回归结果重新拟合
model2<-glm(formula = affairs ~ gender + age + yearsmarried + religiousness + rating, family = binomial(), data = mydata)
#查看模型结果
summary(model2)$coefficients
image.png
性别、年龄、婚龄、宗教信仰程、对婚姻自我评分
性别不存在统一意义
模型比较方法1
#这个方法需要是嵌套模型，模型变量之间是包含关系
anova(model,model2,test="Chisq")
image.png
模型比较方法2
#任何模型
as.matrix(AIC(model,model2))
image.png
生成一个新的数据testdata,评价男性或女性对婚姻自我评分rating取不同取值时，婚外情概率
婚姻自我评分rating(5 分制，1表示非常不幸福，5表示非常幸福)
testdata <- data.frame(rating =rep(c(1, 2, 3, 4, 5),2),
gender=rep(c("male","female"),each = 5),
age = mean(mydata$age),
yearsmarried = mean(mydata$yearsmarried),
religiousness = mean(mydata$religiousness))
testdata$prob <- predict(model2, newdata=testdata, type="response")
testdata<-as.matrix(testdata)
testdata
image.png
结果显示男性对婚姻非常不满意的出轨率是0.5818455
男性对婚姻非常满意的出轨率是0.1767546
其他结果自行读取
生成一个新的数据testdata2,评价男性或女性宗教信仰程度取不同取值时，婚外情概率
宗教信仰程度（5分制，1分表示反对，5分表示非常信仰）
testdata2 <- data.frame(rating = mean(mydata$rating),
gender=rep(c("male","female"),each = 5),
age = mean(mydata$age),
yearsmarried = mean(mydata$yearsmarried),
religiousness =rep(c(1, 2, 3, 4, 5),2))
testdata2$prob <- predict(model2, newdata=testdata2, type="response")
testdata2<-as.matrix(testdata2)
testdata2
image.png
结果显示男性宗教信仰程度非常不信仰的出轨率是0.4141053
男性宗教信仰程度非常信仰的出轨率是0.1603550
其他结果自行读取
清除工作空间
rm(list=ls())
下载.png
