Rを使った分析の準備

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

1. 線形回帰と非線形回帰の違い

  • 選挙での「候補者の成績」と「選挙で使うお金の額」の関係を調べたいとする
    →　つまり、お金を使う候補者ほど選挙結果がいいのか？という疑問

  • 説明したい対象（＝応答変数）としては 2 つ想定できる：
    ①得票率（もしくは得票数）　　　　　・・・連続変数 => 回帰分析
    ②当落・・・当選したら1 、落選したら 0 ・・・2値変数 => ロジスティック回帰分析

  • 回帰分析とロジスティック回帰分析の違いを視覚的に示してみる

回帰分析 ロジスティック回帰分析
応答変数 得票率(%) 当落(0 or 1)
説明変数 選挙費用（百万円） 選挙費用（百万円）

主要な違い

・回帰分析の応答変数は得票率＝「連続変数」
・ロジスティック回帰分析の応答変数は選挙での当落＝ 2 値変数 (1 or 0)
・それぞれの図に回帰直線を図に書き入れてみる

・回帰分析では、選挙費用を使うほど得票率が大きくなっている => 問題なし
・ロジスティック回帰分析でも、選挙費用を使うほど得票率が大きくなっている
→　しかし、3つ問題点がある
1. 当選確率（0〜1）が 1 (100%) を超えてしまうこと
2. 当選確率（0〜1）がマイナスの値をとってしまうこと
3. 当てはまり感がない（むりやり当てはめている感が強い）

2. ロジスティック回帰分析とは

  • 「ロジスティック回帰分析」とは、2値変数をロジット変換した回帰分析

ロジット変換とは

  • 2値を表す散布図を「率」に変換すると S 字になる（左の図 → 中央の図）
  • 率をロジット変換すると直線になる（中央の図 → 右の図）
  • ロジスティック回帰分析ではこの数学の特徴を使う

  • 選挙での「当落」（左側の図）の代わりに「当選率」を考えてみる（中央の図）
  • 選挙費用ごとにどれくらいの割合が当選しているかを棒グラフで表してみる
  • 中央の「当選率」の棒グラフを見ると・・・

・選挙費用が少ない200万円以下の人はほとんど当選していない
・選挙費用が1000万円の候補者は50%程当選している
・選挙費用が1600万円の候補者は75%以上が当選している
・選挙費用ごとの「当選率」の対数をとってロジット変換する
→　Ｓ字が直線になる
→　しかし、実際のデータを使うと完璧な直線ではないものの、かなり直線ぽくなる
→　この曲線を使ってロジスティック回帰分析を推定する

  • 2値変数をロジットに変換してから回帰分析を行う→ロジスティック回帰分析と呼ばれる

・ロジスティック回帰分析では「割合」ではなく「オッズ」で結果が示される

Logit Model の基礎知識：オッズ

オッズ (Odds) とは

  • 例えば、選挙費用を800万円使う候補者は100人いる
  • そのうち当選したのは60人で落選したのは40人
  • この場合の当選者の「割合」と「オッズ」は次のように計算できる

  • 「割合」は 1 を超えることはない
  • しかし、「オッズ」は 1 を超えることはある
  • 「割合」も「オッズ」もどちらも起こりやすさを示す指数

オッズのまとめ

  • ある事象が起こる確率 \(\mathrm{p}\) と起こらない確率 \(\mathrm{1 - p}\) との比

p: ある事象が起こる確率
1-p: ある事象が起こらない確率

\[Odds = \frac{p}{1-p}\]

  • 例えば、p = 0.01（起こる確率が1%）の場合の Odds を計算してみる
  • p = 0.01, 1-p = 1-0.01= 0.99 を Odds を求める式に代入する

\[Odds = \frac{p}{1-p} = \frac{0.01}{0.99} = \frac{1}{99} = 0.01\]

Odds の特徴・オッズが大きいほど、その事象 (p) が起こりやすい
・オッズの最小値は 0
・オッズ = 1 は「その事象 (p) が起こる確率が 50%」という意味
・オッズの最大値は無限大 (∞)

Logit Model の基礎知識：ロジット変換

ロジット変換とは

・ロジット変換とは 2 値変数 (0 or 1) のオッズの対数をとること
・Odds の下限は 0 なので、説明変数としては扱いにくい
→ Odds [= p / (1 - p)] を対数変換し「対数オッズ」(Log Odds)を計算する
→　これが logit

\[LogOdds =log\frac{p}{1-p} = logit\]

  • ロジット (logit) は対数オッズ (Log Odds) として解釈できる

Log Odds (= logit) の特徴 ・最小値も最大値もどちらも無限大 (∞)
・Log Odds = 0 が意味するのは「その事象 (p) が起こる確率が 50%」という意味

3. ロジスティック回帰分析の手順

ロジスティック回帰分析の手順は次のとおりである

  1. 対立仮説と帰無仮説を設定する
  2. 説明変数と応答変数の散布図を表示する
  3. ロジスティック回帰式を求める
  4. 回帰係数の有意性を検定する
  5. 推定結果の意味を解釈する
  6. ロジスティック回帰モデルの当てはまり具合を評価する

3.1 対立仮説と帰無仮説を設定する

  • ここで検証したい仮説（＝対立仮説）は次のとおり

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

  • 従って、否定されるために設定する帰無仮説は次のようになる

帰無仮説 ・選挙費の額は、小選挙区での当選確率とは関係がない

  • ロジスティック回帰分析における検定では、重回帰分析における検定と同様、得られた p 値 が有意水準よりも小さいときに帰無仮説を棄却し、対立仮説を受け容れる

  • 仮説1を検証するために次のモデルを考える

Model 1

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

データの準備 (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 が消えていることがわかる　　

3.2 説明変数と応答変数の散布図を表示する

  • 分析で使うデータ (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
  • 当選回数と当落の散布図を描く.
  • 2つの変数の関係を図示してみる
  • 文字化けを避けるため、Macユーザは次の2行を実行する
theme_set(theme_gray(base_size = 10, 
                     base_family = "HiraginoSans-W3"))
  • wlsmd が 2 値変数 (0 or 1) なので、jitter()関数を使ってデータを散らして表示させる
p <- ggplot(hr21, aes(x = expm, y = wlsmd)) + 
  geom_jitter(size = 1,        # データを散らして表示させる指定
              alpha = 1/3,
              width = 0,
              height = 0.05) +
  labs(x = "選挙費用（100万円）",
       y = "小選挙区での当落")

plot(p)

  • この図にむりやり通常の回帰直線を当てはめてみるとこのようになる
p + geom_smooth(method = "lm", se = FALSE) +
    annotate("label", 
           label = "当落 (0 or 1) = a + b選挙費用", 
           x = 6, y = 0.75,
           size = 5, 
           colour = "blue", 
           family = "HiraginoSans-W3")

  • 応答変数が2値 (o or 1) なのに、むりやり回帰直線を当てはめた時の問題点
    → 応答変数（当選確率）がマイナスの値になったり、100%を超える場合がある

解決策

  • 当選確率を「ロジット変換」する
  • 当選確率 \(P\) は 0 と 1 の間の値をとる 2 値変数
  • \(\frac{P}{1-P}\) はオッズ（コラムLogit Model の基礎知識：オッズを参照）
  • オッズ＝当選する確率 \(\mathrm{p}\) と落選する確率 \(\mathrm{1 - p}\) との比
  • 「当選確率\(P\)」と「選挙費用」の関係は次の式で表現できる

\[\frac{P}{1-p} = a + 選挙費用b + 当選回数c\]

  • 回帰式の左辺のオッズ \(\frac{p}{1-P}\)は曲線を表現している
  • 回帰式の右辺 \(a + 選挙費用b\) は直線を表現している
    →　表現が一致していない
  • 回帰式の左辺をする＝オッズ \(\frac{p}{1-P}\)の対数をとる
  • 当選確率 \(P\) は 0 から 1 の範囲しか取れないが、オッズの対数をとると全ての実数をとれる
    →　曲線の関係を直線の関係に変更でき、推定が可能になる ＝「ロジット変換」

\[log\frac{P}{1-p} = a + 選挙費用b + 当選回数c\]

オッズの対数をとった曲線

p3 <- ggplot(hr21, aes(x = expm, y = wlsmd)) + 
  geom_jitter(size = 1,
              alpha = 1/3,
              width = 0,
              height = 0.05) +
  geom_smooth(method = "glm", 
    color = "red",
    method.args = list(family = binomial(link = "logit"))) +
  labs(x = "選挙費用（100万円）",
       y = "小選挙区での当選確率")
print(p3)

3.3 ロジスティック回帰式を求める

ロジスティック回帰式 (Model 1)

  • ここでは2021年の衆院選データを用い、選挙費用（expm: 単位100万円）で衆院選小選挙区での当落 (wlsmd) を説明するモデルを考える

Model 1

  • このモデルを式で表すとこのようになる

\[Pr(当選)=logit^{−1}(𝛼+\beta_1選挙費用 + \beta_2当選回数)\]

  • Rでロジスティック回帰分析を行うには、一般化線形モデル (Generalized Linear Models: GLM) を当てはめるための関数 glm() を使う
  • glm()はロジスティック回帰以外でも頻繁に使う関数
  • この関数をロジスティック回帰分析に使うには、引数 family = binomial(link = "logit") を指定
  • logit とは「オッズの対数 」= Log Odds のこと
model_1 <- glm(wlsmd ~ expm + previous, 
            data = hr21, 
            family = binomial(link = "logit")) # 係数を「オッズの対数」に指定
  • tidy() を使って、推定結果を確認する
tidy(model_1, conf.int = TRUE)
# A tibble: 3 × 7
  term        estimate std.error statistic  p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 (Intercept)   -2.65     0.194     -13.7  9.90e-43   -3.04     -2.28 
2 expm           0.170    0.0235      7.23 4.86e-13    0.125     0.217
3 previous       0.342    0.0355      9.62 6.38e-22    0.274     0.413

→　表示される回帰式の係数 estimate は 「オッズの対数」 = Log Odds

  • Log Odds は解釈しにくい → 確率に変換する

  • 当選確率を予測するためのロジスティック回帰式は次のとおり

\[Pr(当選) = logit^{-1}(\alpha + \beta_1選挙費用 + \beta_2当選回数)\]

\[= \frac{1}{1+exp(-[\alpha + \beta選挙費用])}\] \[= \frac{1}{1+exp(-[-2.65 + 0.170[expm] + 0.342[previous]])}\]

3.4 回帰係数の有意性を検定する

  • tidy() を使って、推定結果を再度確認してみる
tidy(model_1, conf.int = TRUE)
# A tibble: 3 × 7
  term        estimate std.error statistic  p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 (Intercept)   -2.65     0.194     -13.7  9.90e-43   -3.04     -2.28 
2 expm           0.170    0.0235      7.23 4.86e-13    0.125     0.217
3 previous       0.342    0.0355      9.62 6.38e-22    0.274     0.413

  • expm の p.value は 1.46e-24
    →　有意水準を 0.05 とすれば、選挙費用 (expm) の p値は 1.61e- 5 なので、「影響がない」という帰無仮説を棄却する
    → 選挙費用が当選確率に影響すると判断する

  • previous の p.value は 4.29e-14
    →　有意水準を 0.05 とすれば、過去の当選回数 (previous) の p値は 4.29e-14 なので、「影響がない」という帰無仮説を棄却する
    → 過去の当選回数が当選確率に影響すると判断する

ロジスティック回帰分析の係数

→　表示される回帰式の係数 estimate は 「オッズの対数」 = Log Odds

  • Log Odds は解釈しにくい → 確率に変換する

  • 当選確率を予測するためのロジスティック回帰式は次のとおり

\[Pr(当選) = logit^{-1}(\alpha + \beta_1選挙費用 = \beta_2当選回数)\]

\[= \frac{1}{1+exp(-(\alpha + \beta選挙費用)}\] \[= \frac{1}{1+exp(-(-2.65 + 0.170expm + 0.342previous))}\]

  • 例えば、選挙費用 (expm) の係数の推定値は約 0.170
  • これが重回帰分析の結果であれば、「選挙費用が 1 単位（100 万円）増えるごとに、当選確率が 0.170 ポイント大きくなる」という意味
  • しかし、ロジスティック回帰分析の係数は、重回帰分析の係数と同じように解釈できない

予測当選確率と限界効果

予測当選確率

  • 選挙費用の値ごとに候補者が当選する確率の予測値が当選予測確率: 0 〜 1

限界効果

  • ロジスティック回帰分析の係数は、重回帰分析の係数と同じように解釈できない

  • ロジスティック曲線は直線ではない → 「曲線」
    → 説明変数 1単位の変化が応答変数に与える影響（＝傾き）は説明変数の値によって異なる

  • 説明変数が応答変数に与える影響を考えるためには
    →　説明変数の値を特定し、「特定された値における影響」を計算する必要あり
    →　限界効果 (Marginal Effects)

  • 上図の例では、xの値 1, 5, 10, 15, 20 それぞれの時点での直線の傾きが限界効果

予測当選確率

  • 当選回数 (previous) が3回の候補者に関して次の二つのケースを考えてみる
    例1: 選挙費用を0 円から 100 万円に変化させた場合の予測当選確率
    例2: 選挙費用を100 万円から 200 万円に変化させた場合の予測当選確率

例1: 選挙費用 (0 円 → 100 万円)

  • ここでは当選回数が 3 回の候補者が、選挙費用を0 円から 100 万円に変化させると、応答変数にどのような変化があるか確かめる
    →それぞれの予測当選確率を計算して差を取る
  • 当選回数が 3 回で選挙費用が 0 円の候補者の予測当選確率は、次の回帰式に、previous = 3、expm = 0 を代入すれば求められる

\[Pr（当選） = \frac{1}{1+exp(-(-1.98 + 0.0735expm + 0.285previous))}\]

  • この値を R で求めてみる
  • 当選回数が 3 回で選挙費用が 0 円の候補者の予測当選確率は
p_0 <- predict(model_1, type = "response", #予測当選確率を表示したい時の指定  
               newdata = data_frame(previous = 3, expm = 0))
p_0
        1 
0.1642294
  • 0.1642294 (= 約16.4%) である
  • 当選回数が 3 回で選挙費用が 100 万円の候補者の予測当選確率は 　　
p_1 <- predict(model_1, type = "response",
               newdata = data_frame(previous = 3, expm = 1))
p_1
        1 
0.1889165
  • 0.1889165 (= 約18.9%)だとわかる
  • 当選回数が 3 回の候補者が選挙費用を 0 円から 100 万円に変化させたときの予測当選確率の変化は、上で求めた二つの確率の差 (p_1 - p_0)
p_1 - p_0
         1 
0.02468706
  • 当選回数が 3 回の候補者が選挙費用を 0 円から 100 万円に増やすと、予測当選確率は約 2.47 %ポイント増える

例2: 選挙費用 (100万円 → 200万円)

  • ここでは当選回数が 3 回の候補者が、選挙費用を100万円から 200万円に変化させると、応答変数にどのような変化があるか確かめる
  • 当選回数が 3 回で選挙費用が 200万円の候補者の予測当選確率は 0.2163539（約21.6%ポイント）
p_2 <- predict(model_1, type = "response",
               newdata = data_frame(previous = 3, expm = 2))
p_2
        1 
0.2163539
  • 当選回数が 3 回の候補者が選挙費用を 100 万円から 200 万円に変化させたときの予測当選確率の変化は、上で求めた二つの確率の差 (p_2 - p_1) なので
p_2 - p_1
         1 
0.02743738
  • 当選回数が 3 回の候補者が選挙費用を 100 万円から 200万円に増やすと、予測当選確率は 0.02743738 (= 約2.7% ポイント) 増える

まとめ ・選挙費用を 0 円から 100 万円に増やすと、予測当選確率は 2.47%ポイント上昇する
・選挙費用を 100 万円から 200 万円に増やすと、予測当選確率は 2.7%ポイント上昇する
・選挙費用の0 円から 100 万円への変化も、100 万円から 200 万円への変化も、選挙費用の増分としては同じ
・しかし、それらの変化が応答変数（すなわち当選確率）に与える影響は異なる
・0 円から 100 万円への 100 万円の変化より、100 万円から 200 万円への 100 万円の変化のほうが、応答変数に与える影響が大きい
・実際に何%ポイント予測当選確率が上昇するかを、具体的な予測当選確率を計算せずに係数の推定値(expm = 0.170) から読み取るのは難しい

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

・この推定値は「オッズの対数 」= Log Odds = logit
→ いくつか特定した expm の値で予測確率を計算し、その差を比較する作業が必要

限界効果

margins を使った限界効果の計算

  • 上の例のように、予測当選確率を複数計算し、その差をとるという作業を繰り返すのは面倒
    → 特定の値における限界効果を直接計算するための marginsパッケージを使うのが便利
  • margins()関数を使うと、特定の選挙費用額ごとの限界効果を求めることができる

選挙費用の限界効果 (previous = 3, expm = 400万円)

  • 例えば、当選回数が 3 回で選挙費用が 400 万円のとき、選挙費用 1 単位（= 100 万円）の増加が当選確率に与える影響の大きさ（＝限界効果）は
margins(model_1, variables = "expm",
        at = list(previous = 3, expm = 4))
 at(previous) at(expm)    expm
            3        4 0.03424
  • 選挙費用が当選確率に与える影響（限界効果）は 0.03424 (約3.4%ポイント) だとわかる

選挙費用の限界効果（previous = 3, expm: 0 〜 2,800万円)

  • どの値での限界効果を求めるかは、引数 atに変数のリストを渡すことで指定する
  • 複数の値、例えば、選挙費用が 0円 (previous = 0) から 2,800万円 (previous = 3) までの限界効果を同時に表示してみる
margins(model_1, variables = "expm",
at = list(previous = 3, expm = c(0:28))) #選挙費用を0円から2,800万円まで指定  
 at(previous) at(expm)     expm
            3        0 0.023337
            3        1 0.026052
            3        2 0.028827
            3        3 0.031585
            3        4 0.034239
            3        5 0.036685
            3        6 0.038820
            3        7 0.040540
            3        8 0.041758
            3        9 0.042406
            3       10 0.042449
            3       11 0.041884
            3       12 0.040742
            3       13 0.039087
            3       14 0.037004
            3       15 0.034595
            3       16 0.031964
            3       17 0.029215
            3       18 0.026437
            3       19 0.023710
            3       20 0.021093
            3       21 0.018631
            3       22 0.016352
            3       23 0.014273
            3       24 0.012397
            3       25 0.010723
            3       26 0.009241
            3       27 0.007938
            3       28 0.006801
  • ロジスティック回帰分析の係数の限界効果は、説明変数の値によって変化する
    → 各説明変数の係数を見ただけで結果を理解するのは困難
    → 説明変数が応答変数に影響を与える様子をいくつかの条件について図示する
    → 説明変数の値に応じた限界効果 (ME:Marginal Effects) を図示する必要がある

3.5 推定結果の意味を解釈する

仮説 1 の結果の表示方法 (Model_1)

「予測当選確率」と「限界効果」は異なる

  • 「予測当選確率」は特定の選挙費用を使った候補者が当選する確率

  • 「限界効果」は選挙費用が予測当選確率に与える影響（＝傾き）

  • 予測当選確率と限界効果は margins::cplot() を利用して図示する
    ・Y軸に予測当選確率を表示したい時は what = "prediction" と指定
    ・Y軸に限界効果を表示したい時は what = "effect" と指定

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

hr21_me <- cplot(model_1, 
                  x = "expm", 
                  what = "prediction") %>%
  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総選挙")

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

hr21_pred <- cplot(model_1,  
                  x = "expm",    # x軸に据える変数
                  dx = "expm",   # 説明変数
                  what = "effect") %>%
  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 = "選挙費用", 
     y = "選挙費用の限界効果", 
     title = "選挙費用の限界効果:2021総選挙")

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

hr21_me + hr21_pred

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

左側の図からわかること

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

右側の図からわかること

実質的有意性:
・選挙費用の限界効果は、約1000万円の時が最大
・選挙費用がそれより多くても少なくても効果が小さい
→　選挙費用が約1300万円の時、当選確率に対する影響力が最大
統計的有意性: 全ての選挙費用の範囲で統計的に有意

  • 選挙費用の限界効果の最小値と最大値を確認したければ、次のように入力する
  • この限界効果は、モデルに含まれているコントロール変数 (previous) の値を平均値に固定した値である
margins(model_1, variables = "expm",
  at = list(expm = c(0:28))) #選挙費用を0円から2,800万円まで指定  
 at(expm)    expm
        0 0.01911
        1 0.02082
        2 0.02257
        3 0.02435
        4 0.02614
        5 0.02789
        6 0.02959
        7 0.03119
        8 0.03264
        9 0.03391
       10 0.03494
       11 0.03569
       12 0.03611
       13 0.03618
       14 0.03587
       15 0.03518
       16 0.03413
       17 0.03275
       18 0.03107
       19 0.02916
       20 0.02708
       21 0.02489
       22 0.02266
       23 0.02045
       24 0.01831
       25 0.01626
       26 0.01435
       27 0.01259
       28 0.01099
  • expm = 13 の時に限界効果が最大値 (0.03618) であることが確認できる

当選回数ごとの予測当選確率

予測当選確率

  • 選挙費用の額に関わりなく、選挙費用が予測当選確率に正の影響を与えていることはわかった
  • それは、候補者の当選回数が異なってもいえるのか？
  • 候補者の当選回数ごとに選挙費用が予測当選確率に与える影響をチェックする
  • 候補者の当選回数 (previous) の記述統計を確認する
table(hr21$previous)

  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  16  17 
390  79  55 106  51  31  33  30  24  16   7   3   5   2   1   1   1 
hist(hr21$previous)

  • 候補者の大多数の当選回数が 0 だとわかる

  • 17回の当選回数を0から2回ずつ区切り、0, 2, 4, …, 16まで分割する
    →　previous = seq(0, 16, by = 2)

  • 候補者の選挙費用 (expm) の記述統計を確認する

summary(hr21$expm)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
 0.009319  3.435080  5.899882  6.434585  8.692808 27.443685 
hist(hr21$expm)

- 選挙費用は0円から2700万円強まで分布
→ 500万円ずつ分割する →　expm = seq(0, 20, by = 5))
- 当選回数ごとの予測当選確率を可視化してみる

df_pre <- expand.grid(previous = seq(0, 17, by = 2), 
  expm = seq(0, 20, by = 5)) %>% 
  as_data_frame() 
pred <- predict(model_1, 
  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 = "当選確率の予測値") 
print(plt_prob)

  • 当選回数が少ない候補者は、選挙費用が当選確率の予測値に強い正の影響を与えている
  • 当選回数が増えるにつれて、選挙費用が当選確率の予測値に与える影響が小さくなる
    →　当選回数16回の候補者は、選挙費用を使わなくても使ってもほぼ当選している

限界効果

  • 次に、当選回数が変わっても、選挙費用が当選確率に予測値に与える影響が統計的に有意なのかどうかをチェック

  • 当選回数 (previous) が増えるにつれて、選挙費用 (expm) の影響がどう変化するかを指定
    → x = "previous", dx = "expm"

plt_previous <- cplot(model_1, 
                  x = "previous", # x軸に据える変数
                  dx = "expm", 　 # 説明変数
                  what = "effect") %>%
  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") + 
labs(x = "当選回数", 
     y = "選挙費用の限界効果", 
     title = "当選回数ごとの選挙費用の限界効果:2021総選挙")
print(plt_previous)

  • previous = 15までは、95%信頼区間が y = 0 の赤い点線に触れていない
    → 当選回数が変わっても、previous = 15 までは、選挙費用が当選確率に予測値に与える影響が統計的に有意

当選回数ごとの予測当選確率：まとめ

plt_prob + plt_previous

左側の図からわかること
  • 当選回数が少ない候補者は、選挙費用が当選確率の予測値に強い正の影響を与えている
  • 当選回数が増えるにつれて、選挙費用が当選確率の予測値に与える影響が小さくなる
    →　当選回数16回の候補者は、選挙費用を使わなくても使ってもほぼ当選している
右側の図からわかること
  • このことは当選回数0回から15回までの候補者に当てはまり、統計的に有意
  • 選挙費用が当選確率の予測値に与える影響は当選回数が5回の候補者が最も大きい

限界効果の統計的有意性 ・限界効果は、説明変数 1 単位の増加が応答変数を何単位増加させるかを示す
・ここでの説明変数は「選挙費用」

注意：「限界効果」と「予測当選確率」とは異なる
  • 「限界効果」は選挙費用が予測当選確率に与える影響（＝傾き）
  • 「予測当選確率」は特定の選挙費用を使った候補者が当選する確率

→ 限界効果は説明変数の値によって変化する
→ 説明変数の値に応じた限界効果を示す必要がある
・図中の95% 信頼区間を見ると、有権者数によって信頼区間の幅が異なる
← 限界効果の標準誤差が、説明変数である選挙費用の値に応じて変わるため

- 主な説明変数が応答変数に与える影響が統計的に有意かどうかの判断

→ 説明変数の値によって標準誤差（信頼区間）が変化する様子も示す必要がある
・上の事例では、観察された有権者数の範囲で95%信頼区間全体が 0 より大きい範囲にある
→ 選挙費用が選挙の当落に与える影響は、選挙費用額にかかわらず統計的に有意
(観察された選挙費用の範囲で95%信頼区間全体が 0 より小さい範囲にある時も統計的に有意)

・上図のように、主な説明変数が応答変数に与える影響（限界効果）の符号が、説明変数の値によって変わりうる
・説明費用の値が 1 より小さければ限界効果の値はマイナス、それより大きければ限界効果の値はプラス

・上図に示されているように、ロジスティック回帰分析では、主な説明変数の効果が説明変数の値によって変わるだけでなく、統計的に有意な範囲と有意でない範囲の両方をもつことがあり得る
・説明変数の値が 1 より小さい、もしくは 2 より大きければ統計的に有意
・説明変数の値が 1 と 2 の間では統計的に有意ではない
・説明変数がどの範囲の値をとると限界効果が統計的に有意になるかどうかを、回帰分析の係数の推定値を見ただけで判断することは非常に難しい
→ 限界効果を図示してはじめて、限界効果がどの範囲でどのような符号をもつか、どの範囲で統計的に有意かが明らかになる

参考（オッズ比を使った分析結果の解釈）

  • ここでは選挙費用と選挙の当落の関係を「確率」を使って分析した
  • しかし、オッズ比を使った分析も可能なので、簡単に紹介する

Logit Model の基礎知識：オッズ比 ・ある事象の起こりやすさを 2 つの群で比較して示す尺度

第 1 群の確率- - - \(\mathrm{p}\)
第 2 群の確率- - - \(\mathrm{q}\)

この場合の Odds Ratio は次の式で表すことができる　

\[OddsRatio=\frac{Odds[p]}{Odds[q]}\ =\frac{p/(1-p)}{q/(1-q)}\ = \frac{p(1-q)}{(1-p)q}\]

  • 例）
    肺がんになる確率に関する Odds 比を計算してみる
    第 1 群：肺がん患者100人を調査　 ⇒ 80人が喫煙者、20人が非喫煙者 (肺がんの確率 p = 0.8)
    第 2 群：健康な100人を調査　 ⇒ 20人が喫煙者、80人が非喫煙者 (肺がんの確率 (1 - p) = q = 0.2)

  • 下の式に p = 0.8, q = 0.2 を代入して Odds Ratio を求めると
    \[OddsRatio=\frac{Odds[p]}{Odds[q]}\ =\frac{p/(1-p)}{q/(1-q)}\ = \frac{p(1-q)}{(1-p)q} = \frac{0.8(1-0.2)}{(1-0.8)0.2} = \frac{4.0}{0.25} = 16\]

  • 下の表からも、p = 0.8 の時の Odds Ratio = 4、p = 0.2 の時の Odds Ratio = 0.25 が確認できる

  • Odds Ratio = 1 ⇒ 事象の起こりやすさが両群で同じ
  • Odds Ratio > 1 ⇒ 事象が第 1 群で起こりやすい
  • Odds Ratio < 1 ⇒ 事象が第 2 群で起こりやすい

例）Odds Ratio = 16 が意味するのは「第 2 群と比べると、第 1 群では喫煙者が肺がんである可能性が非喫煙者の16倍」

  • model_1 の分析結果は次のとおり
tidy(model_1)
# A tibble: 3 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)   -2.65     0.194     -13.7  9.90e-43
2 expm           0.170    0.0235      7.23 4.86e-13
3 previous       0.342    0.0355      9.62 6.38e-22
  • 表示される回帰式の係数 estimate は 「オッズの対数」 = Log Odds
  • Log Odds は解釈しにくい → exp()関数を使ってオッズ比に変換する
tidy(model_1, 
  exponentiate = TRUE, 
  conf.int = TRUE)
# A tibble: 3 × 7
  term        estimate std.error statistic  p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 (Intercept)   0.0705    0.194     -13.7  9.90e-43   0.0476     0.102
2 expm          1.19      0.0235      7.23 4.86e-13   1.13       1.24 
3 previous      1.41      0.0355      9.62 6.38e-22   1.32       1.51
  • expm のオッズ比は 1.19
  • Odds Ratio > 1 ⇒ 事象が第 1 群で起こりやすい（つまり当選しがち）
    →　選挙費用 (expm) が100万円増えるにつれて、当選するオッズが 約1.9倍になる
  • previous のオッズ比は 1.4
    →　当選回数 (previous) が1回増えるにつれて、当選するオッズが 約1.4倍になる
  • いずれも統計的に有意 (p.value が0.05以下)

注：もし係数が 1 以下なら、「当選するオッズが下がる」

3.6 ロジスティック回帰モデルの当てはまり具合を評価する

  • ロジスティック回帰分析の当てはまり具合を評価する方法は次の 3 つ:　　

(1) 予測の的中率
(2) ROC 曲線と AUC
(3) 赤池情報量基準: AID

  • 予測の的中率

  • ここで説明しようとしているのは次のモデル

  • 応答変数は二値変数で、各候補者が「当選」か「落選」か (wlsmd = 1 or 0)

model_1 <- glm(wlsmd ~ previous + expm, 
               data = hr21,
               family = binomial(link = "logit"))
  • 上で推定したロジスティック回帰式は、当落をどの程度正確に予測しているかを調べる
  • 「予測値に基づく当落」と「実際の当落」をクロス集計表にする
  • 当選確率が 0.5 以上という予測 = 「当選」と予測したと考える
  • model_1 の予測値を fitted()関数を使って取り出し
    　→「予測値に基づく当落」と「実際の当落」をクロス集計表にする
Pred <- (fitted(model_1) >= 0.5) %>%
  factor(levels = c(FALSE, TRUE),　
         labels = c("落選予測", "当選予測"))

table(Pred, hr21$wlsmd) %>% addmargins()
          
Pred         0   1 Sum
  落選予測 476 120 596
  当選予測  70 169 239
  Sum      546 289 835

【予想値に基づく当落】

・実際に落選した 546 人のうち、落選と予測されたのは 476 人
・残りの 70 人については当選という誤った予測
・実際に当選した 289 人のうち、169 人については当選という正しい予測
・実際に当選した 298 人のうち 120 人については落選という誤った予測
・全体としては、835 人中 645 人 (476 人 + 169 人) については正しい予測
・残りの 190 人については誤った予測
→　従って、このモデルの的中率は、645/835 (約77%)

  • この「77%」という的中率をどう評価すべきか？
  • ロジスティック回帰によって予測の的中率が 0 % から77% に上がったわけではない

【実際の当落 = 的中率】

- ここでは 835 人の候補者中、実際には 289 人が当選
→ 説明変数を何も加えず「全員当選」という予測をすれば
→ 予測の精度は 289/835 (約34%)

予測の的中率 model_1 は「的中率」を 34% から 77% へ 43 ポイント上げた

ロジスティック回帰の予測精度が高いといえるかどうかは、説明変数をいっさい使わなくても得られる「的中率」と比較して評価する

  • ROC 曲線と AUC
  1. 目視によるモデルの当てはまり具合の評価方法
    受信者操作特性（receiver operating characteristic：ROC）曲線

ROC 曲線の描き方
ROC 線の横軸には、偽陽性率（false positive rate：FPR）と呼ばれるものを使う
- 偽陽性とは「本当は陰性なのに誤って陽性と判断されること」　　
- ここでは「本当は落選したのに当選と予測すること」= 偽陽性
- 当落の境界線を当選確率 0.5 に設定すると、偽陽性率は 55/678 （下図参照）

Pred <- (fitted(model_1) >= 0.5) %>%
  factor(levels = c(FALSE, TRUE),　labels = c("落選予測", "当選予測"))

table(Pred, hr21$wlsmd) %>% addmargins()
          
Pred         0   1 Sum
  落選予測 476 120 596
  当選予測  70 169 239
  Sum      546 289 835
  • 偽陽性は誤った判断なので、偽陽性率は小さくなることが望ましい
  • 縦軸には、真陽性率（true positive rate：TPR）を使う
  • 真陽性率は感度（sensitivity）とも呼ばれる
  • 真陽性とは「本当は陽性のときに陽性であると正しく判断されること」
  • ここでは「実際に当選した候補者を当選すると予測すること」= 真陽性　　
  • 当落の境界線を当選確率 0.5 に設定した場合の感度は 169/289（上の表を参照）
  • 感度は正しい判断の確率を表すので、大きいほうが望ましい

- 当てはまりがよいモデルのROC 曲線:
→ 点(0, 0) から点(0,1) の近くに進み、そこから点(1, 1) に向かって進む曲線
- 当てはまりの悪いモデルのROC 曲線:
→ ROC 曲線が 45 度線の近くを通過する　　

  • ここで推定している二つのモデルは次のとおり
model_1 <- glm(wlsmd ~ expm, data = hr21,
               family = binomial(link = "logit"))
model_2 <- glm(wlsmd ~ expm + previous, data = hr21,
               family = binomial(link = "logit"))
  • 交差項を含まない 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_2) のROC 曲線は 45 度線（点線）から点(0, 1) のほうに離れており、モデルの当てはまりがいい

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

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.7825171
auc2 <- performance(pr2, measure = "auc")
auc2@y.values[[1]]   # model_2 のAUC
[1] 0.853502

AUC による当てはまり AUCを見る限り、model_2 の当てはまりの方がよい

赤池情報量基準 (AID)

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

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

4. ロジスティック回帰分析の実例

4.1 出場順番と優勝との関係（M1 決勝）

M1グランプリについての解説はこちらを参照：[M-1グランプリ2022](https://ja.wikipedia.org/wiki/M-1%E3%82%B0%E3%83%A9%E3%83%B3%E3%83%97%E3%83%AA2022)

M1グランプリについての解説はこちらを参照：M-1グランプリ2022

  • この分析で使うデータ
変数名 内容 備考
ID
No 回
Year 開催年
Rank 最終順位
Name コンビ名
Company 所属事務所名
Entry_No エントリー番号
Since コンビ結成年
No_Finals これまでの決勝進出回数 当該年を除く
Catchphrase キャッチフレーズ
Final 最終決戦進出（0：進出せず、1:進出） 従属変数
Order10 決勝における出場順番 独立変数
Order3 最終決戦における出場順番 最終決戦に進出しなかったコンビは欠損値
Score10 決勝における評価
Score3 最終決戦における評価 最終決戦に進出しなかったコンビは欠損値
No_Reviewer 審査委員数

・この分析で使うデータはここからダウンロード： M1_Grand_Pix_2022.csv
・出典：https://github.com/JaehyunSong/M-1_Grand_Pix

データの準備

  • ファイル (M1_Grand_Pix_2022.csv) をダウンロードする
  • Rプロジェクトフォルダ内に data という名称のフォルダを作成
  • data フォルダ内に M1_Grand_Pix_2022.csvを置く
  • ファイル (M1_Grand_Pix_2022.csv) を読み込み df_m1 と名前を付ける
df_m1 <- read_csv("data/M1_Grand_Pix_2022.csv")
  • 分析しやすいように Zombie と Winner という二つの変数を作る
df_m1 <- df_m1 %>%
  mutate(Zombie = if_else(Catchphrase == "（敗者復活）", 1, 0), # 敗者復活したかどうか（0:しない、1:した）　　
         Winner = if_else(Rank == 1, 1, 0)) # 最終決戦で優勝したどうか（0: 2 位もしくは 3 位、1: 優勝）
  • df_m1 の中身ははこんな感じ
DT::datatable(df_m1)
  • ロジステック回帰分析で使う変数一覧
変数名 内容 備考
Final 最終決戦でトップ 3 に残れたどうか（0: 残れず、1: 残れた） 従属変数
Order10 決勝における出場順番 (1 〜 10） 独立変数
Since コンビ結成年
No_Finals これまでの決勝進出回数 当該年を除く
Zombie 敗者復活したコンビなら 1、それ以外は 0
  • データの記述統計を表示
stargazer(as.data.frame(df_m1),
          type = "html")
Statistic N Mean St. Dev. Min Max
No 169 85.000 48.930 1 169
Year 169 2,011.509 7.071 2,001 2,022
Rank 169 5.189 2.719 1 10
Entry_No 169 2,610.000 1,477.452 6 4,948
Since 169 2,003.609 6.336 1,992 2,020
No_Finals 169 0.905 1.346 0 8
Final 169 0.314 0.465 0 1
Order10 169 5.207 2.725 1 10
Order3 53 1.981 0.820 1 3
Score10 169 628.941 72.934 436 834
Score3 53 2.377 2.177 0 7
No_Reviewer 169 7.000 0.655 5 9
Zombie 169 0.101 0.302 0 1
Winner 169 0.107 0.309 0 1

ロジステック回帰分析の実行

model_1 <- glm(Final ~ Order10 + Since + No_Finals + Zombie, 
            data = df_m1, 
            family = binomial("logit"))
  • 分析結果を表示する　　
stargazer(model_1,
          type = "html")
Dependent variable:
Final
Order10 0.243***
(0.072)
Since -0.030
(0.029)
No_Finals 0.345**
(0.136)
Zombie 0.126
(0.559)
Constant 58.151
(57.454)
Observations 169
Log Likelihood -94.747
Akaike Inf. Crit. 199.494
Note: p<0.1; p<0.05; p<0.01 
plot_summs(model_1)

最終決戦に選ばれる確率

  • Order10 の係数 0.243 を見ると、「決勝における出場順番 (Order10)」と「最終決戦進出 (Final)」に正の関係があることがわかる
  • しかし、この係数係数 0.243 は logit = Log Odds を表している
    → このままでは解釈が難しい
    → 出場順番が 1 から 10 まで変化する（= 出場が遅くなる）場合それぞれにおける最終決戦に選ばれる確率を計算してみる
model_1 %>% 
  prediction(at = list(Order10 = 1:10)) %>%
  summary() %>%
  rename("Order" = "at(Order10)") 
 Order Prediction      SE     z         p   lower  upper
     1     0.1400 0.04326 3.236 1.213e-03 0.05519 0.2248
     2     0.1700 0.04250 3.999 6.352e-05 0.08668 0.2533
     3     0.2048 0.04042 5.066 4.057e-07 0.12557 0.2840
     4     0.2446 0.03763 6.499 8.110e-11 0.17081 0.3183
     5     0.2892 0.03576 8.087 6.096e-16 0.21913 0.3593
     6     0.3384 0.03738 9.052 1.399e-19 0.26511 0.4116
     7     0.3913 0.04417 8.860 7.991e-19 0.30477 0.4779
     8     0.4471 0.05524 8.094 5.763e-16 0.33884 0.5554
     9     0.5044 0.06834 7.380 1.581e-13 0.37044 0.6383
    10     0.5618 0.08133 6.907 4.941e-12 0.40236 0.7212
  • この結果を可視化すると次のようになる
fig_1 <- cplot(model_1, 
                  x = "Order10", 
                  what = "prediction") %>%
  as_data_frame() %>%
  ggplot(aes(x = xvals, y = yvals, ymin = lower, ymax = upper)) +
  geom_ribbon(fill = "gray") +
  geom_line() +
labs(x = "出場順番", 
     y = "最終決戦に出場する確率の予測値", 
     title = "Fig.1:出場順と最終決戦出場予測値") +
  theme_bw(base_family = "HiraKakuProN-W3") +
  scale_x_continuous(breaks = seq(0, 10, by = 1),
                     labels = seq(0, 10, by = 1))

  • 「出場順番」が変化するにつれて「出場順番」が「最終決戦に選ばれる確率」に与える影響力（＝傾き＝限界効果）は大きくなってる
  • 出場順番が 1 から 10 まで変化するにつれて（= 出場が遅くなるほど）
    →　最終決戦に選ばれる確率が大きくなっていることがわかる
  • 出場順番が 1 番の場合　→　Order が 1 の Prediction の値は 0.1400
    →　最終決戦に選ばれる確率は約 0.14 (= 14%)
  • 95% 信頼区間を見るとその確率のブレ幅は 5.5% から 22% の間
  • 出場順番が 10 番の場合　→　Order が 10 の Prediction の値は 0.5618
    →　最終決戦に選ばれる確率は約 56%
  • その確率のブレ幅は 40% から 72% の間
ロジステック回帰分析では、説明変数（この場合だと Order10）の値ごとに傾きが異なる

→　Order10 の値ごとの傾き（＝限界効果）を調べる必要がある

限界効果

  • 「出場順番」が「最終決戦に選ばれる確率」に与える影響力（＝傾き＝限界効果）は「出場順番」の値によって異なる
    →　「出場順番」ごとの「限界効果」をチェック
  • 横軸が「出場順番」(Order10)
  • 縦軸が「限界効果」(Marginal effect of Order10)
fig_2 <- cplot(model_1, 
                  x = "Order10",   # x軸に据える変数
                  dx = "Order10",  # 説明変数
                  what = "effect") %>%
  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, 0.1) +
labs(x = "出場順番", 
     y = "出場順番の平均限界効果", 
     title = "Fig.2: 出場順番の平均限界効果") +
  theme_bw(base_family = "HiraKakuProN-W3") +
  scale_x_continuous(breaks = seq(0, 10, by = 1),
                     labels = seq(0, 10, by = 1))

  • 「最終決戦に選ばれる確率」と「限界効果」の結果を一つのグラフにまとめて表示する
fig_1 + fig_2

左側の fig_1

  • 出場順番が 1 から 10 まで変化するにつれて（= 出場が遅くなるほど）
    →　最終決戦に選ばれる確率が大きくなっている

右側の fig_2

  • 95%信頼区間が縦軸 = 0 を交わっていない
    →　「出場順番」(Order10) が 1 から 10 のいずれにおいても、統計的に有意

棒グラフと線グラフで可視化

model_1 <- glm(Final ~ Order10 + Since + No_Finals + Zombie, 
            data = df_m1, 
            family = binomial("logit"))
  • df_m1 を使って、出場順番 (Order10) ごとに最終決戦進出に進出する確率を計算する
  • 結果を Mean_final というデータフレームに格納する
bar_df_m1 <- df_m1 %>%
  group_by(Order10) %>%
  summarise(Mean_Final = mean(Final))
bar_df_m1
# A tibble: 10 × 2
   Order10 Mean_Final
     <dbl>      <dbl>
 1       1      0.111
 2       2      0.167
 3       3      0.111
 4       4      0.389
 5       5      0.444
 6       6      0.389
 7       7      0.222
 8       8      0.278
 9       9      0.611
10      10      0.571
  • 出場順番 (Order10) ごとに最終決戦進出に進出する確率を棒グラフで可視化する
  • model_1 の推定結果を線グラフで可視化する
model_1 %>% 
  prediction(at = list(Order10 = 1:10)) %>%
  summary() %>%
  rename("Order" = "at(Order10)") %>%
  ggplot() +
  geom_bar(data = bar_df_m1,
           aes(x = Order10, 
               y = Mean_Final), 
           stat = "Identity", 
           fill = "gray70") +
  geom_pointrange(aes(x = Order, 
                      y = Prediction, 
                      ymin = lower,  # 進出確率の最小値
                      ymax = upper), # 進出確率の最大値
                  size = 1.2) +
  geom_line(aes(x = Order, 
                y = Prediction), 
            size = 1.2) +
  scale_x_continuous(breaks = 1:10, labels = 1:10) +
  coord_cartesian(ylim = c(0, 0.75)) +
  labs(x = "出場順番",
       y = "最終決戦へ進出する確率",
       title = "M-1グランプリ (2001-2022)") +
  theme_minimal(base_family = "HiraKakuProN-W3") +
  theme(panel.grid.minor = element_blank(),
        text = element_text(size = 16))

グラフの解釈：

  • 棒グラフは単純にそれぞれの出場順番の出場者が最終決戦に出場した確率を表す
bar_df_m1 <- df_m1 %>%
  group_by(Order10) %>%
  summarise(Mean_Final = mean(Final))
bar_df_m1
# A tibble: 10 × 2
   Order10 Mean_Final
     <dbl>      <dbl>
 1       1      0.111
 2       2      0.167
 3       3      0.111
 4       4      0.389
 5       5      0.444
 6       6      0.389
 7       7      0.222
 8       8      0.278
 9       9      0.611
10      10      0.571
  • 線グラフはロジステック回帰分析 (model_1) による推定結果
model_1 %>% 
  prediction(at = list(Order10 = 1:10)) %>%
  summary() %>%
  rename("Order" = "at(Order10)")
 Order Prediction      SE     z         p   lower  upper
     1     0.1400 0.04326 3.236 1.213e-03 0.05519 0.2248
     2     0.1700 0.04250 3.999 6.352e-05 0.08668 0.2533
     3     0.2048 0.04042 5.066 4.057e-07 0.12557 0.2840
     4     0.2446 0.03763 6.499 8.110e-11 0.17081 0.3183
     5     0.2892 0.03576 8.087 6.096e-16 0.21913 0.3593
     6     0.3384 0.03738 9.052 1.399e-19 0.26511 0.4116
     7     0.3913 0.04417 8.860 7.991e-19 0.30477 0.4779
     8     0.4471 0.05524 8.094 5.763e-16 0.33884 0.5554
     9     0.5044 0.06834 7.380 1.581e-13 0.37044 0.6383
    10     0.5618 0.08133 6.907 4.941e-12 0.40236 0.7212

4.2 出場順番と優勝との関係（M1 最終決戦）

  • 「決勝」において 10 組から 3 組だけが「最終決戦」に選ばれる

  • ここでは「最終決戦」に選ばれた 3 組の中における出場順番と優勝との関係を分析する

  • ロジステック回帰分析で使う変数一覧

変数名 内容 備考
Winner 最終決戦で優勝したかどうか（0：2 位もしくは 3 位、1 = 優勝） 従属変数
Order3 最終決戦における出場順番 独立変数
Since コンビ結成年
No_Finals これまでの決勝進出回数 当該年を除く
Zombie 敗者復活したコンビなら 1、それ以外は 0

ロジステック回帰分析の実行

model_2 <- glm(Winner ~ Order3 + Since + No_Finals + Zombie, 
            data = df_m1, 
            family = binomial("logit"))
  • 分析結果を表示する　　
stargazer(model_2,
          type = "html")
Dependent variable:
Winner
Order3 0.775*
(0.411)
Since -0.090
(0.057)
No_Finals -0.288
(0.221)
Zombie -1.296
(0.981)
Constant 179.137
(113.774)
Observations 53
Log Likelihood -30.299
Akaike Inf. Crit. 70.598
Note: p<0.1; p<0.05; p<0.01 
plot_summs(model_2)

優勝する確率

  • Order3 の係数 0.775 を見ると、「決勝における出場順番 (Order3)」と「最終決戦で優勝したかどうか (Final)」に正の関係があることがわかる
  • しかし、この係数係数 0.775 は logit = Log Odds を表している
    → このままでは解釈が難しい
    → 出場順番が 1 から 3 まで変化する（= 出場が遅くなる）場合それぞれにおける優勝する確率を計算してみる
model_2 %>% 
  prediction(at = list(Order3 = 1:3)) %>%
  summary() %>%
  rename("Order" = "at(Order3)") 
 Order Prediction      SE     z         p   lower  upper
     1     0.2011 0.07818 2.573 1.009e-02 0.04791 0.3544
     2     0.3385 0.06355 5.326 1.003e-07 0.21393 0.4630
     3     0.5066 0.10855 4.666 3.066e-06 0.29379 0.7193
  • この結果を可視化すると次のようになる
fig_3 <- cplot(model_2, 
                  x = "Order3", 
                  what = "prediction") %>%
  as_data_frame() %>%
  ggplot(aes(x = xvals, y = yvals, ymin = lower, ymax = upper)) +
  geom_ribbon(fill = "gray") +
  geom_line() +
labs(x = "出場順番", 
     y = "優勝する確率の予測値", 
     title = "Fig.3:出場順と優勝予測値") +
  theme_bw(base_family = "HiraKakuProN-W3") +
  scale_x_continuous(breaks = seq(0, 3, by = 1),
                     labels = seq(0, 3, by = 1))

  • 「出場順番」が変化するにつれて「出場順番」が「優勝する確率」に与える影響力（＝傾き＝限界効果）は大きくなってる
  • 出場順番が 1 から 3 まで変化するにつれて（= 出場が遅くなるほど）
    →　優勝する確率が大きくなっていることがわかる
  • 出場順番が 1 番の場合　→　Order が 1 の Prediction の値は 0.2011
    →　最終決戦に選ばれる確率は約 0.2011 (= 20%)
  • 95% 信頼区間を見るとその確率のブレ幅は 4.8% から 35% の間
  • 出場順番が 3 番の場合　→　Order が 3 の Prediction の値は 0.5066
    →　最終決戦に選ばれる確率は約 51%
  • その確率のブレ幅は 29% から 72% の間
ロジステック回帰分析では、説明変数（この場合だと Order3）の値ごとに傾きが異なる

→　Order3 の値ごとの傾き（＝限界効果）を調べる必要がある

限界効果

  • 「出場順番」が「最終決戦に選ばれる確率」に与える影響力（＝傾き＝限界効果）は「出場順番」の値によって異なる
    →　「出場順番」ごとの「限界効果」をチェック
  • 横軸が「出場順番」(Order3)
  • 縦軸が「限界効果」(Marginal effect of Order3)
fig_4 <- cplot(model_2, 
                  x = "Order3",    # x軸に据える変数
                  dx = "Order3",   # 説明変数
                  what = "effect") %>%
  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.1, 0.36) +
labs(x = "出場順番", 
     y = "出場順番の平均限界効果", 
     title = "Fig.4: 出場順番の平均限界効果") +
  theme_bw(base_family = "HiraKakuProN-W3") +
  scale_x_continuous(breaks = seq(0, 3, by = 1),
                     labels = seq(0, 3, by = 1))

  • 「優勝する確率」と「限界効果」の結果を一つのグラフにまとめて表示する
fig_3 + fig_4

左側の fig_3

  • 出場順番が 1 から 3 まで変化するにつれて（= 出場が遅くなるほど）
    →　優勝する確率が大きくなっている

右側の fig_2

  • 95%信頼区間が縦軸 = 0 と一部交わっている（出場順番が 2 と 3 の場合）
    →　「出場順番」(Order3) が 1 の場合のみ統計的に有意

棒グラフと線グラフで可視化

model_2 <- glm(Winner ~ Order3 + Since + No_Finals + Zombie, 
            data = df_m1, 
            family = binomial("logit"))
  • df_m1 を使って、出場順番 (Order3) ごとに最終決戦進出に進出する確率を計算する
  • 結果を Mean_final というデータフレームに格納する
bar_df_m1_w <- df_m1 %>%
  group_by(Order3) %>%
  drop_na() |> 
  summarise(Mean_Winner = mean(Winner))
bar_df_m1_w
# A tibble: 3 × 2
  Order3 Mean_Winner
   <dbl>       <dbl>
1      1       0.222
2      2       0.333
3      3       0.471
  • 出場順番 (Order3) ごとに優勝する確率を棒グラフで可視化する
  • model_2 の推定結果を線グラフで可視化する
model_2 %>% 
  prediction(at = list(Order3 = 1:3)) %>%
  summary() %>%
  rename("Order" = "at(Order3)") %>%
  ggplot() +
  geom_bar(data = bar_df_m1_w,
           aes(x = Order3, 
               y = Mean_Winner), 
           stat = "Identity", 
           fill = "gray70") +
  geom_pointrange(aes(x = Order, 
                      y = Prediction, 
                      ymin = lower,  # 進出確率の最小値
                      ymax = upper), # 進出確率の最大値
                  size = 1.2) +
  geom_line(aes(x = Order, 
                y = Prediction), 
            size = 1.2) +
  scale_x_continuous(breaks = 1:3, labels = 1:3) +
  coord_cartesian(ylim = c(0, 0.75)) +
  labs(x = "出場順番",
       y = "優勝する確率",
       title = "M-1グランプリ (2001-2022)：最終決戦") +
  theme_minimal(base_family = "HiraKakuProN-W3") +
  theme(panel.grid.minor = element_blank(),
        text = element_text(size = 16))

グラフの解釈：

  • 棒グラフは単純にそれぞれの出場順番の出場者が最終決戦に出場した確率を表す
bar_df_m1_w <- df_m1 %>%
  group_by(Order3) %>%
  drop_na() |> 
  summarise(Mean_Winner = mean(Winner))
bar_df_m1_w
# A tibble: 3 × 2
  Order3 Mean_Winner
   <dbl>       <dbl>
1      1       0.222
2      2       0.333
3      3       0.471
  • 線グラフはロジステック回帰分析 (model_2) による推定結果
model_2 %>% 
  prediction(at = list(Order3 = 1:3)) %>%
  summary() %>%
  rename("Order" = "at(Order3)")
 Order Prediction      SE     z         p   lower  upper
     1     0.2011 0.07818 2.573 1.009e-02 0.04791 0.3544
     2     0.3385 0.06355 5.326 1.003e-07 0.21393 0.4630
     3     0.5066 0.10855 4.666 3.066e-06 0.29379 0.7193

結論

・2001年から2022年までのM1グランプリにおける「出場順番」と「優勝」の関係は次の世にまとめることができる

結論・M1決勝では「出場順番」が遅くなるにつれて最終決戦に選ばれる確率大きくなる
・1 番目に出場する場合、約14%（統計的に有意）
・5 番目に出場する場合、約29%（統計的に有意）
・・・
・10 番目に出場する場合、約56%（統計的に有意）

・M1最終決戦では「出場順番」が遅くなるにつれて優勝する確率大きくなる「傾向にある」
・1 番目に出場する場合、約20%（統計的に有意）
・2 番目に出場する場合、約34%
・3 番目に出場する場合、約51%

5. Excercise

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

Q1.当落（wlsmd）を縦軸に、選挙費用（expm）を横軸にとった散布図（むりやり通常の回帰直線を当てはめた直線）を描きなさい。

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

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

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

Q5. model_1 と model_2 の分析結果を stargazer()関数を使って同一の表に示しなさい。

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

Q7. margins::cplot()関数を使って model_2 の推定結果を表示し、実質的有意性と統計的有意性に関して結論を述べなさい。分析結果は patchworkパッケージを使って表示させること。

Q8. model_2 の予測の的中率を上げるのか下げるのか、解説しなさい。

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

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

参考文献