必要なパッケージを読み込む。
::p_load(tidyverse,
pacman
broom,
irr,
gmodels,
caret,
rsample, ISLR)
次に、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),
LDP = ifelse (party == "LDP", 1, 0))
ロジスティック回帰 (2) で実施した内容の繰り返しだが、 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
である。
混同行列を作る。
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\)は、
<- (641 + 168) / 989) (p_o
## [1] 0.817998
である。
ランダムに分類を行ったときに分類が正しくなる確率\(p_e\)は、
<- (773/989) * (689/989) + (216/989) * (300 / 989)) (p_e
## [1] 0.6107598
である。
コーエンの\(\kappa\)は \[ \kappa = \frac{p_o - p_e}{1 - p_e} \] なので、
- p_e) / (1 - p_e) (p_o
## [1] 0.5324172
である。
この値は、irrパッケージのkappa2()
で求めることができる。1列目が観測された応答変数、2列目が予測値であるデータフレームを渡す必要がある。
<- tibble(observed = HR2005$wlsmd,
df 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については過去の資料を参照されたい。
\(k\)分割交差検証を\(k=10\)で行ってみよう。
caretパッケージの createFolds()
で10個のフォールドを作る。
<- createFolds(HR2005$wlsmd, k = 10) HR2005_folds
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\)を計算する関数を作る。
<- function(fold) {
logistic_cv <- HR2005[-fold, ]
df_train <- HR2005[fold, ]
df_test <- glm(wlsmd ~ previous,
fit data = df_train,
family = binomial(link = "logit"))
<- predict(fit, newdata = df_test, type = "response")
pred <- tibble(observed = df_test$wlsmd,
df_k 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つを使う。
expm
previous
LDP
自民党ダミーを含まないモデル1と、自民党ダミーを含むモデル2を比べる。
<- formula(voteshare ~ expm_c + previous)
f1 <- formula(voteshare ~ expm_c + previous + LDP) f2
性能評価のために、データを訓練用とテスト用に分ける。
<- HR2005 %>%
HR05 filter(!is.na(expm)) %>%
mutate(expm_c = expm - mean(expm))
<- initial_split(HR05, prop = 2/3)
HR05_split <- training(HR05_split)
HR05_train <- testing(HR05_split) HR05_test
訓練用データで学習を行う。
<- lm(f1,
ols1 data = HR05_train)
<- lm(f2,
ols2 data = HR05_train)
自由度調整済み決定係数は、
cbind(summary(ols1)$adj.r.squared,
summary(ols2)$adj.r.squared)
## [,1] [,2]
## [1,] 0.6245738 0.6878217
であり、モデル2のほうが当てはまりがよい。
訓練データでRMSEを計算する。
<- (HR05_train$voteshare - ols1$fitted.values)^2 %>%
rmse1 sum() %>%
sqrt()
<- (HR05_train$voteshare - ols2$fitted.values)^2 %>%
rmse2 sum() %>%
sqrt()
cbind(rmse1, rmse2)
## rmse1 rmse2
## [1,] 298.3754 271.875
(当てはまり良いので当然だが)RMSEもモデル2のほうが小さい。
テストデータでRMSEを計算する。
<- predict(ols1, newdata = HR05_test)
yhat1 <- predict(ols2, newdata = HR05_test)
yhat2 <- (HR05_test$voteshare - yhat1)^2 %>%
rmse1t sum() %>%
sqrt()
<- (HR05_test$voteshare - yhat2)^2 %>%
rmse2t sum() %>%
sqrt()
cbind(rmse1t, rmse2t)
## rmse1t rmse2t
## [1,] 195.5505 180.4921
やはりモデル2のほうが性能が良い。
\(k\)分割交差検証を\(k=10\)で行ってみよう。
フォールドを作り直す。
<- createFolds(HR05$wlsmd, k = 10) HR05_folds
作ったそれぞれのフォールドをテストデータに、残りを学習データとして用いてRMSEを計算する。 まず、1つのフォールドについてRMSEを計算する関数を作る。
<- function(fold, formula) {
ols_cv <- HR05[-fold, ]
df_train <- HR05[fold, ]
df_test <- lm(formula,
fit data = df_train)
<- predict(fit, newdata = df_test)
yhat $voteshare - yhat)^2 %>%
(df_testsum() %>%
sqrt()
}
10個のフォールドのそれぞれでこの関数を使って\(\kappa\)の値を求め、10個の\(\kappa\)の平均値を求める。
<- HR05_folds %>%
m1 map(.f = ols_cv,
formula = f1) %>%
unlist() %>%
mean()
<- HR05_folds %>%
m2 map(.f = ols_cv,
formula = f2) %>%
unlist() %>%
mean()
cbind(m1, m2)
## m1 m2
## [1,] 113.055 103.3513
やはり、モデル2のほうが予測性能が良さそうだ。