情報科学3
モデルの性能評価
矢内 勇生
2021-11-12
準備
パッケージの読み込み
必要なパッケージを読み込む。
pacman::p_load(tidyverse,
broom,
irr,
gmodels,
caret,
rsample,
ISLR)次に、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),
LDP = ifelse (party == "LDP", 1, 0))分類の性能評価
ロジスティック回帰による分類
ロジスティック回帰 (2) で実施した内容の繰り返しだが、 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である。
分類性能の評価
混同行列を作る。
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 |
## -------------|-----------|-----------|-----------|
##
## この予測の正解率\(p_o\)は、
(p_o <- (641 + 168) / 989)## [1] 0.817998である。
ランダムに分類を行ったときに分類が正しくなる確率\(p_e\)は、
(p_e <- (773/989) * (689/989) + (216/989) * (300 / 989))## [1] 0.6107598である。
コーエンの\(\kappa\)は \[ \kappa = \frac{p_o - p_e}{1 - p_e} \] なので、
(p_o - p_e) / (1 - p_e)## [1] 0.5324172である。
この値は、irrパッケージのkappa2()で求めることができる。1列目が観測された応答変数、2列目が予測値であるデータフレームを渡す必要がある。
df <- tibble(observed = HR2005$wlsmd,
predicted = y_hat,)
kappa2(df)## Cohen's Kappa for 2 Raters (Weights: unweighted)
##
## Subjects = 989
## Raters = 2
## Kappa = 0.532
##
## z = 17.2
## p-value = 0混同行列に表示されている値を利用して計算した結果と同じ\(\kappa\)が得られた。 \(\kappa\)の値のみ取り出したいときは、次のようにする。
kappa2(df)$value## [1] 0.5324172正解率は0.82でそれなりに良さそうに見えるが、ランダムに予測してもそこそこ当たることを考慮に入れると、このモデルの予測性能はあまり高くないということがわかる。
分類性能を測るその他の指標やROC曲線、AUCについては過去の資料を参照されたい。
10分割交差検証
\(k\)分割交差検証を\(k=10\)で行ってみよう。
caretパッケージの createFolds() で10個のフォールドを作る。
HR2005_folds <- createFolds(HR2005$wlsmd, k = 10)glimpse(HR2005_folds)## List of 10
## $ Fold01: int [1:99] 2 13 18 28 61 69 86 87 110 121 ...
## $ Fold02: int [1:99] 20 65 78 98 114 116 117 122 125 129 ...
## $ Fold03: int [1:99] 6 7 32 35 44 49 70 73 75 80 ...
## $ Fold04: int [1:99] 3 5 12 24 47 63 71 72 79 81 ...
## $ Fold05: int [1:99] 21 36 42 66 77 90 92 93 103 107 ...
## $ Fold06: int [1:99] 1 16 22 23 38 40 50 59 67 76 ...
## $ Fold07: int [1:99] 8 30 34 41 48 52 53 55 56 57 ...
## $ Fold08: int [1:98] 11 15 25 29 39 46 51 60 62 95 ...
## $ Fold09: int [1:99] 9 14 19 31 37 43 45 64 74 102 ...
## $ Fold10: int [1:99] 4 10 17 26 27 33 54 89 99 100 ...10個のフォールドができている。このリストのそれぞれの要素に含まれる数値は、データセットの行番号で、どの行がどのfold に割り当てられたかを示している。
作ったそれぞれのフォールドをテストデータに、残りを学習データとして用いてコーエンの\(\kappa\)を計算する。 まず、1つのフォールドについて\(\kappa\)を計算する関数を作る。
logistic_cv <- function(fold) {
df_train <- HR2005[-fold, ]
df_test <- HR2005[fold, ]
fit <- glm(wlsmd ~ previous,
data = df_train,
family = binomial(link = "logit"))
pred <- predict(fit, newdata = df_test, type = "response")
df_k <- tibble(observed = df_test$wlsmd,
predicted = ifelse(pred > 0.5, 1, 0))
kappa2(df_k)$value
}10個のフォールドのそれぞれでこの関数を使って\(\kappa\)の値を求め、10個の\(\kappa\)の平均値を求める。
HR2005_folds %>%
map(.f = logistic_cv) %>%
unlist() %>%
mean()## [1] 0.5303089やはり、あまり性能は良くないようだ。
回帰の性能評価
2005年の衆院選データについて、候補者の小選挙区における得票率 voteshare を応答変数とする回帰分析を行う。 説明変数として、以下の3つを使う。
- 選挙費用（100万円）：
expm - 過去の当選回数：
previous - 自民党ダミー:
LDP
自民党ダミーを含まないモデル1と、自民党ダミーを含むモデル2を比べる。
f1 <- formula(voteshare ~ expm_c + previous)
f2 <- formula(voteshare ~ expm_c + previous + LDP)性能評価のために、データを訓練用とテスト用に分ける。
HR05 <- HR2005 %>%
filter(!is.na(expm)) %>%
mutate(expm_c = expm - mean(expm))
HR05_split <- initial_split(HR05, prop = 2/3)
HR05_train <- training(HR05_split)
HR05_test <- testing(HR05_split)訓練用データで学習を行う。
ols1 <- lm(f1,
data = HR05_train)
ols2 <- lm(f2,
data = HR05_train)自由度調整済み決定係数は、
cbind(summary(ols1)$adj.r.squared,
summary(ols2)$adj.r.squared)## [,1] [,2]
## [1,] 0.6245738 0.6878217であり、モデル2のほうが当てはまりがよい。
訓練データでRMSEを計算する。
rmse1 <- (HR05_train$voteshare - ols1$fitted.values)^2 %>%
sum() %>%
sqrt()
rmse2 <- (HR05_train$voteshare - ols2$fitted.values)^2 %>%
sum() %>%
sqrt()
cbind(rmse1, rmse2)## rmse1 rmse2
## [1,] 298.3754 271.875（当てはまり良いので当然だが）RMSEもモデル2のほうが小さい。
テストデータでRMSEを計算する。
yhat1 <- predict(ols1, newdata = HR05_test)
yhat2 <- predict(ols2, newdata = HR05_test)
rmse1t <- (HR05_test$voteshare - yhat1)^2 %>%
sum() %>%
sqrt()
rmse2t <- (HR05_test$voteshare - yhat2)^2 %>%
sum() %>%
sqrt()
cbind(rmse1t, rmse2t)## rmse1t rmse2t
## [1,] 195.5505 180.4921やはりモデル2のほうが性能が良い。
10分割交差検証
\(k\)分割交差検証を\(k=10\)で行ってみよう。
フォールドを作り直す。
HR05_folds <- createFolds(HR05$wlsmd, k = 10)作ったそれぞれのフォールドをテストデータに、残りを学習データとして用いてRMSEを計算する。 まず、1つのフォールドについてRMSEを計算する関数を作る。
ols_cv <- function(fold, formula) {
df_train <- HR05[-fold, ]
df_test <- HR05[fold, ]
fit <- lm(formula,
data = df_train)
yhat <- predict(fit, newdata = df_test)
(df_test$voteshare - yhat)^2 %>%
sum() %>%
sqrt()
}10個のフォールドのそれぞれでこの関数を使って\(\kappa\)の値を求め、10個の\(\kappa\)の平均値を求める。
m1 <- HR05_folds %>%
map(.f = ols_cv,
formula = f1) %>%
unlist() %>%
mean()
m2 <- HR05_folds %>%
map(.f = ols_cv,
formula = f2) %>%
unlist() %>%
mean()
cbind(m1, m2)## m1 m2
## [1,] 113.055 103.3513やはり、モデル2のほうが予測性能が良さそうだ。