Rを使った分析の準備
- ここで使うRのパッケージは次のとおり。

library(corrplot)
library(jtools)
library(margins)
library(ROCR)
library(patchwork)
library(prediction)
library(stargazer)
library(tidyverse)

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

  • 線形回帰での応答変数は「得票率」のような連続変数
  • 非線形回帰での応答変数は「当落」のようなカテゴリ変数（ダミー変数）
  • ダミー変数は 2 通りの値 (例えば 0 = 落選、 1 = 当選) をとる
  • このようなダミー変数は「binary変数」とも呼ばれる
分析方法 応答変数 推定方法 確率モデル
回帰分析 連続変数 OLS（最小二乗法） 線形確率モデル
ロジスティック回帰分析 質的変数 MLE（最尤法） 非線形確率モデル
プロビット回帰分析 質的変数 MLE（最尤法） 非線形確率モデル
  • 線形回帰と非線形回帰（ここではロジスティック回帰）の違いを視覚的に示してみる

まとめ線形モデルでは「直線」を使う
→ 説明変数の値にかかわらず、常に影響力（= 傾き）は一定
→ 「説明変数が 1 変化した時に確率がどれだけ変化するか（傾き）」は変わらない

非線形モデルでは「曲線」を使う
→ 説明変数の値しだいで、影響力（= 傾き）が異なる
「説明変数 x が 1 変化した時に確率がどれだけ変化するか（傾き）」 は、その x の値によって異なる」 → 「説明変数の値ごとの傾き」を求める必要あり

2. 線形確率モデルの問題点

  • 応答変数がダミー変数（`0 = 落選、 1 = 当選) を線形確率モデル（回帰分析）で分析すると想定
  • 線形確率モデルは「直感的に」理解しやすい
  • しかし次の 3 つの問題点がある
問題1: 不均一分散になる
問題2: Y の予測値（確率）が 0 と 1 を超えてしまう
問題3: 説明変数の値によって傾きが異なる→係数の推定にはバイアスが発生（深刻な問題）
  • これらの問題点に対するた対処法は次のとおりである　　　
対処法
問題1: 頑強な標準誤差や GLS を使えば正しく推定可能→深刻な問題ではない
問題2: 説明変数が極端でない値をとる限りは、深刻な問題ではない
問題3: 非線形モデル (「プロビットモデル」や「ロジットモデル」 を使って推定する

Logit Model と Probit Model の違い

  • 線形確率モデルの 3 つの問題点を解決するための方法：
    線形確率モデルが想定する関係
    「説明変数 x が 1 増える　→　確率 y が β1 増える」という次の単純比例（線形）の関係

\[y = β_0 + β_1x + ε\]

この関係を次のような非線形確率モデルに変形する

\[y = F(β_0 + β_1x) + ε\]

つまり「説明変数 x がどのような値でも、確率 y が 0 から 1 の間に収まる関係」に、累積分布関数 F を使って変形する

  • この変形には様々な分布を使うことが可能だが、実際に使われるのは次の二つの方法
変換方法 使用する分布関数 特徴 使用分野
Logit Model ロジスティク分布 計算が簡素で解釈が容易 社会- 心理- 医学系
Probit Model 正規分布 計算式が複雑だが理論的に正当 経済学系

  • ここで Probit Model の「理論的に正当」という意味は、私たちの身長など世の中の多くは正規分布で近似できるため、正規分布を使うことにより合理性があるということ

  • 線形モデルにおける OLS と同様、これら非線形モデルについても、データに最も適合するようなパラメータ（係数）を探すことが可能

  • Logit Model と Probit Model の推定結果はほとんど同じ →　どちらを使っても良い

  • ここでは計算が簡素で解釈が容易な Logit Model を使って解説

3. Logit Model の基礎知識

3.1 オッズ (Odds):

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

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

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

  • 例えば、p = 0.01 の場合の 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%」という意味
- オッズの最大値は無限大 (∞)

3.2 オッズ比 (Odds Ratio):

  • ある事象の起こりやすさを 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の特徴 - Odds Ratio = 1 ⇒ 事象の起こりやすさが両群で同じ
- Odds Ratio > 1 ⇒ 事象が第 1 群で起こりやすい
- Odds Ratio < 1 ⇒ 事象が第 2 群で起こりやすい

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

3.3 Log Odds (= logit)

  • 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%」という意味

4. ロジスティック関数

  • 「選挙で当選したか落選したか」を変数で表す
  • 当選なら 1、落選なら 0 という値をとる変数を考える
  • このような変数は、とり得る値が 2 通りしかない = 2 値変数（binary variable)
  • 2 値変数 y が応答変数、説明変数を x とする単回帰分析の回帰直線は次の式で表すことができる

\[\hat{y} =　b_0 + b_1x \]

  • \(b_0\)と\(b_1\)は様々な値をとり得る
  • 応答変数が 2 値変数のとき、通常の線形回帰はおかしな予測をしてしまう
  • 応答変数として二値変数を扱う時 → 2 値変数に適したモデルを考える必要あり
  • 2 値変数の特徴:
  • 0 と1 という二つの結果のいずれかしか起こらない → ベルヌーイ試行

→第 i 番目の \(Y_i（i = 1, 2, …, n）\)は、1 をとる確率（成功確率）が \(π_i\) のベルヌーイ分布に従う

  • \(Y_i\) が 1 をとる確率- - - \(Pr(Y_i=1) = π_i\)
  • \(Y_i\) が 0 をとる確率- - - \(Pr(Y_i=0) = 1 - π_i\)
  • \(Y_i\) の期待値- - - \(E(Y_i) = π_i\)
  • \(Y_i\) が 0 になるか1 になるかを決めるメカニズム

→ 1 つの母数 \(E(Y_i) = π_i\) を使って表すことができる
- \(Y_i\) の値そのものではなく、 \(Y_i\) が 1 になる確率 \(E(Y_i) = π_i\) を予測するモデルを考える
- 確率は 0 以上 1 以下の値をとる
→確率を表すためには 0 以上 1 以下の値しかとらないモデルを考える必要がある

→ロジスティック関数は連続値 x を(0, 1) の範囲の値に変換する関数
→ロジスティック関数はロジットの逆関数 : \(logit^{-1}(x)\)

\[logistic(x) = logit^{-1}(x) = \frac{exp(x)}{1+exp(x)}= \frac{1}{1+exp(-x)}\]

  • 説明変数が 1 つ（\(x_1\)）のときに \(Y_i\) が 1 をとる確率は次のように表すことができる

\[Pr(Y_i = 1) = logistic(b_0 + b_1x_1) = \frac{1}{1+exp(-[b_0 + b_1x_1])}\]

  • \(b_0 = b_1x_1\) が含まれている項は線形関数- - - 線形予測子 (linear predictor)　　

  • ロジスティック回帰分析では、線形予測子の \(b_0\) や \(b_1\) の値を推定する

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

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

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

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

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

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

  • 仮説 2：当選回数によって、選挙費が当選確率に与える影響力は異なる
  • 従って、否定されるために設定する帰無仮説は次のようになる

帰無仮説 - 帰無仮説 1：選挙費の額は、小選挙区での当選確率とは関係がない
- 帰無仮説 2：当選回数によって、選挙費が当選確率に与える影響力は変わらない

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

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

Model 1（仮説 1 を検証するためのモデル）

Model 2（仮説 2 を検証するためのモデル）

5.2 説明変数と応答変数の散布図

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

データの準備 (hr-data.Rds)

  • 衆議院議員総選挙の得票データ hr96-17.Rds をダウンロード
download.file(url = "https://git.io/fACk6",
              destfile = "data/hr-data.Rds")

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

df <- read_rds("data/hr-data.Rds")

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

names(df)
 [1] "party"      "party_code" "year"       "ku"         "kun"       
 [6] "name"       "status"     "previous"   "wl"         "voteshare" 
[11] "age"        "nocand"     "rank"       "vote"       "eligible"  
[16] "turnout"    "exp"        "expm"       "vs"         "exppv"     
[21] "smd"        "party_jpn" 
head(df)
# A tibble: 6 × 22
  party party_code  year ku       kun name    status previ…¹ wl    votes…²   age
  <chr>      <int> <int> <chr>  <int> <chr>   <fct>    <int> <fct>   <dbl> <int>
1 LDP            1  1996 aichi      6 ITO, K… 新人         0 落選     26.1    51
2 LDP            1  1996 aichi      7 NIWA, … 新人         0 落選     25.9    49
3 LDP            1  1996 aichi      9 YOSHIK… 新人         0 落選     24.8    73
4 LDP            1  1996 aichi     10 MORI, … 新人         0 落選     30.2    50
5 LDP            1  1996 aomori     4 TSUSHI… 新人         0 落選     36.6    42
6 LDP            1  1996 chiba      3 MATSUN… 新人         0 落選     34.5    34
# … with 11 more variables: nocand <int>, rank <int>, vote <int>,
#   eligible <int>, turnout <dbl>, exp <int>, expm <dbl>, vs <dbl>,
#   exppv <dbl>, smd <fct>, party_jpn <chr>, and abbreviated variable names
#   ¹​previous, ²​voteshare
  • select() 関数を使って year, smd, previous, expm という 4 つの変数だけを選ぶ
  • filter() 関数を使って 2005年衆院選だけのデータを残す
  • wl の値を（落選=0, 当選=1) に変換した wlsmd という変数を新たにつくる
df1 <- df %>%
  select(year, smd, previous, expm) %>%
  filter(year == 2005) %>% 
  mutate(wlsmd = ifelse(smd == "当選", 1, 0))

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

df1
# A tibble: 989 × 5
    year smd   previous  expm wlsmd
   <int> <fct>    <int> <dbl> <dbl>
 1  2005 落選         0  5.29     0
 2  2005 落選         0 10.9      0
 3  2005 落選         0 10.4      0
 4  2005 落選         0 10.7      0
 5  2005 落選         0 11.8      0
 6  2005 落選         0 10.3      0
 7  2005 落選         0  8.22     0
 8  2005 落選         0 21.4      0
 9  2005 落選         0  7.95     0
10  2005 落選         0  5.52     0
# … with 979 more rows

データ：

変数名 詳細
year 衆院選が行われた年
wlsmd 小選挙区での当落ダミー(当選 = 1, 落選 = 0)
previous 当選回数
expm 候補者の選挙費用（百万円）
smd 小選挙区での当落結果（「当選」「落選」）
  • データフレーム df1 の記述統計を表示させる
summary(df1)
      year        smd         previous           expm              wlsmd       
 Min.   :2005   落選:689   Min.   : 0.000   Min.   : 0.06271   Min.   :0.0000  
 1st Qu.:2005   当選:300   1st Qu.: 0.000   1st Qu.: 2.91743   1st Qu.:0.0000  
 Median :2005              Median : 1.000   Median : 7.69602   Median :0.0000  
 Mean   :2005              Mean   : 1.975   Mean   : 8.14224   Mean   :0.3033  
 3rd Qu.:2005              3rd Qu.: 3.000   3rd Qu.:11.82280   3rd Qu.:1.0000  
 Max.   :2005              Max.   :16.000   Max.   :24.64971   Max.   :1.0000  
                                            NA's   :4
  • データフレーム df1 の expm に欠損値 (missing data = NA's) が4個あることがわかる
  • na.omit() を使って欠測のない観測だけを残す
df1 <- na.omit(df1)
  • データフレーム df1 の記述統計を表示して確認する
summary(df1)
      year        smd         previous          expm              wlsmd       
 Min.   :2005   落選:687   Min.   : 0.00   Min.   : 0.06271   Min.   :0.0000  
 1st Qu.:2005   当選:298   1st Qu.: 0.00   1st Qu.: 2.91743   1st Qu.:0.0000  
 Median :2005              Median : 1.00   Median : 7.69602   Median :0.0000  
 Mean   :2005              Mean   : 1.98   Mean   : 8.14224   Mean   :0.3025  
 3rd Qu.:2005              3rd Qu.: 3.00   3rd Qu.:11.82280   3rd Qu.:1.0000  
 Max.   :2005              Max.   :16.00   Max.   :24.64971   Max.   :1.0000

  • 4個の NA's が消えていることがわかる　　

  • 当選回数と当落の散布図

  • Macユーザは次の2行を ggplot2 のテーマを設定する

theme_set(theme_gray(base_size = 10, 
                     base_family = "HiraginoSans-W3"))
figure_1 <- ggplot(df1, aes(x = previous, y = wlsmd)) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(x = "当選回数", y = "選挙での当落", title = "回帰直線（当選回数ー選挙での当落）") +
  geom_hline(yintercept = c(0, 1), color = "gray") +
  geom_jitter(width = 0.05, height = 0.05)
　　　　　　#jitterで重複したデータを散らす

print(figure_1)

  • 選挙費用と当落の散布図
figure_2 <- ggplot(df1, aes(x = expm, y = wlsmd)) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(x = "選挙費用（100万円）", y = "選挙での当落", title = "回帰直線（選挙費用ー選挙での当落）") +
  geom_hline(yintercept = c(0, 1), color = "gray") +
  geom_jitter(width = 0.05, height = 0.05)

print(figure_2)

  • 変数間の相関係数を調べる

  • cor() をデータフレームに適用すると、デー タフレーム内にある変数のすべてのペアについて、相関係数を示してくれる

  • smd と year を除いて相関係数を計算すると

df1  %>%
  select(-smd, -year) %>%
  cor()
          previous      expm     wlsmd
previous 1.0000000 0.5549922 0.6272670
expm     0.5549922 1.0000000 0.5348002
wlsmd    0.6272670 0.5348002 1.0000000

結果
- 小選挙区の当落(wlsmd)と当選回数(previous)には中程度の正の相関（0.63）
- 小選挙区の当落(wlsmd)と選挙費用(expm)には中程度の強さの正の相関（0.53）
- 当選回数(previous)と選挙費用(expm)には中程度の正の相関（0.55）

  • corrplotパッケージを使って、結果を可視化して示すこともできる
library(corrplot)
  • df1 から smd と year を削除し、vis.cor というデータフレームを作る
vis.cor <- df1 %>% 
  select(-smd, -year)
  • previous, expm, wlsmd 三つの変数の相関係数を可視化する
correlations <- cor(vis.cor[,1:3])
corrplot(correlations, method="circle")

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

  • 当選確率を予測するため、ロジスティック回帰式を推定する

\[Pr(wlsmd_i = 1) = logistic(b_0 + b_1[previous]_i + b_2[expm]_i)\]

\[= \frac{1}{1+exp(-[b_0 + b_1[previous]_i + b_2[expm]_i])}\]

  • ロジスティック回帰のリンク関数を logit と指定
    →　表示される回帰式の係数は logit (Log Odds) 　　
model_1 <- glm(wlsmd ~ previous + expm, data = df1,
               family = binomial(link = "logit"))

summary(model_1)

Call:
glm(formula = wlsmd ~ previous + expm, family = binomial(link = "logit"), 
    data = df1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0941  -0.5083  -0.2706   0.3388   2.2611  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.66819    0.23864 -15.371  < 2e-16 ***
previous     0.56727    0.05002  11.340  < 2e-16 ***
expm         0.16240    0.02075   7.827 5.01e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1207.61  on 984  degrees of freedom
Residual deviance:  707.83  on 982  degrees of freedom
AIC: 713.83

Number of Fisher Scoring iterations: 5
  • 候補者が当選する確率Pr(wlsmdi= 1) の予測値は次の式で表すことができる

\[\hat{pi} = \frac{1}{1+exp(-[-3.67 + 0.57previous_i + 0.16expm_i])}\]

  • z値・・・Wald statistic

  • 帰無仮説:「推定値は平均値 0、標準偏差 1 の正規分布である」

  • ロジスティック回帰式の係数 (coefficients)は推定値（Estimate) = logit (Log Odds)

  • 当選回数（previous）と選挙費用 (expm) はどちらも選挙の当落に正の影響を与えている

  • しかし、これらの推定値から得られる知見は次のとおり

  • 回帰式の expm の係数: 0.16240

  • 回帰式の係数は logit (Log Odds)
    → Odds と当選確率を計算できる
    → 当選回数が平均値の候補者が選挙費用を 1 円使った時の Odds は

exp(0.16240)
[1] 1.176331
  • 当選回数が平均値の候補者が選挙費用を 100万円使った時、当選確率は
1/(1+exp(-0.16240))
[1] 0.540511

  • 約54%である

  • 回帰式の previous の係数: 0.56727

→ 当選回数が 1回多い候補者の Odds は次の式で計算できる（選挙費用を平均値に固定）

exp(0.56727)
[1] 1.763446
  • 選挙費用が平均値の候補者の当選回数が一回多ければ、当選確率は
1/(1+exp(-0.56727))
[1] 0.638133

  • 約64%である

  • 回帰式の Intercept の係数: -3.44865

→ 選挙費用が0円、当選回数が 0 回の候補者の Odds は次の式で計算できる

exp(-3.66819)
[1] 0.02552262

→ 選挙費用が0円、当選回数が 0 回の候補者の当選確率は

1/(1+exp(-(-3.66819)))
[1] 0.02488743
  • 約2.5%なので、当選する可能性は低いことがわかる
  • この当選確率は predict()関数を使っても計算できる
predict(model_1, type = "response",
        newdata = data_frame(previous = 0, expm = 0))
         1 
0.02488734

  • 以上の結果は、仮説を検証する上であまり有益とはいえない

  • 後の節ではこの結果を可視化して、回帰係数の有意性を検定する

5.4 モデルの「当てはまり」

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

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

  • 予測の的中率

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

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

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

table(Pred, df1$smd) %>% addmargins()
          
Pred       落選 当選 Sum
  落選予測  632  100 732
  当選予測   55  198 253
  Sum       687  298 985

【予想値に基づく当落】

- 実際に落選した 687 人のうち、落選と予測されたのは 632 人
- 残りの 55 人については当選という誤った予測
- 実際に当選した 298 人のうち、198 人については当選という正しい予測
- 実際に当選した 298 人のうち 100 人については落選という誤った予測
- 全体としては、985 人中 830 人 (632 人 + 198 人) については正しい予測
- 残りの 155 人については誤った予測
→　従って、このモデルの的中率は、830/985 (約84%)

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

【実際の当落 = 的中率】

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

予測の的中率 ✔ model_1 は「的中率」を 30% から 84% へ 54 ポイント上げた

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

  • 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, df1$smd) %>% addmargins()
          
Pred       落選 当選 Sum
  落選予測  632  100 732
  当選予測   55  198 253
  Sum       687  298 985
  • 偽陽性は誤った判断なので、偽陽性率は小さくなることが望ましい
  • 縦軸には、真陽性率（true positive rate：TPR）を使う
  • 真陽性率は感度（sensitivity）とも呼ばれる
  • 真陽性とは「本当は陽性のときに陽性であると正しく判断されること」
  • ここでは「実際に当選した候補者を当選すると予測すること」= 真陽性　　
  • 当落の境界線を当選確率 0.5 に設定した場合の感度は 198/298（上の表を参照）
  • 感度は正しい判断の確率を表すので、大きいほうが望ましい

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

  • ここで推定している二つのモデルは次のとおり
model_1 <- glm(wlsmd ~ previous + expm, data = df1,
               family = binomial(link = "logit"))
model_2 <- glm(wlsmd ~ expm*previous, data = df1,
               family = binomial(link = "logit"))
  • 交差項を含まない model_1 と交差項を含む model_2 のどちらの当てはまりがよいかを、ROC曲線を描いて考える
  • ROC曲線を描くために ROCR パッケージを使う
library(ROCR)
  • prediction() という名前の関数は ROCR パッケージだけでなく margins パッケージにもある
    →どちらの関数を使うか明記する　　
pi1 <- predict(model_1, type = "response")
pi2 <- predict(model_2, type = "response")
pr1 <- ROCR::prediction(pi1, labels = df1$smd == "当選")
pr2 <- ROCR::prediction(pi2, labels = df1$smd == "当選")
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曲線による当てはまり ✔ 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.9145834
auc2 <- performance(pr2, measure = "auc")
auc2@y.values[[1]]   # model_2 のAUC
[1] 0.9108907

AUC による当てはまり ✔ AUCがほとんど差はないが、わずかに model_1 の当てはまりの方がよい

赤池情報量基準 (AID)

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

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

  • 当選回数と当選確率のプロット

  • 横軸 = 当選回数（previous）、縦軸 = 当選確率とした散布図に、ロジスティック回帰分析で推定された曲線を上書きする

  • このとき、横軸に使わないもう一つの説明変数（つまり選挙費用: expm）は、平均値に固定する

  • 手順は次のとおり:
    data_frame()関数を使って予測に利用するデータフレーム (pred_prev) を作る

pred_prev <- data_frame(previous = min(df1$previous):max(df1$previous),
                      expm = mean(df1$expm))
  • predict()関数を使って予測値を計算し、作ったデータフレーム (pred_prev) に加える
  • type = "response" と指定することで、確率（ここでは、当選確率）が計算される
pred_prev$fit <- predict(model_1, 
                         type = "response", 
                         newdata = pred_prev)
  • df1 で散布図（丸い点［pch = 16］）を描いた上に、上で作ったデータフレーム (pred_prev) で当選確率の予測値（正方形の点［pch = 18］）を上書きする
  • 当選回数は連続変数でなく離散変数
    →各当選回数に対して計算された予測値を線で結び、曲線のようなものを描く
figure_3 <- ggplot(df1, aes(x = previous)) +
  geom_hline(yintercept = c(0, 1), color = "gray") +
  geom_jitter(aes(y = wlsmd), pch = 16, size = 1, width = 0.05, height = 0.05) +
  geom_line(data = pred_prev, aes(y = fit), color = "red") +
  geom_point(data = pred_prev, aes(y = fit), pch = 18, size = 4, color = "red") +
  labs(x = "当選回数", y = "当選確率", title = "当選の予測値（当選回数ー当選確率）")

print(figure_3)

  • 当選回数が増えるにつれて、当選確率も大きいことがわかる

  • 選挙費用と当選確率のプロット

  • 横軸 = 選挙費用（expm）、縦軸 = 当選確率とした散布図に、ロジスティック回帰分析で推定された曲線を上書きする

  • 当選回数は平均値に固定

  • 当選回数とは違い、選挙費用は連続変数
    =→ 特定の選挙費用に対する点は描かずに、関数を表す曲線を描く

pred_expm <- data_frame(expm = seq(0, max(df1$expm), length.out = 100),
                        previous = mean(df1$previous))
pred_expm$fit <- predict(model_1, type = "response", newdata = pred_expm)

figure_4 <- ggplot(df1, aes(x = expm)) +
  geom_hline(yintercept = c(0, 1), color = "gray") +
  geom_jitter(aes(y = wlsmd), pch = 16, size = 1, width = 0.05, height = 0.05) +
  geom_line(data = pred_expm, aes(y = fit), color = "red") +
  labs(x = "選挙費用（100万円）", y = "当選確率", title = "当選の予測値（選挙費用ー当選確率）")

print(figure_4)

  • 選挙費用が増えるにつれて、当選確率も大きいことがわかる

  • 「回帰直線」と「当選確率の予測値」の違い

library(patchwork)

「当選回数」と「当選確率」

figure_1 + figure_3

「選挙費用」と「当選確率」

figure_2 + figure_4

- 回帰「直線」→どの選挙費用額でも、選挙費用が選挙当落に与える「傾き」は同じ
- ロジスティック「曲線」→選挙費用額の値によって、選挙費用が選挙当落に与える「傾き」は異なる
- ロジスティック曲線- - - 「当選確率」は 0 以上 1 以下の範囲内

  • ロジスティック回帰分析モデルは予測値が 0 〜 1 に収まるが、非線形モデルなので、推定されたパラメータそのものの解釈は困難

→ 説明変数 1 単位の変化が、応答変数の確率の変化にどの程度影響を及ぼすか（= 限界効果）を計算する必要がある

限界効果 (ME: Marginal Effects)- - - 説明変数の変化の条件付き確率に対する影響

5.5 回帰係数の有意性検定

  • ロジスティック回帰分析の結果は次のとおり
summary(model_1)

Call:
glm(formula = wlsmd ~ previous + expm, family = binomial(link = "logit"), 
    data = df1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0941  -0.5083  -0.2706   0.3388   2.2611  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.66819    0.23864 -15.371  < 2e-16 ***
previous     0.56727    0.05002  11.340  < 2e-16 ***
expm         0.16240    0.02075   7.827 5.01e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1207.61  on 984  degrees of freedom
Residual deviance:  707.83  on 982  degrees of freedom
AIC: 713.83

Number of Fisher Scoring iterations: 5
  • stargazer()を使って、見やすく表示できる
stargazer(model_1, type = "text", single.row=TRUE,
          ci = TRUE) # print 95% CIs

=============================================
                      Dependent variable:    
                  ---------------------------
                             wlsmd           
---------------------------------------------
previous            0.567*** (0.469, 0.665)  
expm                0.162*** (0.122, 0.203)  
Constant          -3.668*** (-4.136, -3.200) 
---------------------------------------------
Observations                  985            
Log Likelihood             -353.915          
Akaike Inf. Crit.           713.830          
=============================================
Note:             *p<0.1; **p<0.05; ***p<0.01
  • この結果からわかること
  • 有意水準を 0.05 に設定すると：
  1. 当選回数（previous）のp値 = 2e-16
    → 当選回数は「影響がない」という帰無仮説を棄却
  2. 選挙費用（expm）の p値 = 2e-16
    → 選挙費用は「影響がない」という帰無仮説を棄却
    → 当選回数と選挙費用は、どちらも当選確率に影響する（統計的に有意）
  • jtoolsパッケージを使うと、統計的有意性を視覚的に確認できる
library("jtools")
plot_summs(model_1)

  • 但し、選挙費用が当選確率に与える「影響の大きさ」は確認できない データ解析において重要なこと: - 統計的有意性 (statistical significance) → 限界効果を計算して図示する - 実質的有意性 (substantial significance) → 限界効果を計算して図示する

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

  • 限界効果と予測確率
    限界効果 (ME) と平均限界効果 (AME)

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

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

  • 限界効果 (Marginal Effects) - - - 説明変数の変化の条件付き確率に対する影響
    → 説明変数の係数だけではなく説明変数の値にも依存 → 平均効果を計算するのが一般的
  • 平均限界効果 (Average Marginal Effects) - - - 説明変数の各値ごとに限界効果を計算し、その平均を取る

→説明変数が応答変数に与える影響を考えるためには、説明変数(x)の値を特定し、特定された値における影響（＝傾き＝平均限界効果）を計算する必要がある
- 上図の例では、xの値が -4, -2, 0, 2 それぞれの時点での傾きが平均限界効果

予測確率

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

例1:
選挙費用 (0 円→100 万円) →予測当選確率:1.86%

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

\[\hat{pi} = \frac{1}{1+exp(-[-3.66819 + 0.56727previous_i + 0.16240expm_i])}\]

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

例2

選挙費用 (100万円→200万円) →予測当選確率: 2.1%

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

→ 当選回数が 3 回の候補者が選挙費用を0 円から 100 万円へ 100 万円変化する (当選確率は約1.8 %ポイント↑）より、100 万円から200 万円へ 100 万円変化する (当選確率は約2.1%ポイント↑）ほうが、応答変数に与える影響が大き い

  • 0 円から100 万円への変化も、100 万円から200 万円への変化も、選挙費用の増分としては同じ
  • しかし、それらの変化が応答変数（= 予測当選確率）に与える影響は異なる
  • 実際に何パーセンテージ- ポイント予測当選確率が上昇するかを、具体的な予測当選確率を計算せずに係数の推定値 (expm = 0.162) から読み取るのは難しい
  • Model 1 のロジスティック回帰分析結果は次のとおり
stargazer(model_1, type = "html")
Dependent variable:
wlsmd
previous 0.567***
(0.050)
expm 0.162***
(0.021)
Constant -3.668***
(0.239)
Observations 985
Log Likelihood -353.915
Akaike Inf. Crit. 713.830
Note: p<0.1; p<0.05; p<0.01
  • 上の例のように、予測当選確率を複数計算し、その差をとるという作業を繰り返すのは面倒
    → 特定の値における限界効果を直接計算するための marginsパッケージを使うのが便利

限界効果：説明変数が（特定の値において）応答変数に与える影響力の強さ

平均限界効果の計算

  • margins()関数を使うと、特定の選挙費用額ごとの平均限界効果 (AME: Average Marginal Effects) を求めることができる
library(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.02707
  • 選挙費用が当選確率に与える影響（平均限界効果）は 0.02707 (約2.7%ポイント) だとわかる

選挙費用の平均限界効果（previous = 3, expm: 100〜1000万円)

  • どの値での平均限界効果を求めるかは、引数 atに変数のリストを渡すことで指定する
  • 複数の値、例えば、選挙費用が 0円 (previous = 0) から 1000万円 (previous = 10) までの平均限界効果を同時に表示してみる
margins(model_1, variables = "expm",
at = list(previous = 3, expm = c(0:28))) #選挙費用を0円から1000万円まで指定  
 at(previous) at(expm)    expm
            3        0 0.01749
            3        1 0.01971
            3        2 0.02207
            3        3 0.02454
            3        4 0.02707
            3        5 0.02959
            3        6 0.03205
            3        7 0.03434
            3        8 0.03640
            3        9 0.03812
            3       10 0.03943
            3       11 0.04027
            3       12 0.04060
            3       13 0.04039
            3       14 0.03966
            3       15 0.03844
            3       16 0.03680
            3       17 0.03481
            3       18 0.03256
            3       19 0.03014
            3       20 0.02762
            3       21 0.02509
            3       22 0.02260
            3       23 0.02021
            3       24 0.01796
            3       25 0.01586
            3       26 0.01394
            3       27 0.01220
            3       28 0.01063
  • ロジスティック回帰分析の係数の限界効果は、説明変数の値によって変化する
    → 各説明変数の係数を見ただけで結果を理解するのは困難
    → 説明変数が応答変数に影響を与える様子をいくつかの条件について図示する
    → 説明変数の値に応じた平均限界効果 (AME: Average Marginal Effects) を図示する必要がある

6. 仮説 1 と仮説 2 の検証結果

6.1 仮説 1 の検証結果（Model 1）

cplot()を使った「平均限界効果」と「予測確率」の表示

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

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

figure_5 <- 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 = "Fig.5:選挙費用と予測当選確率の予測値")

print(figure_5)

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

figure_6 <- cplot(model_1, 
                  x = "expm", 
                  dx = "expm", #調整変数を dx で指定
                  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 = "Fig.6: 選挙費用の平均限界効果（選挙費用ごと）")

print(figure_6)

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

library(patchwork)
figure_5 + figure_6

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

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

実質的有意性 1:
- 選挙費用が大きくなるにつれて、候補者の当選確率は大きくなる → See Fig.5
- 選挙費用が小さいとき、予測当選確率は緩やかに上昇
- 選挙費用が中程度になると曲線の傾きが少しずつ急になる
- 選挙費用が大きくなると再び傾きが小さくなる
統計的有意性: 全ての選挙費用の範囲で統計的に有意 → See Fig.6　　

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

6.2 仮説 2 の検証結果（Model 2）

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

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

ロジスティック回帰式を推定
- 当選確率を予測するため、ロジスティック回帰式を推定する

model_2 <- glm(wlsmd ~ expm*previous, data = df1,
               family = binomial(link = "logit"))

summary(model_2)

Call:
glm(formula = wlsmd ~ expm * previous, family = binomial(link = "logit"), 
    data = df1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.1103  -0.5052  -0.2121   0.4729   2.3289  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -4.316174   0.336556 -12.825  < 2e-16 ***
expm           0.228433   0.029646   7.705 1.31e-14 ***
previous       0.918340   0.117638   7.806 5.88e-15 ***
expm:previous -0.032462   0.009126  -3.557 0.000375 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1207.61  on 984  degrees of freedom
Residual deviance:  696.64  on 981  degrees of freedom
AIC: 704.64

Number of Fisher Scoring iterations: 6
  • 候補者が当選する確率Pr(wlsmdi= 1) の予測値は次の式で表すことができる

\[\hat{pi} = \frac{1}{1+exp(-[-4.32 + 0.92previous_i + 0.23expm_i - 0.03expm:previous_i])}\]

  • 交差項 expm:previous が統計的に有意 (p value = 0.000375)
    → 当選回数によって、選挙費が当選確率に与える影響力は異なる
    → 具体的にどのように異なるのか、予測当選確率と平均限界効果を確かめる

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

figure_7 <- cplot(model_2, 
                  x = "expm", # 説明変数を 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 = "Fig.7:選挙費用と予測当選確率の予測値")

print(figure_7)

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

figure_8 <- cplot(model_2, 
                  dx = "expm", 　　　　# 説明変数を expm に指定
                  x = "expm",  　　　　# 調整変数を 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 = "Fig.8: 選挙費用の平均限界効果（選挙費用ごと）")

print(figure_8)

当選回数（previous）が選挙費用（expm）の限界効果をどのように変化させるか図示

「当選回数」ごとの「選挙費用の限界効果」」(Figure 9)

figure_9 <- cplot(model_2, 
                  dx = "expm",     # 説明変数を expm に指定
                  x = "previous",  # 調整変数を previous に指定
                 what = "effect",  # 限界効果を指定
                 draw = FALSE) %>% 
    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 = "Fig.9: 選挙費用の平均限界効果（当選回数ごと）")

print(figure_9)

「当選回数」ごとの「選挙費用の予測当選確率」」(Figure 10)

df_pre <- expand.grid(expm = seq(0, 25, by = 5),
                      previous = seq(0, 16, by = 2)) %>%
  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) #95%信頼区間を表示
df_pre$upper <- with(pred, fit + 2 * se.fit) #95%信頼区間を表示
df_pre <- df_pre %>%
  mutate(lower = ifelse(lower < 0, 0, lower),
         upper = ifelse(upper > 1, 1, upper))

figure_10 <- ggplot(df_pre, aes(x = expm, y = fit)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = "gray") +
  geom_line() +
  facet_wrap(. ~ previous) +
  labs(x = "選挙費用", y = "予測当選確率の予測値", title = "Fig.10: 当選回数ごとの予測当選確率（説明変数 = 選挙費用）")

print(figure_10)

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

library(patchwork)
figure_7 + figure_10

figure_8 + figure_9

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

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

実質的有意性 1: 
- 当選回数ごとの選挙費用の平均限界効果は、当選回数が約 2〜3 回の時が最大 → See Fig.9 
- 当選回数が 0 回から約 2〜3 回まで増えるにつれ、選挙費用の平均限界効果は大きくなる →　See Fig.9   - 当選回数が約 2〜3 回を超えると、選挙費用の平均限界効果は小さくなる →　See Fig.9

統計的有意性:  当選回数が 0 回から 5 回で統計的に有意 →　See Fig.9  
- 当選回数が 4 回の候補者は、選挙費用を使っても当選確率がそれほど高くならない→ See Fig.10
- 当選回数が 6 回以上の候補者は、選挙費用にかかわらず当選確率が高く、選挙費用の影響がほぼなくなる→ See Fig.10

統計的有意性: 当選回数が 0 回から 5 回までは選挙費用の平均限界効果が有意だが、当選回数が 5 回以上になると有意ではなくなる（つまり、選挙費用が当選確率と関係なくなる）→ See Fig.9

実質的有意性 2:
- 選挙費用が大きくなるにつれて、候補者の当選確率は大きくなる → See Fig.7 & Fig.8
- 選挙費用が小さいとき、予測当選確率は緩やかに上昇→ See Fig.7
- 選挙費用が中程度になると曲線の傾きが少しずつ急になる→ See Fig.7
- 選挙費用が大きくなると再び傾きが小さくなる→ See Fig.7
統計的有意性: 全ての選挙費用の範囲が統計的に有意 See Fig.8　　

実質的有意性 3:
- 選挙費用の平均限界効果は、約1700-1800万円の時が最大 → See Fig.8
- 選挙費用がそれより多くても少なくても効果が小さい
- 選挙費用が増えるにつれて「選挙費用が当落に与える影響力（選挙費用の平均限界効果）」が徐々に大きくなっている
- 選挙費用が約1700-1800万円の時、当選確率に対する影響力が最大
- しかし、その影響力は、選挙費用が約1700-1800万円を越えると小さくなる
統計的有意性: 全ての選挙費用の範囲で統計的に有意 　　

6.3 限界効果の統計的有意性

  • 限界効果は、説明変数 1 単位の増加が応答変数を何単位増加させるかを示す

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

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

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

  • 上図のように、主な説明変数が応答変数に与える影響（限界効果）の符号が、調整変数の値によって変わりうる

  • 調整変数の値が -1.5 くらいより小さければ限界効果の値はマイナス、それより大きければ限界効果の値はプラス

  • 上図に示されているように、交差項を含む回帰分析では、主な説明変数の効果が調整変数の値によって変わるだけでなく、統計的に有意な範囲と有意でない範囲の両方をもつことがあり得る

  • 調整変数の値が -2 より小さい、もしくは -1 より大きければ統計的に有意

  • 調整変数の値が -2 と -1 の間では統計的に有意ではない

  • 調整変数がどの範囲の値をとると限界効果が統計的に有意になるかどうかを、回帰分析の係数の推定値を見ただけで判断することは非常に難しい
    → 限界効果を図示してはじめて、限界効果がどの範囲でどのような符号をもつか、どの範囲で統計的に有意かが明らかになる

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

7.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", 
                  dx = "Order10", #調整変数を dx で指定
                  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

7.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", 
                  dx = "Order3", #調整変数を dx で指定
                  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%

8. Excercise

自民党が政権交代を果たした2012年総選挙データを使って、立候補者が小選挙区で当選したか否か（smd）を応答変数、過去の当選回数（previous）と選挙費用（expm）を説明変数としたロジスティック回帰分析を実行し、以下の各問に答えなさい。衆院選選挙データ（hr-data.Rds）を使うこと。Q3, Q7, Q9-Q14に関しては、統計的有意性と実質的有意性を考慮すること。

Q1. 当落（smd）を縦軸に、選挙費用（expm）を横軸にとった散布図を描きなさい。

Q2. 当落（smd）を縦軸に、過去の当選回数（previous）を横軸にとった散布図を描きなさい。

Q3. 立候補者が小選挙区で当選したか否か（smd）を応答変数、過去の当選回数（previous）と選挙費用（expm）を説明変数としたモデルにmodel_4 という名前を付けてロジスティック回帰分析し、その結果を推定し解釈しなさい。

Q4. model_4 の結果から、候補者が当選する確率Pr(smd= 1) の予測値の式を書きなさい。

Q5. 上の式を可視化し、「当選確率」を縦軸、「過去の当選回数（previous）」を横軸にしたグラフを描きなさい。

Q6. 上の式を可視化し、「当選確率」を縦軸、「選挙費用（expm）」を横軸にしたグラフを描きなさい。

Q7. 立候補者が小選挙区で当選したか否か（smd）を応答変数、過去の当選回数（previous）と選挙費用（expm）を説明変数とした交差項（previous:expm）を含むモデルにmodel_5 という名前を付けてロジスティック回帰分析し、その結果を推定し解釈しなさい。

Q8. ROC 曲線とAUC を使って、model_4 とmodel_5 それぞれのモデルの当てはまりのよさを評価しなさい。

Q9. margins パッケージを使って、選挙費用（expm）が選挙費用（expm）の限界効果をどのように変化させるか図示しなさい。

Q10. margins パッケージを使って、当選回数（previous）が選挙費用（expm）の限界効果をどのように変化させる か図示しなさい。

Q11. margins パッケージを使って、選挙費用（expm）が当選回数（previous）の限界効果をどのように変化させるか図示しなさい。

Q12. margins パッケージを使って、当選回数（previous）が当選回数（previous）の限界効果をどのように変化させるか図示しなさい。

Q13. 横軸を選挙費用（expm）、縦軸を「当選確率の予測値」に指定し、当選回数（previous）が0 回から16 回まで増えることによって当選確率の予測値がどのように変化するか図示しなさい。

Q14. 横軸を選挙費用（expm）、縦軸を「当選確率の予測値」に指定し、当選回数（previous）が増えることによって当選確率の予測値がどのように変化するかということに関して、2005 年（このセッションで演習した結果）と2012年の総選挙を比較しその違いを説明しなさい。

参考文献