分類（ロジスティック回帰２）
東京国際大学 データサイエンス教育研究所 竹田 恒
2024-09-15
次のデータを用いて新車購入の確率を予測する。
- NAデータは含まない。
- 新車購入の有無は，既に2値（バイナリ）化されている。
1 データ
- お客様番号
- 性別（男性，女性）
- 年齢
- 年収（US＄）
- 購入判断（1：購入，0：購入せず）
options(digits = 2) # 表示の有効桁数２桁に設定
d0 <- read.csv(file = 'https://stats.dip.jp/01_ds/data/car_data_jp.csv')
summary(d0) # 要約統計量## お客様番号 性別 年齢 年収 購入判断
## Min. : 1 Length:1000 Min. :18 Min. : 15000 Min. :0.0
## 1st Qu.: 251 Class :character 1st Qu.:32 1st Qu.: 46375 1st Qu.:0.0
## Median : 500 Mode :character Median :40 Median : 72000 Median :0.0
## Mean : 500 Mean :40 Mean : 72689 Mean :0.4
## 3rd Qu.: 750 3rd Qu.:48 3rd Qu.: 90000 3rd Qu.:1.0
## Max. :1000 Max. :63 Max. :152500 Max. :1.0head(d0) # 先頭の5レコード表示## お客様番号 性別 年齢 年収 購入判断
## 1 385 男性 35 20000 0
## 2 681 男性 40 43500 0
## 3 353 男性 49 74000 0
## 4 895 男性 40 107500 1
## 5 661 男性 25 79000 0
## 6 846 女性 47 33500 1tail(d0) # 末尾の5レコード表示## お客様番号 性別 年齢 年収 購入判断
## 995 951 女性 53 104000 1
## 996 863 男性 38 59000 0
## 997 800 女性 47 23500 0
## 998 407 女性 28 138500 1
## 999 299 女性 48 134000 1
## 1000 687 女性 44 73500 0(n <- nrow(d0)) # 全レコードサイズ## [1] 1000# 乱数シード（乱数の種）：
# 無作為標本抽出でプログラム実行のたびに同じ標本が抽出される。
# 数字を変えると異なる組み合わせで毎回同じ標本が抽出される。
set.seed(5)
# 訓練データと試験データに分割
# 層化無作為標本抽出
# prop：訓練データの分割割合
# strata：設定した変数の割合がそれぞれのデータセットで均等になる（層化分割）。
library(rsample)
d.trte <- initial_split(d0, prop = 4/5, strata = 購入判断)
d.trte## <Training/Testing/Total>
## <799/201/1000>d.tr <- training(d.trte) # 訓練データ
d.te <- testing (d.trte) # 試験データ2 ロジスティック回帰モデル
\[ \begin{align} \hat{p}&=\frac{1}{1+\exp\{-(\beta_0+\beta_1 x_{i}+\cdots+\beta_k x_{k})\}}\\ &ここで，\\ &\quad \hat{p}：確率の推定値\\ &\quad x_{k}： 説明変数\\ &\quad \beta_k： 偏回帰係数\\ \end{align} \]
3 ロジスティック回帰分析
3.1 フィッティング
fit.all <- glm(購入判断 ~ ., data = d.tr, family = 'binomial')
summary(fit.all)##
## Call:
## glm(formula = 購入判断 ~ ., family = "binomial", data = d.tr)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.24e+01 9.34e-01 -13.25 <2e-16 ***
## お客様番号 1.30e-04 3.56e-04 0.37 0.71
## 性別男性 2.49e-01 2.10e-01 1.19 0.23
## 年齢 2.22e-01 1.72e-02 12.93 <2e-16 ***
## 年収 3.43e-05 3.68e-06 9.32 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1076.60 on 798 degrees of freedom
## Residual deviance: 584.57 on 794 degrees of freedom
## AIC: 594.6
##
## Number of Fisher Scoring iterations: 6#（注意）
# 回帰係数の統計的有意性の判定には，
# 二値データを扱うロジスティック回帰分析ではｔ値ではなくてｚ値を使用する。お客様番号と性別の回帰係数は有意水準5％で有意でない。 これらは購入判断に寄与しない説明変数だが， モデルを予測用に使うときは，そのまま残しておいてもよい。
3.2 モデル改良
有意だった説明変数だけでフィッティング
fit <- glm(購入判断 ~ 年齢 + 年収, data = d.tr, family = 'binomial')
summary(fit)##
## Call:
## glm(formula = 購入判断 ~ 年齢 + 年収, family = "binomial",
## data = d.tr)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.20e+01 8.74e-01 -13.78 <2e-16 ***
## 年齢 2.20e-01 1.69e-02 12.99 <2e-16 ***
## 年収 3.40e-05 3.67e-06 9.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1076.60 on 798 degrees of freedom
## Residual deviance: 586.09 on 796 degrees of freedom
## AIC: 592.1
##
## Number of Fisher Scoring iterations: 63.3 予測
年収80,000US＄の45歳の人の新車購入確率を求める。
d.new <- data.frame(年齢 = 45, 年収 = 80000)
p.hat <- predict(fit, type = 'response', newdata = d.new) # 予測確率
sprintf('新車購入確率：%2.1f％', p.hat * 100)## [1] "新車購入確率：63.5％"3.4 混同行列（confusion matrix）による精度評価
混同行列を作成して精度を計算する。
\[ Accuracy = \frac{TP+TN}{TP+FP+TN+FN} \] ここで，\(Accuracy\)は精度
TP: True positive；真陽性数
TN: True negative；真陰性数
FP: False positive；偽陽性数
FN: False negative；偽陰性数
この指標は2値クラスの要素がいくつ正確に分類されたかを示す。 2値クラスの要素数が大きく異なる（アンバランスな）場合は，この指標を評価には用いては行けない。 なぜなら，大多数側に単純に分類することで高精度が得られるためである。
p.hat <- predict(fit, type = 'response', newdata = d.te)
threshold <- 0.5 # 閾値（しきいち），カットオフ値（cut-off）とも呼ばれる。
is.pred <- p.hat > threshold
is.ref <- d.te$購入判断 == 1
table(予測値 = is.pred, 真値 = is.ref)## 真値
## 予測値 FALSE TRUE
## FALSE 106 24
## TRUE 14 57is.ok <- is.pred == is.ref
n.ok <- sum(is.ok)
sprintf('新車購入予測精度：%2.1f％', n.ok / nrow(d.te) * 100)## [1] "新車購入予測精度：81.1％"# 混同（こんどう）行列 (confusion matrix) caretパッケージ利用による詳細分析
library(caret)
confusionMatrix(data = as.factor(is.pred),
reference = as.factor(is.ref))## Confusion Matrix and Statistics
##
## Reference
## Prediction FALSE TRUE
## FALSE 106 24
## TRUE 14 57
##
## Accuracy : 0.811
## 95% CI : (0.75, 0.863)
## No Information Rate : 0.597
## P-Value [Acc > NIR] : 7e-11
##
## Kappa : 0.599
##
## Mcnemar's Test P-Value : 0.144
##
## Sensitivity : 0.883
## Specificity : 0.704
## Pos Pred Value : 0.815
## Neg Pred Value : 0.803
## Prevalence : 0.597
## Detection Rate : 0.527
## Detection Prevalence : 0.647
## Balanced Accuracy : 0.794
##
## 'Positive' Class : FALSE
## 3.4.1 ROC曲線
真陽性率（TPR）と偽陽性率（FPR）のトレードオフを表す曲線。 分類器としての性能はAUC（Area under the curve），つまり曲線下側面積の大きさで表される。 精度とは異なり，真陽性と偽陽性の確率値で表されるため，クラスのどちらかが大多数のときに生ずる問題がない。
2値判別の閾値（しきいち）は，確率50％のところに置くことが基本だが， ROC曲線を描くことにより，最適な閾値を求めることができる。 閾値が50%以外のときは，予測値には「確率」ではなく「判別スコア」という用語を使うこと。
3.4.1.1 用語
- ROC (Receiver Operating Characteristic)：受信者動作特性 ※
元々の意味は忘れて良い
- AUC (Area Under Curve) ：ROC曲線の下側の面積（灰色部分）
- Sensitivity (true positive rate) ：感度（真陽性率）
- Specificity (true negative rate) ：特異度（真陰性率）= 1-偽陽性率； 軸の向きに注意
- Threshold ：閾値
library(pROC)
roc1 <- roc(response = d.te$購入判断, predict = p.hat,
of = 'thresholds', thresholds = 'best', print.thres = 'best',
percent = F, plot = T, print.auc = T, grid = T, ci = T, auc.polygon=T)
coords(roc1, 'best') # 最適閾値 (optimal threshold)## threshold specificity sensitivity
## 1 0.31 0.77 0.9c <- coords(roc1) # 閾値 (threshold)ごとの値 (sensitivity, specificity)
library(DT)
datatable(round(c, 3))4 演習課題
次のタイタニック号の乗船者データを用いてロジスティック回帰分析を行え。
また，次の問いに答えよ。
- 運賃は生死を分ける要素と考えて良いか？
- 2等客室にいる45才の女性の生存確率を推定せよ。
- 生死を判定するのに最適な，判別スコアの閾値を求めよ。
4.1 データ
d0 <- read.csv(file = 'https://stats.dip.jp/01_ds/data/titanic_data_jp.csv')
datatable(d0)