title: "第一章-二分類logistic回歸預(yù)測(cè)模型"
author: "大鵬統(tǒng)計(jì)工作室R語(yǔ)言"
date: "2020/4/28"

案例1

二分類logsitic回歸

?以AER包中的數(shù)據(jù)框Affairs為例，我們將通過(guò)探究婚外情的數(shù)據(jù)來(lái)闡述Logistic回歸分析的過(guò)程。

?婚外情數(shù)據(jù)取自1969年（Psychology Today）所做一個(gè)非常有代表性的調(diào)查。

該數(shù)據(jù)從601個(gè)參與者收集了9個(gè)變量，包括一年來(lái)婚外私通的頻率以及參與者性別、年齡、婚齡、是否有小孩、宗教信仰程度（5分制，1分表示反對(duì)，5分表示非常信仰）、學(xué)歷、職業(yè)，還有對(duì)婚姻自我評(píng)分（5 分制，1表示非常不幸福，5表示非常幸福）。

使用軟件版本信息R 3.6.1

R version 3.6.1 (2019-07-05)

統(tǒng)計(jì)分析代做，聯(lián)系微信，1265950003

課程購(gòu)買，聯(lián)系微信，1265950003

下載.png

安裝本小節(jié)需要的R包，如已安裝請(qǐng)忽略

#install.packages("AER")

#調(diào)用AER包中的Affairs數(shù)據(jù),并將Affairs數(shù)據(jù)傳遞給mydata
data(Affairs, package="AER")
mydata<-Affairs
#查看數(shù)據(jù)集mydata中各變量基本統(tǒng)計(jì)信息
summary(mydata)

image.png

#查看因變量affairs的分布
table(mydata$affairs)

image.png

其中0為過(guò)去幾年中沒(méi)有婚外性行為，1、2、3分別為有1、2和3次婚外性行為；當(dāng)婚外性行為的次數(shù)在4-10之間時(shí)，affairs定義為7；當(dāng)婚外性行為頻率較高，為每周，每月或每日時(shí)，affairs定義為12

#將婚外情次數(shù)大于0的修改成1
mydata$affairs[mydata$affairs> 0] <-1
table(mydata$affairs)

image.png

0表示沒(méi)有婚外情，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())
#查看模型結(jié)果
summary(model)$coefficients

image.png

存在統(tǒng)計(jì)學(xué)意義的變量age、yearsmarried、religiousness、rating

年齡、婚齡、宗教信仰程、對(duì)婚姻自我評(píng)分

#第一列Estimate  → coef
coef<-summary(model)$coefficients[,1]
#第二列Std. Error  → se
se<-summary(model)$coefficients[,2]
#第四列Pr(>|z|)  → pvalue
pvalue<-summary(model)$coefficients[,4]

第一種方法計(jì)算OR及其可信區(qū)間

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

存在統(tǒng)計(jì)學(xué)意義的變量age、yearsmarried、religiousness、rating

年齡、婚齡、宗教信仰程、對(duì)婚姻自我評(píng)分

第二種方法計(jì)算OR及其可信區(qū)間

#exp(cbind("OR"=coef(model),confint(model)))
confint(model)
exp(confint(model))

image.png

image.png

#####逐步回歸
step(model,direction ="backward")

image.png

AIC赤池信息準(zhǔn)則

衡量統(tǒng)計(jì)模型擬合優(yōu)良性(Goodness of fit)的一種標(biāo)準(zhǔn)，通常最小AIC值相對(duì)應(yīng)的模型作為最優(yōu)對(duì)象

#根據(jù)逐步回歸結(jié)果重新擬合
model2<-glm(formula = affairs ~ gender + age + yearsmarried + religiousness + rating, family = binomial(), data = mydata)
#查看模型結(jié)果
summary(model2)$coefficients

image.png

性別、年齡、婚齡、宗教信仰程、對(duì)婚姻自我評(píng)分

性別不存在統(tǒng)一意義

模型比較方法1

#這個(gè)方法需要是嵌套模型，模型變量之間是包含關(guān)系  
anova(model,model2,test="Chisq")

image.png

模型比較方法2

#任何模型 
as.matrix(AIC(model,model2))

image.png

生成一個(gè)新的數(shù)據(jù)testdata,評(píng)價(jià)男性或女性對(duì)婚姻自我評(píng)分rating取不同取值時(shí)，婚外情概率

婚姻自我評(píng)分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

結(jié)果顯示男性對(duì)婚姻非常不滿意的出軌率是0.5818455

男性對(duì)婚姻非常滿意的出軌率是0.1767546

其他結(jié)果自行讀取

生成一個(gè)新的數(shù)據(jù)testdata2,評(píng)價(jià)男性或女性宗教信仰程度取不同取值時(shí)，婚外情概率

宗教信仰程度（5分制，1分表示反對(duì)，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

結(jié)果顯示男性宗教信仰程度非常不信仰的出軌率是0.4141053

男性宗教信仰程度非常信仰的出軌率是0.1603550

其他結(jié)果自行讀取

清除工作空間

rm(list=ls())

下載.png