31. ロジスティック回帰分析1
Masahiko Asano
2021-04-16
Rを使った分析の準備
- ここで使うRのパッケージは次のとおり。
library(tidyverse)
library(stargazer)
library(patchwork)
library(ROCR)
library(margins)
library(corrplot)
library(jtools)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 x 22
party party_code year ku kun name status previous wl voteshare age
<chr> <int> <int> <chr> <int> <chr> <fct> <int> <fct> <dbl> <int>
1 LDP 1 1996 aichi 6 ITO,… 新人 0 落選 26.1 51
2 LDP 1 1996 aichi 7 NIWA… 新人 0 落選 25.9 49
3 LDP 1 1996 aichi 9 YOSH… 新人 0 落選 24.8 73
4 LDP 1 1996 aichi 10 MORI… 新人 0 落選 30.2 50
5 LDP 1 1996 aomo… 4 TSUS… 新人 0 落選 36.6 42
6 LDP 1 1996 chiba 3 MATS… 新人 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>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 x 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
- 目視によるモデルの当てはまり具合の評価方法
受信者操作特性（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.9145834auc2 <- performance(pr2, measure = "auc")
auc2@y.values[[1]] # model_2 のAUC[1] 0.9108907AUC による当てはまり ✔ AUCがほとんど差はないが、わずかに model_1 の当てはまりの方がよい
赤池情報量基準 (AID)
- 赤池情報量基準 (
AID: Akaike information criterion) を使う
- 多くの説明変数を使っても(= パラメータ数が増えても)、それに応じた分だけ尤度が増加しなければ、そのモデルは「良いモデル」とは言えない
- この基準で評価するのが
AID
AICが最小のモデルが「最も良いモデル」とされる
AIC(model_1)[1] 713.8303AIC(model_2)[1] 704.6358AID による当てはまり ✔ 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: 5stargazer()を使って、見やすく表示できる
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 に設定すると：
- 当選回数（previous）の
p値 = 2e-16
→ 当選回数は「影響がない」という帰無仮説を棄却
- 選挙費用（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:選挙費用と予測当選確率の予測値") xvals yvals upper lower
1 0.062710 0.07344159 0.1039941 0.04288908
2 1.087168 0.08559711 0.1178147 0.05337949
3 2.111627 0.09954837 0.1331773 0.06591946
4 3.136085 0.11548633 0.1502069 0.08076574
5 4.160543 0.13359742 0.1690458 0.09814904
6 5.185002 0.15405407 0.1898674 0.11824079
7 6.209460 0.17700316 0.2128989 0.14110741
8 7.233918 0.20255236 0.2384498 0.16665487
9 8.258377 0.23075539 0.2669313 0.19457950
10 9.282835 0.26159715 0.2988276 0.22436675
11 10.307293 0.29498038 0.3345676 0.25539317
12 11.331752 0.33071571 0.3743104 0.28712098
13 12.356210 0.36851732 0.4177842 0.31925039
14 13.380668 0.40800571 0.4642999 0.35171157
15 14.405127 0.44871875 0.5128879 0.38454960
16 15.429585 0.49013082 0.5624446 0.41781708
17 16.454043 0.53167880 0.6118416 0.45151601
18 17.478502 0.57279213 0.6600078 0.48557648
19 18.502960 0.61292382 0.7059934 0.51985429
20 19.527418 0.65157878 0.7490187 0.55413883
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:選挙費用と予測当選確率の予測値") xvals yvals upper lower
1 0.062710 0.07671468 0.1104017 0.04302770
2 1.087168 0.08950762 0.1251322 0.05388300
3 2.111627 0.10419315 0.1414822 0.06690406
4 3.136085 0.12096802 0.1595678 0.08236825
5 4.160543 0.14002147 0.1795155 0.10052744
6 5.185002 0.16152453 0.2014746 0.12157443
7 6.209460 0.18561706 0.2256370 0.14559707
8 7.233918 0.21239288 0.2522639 0.17252189
9 8.258377 0.24188398 0.2817090 0.20205899
10 9.282835 0.27404517 0.3144076 0.23368270
11 10.307293 0.30874107 0.3507780 0.26670413
12 11.331752 0.34573766 0.3910180 0.30045730
13 12.356210 0.38470033 0.4349038 0.33449683
14 13.380668 0.42520028 0.4817498 0.36865080
15 14.405127 0.46672973 0.5305385 0.40292097
16 15.429585 0.50872545 0.5800985 0.43735243
17 16.454043 0.55059839 0.6292454 0.47195141
18 17.478502 0.59176633 0.6768802 0.50665241
19 18.502960 0.63168574 0.7220576 0.54131388
20 19.527418 0.66987922 0.7640303 0.57572816
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. 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年の総選挙を比較しその違いを説明しなさい。
- Kieran Healy, DATA VISUALIZATION, Princeton, 2019
- 宋財泫 (Jaehyun Song)・矢内勇生 (Yuki Yanai)「私たちのR: ベストプラクティスの探究」
- 浅野正彦, 矢内勇生.『Rによる計量政治学』オーム社、2018年
- 浅野正彦, 中村公亮.『初めてのRStudio』オーム社、2018年
- Winston Chang, R Graphics Cookbook, O’Reilly Media, 2012.
- Kosuke Imai, Quantitative Social Science: An Introduction, Princeton University Press, 2017