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

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

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

下面開始實際操作

這個數據集來自R語言包AER，如果要用這個數據集需要先安裝這個包

install.packages("AER")

然后使用data()函數獲取這個數據集

data(Affairs,package = "AER")

然后就可以在環境的窗口里看到如下

這個數據集總共有601個觀察值，總共9個變量

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()函數

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

通過summary()函數查看擬合結果

summary(fit.full)

 

根據回歸系數的P值可以看到 性別、是否有孩子、學歷、職業對方程的貢獻都不顯著。去除這些變量重新擬合模型

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

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

 

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

構造一個測試集

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

預測概率

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()

 

從這些結果可以看到，當婚姻評分從1（很不幸福）變為5（非常幸福）時，婚外情概率從0.53降低到了0.15。模型的預測結果和我們的經驗還挺符合的