R語言怎樣做Logistic回歸

當(dāng)因變量為二值型結(jié)果變量，自變量包括連續(xù)型和類別型的數(shù)據(jù)時(shí)，Logistic回歸是一個(gè)非常常用的工具。比如今天的例子中用到的婚外情數(shù)據(jù) “Fair's Affairs”。因變量是是否有過婚外情，自變量有8個(gè)，分別是

• 性別
• 年齡
• 婚齡
• 是否有小孩
• 宗教信仰程度 （5分制，1表示反對，5表示非常信仰）
• 學(xué)歷
• 職業(yè) （逆向編號(hào)的戈登7種分類）這個(gè)是啥意思？）
• 對婚姻的自我評分

因變量y是出軌次數(shù)，我們將其轉(zhuǎn)換成二值型，出軌次數(shù)大于等于1賦值為1，相反賦值為0

下面開始實(shí)際操作

這個(gè)數(shù)據(jù)集來自R語言包AER，如果要用這個(gè)數(shù)據(jù)集需要先安裝這個(gè)包

install.packages("AER")

然后使用data()函數(shù)獲取這個(gè)數(shù)據(jù)集

data(Affairs,package = "AER")

然后就可以在環(huán)境的窗口里看到如下

這個(gè)數(shù)據(jù)集總共有601個(gè)觀察值，總共9個(gè)變量

df<-Affairs
df$ynaffairs<-ifelse(df$affairs>0,1,0)
table(df$ynaffairs)
df$ynaffairs<-factor(df$ynaffairs,
                     levels = c(0,1),
                     labels = c("No","Yes"))table
table(df$ynaffairs)

擬合模型用到的是glm()函數(shù)

fit.full<-glm(ynaffairs~gender+age+yearsmarried+
                children+religiousness+education+occupation+rating,
              data=df,family = binomial())

通過summary()函數(shù)查看擬合結(jié)果

summary(fit.full)

 

根據(jù)回歸系數(shù)的P值可以看到 性別、是否有孩子、學(xué)歷、職業(yè)對方程的貢獻(xiàn)都不顯著。去除這些變量重新擬合模型

fit.reduced<-glm(ynaffairs~age+yearsmarried+
                religiousness+rating,
              data=df,family = binomial())

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

 

可以看到結(jié)果中p值等于0.2108大于0.05，表明四個(gè)變量和9個(gè)變量的模型你和程度沒有差別

構(gòu)造一個(gè)測試集

testdata<-data.frame(rating=c(1,2,3,4,5),
                     age=mean(df$age),
                     yearsmarried=mean(df$yearsmarried),
                     religiousness=mean(df$religiousness)

預(yù)測概率

testdata$prob<-predict(fit.reduced,newdata = testdata,
                       type = "response")

library(ggplot2)
ggplot(testdata,aes(x=rating,y=prob))+
  geom_col(aes(fill=factor(rating)),show.legend = F)+
  geom_label(aes(label=round(prob,2)))+
  theme_bw()

 

從這些結(jié)果可以看到，當(dāng)婚姻評分從1（很不幸福）變?yōu)?（非常幸福）時(shí)，婚外情概率從0.53降低到了0.15。模型的預(yù)測結(jié)果和我們的經(jīng)驗(yàn)還挺符合的