Rを使った分析の準備

  • ここで使うRのパッケージは次のとおり。
library(corrplot)
library(jtools)
library(margins)
library(ROCR)
library(patchwork)
library(prediction)
library(stargazer)
library(tidyverse)

1. Model_1 と Model_2 の比較

  • 交差項を含まないロジスティック回帰分析・・・・Model_1

  • 交差項を含むロジスティック回帰分析・・・・Model_2

  • それぞれのモデルで検証すべき仮説が異なる

  • ここでは2つのモデルを比較しながら、交差項を含むロジスティック回帰分析で何が明らかになるかを解説する

1.1 データの準備 (hr96-24.csv)

  • 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))
  • 衆議院議員総選挙の得票データhr96-24.csv をダウンロード

データのダウンロードが終わったら、データを読み込み hr と名前を付ける。

hr <- read_csv("data/hr96-24.csv", na = ".")

データフレーム hr の中身を表示する

names(hr)
 [1] "year"          "pref"          "ku"            "kun"          
 [5] "wl"            "rank"          "nocand"        "seito"        
 [9] "j_name"        "gender"        "name"          "previous"     
[13] "age"           "exp"           "status"        "vote"         
[17] "voteshare"     "eligible"      "turnout"       "seshu_dummy"  
[21] "jiban_seshu"   "nojiban_seshu"
  • select() 関数を使って year, wl, previous, exp という 4 つの変数だけを選ぶ
  • filter() 関数を使って 2021年衆院選だけのデータを残す
  • wl の中身を確認
unique(hr$wl)
[1] 1 0 2

  • 0・・・小選挙区での落選

  • 1・・・小選挙区での当選

  • 2・・・復活当選

  • wl を使って、選挙結果を示す変数 wlsmd （1 = 小選挙区当選、0 = その他）を作る

  • exp を使って、選挙費用を示す変数 expm（単位は「百万円」）を作る

hr21 <- hr %>%
  select(year, wl, previous, exp) %>%
  filter(year == 2021) %>% 
  mutate(wlsmd = ifelse(wl == 1, 1, 0)) |> 
  mutate(expm = exp / 1000000) |> 
  
  select(year, wlsmd, expm, previous)

データフレーム hr21 の中身を表示する

hr21
# A tibble: 857 × 4
    year wlsmd  expm previous
   <dbl> <dbl> <dbl>    <dbl>
 1  2021     1 13.4         3
 2  2021     0  9.62        2
 3  2021     0 NA           0
 4  2021     1  8.82        8
 5  2021     0 11.4         0
 6  2021     1  7.42        8
 7  2021     0 11.9         3
 8  2021     1 13.0         3
 9  2021     0  7.48        6
10  2021     0  4.61        0
# ℹ 847 more rows

データ：

変数名 詳細
year 衆院選が行われた年
wlsmd 小選挙区での当落ダミー(当選 = 1, 落選 = 0)
previous 当選回数
expm 候補者の選挙費用（百万円）
  • データフレーム hr21 の記述統計を表示させる
summary(hr21)
      year          wlsmd             expm              previous     
 Min.   :2021   Min.   :0.0000   Min.   : 0.009319   Min.   : 0.000  
 1st Qu.:2021   1st Qu.:0.0000   1st Qu.: 3.435080   1st Qu.: 0.000  
 Median :2021   Median :0.0000   Median : 5.899882   Median : 1.000  
 Mean   :2021   Mean   :0.3372   Mean   : 6.434585   Mean   : 2.162  
 3rd Qu.:2021   3rd Qu.:1.0000   3rd Qu.: 8.692808   3rd Qu.: 3.000  
 Max.   :2021   Max.   :1.0000   Max.   :27.443685   Max.   :17.000  
                                 NA's   :22
  • データフレーム hr21 の expm に欠損値 (missing data = NA's) が22個あることがわかる
  • na.omit() を使って欠測のない観測だけを残す
hr21 <- na.omit(hr21)
  • データフレーム hr21 の記述統計を表示して確認する
summary(hr21)
      year          wlsmd             expm              previous     
 Min.   :2021   Min.   :0.0000   Min.   : 0.009319   Min.   : 0.000  
 1st Qu.:2021   1st Qu.:0.0000   1st Qu.: 3.435080   1st Qu.: 0.000  
 Median :2021   Median :0.0000   Median : 5.899882   Median : 1.000  
 Mean   :2021   Mean   :0.3461   Mean   : 6.434585   Mean   : 2.211  
 3rd Qu.:2021   3rd Qu.:1.0000   3rd Qu.: 8.692808   3rd Qu.: 3.000  
 Max.   :2021   Max.   :1.0000   Max.   :27.443685   Max.   :17.000
  • 22個の NA's が消えていることがわかる　　

1.2 Model_1: 交差項を含まないモデル

仮説 1・選挙費を使えば使うほど、当選確率は大きい

  • この仮説を検証するために Model_1 を推定する

  • ダウンロードした 2021年衆院選挙データを使い、候補者の小選挙区における当落 (wlsmd) を縦軸に、選挙費用 (expm) を横軸にした散布図を表示する

  • 分析で使うデータ (hr21) の記述統計を表示

stargazer::stargazer(as.data.frame(hr21), 
  type = "html")
Statistic N Mean St. Dev. Min Max
year 835 2,021.000 0.000 2,021 2,021
wlsmd 835 0.346 0.476 0 1
expm 835 6.435 4.287 0.009 27.444
previous 835 2.211 2.921 0 17

Model 1 の分析結果

  • ここでは2021年の衆院選データを用い、選挙費用（expm: 単位100万円）で衆院選小選挙区での当落 (wlsmd) を説明するモデルを推定する
model_1 <- glm(wlsmd ~ expm + previous, 
            data = hr21, 
            family = binomial(link = "logit")) # 係数を「オッズの対数」に指定
  • この結果を stargazerパッケージを使ってまとめてみる
stargazer::stargazer(model_1,
  type = "html")
Dependent variable:
wlsmd
expm 0.170***
(0.024)
previous 0.342***
(0.036)
Constant -2.652***
(0.194)
Observations 835
Log Likelihood -397.969
Akaike Inf. Crit. 801.939
Note: p<0.1; p<0.05; p<0.01

→　表示される回帰式の係数 estimate は「オッズの対数」: log(odds)
= log-odds = logit と呼ばれる

  • ロジスティック回帰分析の係数 log-oddsは、重回帰分析の係数と同じように解釈できない
log-odds は解釈しにくい → オッズ比や確率に変換する

Model_1 の結果

「選挙費用」と「予測当選確率」(hr21_me)

hr21_pred <- cplot(model_1, 
                  x = "expm", 
                  what = "prediction") %>% # Y軸に予測当選確率を表示させる設定  
  as_data_frame() %>%
  ggplot(aes(x = xvals, 
    y = yvals, 
    ymin = lower, 
    ymax = upper)) +
  geom_ribbon(fill = "gray") +
  geom_line() +
labs(x = "選挙費用（単位：100万円）", 
     y = "予測当選確率の予測値（確率: 0 〜 1）", 
     title = "予測当選確率の予測値:2021総選挙")

「選挙費用」と「選挙費用の限界効果」」(hr21_pred)

hr21_me <- cplot(model_1,  
                  x = "expm",    # x軸に据える変数
                  dx = "expm",   # 説明変数
                  what = "effect") %>%　# Y軸に限界効果を表示させる設定  
  as_data_frame() %>%
  ggplot(aes(x = xvals, 
    y = yvals, 
    ymin = lower, 
    ymax = upper)) +
  geom_ribbon(fill = "gray") +
  geom_line() +
  geom_abline(intercept = 0, slope = 0, linetype = "dashed", color = "red") + 
  ylim(-0.0001, 0.06) +
labs(x = "選挙費用（単位：100万円）", 
     y = "選挙費用の限界効果（確率: 0 〜 1）", 
     title = "選挙費用の限界効果:2021総選挙")

仮説検証のまとめ（仮説 1）:

hr21_pred + hr21_me

仮説検証の結論（仮説 1）:当選回数が平均値 (2.211回)の候補者の場合

左側の図からわかること

・縦軸・・・確率（0 〜 1）
・選挙費用が大きくなるにつれて、候補者の当選確率は大きくなる
・選挙費用が小さいとき、予測当選確率は緩やかに上昇
・選挙費用が1000万円の候補者の予測当選確率は0.4933899（= 49%）
・選挙費用が2000万円の候補者の予測当選確率は0.8420812（= 84%）
・選挙費用が中程度になると曲線の傾きが少しずつ急になる
・選挙費用が大きくなると再び傾きが小さくなる

右側の図からわかること

・縦軸・・・確率（0 〜 1）
・選挙費用の増加は当選確率に対して逓減的な効果をもつ
・0円〜1,100万円では限界効果が急上昇するが、1,100万円を超えると追加的効果は徐々に小さくなる
・選挙費用の限界効果は、1,100万円の時が最大: 0.04250（4.25%ポイント）
・選挙費用がそれより多くても少なくても効果が小さい
・全ての選挙費用の範囲で統計的に有意

1.3 Model_2: 交差項を含むモデル

仮説 2・選挙費が当選に与える影響は、当選回数によって異なる

  • この仮説を検証するために Model_2 を推定する

Model 2

Model 2 の分析結果

model_2 <- glm(wlsmd ~ expm*previous, 
            data = hr21, 
            family = binomial(link = "logit")) # 係数を「オッズの対数」に指定
  • この結果を stargazerパッケージを使ってまとめてみる
  • model_1 と model_2 の結果を並べて表示させてみる
stargazer::stargazer(model_1, model_2,
  type = "html")
Dependent variable:
wlsmd
(1) (2)
expm 0.170*** 0.282***
(0.024) (0.033)
previous 0.342*** 0.689***
(0.036) (0.075)
expm:previous -0.044***
(0.008)
Constant -2.652*** -3.404***
(0.194) (0.263)
Observations 835 835
Log Likelihood -397.969 -382.279
Akaike Inf. Crit. 801.939 772.557
Note: p<0.1; p<0.05; p<0.01
  • カッコの中は標準誤差 (standard error) を表す
  • 表示される回帰式の係数 estimate は「オッズの対数」: log(odds)
    = log-odds = logit と呼ばれる
    →　重回帰分析の係数と同じように解釈できない
  • model_2 の結果（図の右側）が意味しているのは
    ⇒ 調整変数が 0 の時 (previous = 0)の結果
  • つまり、previous = 0の候補者が100万円使うと当選確率に正の影響 (0.282) がある

結果のまとめ（図の右側の結果）

  • いずれの結果も統計的に有意
expmの係数 (0.282): expm が増えると当選確率は増加する
previous の係数 (0.689): previous が多いほど当選確率は増加する
expm:previousの係数 (-0.044): previous が大きくなると expm の効果は弱まる
  • しかし、数値は「オッズの対数」なので、確率 (%) に変換して可視化する必要がある

margins::cplot を使って「予測当選確率」と「限界効果」を可視化する

「選挙費用」と「予測当選確率」(hr21_pred)

  • 注意：チャンクオプションに```{r, echo = TRUE, fig.show = "hide"}と指定する
hr21_pred <- margins::cplot(model_2, 
                  x = "expm", 
                  what = "prediction") %>% # Y軸に予測当選確率を表示させる設定  
  as_data_frame() %>%
  ggplot(aes(x = xvals, 
    y = yvals, 
    ymin = lower, 
    ymax = upper)) +
  geom_ribbon(fill = "gray") +
  geom_line() +
labs(x = "選挙費用（単位：100万円）", 
     y = "予測当選確率の予測値", 
     title = "予測当選確率の予測値:2021総選挙")
  • expm が増えると当選確率は増加している

「選挙費用」と「選挙費用の限界効果」(hr21_me)

  • 注意：チャンクオプションに```{r, echo = TRUE, fig.show = "hide"}と指定する
hr21_me <- margins::cplot(model_2,  
                  x = "previous",    # x軸に据える変数
                  dx = "expm",   # 説明変数
                  what = "effect") %>%　# Y軸に限界効果を表示させる設定  
  as_data_frame() %>%
  ggplot(aes(x = xvals, 
    y = yvals, 
    ymin = lower, 
    ymax = upper)) +
  geom_ribbon(fill = "gray") +
  geom_line() +
  geom_abline(intercept = 0, slope = 0, linetype = "dashed", color = "red") + 
  ylim(-0.05, 0.05) +
labs(x = "当選回数", 
     y = "選挙費用の限界効果（当選回数ごと）", 
     title = "選挙費用の限界効果:2021総選挙")

仮説検証のまとめ（仮説 2）:

hr21_pred + hr21_me

仮説検証の結論（仮説 2）

左側の図からわかること

・選挙費用が大きくなるにつれて、候補者の当選確率は大きくなる
・選挙費用が小さいとき、予測当選確率は緩やかに上昇
・選挙費用が中程度になると曲線の傾きが少しずつ急になる
・選挙費用が大きくなると再び傾きが小さくなる

右側の図からわかること

・選挙費用の限界効果は、当選回数が増える程小さくなる
・当選回数が5回程度まではこの傾向が続き、統計的にも有意
・当選回数が6回〜8回あたりまでは統計的に有意ではない
・当選回数が9回以降は選挙費用が当選回数に負の影響を与える

選挙費用が予測当選確率に与える「影響」を当選回数別に可視化する

  • 綺麗に可視化するためには、expm と previous の記述統計をチェックする必要がある
stargazer::stargazer(as.data.frame(hr21),
  type = "html")
Statistic N Mean St. Dev. Min Max
year 835 2,021.000 0.000 2,021 2,021
wlsmd 835 0.346 0.476 0 1
expm 835 6.435 4.287 0.009 27.444
previous 835 2.211 2.921 0 17
  • 下のスクリプトで当選回数 (preious) と選挙費用 (expm) を適切な大きさに分割する際、上の記述統計が参考になる
df_pre <- expand.grid(previous = seq(0, 17, by = 2), # 当選回数を0から17回まで2回ずつ離して指定
  expm = seq(0, 27, by = 0.1)) |> # 選挙費用を0から1700万円まで200万円ずつ離して指定
  as_data_frame() 
pred <- predict(model_2, 
  type = "response", 
  newdata = df_pre, 
  se.fit = TRUE) 
df_pre$fit <- pred$fit 
df_pre$lower <- with(pred, fit - 2 * se.fit) 
df_pre$upper <- with(pred, fit + 2 * se.fit) 
df_pre <- df_pre %>%

mutate(lower = ifelse(lower < 0, 0, lower), 
  upper = ifelse(upper > 1, 1, upper)) 
plt_prob <- ggplot(df_pre, aes(x = expm, y = fit)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = "gray") +
geom_line() +
facet_wrap(. ~ previous) + 
  labs(x = "選挙費用（100万円）", y = "当選確率の予測値") +
ggtitle("選挙費用の限界効果（当選回数ごと）")

print(plt_prob)

この図からわかること・当選回数が 0 回の候補者は、選挙費用が増えると当選確率高くなっている
・当選回数が 6 回の候補者はお金を使っても使わなくてもほとんど当選確率に変化はなく、高い当選確率を維持している
・当選回数が 8 回を超える候補者だと、選挙費用はむしろマイナスの影響すら与えるようになることが読み取れる

1.4 モデルの当てはまり具合を評価する

model_1 <- glm(wlsmd ~ expm + previous, 
            data = hr21, 
            family = binomial(link = "logit")) # 係数を「オッズの対数」に指定

model_2 <- glm(wlsmd ~ expm*previous, 
            data = hr21, 
            family = binomial(link = "logit")) # 係数を「オッズの対数」に指定

ROC曲線によるモデルの当てはまり具合の評価方法

  • 交差項を含まない model_1 と交差項を含む model_2 のどちらの当てはまりがよいかを、ROC曲線を描いて考える
  • ROC曲線を描くために ROCR パッケージを使う
  • prediction() という名前の関数は ROCR パッケージだけでなく margins パッケージにもある
    →どちらの関数を使うか明記する　　
pi1 <- predict(model_1, type = "response")
pi2 <- predict(model_2, type = "response")
pr1 <- ROCR::prediction(pi1, labels = hr21$wlsmd == "1")
pr2 <- ROCR::prediction(pi2, labels = hr21$wlsmd == "1")
roc1 <- performance(pr1, measure = "tpr", x.measure = "fpr")
roc2 <- performance(pr2, measure = "tpr", x.measure = "fpr")
df_roc <- data_frame(fpr = c(roc1@x.values[[1]], roc2@x.values[[1]]),
                     tpr = c(roc1@y.values[[1]], roc2@y.values[[1]])) %>% 
  mutate(model = rep(c("Model 1", "Model 2"), c(n()/2, n()/2)))
roc <- ggplot(df_roc, aes(x = fpr, y = tpr,
                          color = model, linetype = model)) + 
  geom_line() +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  scale_linetype_discrete(name = "") +
  scale_color_discrete(name = "") +
  coord_fixed() +
  labs(x = "偽陽性率（1 - 特異度）", y = "真陽性率（感度）")
print(roc)

ROC曲線による当てはまり オレンジ色 (model_1) も青色 (model_2) のROC 曲線もどちらも同じくらい 45 度線（点線）から点(0, 1) のほうに離れており、モデルの当てはまりが同じくらい良い

数値指標によるモデルの当てはまり具合の評価方法

AUC

  • 目視だけで当てはまり具合の良し悪しを判断するのには限界がある
    → AUC（area under the curve：ROC 曲線の下側の面積）を使って評価する

  • 上のROCを求めた図中の 0 ≤x ≤1、0 ≤y ≤1 の範囲で ROC 曲線より下の面積を求める

  • すべての ROC 曲線が2 点(0, 0) と(1, 1) を通り、当てはまりのよいモデル ほど(1, 1) の近くを通る

  • ROC 曲線が 3 点(0, 0)、(0, 1)、(1, 1) を通るなら →曲線の下側の面積は 1

  • 当てはまりのよいモデルの AUC → 1 に近くなる

  • 当てはまりの悪いモデルの AUC → 0.5 に近くなる
    　　（←ROC 曲線は 45 度線に近づくため）

  • ROC曲線で見る限り、二つのモデルの当てはまりのよさに大きな差はなさそう

  • 念のため、AUCを計算する

auc1 <- performance(pr1, measure = "auc")
auc1@y.values[[1]]   # model_1 のAUC
[1] 0.853502
auc2 <- performance(pr2, measure = "auc")
auc2@y.values[[1]]   # model_2 のAUC
[1] 0.8552955

AUC による当てはまり ・AUCを見る限り、ほんのわずかだけmodel_2 の当てはまりの方がよい
・ほとんど変わらない

赤池情報量基準 (AID)

  • 赤池情報量基準 (AID: Akaike information criterion) を使う
  • 多くの説明変数を使っても(= パラメータ数が増えても)、それに応じた分だけ尤度が増加しなければ、そのモデルは「良いモデル」とは言えない
  • この基準で評価するのが AID
  • AIC が最小のモデルが「最も良いモデル」とされる
AIC(model_1)
[1] 801.9386
AIC(model_2)
[1] 772.5572

AID による当てはまり model_2 の当てはまりの方がよい

2. Excercise

  • 2021年総選挙データを使って、立候補者が小選挙区で当選したか否か（wlsmd）を応答変数、過去の当選回数（previous）と選挙費用（expm）を説明変数としたロジスティック回帰分析を実行し、以下の各問に答えなさい。
  • 分析で使うデータは衆院選選挙データは hr96-24.csv からダウンロード。

Q1.当選予測確率を縦軸に、選挙費用（expm）を横軸にとった散布図（wlsmdのオッズをとった曲線）を描きなさい。

Q2.立候補者が小選挙区で当選したか否か（wlsmd）を応答変数、選挙費用（expm）を説明変数としたモデルに model_1 という名前を付けてロジスティック回帰分析し、stargazer()関数を使ってその結果を表示しなさい。

Q3.立候補者が小選挙区で当選したか否か（wlsmd）を応答変数、選挙費用（expm）を説明変数、そして過去の当選回数（previous）を調整変数として expm:previous という交差項をモデルに加え model_2 という名前を付けてロジスティック回帰分析し、stargazer()関数を使って model_1 と model_2 の結果を並べて表示させなさい。

Q4. model_2 の分析結果の expm と previous の係数が何を意味しているのかを説明し、この結果からわかることを述べなさい。

Q5. margins::cplot()関数を使って model_2 の推定結果を表示し結論を述べなさい。分析結果は patchworkパッケージを使って表示させること。

Q6. model_1 と model_2 それぞれの ROC 曲線を比較して示し、どちらがよりよいモデルかこたえなさい。

Q7. model_2 において、当選回数ごとの予測当選確率とその統計的有意性を margins::cplot() 関数を使って推定しなさい。

Q8. 選挙費用が予測当選確率に与える「影響」を当選回数別に可視化して示し、その傾向を簡単に説明しなさい。

Q9. model_1 と model_2 それぞれの ROC 曲線を比較して示し、どちらがよりよいモデルかこたえなさい。

Q10. model_1 と model_2 それぞれの 数値指標によるモデルの当てはまり具合の評価方法 (AUC と AID) を比較して示し、どちらがよりよいモデルかこたえなさい。

参考文献