必要なパッケージを読み込む。
::p_load(tidyverse,
pacman
broom,
gmodels,
caret, ROCR)
次に、ggplot2のテーマとフォントの設定を行う。自分好みの設定がある場合は自由に変えてよい。
if (.Platform$OS.type == "windows") {
if (require(fontregisterer)) {
<- "Yu Gothic"
my_font else {
} <- "Japan1"
my_font
}else if (capabilities("aqua")) {
} <- "HiraginoSans-W3"
my_font else {
} <- "IPAexGothic"
my_font
}
theme_set(theme_gray(base_size = 9,
base_family = my_font))
前回同様、説明のために『Rによる計量政治学』(浅野正彦, 矢内勇生. 2018)で使用されているデータ(hr-data.csv) を使う。
<- read_csv("data/hr-data.csv") HR
小選挙区での当落を表すダミー変数(これを二値の応答変数として扱う)を作る。
<- HR %>%
HR mutate(wlsmd = ifelse(wl == "当選", 1, 0))
2005年のデータを取り出し、候補者の小選挙区における当落 (wlsmd) を縦軸に、それまでの当選回数 (previous) を横軸にした散布図を表示してみよう。
<- HR %>%
HR2005 filter(year == 2005)
<- ggplot(HR2005, aes(previous, wlsmd)) +
p geom_jitter(size = 1,
alpha = 1/3,
height = 0.02,
width = 0.2) +
labs(x = "過去の当選回数",
y = "小選挙区での当選")
plot(p)
このデータを使い、当落を応答変数、過去の当選回数を説明変数として、ロジスティック回帰分析を行おう。
<- glm(wlsmd ~ previous,
fit 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 +
p_fit geom_smooth(method = "glm",
method.args = list(family = binomial(link = "logit")))
plot(p_fit)
次に、成功確率(当選確率)の予測値を求める。
<- predict(fit, type = "response") pred
予測確率が0.5以上の候補者を「当選 (1)」に、0.5未満の候補者を「落選 (0)」に分類する。
<- ifelse(pred >= 0.5, 1, 0) y_hat
分類結果は、
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 |
## -------------|-----------|-----------|-----------|
##
##
この混同行列から、以下のことがわかる。
また、
これらの値は、自分で計算することもできる。
<- sum(HR2005$wlsmd == 0 & y_hat == 0)) # 真陰性数 (TN
## [1] 641
<- sum(HR2005$wlsmd == 0 & y_hat == 1)) # 偽陽性数 (FP
## [1] 48
<- sum(HR2005$wlsmd == 1 & y_hat == 0)) # 偽陰性数 (FN
## [1] 132
<- sum(HR2005$wlsmd == 1 & y_hat == 1)) # 真陽性数 (TP
## [1] 168
<- mean(HR2005$wlsmd == 0 & y_hat == 0)) # 真陰性率 (TNR
## [1] 0.6481294
<- mean(HR2005$wlsmd == 0 & y_hat == 1)) # 偽陽性率 (FPR
## [1] 0.04853387
<- mean(HR2005$wlsmd == 1 & y_hat == 0)) # 偽陰性率 (FNR
## [1] 0.1334681
<- mean(HR2005$wlsmd == 1 & y_hat == 1)) # 真陽性率 (TPR
## [1] 0.1698686
これらの値を元に、いくつかの評価指標を確認してみよう。 まず、正解率(正確度)は、
+ TNR TPR
## [1] 0.817998
#(TP + TN) / (TP + TN + FP + FN) # これと同じ
である。正解率はそれなり高そうである。(当選者よりも落選者のほうが多いので)全員を落選に分類すれば、
mean(HR2005$wlsmd == 0)
## [1] 0.6966633
すなわち約70%の正解率になる。ロジスティック回帰を使うことによって正解率が約12ポイント上昇している。
次に、適合率は、
/ (TP + FP) TP
## [1] 0.7777778
であり、「当選」と分類されたが実際には落選する候補者が22%ほどいることがわかる。
真陽性率(感度)は、
/ (TP + FN) TP
## [1] 0.56
であり、実際に当選した候補者のうち、44%が誤って落選に分類されていることがわかる。つまり、感度は高くない。
偽陽性率は、
/ (TN + FP) 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") # 上で実行済み
<- prediction(pred, labels = HR2005$wlsmd == 1)
pr <- performance(pr, measure = "tpr", x.measure = "fpr")
roc
<- tibble(fpr = roc@x.values[[1]],
df_roc tpr = roc@y.values[[1]])
<- ggplot(df_roc, aes(x = fpr, y = tpr)) +
ROC geom_line(color = "dodgerblue") +
geom_abline(color = "gray",
linetype = "dashed") +
labs(x = "偽陽性率(1 - 特異度)",
y = "真陽性率(感度)") +
coord_fixed()
plot(ROC)
モデルがない場合に比べれば、それなりに役立つ予測ができそうだということがわかる。
AUCは、
<- performance(pr, measure = "auc")
auc @y.values[[1]] auc
## [1] 0.8965094
である。慣例に従えば、「良い」モデルだと言えそうだ。
浅野・矢内 (2018) の衆院選データには、1996 年から2017年までの8回の衆院選データが含まれている。 1996年から2009年までの5回分のデータを使って学習を行い、2012年の小選挙区で当落を予測(分類)してみよう(2014年と2017年の選挙費用データが欠測なので、これらの年は使わない)。
以下のモデルで選挙費用を利用するので、選挙費用が欠測になっている観測値をあらかじめ取り除く。また、自民党 (LDP) に所属 することを示すダミー変を数を作る。
<- HR %>%
HR_sub filter(!is.na(expm)) %>%
mutate(LDP = ifelse(party == "LDP", 1, 0))
訓練データ(2009年まで)と検証データ(2012年)を用意する。
<- HR_sub %>%
HR_train filter(year < 2012)
<- HR_sub %>%
HR_test 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}\]
これらのモデルを使い、訓練データでロジスティック回帰の係数を推定する。
<- alist(
models model1 = wlsmd ~ expm,
model2 = wlsmd ~ previous,
model3 = wlsmd ~ expm + previous,
model4 = wlsmd ~ expm + previous + LDP) %>%
map(formula)
<- models %>%
fits map(.f = glm,
data = HR_train,
family = binomial(link = "logit"))
推定された係数を利用して、検証データ(2012年)の小選挙区での当落の予測確率を計算する。
<- fits %>%
preds map(.f = predict,
newdata = HR_test,
type = "response")
予測確率を利用して分類を行う。
<- ifelse(preds$model1 >= 0.5, 1, 0)
y_hat_1 <- ifelse(preds$model2 >= 0.5, 1, 0)
y_hat_2 <- ifelse(preds$model3 >= 0.5, 1, 0)
y_hat_3 <- ifelse(preds$model4 >= 0.5, 1, 0) y_hat_4
モデル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曲線を表示する。
<- preds %>%
prs map(.f = prediction,
labels = HR_test$wlsmd == 1)
<- prs %>%
rocs map(.f = performance,
measure = "tpr",
x.measure = "fpr")
<- list(
fprs 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]])
<- list(
tprs 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]])
<- sapply(fprs, length)
roc_lengths
<- tibble(
df_rocs fpr = unlist(fprs),
tpr = unlist(tprs),
model = factor(rep(1:4, times = roc_lengths))
)
<- ggplot(df_rocs, aes(x = fpr, y = tpr, color = model)) +
ROCs geom_line() +
geom_abline(color = "gray",
linetype = "dashed") +
labs(x = "偽陽性率(1 - 特異度)",
y = "真陽性率(感度)") +
scale_color_brewer(palette = "Set2",
name = "モデル") +
coord_fixed()
plot(ROCs)
モデル4の分類性能が最も高そうだということがわかる。
AUCを計算してみよう。
<- prs %>%
aucs map(.f = performance,
measure = "auc")
cbind(aucs$model1@y.values[[1]],
$model2@y.values[[1]],
aucs$model3@y.values[[1]],
aucs$model4@y.values[[1]]) aucs
## [,1] [,2] [,3] [,4]
## [1,] 0.8174505 0.8582005 0.8859449 0.9201476
モデル4のAUCが0.9を超えており、「とても良い」分類ができていると考えられる。 ただし、このモデルが2012年以外の選挙をうまく予測できるとは限らない。
実習課題