データプロセシング

logistic regressionに使えるようにデータを調整します。

データセットをtrainigとtestに分割します。割合は、８：２とします。

# trainとtestの列を作る
splits <- c(rep("train", 0.8*nrow(df)),
            rep("test", 0.2*nrow(df)))
# 元のデータセットに1列加える
# sampleでランダムにサンプルされる
df <- df %>% mutate(splits=sample(splits))

#データセットを分割
default_train <- df %>% 
  filter(splits=="train") %>% 
  select(-splits)

default_test <- df %>% 
  filter(splits=="test") %>% 
  select(-splits)

logistic regressionを実施

glmでロジステック回帰が使えます。 使い方は以下の通りです。

glm(式, family = binomial, data = トレーニングデータ)

defaultがyで、その他がxです。defaultを予測します。 defaultがfactorで、factor 1がNo、factor 2がYesです。

glmは、factorのレベルが大きい方に所属する場合のlog oddsをだしてくれます。

levels(df$default)[1]

## [1] "No"

levels(df$default)[2]

## [1] "Yes"

logistic <- glm(default ~ .,
                family = binomial,
                data = default_train)

予測

モデルができたので、predictを使って予測したいと思います。

predict（モデル,　newdata = 予測したいデータ）

log_prediction <- predict(logistic,
                newdata = default_test)

予測値を見てみましょう。log oddsなので解釈がむずかしいです。

head(log_prediction)

##          1          2          3          4          5          6 
##  -4.572440  -6.215661 -10.587827  -9.267444  -7.116755  -8.357313

各係数×各変数の合計が、logoddsです。 log odds = ax1 + ax2 + ax3 …

log oddsを、確率に変換します。 確率 = 1 / ( 1 + e^(-logodds) )

probability <- 1/(1+exp(-log_prediction))

head(probability)

##            1            2            3            4            5            6 
## 1.022705e-02 1.993911e-03 2.522052e-05 9.444074e-05 8.107379e-04 2.346190e-04

1.022e-02は、1.022×10^-2を示しますが、わかりずらいです。指数表示をやめるおまじない、scipen=100を入れます。

options(scipen=100)
head(probability)

##             1             2             3             4             5 
## 0.01022704628 0.00199391077 0.00002522052 0.00009444074 0.00081073786 
##             6 
## 0.00023461896

実は、こんなにめんどくさい変換をしなくても自動で確率に変換する機能があります。 predictの中に、type = “response”を入れるだけです。

log_prediction <- predict(logistic,
                newdata = default_test,
                type = "response" )
head(probability)

##             1             2             3             4             5 
## 0.01022704628 0.00199391077 0.00002522052 0.00009444074 0.00081073786 
##             6 
## 0.00023461896

精度

予想の精度を見てみます。

まず、予想確率をカットオフ値0.5として、Yes、Noに分類したいと思います。 0.5以上はYes、それ未満はNoです。

pred_class <- ifelse(log_prediction>=0.5,"Yes","No")

confusion matrixを作ります。

table(True_class = default_test$default,
      prediction = pred_class)

##           prediction
## True_class   No  Yes
##        No  1920    6
##        Yes   52   22

TP = 22
FP = 6
TN = 1920
FN = 52

accuracy = (TP + TN)/(TP + TN + FP + FN)
f1score = TP/(TP+0.5*(FP + FN))
sensitivity = TP/(TP + FN)
specificity = TN/(TN + FP)
false_positive_rate = FP/(FP + TN)
positive_predictive_value = TP/(TP + FP)
negative_predictive_value = TN/(TN + FN)

cat(paste("accuracy:",accuracy,"\n",
          "f1score:",f1score,"\n",
          "sensitivity:",sensitivity,"\n",
          "specificity:",specificity,"\n",
          "false_positive_rate:",false_positive_rate,"\n",
          "positive_predictive_value:",positive_predictive_value,"\n",
          "negative_predictive_value:",negative_predictive_value)) 

## accuracy: 0.971 
##  f1score: 0.431372549019608 
##  sensitivity: 0.297297297297297 
##  specificity: 0.996884735202492 
##  false_positive_rate: 0.00311526479750779 
##  positive_predictive_value: 0.785714285714286 
##  negative_predictive_value: 0.973630831643002

係数の解釈

coefficientを見てみます。

coef(logistic)

##      (Intercept)       studentYes          balance           income 
## -10.733745051171  -0.661322126988   0.005591702312   0.000004984331

studentだけYesがついています。これは、Student変数がNoかyesかのカテゴリだったため、自動的にダミー変数にされて計算されています。Noが基準となるため、yesの場合だけ１となって計算されています。

Noのときのdefaultのオッズと、Yesの場合のdefaultのオッズがどれくらいの比か、というのを計算することができます。

exp(cbind(coef(logistic),confint(logistic)))

##                                    2.5 %        97.5 %
## (Intercept) 0.00002179685 0.000007300467 0.00006129732
## studentYes  0.51616844293 0.305089474609 0.87570499644
## balance     1.00560736506 1.005122689800 1.00612469199
## income      1.00000498434 0.999986923707 1.00002307840

student Yesのよこの数字が0.51となっています。これは、student Noのときを基準としたときのstudent Yesオッズ比です。student yesの場合は、student noのときに比べて、0.51のオッズ比です。比が１より低いので、yesの場合のほうがdefaultになる確率は低いということが分かります。

各係数と各変数の合計が、logoddsです。

log odds = a1 + a2 × X1 + a3 × X2 + …

確率 = 1 / ( 1 + e^(-log odds) )です。

例えば、以下の場合はdefaultになる可能性はいくつでしょうか。

studentがYes balanceが5000 incomeが20000

logodds <- coef(logistic)[1] + 1*coef(logistic)[2] + 5000*coef(logistic)[3] + 20000*coef(logistic)[4]

1 / (1 + exp(-logodds))

## (Intercept) 
##   0.9999999

logistic regressionと、linear discriminant analysisの使い分け。

logistic regression：2カテゴリー

LDA：3以上カテゴリー、サンプル数が少ない時

LDAでカテゴリー予想 in R