情報科学3
ロジスティック回帰 (2)
矢内 勇生
2021-10-28（改訂: 2021-11-11）
準備
必要なパッケージを読み込む。
pacman::p_load(tidyverse,
broom,
gmodels,
caret,
ROCR)次に、ggplot2のテーマとフォントの設定を行う。自分好みの設定がある場合は自由に変えてよい。
if (.Platform$OS.type == "windows") {
if (require(fontregisterer)) {
my_font <- "Yu Gothic"
} else {
my_font <- "Japan1"
}
} else if (capabilities("aqua")) {
my_font <- "HiraginoSans-W3"
} else {
my_font <- "IPAexGothic"
}
theme_set(theme_gray(base_size = 9,
base_family = my_font))前回同様、説明のために『Rによる計量政治学』（浅野正彦, 矢内勇生. 2018）で使用されているデータ（hr-data.csv） を使う。
HR <- read_csv("data/hr-data.csv")小選挙区での当落を表すダミー変数（これを二値の応答変数として扱う）を作る。
HR <- HR %>%
mutate(wlsmd = ifelse(wl == "当選", 1, 0))ロジスティック回帰
例：過去の当選回数で当選確率を説明する
2005年のデータを取り出し、候補者の小選挙区における当落 (wlsmd) を縦軸に、それまでの当選回数 (previous) を横軸にした散布図を表示してみよう。
HR2005 <- HR %>%
filter(year == 2005)
p <- ggplot(HR2005, aes(previous, wlsmd)) +
geom_jitter(size = 1,
alpha = 1/3,
height = 0.02,
width = 0.2) +
labs(x = "過去の当選回数",
y = "小選挙区での当選")
plot(p)
このデータを使い、当落を応答変数、過去の当選回数を説明変数として、ロジスティック回帰分析を行おう。
fit <- glm(wlsmd ~ previous,
data = HR2005,
family = binomial(link = "logit"))
tidy(fit, conf.int = TRUE) %>%
select(term, estimate, std.error, conf.low, conf.high)## # A tibble: 2 × 5
## term estimate std.error conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -2.34 0.133 -2.61 -2.09
## 2 previous 0.713 0.0475 0.623 0.809この結果を図示しよう。
p_fit <- p +
geom_smooth(method = "glm",
method.args = list(family = binomial(link = "logit")))
plot(p_fit)
次に、成功確率（当選確率）の予測値を求める。
pred <- predict(fit, type = "response")予測確率が0.5以上の候補者を「当選 (1)」に、0.5未満の候補者を「落選 (0)」に分類する。
y_hat <- ifelse(pred >= 0.5, 1, 0)分類結果は、
table(y_hat)## y_hat
## 0 1
## 773 216である。
分類性能の評価
この分類の性能を確認するために、混同行列 (confusion matrix) を作る。 gmodels::CrossTable() で作れる。
CrossTable(x = HR2005$wlsmd,
y = y_hat,
dnn = c("実際の当落", "当落の予測"),
prop.r = FALSE, # 行%は表示しない
prop.c = FALSE, # 列%は表示しない
prop.t = TRUE, # 合計%は表示する
prop.chisq = FALSE) # chi^2値は表示しない##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 989
##
##
## | 当落の予測
## 実際の当落 | 0 | 1 | Row Total |
## -------------|-----------|-----------|-----------|
## 0 | 641 | 48 | 689 |
## | 0.648 | 0.049 | |
## -------------|-----------|-----------|-----------|
## 1 | 132 | 168 | 300 |
## | 0.133 | 0.170 | |
## -------------|-----------|-----------|-----------|
## Column Total | 773 | 216 | 989 |
## -------------|-----------|-----------|-----------|
##
## この混同行列から、以下のことがわかる。
- 真陰性の個体数は641
- 偽陽性の個体数は48
- 偽陰性の個体数は132
- 真陽性の個体数は168
また、
- 真陰性率は0.648
- 偽陽性率は0.049
- 偽陰性率は0.133
- 真陽性率は0.170
これらの値は、自分で計算することもできる。
(TN <- sum(HR2005$wlsmd == 0 & y_hat == 0)) # 真陰性数## [1] 641(FP <- sum(HR2005$wlsmd == 0 & y_hat == 1)) # 偽陽性数## [1] 48(FN <- sum(HR2005$wlsmd == 1 & y_hat == 0)) # 偽陰性数## [1] 132(TP <- sum(HR2005$wlsmd == 1 & y_hat == 1)) # 真陽性数## [1] 168(TNR <- mean(HR2005$wlsmd == 0 & y_hat == 0)) # 真陰性率## [1] 0.6481294(FPR <- mean(HR2005$wlsmd == 0 & y_hat == 1)) # 偽陽性率## [1] 0.04853387(FNR <- mean(HR2005$wlsmd == 1 & y_hat == 0)) # 偽陰性率## [1] 0.1334681(TPR <- mean(HR2005$wlsmd == 1 & y_hat == 1)) # 真陽性率## [1] 0.1698686これらの値を元に、いくつかの評価指標を確認してみよう。 まず、正解率（正確度）は、
TPR + TNR## [1] 0.817998#(TP + TN) / (TP + TN + FP + FN) # これと同じである。正解率はそれなり高そうである。（当選者よりも落選者のほうが多いので）全員を落選に分類すれば、
mean(HR2005$wlsmd == 0)## [1] 0.6966633すなわち約70％の正解率になる。ロジスティック回帰を使うことによって正解率が約12ポイント上昇している。
次に、適合率は、
TP / (TP + FP)## [1] 0.7777778であり、「当選」と分類されたが実際には落選する候補者が22%ほどいることがわかる。
真陽性率（感度）は、
TP / (TP + FN)## [1] 0.56であり、実際に当選した候補者のうち、44%が誤って落選に分類されていることがわかる。つまり、感度は高くない。
偽陽性率は、
FP / (TN + FP)## [1] 0.06966618であり、落選した候補者のうち誤って当選に分類されるのは7%未満である。特異度は、
1 - FP / (TN + FP)## [1] 0.9303338であり、高いと言えるだろう。
感度 (sensitivity) と特異度 (specificity) は、caret パッケージのconfusionMatrix() で混同行列を作ると自動的に計算される。 ただし、この関数を使うためにはラベルを factor 型にする必要がある。さらに、positive 引数で陽性のラベルを指定する。実際の値が列に、予測が行に表示されるので注意されたい。
confusionMatrix(reference = as.factor(HR2005$wlsmd),
data = as.factor(y_hat),
positive = "1")## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 641 132
## 1 48 168
##
## Accuracy : 0.818
## 95% CI : (0.7925, 0.8416)
## No Information Rate : 0.6967
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.5324
##
## Mcnemar's Test P-Value : 6.153e-10
##
## Sensitivity : 0.5600
## Specificity : 0.9303
## Pos Pred Value : 0.7778
## Neg Pred Value : 0.8292
## Prevalence : 0.3033
## Detection Rate : 0.1699
## Detection Prevalence : 0.2184
## Balanced Accuracy : 0.7452
##
## 'Positive' Class : 1
## F値は
(2 * TP) / (2 * TP + FP + FN)## [1] 0.6511628である。
最後に、ROC曲線を描いてみよう。ROCR パッケージを利用する。
#pred <- predict(fit, type = "response") # 上で実行済み
pr <- prediction(pred, labels = HR2005$wlsmd == 1)
roc <- performance(pr, measure = "tpr", x.measure = "fpr")
df_roc <- tibble(fpr = roc@x.values[[1]],
tpr = roc@y.values[[1]])
ROC <- ggplot(df_roc, aes(x = fpr, y = tpr)) +
geom_line(color = "dodgerblue") +
geom_abline(color = "gray",
linetype = "dashed") +
labs(x = "偽陽性率（1 - 特異度）",
y = "真陽性率（感度）") +
coord_fixed()
plot(ROC)
モデルがない場合に比べれば、それなりに役立つ予測ができそうだということがわかる。
AUCは、
auc <- performance(pr, measure = "auc")
auc@y.values[[1]]## [1] 0.8965094である。慣例に従えば、「良い」モデルだと言えそうだ。
ロジスティック回帰による分類
浅野・矢内 (2018) の衆院選データには、1996 年から2017年までの8回の衆院選データが含まれている。 1996年から2009年までの5回分のデータを使って学習を行い、2012年の小選挙区で当落を予測（分類）してみよう（2014年と2017年の選挙費用データが欠測なので、これらの年は使わない）。
データの準備
以下のモデルで選挙費用を利用するので、選挙費用が欠測になっている観測値をあらかじめ取り除く。また、自民党 (LDP) に所属 することを示すダミー変を数を作る。
HR_sub <- HR %>%
filter(!is.na(expm)) %>%
mutate(LDP = ifelse(party == "LDP", 1, 0))訓練データ（2009年まで）と検証データ（2012年）を用意する。
HR_train <- HR_sub %>%
filter(year < 2012)
HR_test <- HR_sub %>%
filter(year == 2012)モデル
以下の4つのロジスティック回帰モデルを考える。 \[\begin{align} \mathrm{logit}(\theta)_n & = \beta_0 + \beta_1 \mathrm{選挙費用}_n \\ \mathrm{logit}(\theta)_n & = \gamma_0 + \gamma_1 \mathrm{過去の当選回数}_n \\ \mathrm{logit}(\theta)_n & = \tau_0 + \tau_1 \mathrm{選挙費用}_n + \tau_2 \mathrm{過去の当選回数}_n\\ \mathrm{logit}(\theta)_n & = \lambda_0 + \lambda_1 \mathrm{選挙費用}_n + \lambda_2 \mathrm{過去の当選回数}_n + \lambda_3 \mathrm{自民党ダミー}_n\\ \end{align}\]
これらのモデルを使い、訓練データでロジスティック回帰の係数を推定する。
models <- alist(
model1 = wlsmd ~ expm,
model2 = wlsmd ~ previous,
model3 = wlsmd ~ expm + previous,
model4 = wlsmd ~ expm + previous + LDP) %>%
map(formula)
fits <- models %>%
map(.f = glm,
data = HR_train,
family = binomial(link = "logit"))推定された係数を利用して、検証データ（2012年）の小選挙区での当落の予測確率を計算する。
preds <- fits %>%
map(.f = predict,
newdata = HR_test,
type = "response")予測確率を利用して分類を行う。
y_hat_1 <- ifelse(preds$model1 >= 0.5, 1, 0)
y_hat_2 <- ifelse(preds$model2 >= 0.5, 1, 0)
y_hat_3 <- ifelse(preds$model3 >= 0.5, 1, 0)
y_hat_4 <- ifelse(preds$model4 >= 0.5, 1, 0)モデル1の混同行列を表示する。
CrossTable(x = HR_test$wlsmd,
y = y_hat_1,
dnn = c("実際の当落", "当落の予測"),
prop.r = FALSE,
prop.c = FALSE,
prop.t = TRUE,
prop.chisq = FALSE)##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 1280
##
##
## | 当落の予測
## 実際の当落 | 0 | 1 | Row Total |
## -------------|-----------|-----------|-----------|
## 0 | 967 | 16 | 983 |
## | 0.755 | 0.013 | |
## -------------|-----------|-----------|-----------|
## 1 | 262 | 35 | 297 |
## | 0.205 | 0.027 | |
## -------------|-----------|-----------|-----------|
## Column Total | 1229 | 51 | 1280 |
## -------------|-----------|-----------|-----------|
##
## モデル2の混同行列。
CrossTable(x = HR_test$wlsmd,
y = y_hat_2,
dnn = c("実際の当落", "当落の予測"),
prop.r = FALSE,
prop.c = FALSE,
prop.t = TRUE,
prop.chisq = FALSE)##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 1280
##
##
## | 当落の予測
## 実際の当落 | 0 | 1 | Row Total |
## -------------|-----------|-----------|-----------|
## 0 | 915 | 68 | 983 |
## | 0.715 | 0.053 | |
## -------------|-----------|-----------|-----------|
## 1 | 153 | 144 | 297 |
## | 0.120 | 0.113 | |
## -------------|-----------|-----------|-----------|
## Column Total | 1068 | 212 | 1280 |
## -------------|-----------|-----------|-----------|
##
## モデル3の混同行列。
CrossTable(x = HR_test$wlsmd,
y = y_hat_3,
dnn = c("実際の当落", "当落の予測"),
prop.r = FALSE,
prop.c = FALSE,
prop.t = TRUE,
prop.chisq = FALSE)##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 1280
##
##
## | 当落の予測
## 実際の当落 | 0 | 1 | Row Total |
## -------------|-----------|-----------|-----------|
## 0 | 942 | 41 | 983 |
## | 0.736 | 0.032 | |
## -------------|-----------|-----------|-----------|
## 1 | 188 | 109 | 297 |
## | 0.147 | 0.085 | |
## -------------|-----------|-----------|-----------|
## Column Total | 1130 | 150 | 1280 |
## -------------|-----------|-----------|-----------|
##
## モデル4の混同行列。
CrossTable(x = HR_test$wlsmd,
y = y_hat_4,
dnn = c("実際の当落", "当落の予測"),
prop.r = FALSE,
prop.c = FALSE,
prop.t = TRUE,
prop.chisq = FALSE)##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 1280
##
##
## | 当落の予測
## 実際の当落 | 0 | 1 | Row Total |
## -------------|-----------|-----------|-----------|
## 0 | 946 | 37 | 983 |
## | 0.739 | 0.029 | |
## -------------|-----------|-----------|-----------|
## 1 | 168 | 129 | 297 |
## | 0.131 | 0.101 | |
## -------------|-----------|-----------|-----------|
## Column Total | 1114 | 166 | 1280 |
## -------------|-----------|-----------|-----------|
##
## それぞれのモデルについて、ROC曲線を表示する。
prs <- preds %>%
map(.f = prediction,
labels = HR_test$wlsmd == 1)
rocs <- prs %>%
map(.f = performance,
measure = "tpr",
x.measure = "fpr")
fprs <- list(
fpr_1 = rocs$model1@x.values[[1]],
fpr_2 = rocs$model2@x.values[[1]],
fpr_3 = rocs$model3@x.values[[1]],
fpr_4 = rocs$model4@x.values[[1]])
tprs <- list(
tpr_1 = rocs$model1@y.values[[1]],
tpr_2 = rocs$model2@y.values[[1]],
tpr_3 = rocs$model3@y.values[[1]],
tpr_4 = rocs$model4@y.values[[1]])
roc_lengths <- sapply(fprs, length)
df_rocs <- tibble(
fpr = unlist(fprs),
tpr = unlist(tprs),
model = factor(rep(1:4, times = roc_lengths))
)
ROCs <- ggplot(df_rocs, aes(x = fpr, y = tpr, color = model)) +
geom_line() +
geom_abline(color = "gray",
linetype = "dashed") +
labs(x = "偽陽性率（1 - 特異度）",
y = "真陽性率（感度）") +
scale_color_brewer(palette = "Set2",
name = "モデル") +
coord_fixed()
plot(ROCs)
モデル4の分類性能が最も高そうだということがわかる。
AUCを計算してみよう。
aucs <- prs %>%
map(.f = performance,
measure = "auc")
cbind(aucs$model1@y.values[[1]],
aucs$model2@y.values[[1]],
aucs$model3@y.values[[1]],
aucs$model4@y.values[[1]])## [,1] [,2] [,3] [,4]
## [1,] 0.8174505 0.8582005 0.8859449 0.9201476モデル4のAUCが0.9を超えており、「とても良い」分類ができていると考えられる。 ただし、このモデルが2012年以外の選挙をうまく予測できるとは限らない。
実習課題
- それぞれのモデルについて、以下の評価指標を計算しなさい。
- 正解率（正確度）
- 適合率
- 真陽性率（感度）
- 特異度 (1 - 偽陽性率)
- F値（F1スコア）