Rを使った分析の準備
library(corrplot)
library(jtools)
library(margins)
library(ROCR)
library(patchwork)
library(prediction)
library(stargazer)
library(tidyverse)
選挙での「候補者の成績」と「選挙で使うお金の額」の関係を調べたいとする
→ つまり、お金を使う候補者ほど選挙結果がいいのか?という疑問
説明したい対象(=応答変数)としては 2 つ想定できる:
①得票率(もしくは得票数) ・・・連続変数 => 回帰分析
②当落・・・当選したら1 、落選したら 0 ・・・2値変数 => ロジスティック回帰分析
回帰分析とロジスティック回帰分析の違いを視覚的に示してみる
回帰分析 | ロジスティック回帰分析 | |
---|---|---|
応答変数 | 得票率(%) | 当落(0 or 1) |
説明変数 | 選挙費用(百万円) | 選挙費用(百万円) |
・回帰分析の応答変数は得票率=「連続変数」
・ロジスティック回帰分析の応答変数は選挙での当落= 2 値変数 (1 or
0)
・それぞれの図に回帰直線を図に書き入れてみる
・回帰分析では、選挙費用を使うほど得票率が大きくなっている =>
問題なし
・ロジスティック回帰分析でも、選挙費用を使うほど得票率が大きくなっている
→ しかし、3つ問題点がある
1. 当選確率(0〜1)が 1 (100%) を超えてしまうこと
2. 当選確率(0〜1)がマイナスの値をとってしまうこと
3. 当てはまり感がない(むりやり当てはめている感が強い)
・選挙費用が少ない200万円以下の人はほとんど当選していない
・選挙費用が1000万円の候補者は50%程当選している
・選挙費用が1600万円の候補者は75%以上が当選している
・選挙費用ごとの「当選率」の対数をとってロジット変換する
→ S字が直線になる
→ しかし、実際のデータを使うと完璧な直線ではないものの、かなり直線ぽくなる
→ この曲線を使ってロジスティック回帰分析を推定する
・ロジスティック回帰分析では「割合」ではなく「オッズ」で結果が示される
Logit Model
の基礎知識:オッズ
Odds
) とは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%」という意味
ロジスティック回帰分析の手順は次のとおりである
仮説1・選挙費を使えば使うほど、小選挙区での当選確率は大きい
帰無仮説
・選挙費の額は、小選挙区での当選確率とは関係がない
ロジスティック回帰分析における検定では、重回帰分析における検定と同様、得られた
p 値
が有意水準よりも小さいときに帰無仮説を棄却し、対立仮説を受け容れる
仮説1を検証するために次のモデルを考える
wlsmd
) を縦軸に、選挙費用 (expm
)
を横軸にした散布図を表示するデータの準備 (hr96-24.csv
)
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
の中身を表示する
[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
の中身を確認[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
の中身を表示する
# 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 |
候補者の選挙費用(百万円) |
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
の記述統計を表示して確認する 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
NA's
が消えていることがわかる hr21
) の記述統計を表示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 |
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")
Logit Model の基礎知識:オッズ
を参照)\[\frac{P}{1-p} = a + 選挙費用b + 当選回数c\]
\[log\frac{P}{1-p} = a + 選挙費用b + 当選回数c\]
Model 1
)expm
:
単位100万円)で衆院選小選挙区での当落 (wlsmd
)
を説明するモデルを考えるModel 1
\[Pr(当選)=logit^{−1}(𝛼+\beta_1選挙費用 + \beta_2当選回数)\]
glm()
を使うglm()
はロジスティック回帰以外でも頻繁に使う関数family = binomial(link = "logit")
を指定logit
とは「オッズの対数
」= Log Odds
のことmodel_1 <- glm(wlsmd ~ expm + previous,
data = hr21,
family = binomial(link = "logit")) # 係数を「オッズの対数」に指定
tidy()
を使って、推定結果を確認する# 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]])}\]
tidy()
を使って、推定結果を再度確認してみる# 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予測当選確率と限界効果
0 〜 1
ロジスティック回帰分析の係数は、重回帰分析の係数と同じように解釈できない
ロジスティック曲線は直線ではない → 「曲線」
→ 説明変数
1単位の変化が応答変数に与える影響(=傾き)は説明変数の値によって異なる
説明変数が応答変数に与える影響を考えるためには
→ 説明変数の値を特定し、「特定された値における影響」を計算する必要あり
→ 限界効果 (Marginal
Effects)
上図の例では、x
の値 1, 5, 10, 15, 20
それぞれの時点での直線の傾きが限界効果
previous
)
が3回の候補者に関して次の二つのケースを考えてみる\[Pr(当選) = \frac{1}{1+exp(-(-1.98 + 0.0735expm + 0.285previous))}\]
p_0 <- predict(model_1, type = "response", #予測当選確率を表示したい時の指定
newdata = data_frame(previous = 3, expm = 0))
p_0
1
0.1642294
1
0.1889165
(p_1 - p_0)
1
0.02468706
1
0.2163539
(p_2 - p_1)
なので 1
0.02743738
まとめ ・選挙費用を 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万円)
at(previous) at(expm) expm
3 4 0.03424
選挙費用の限界効果(previous = 3, expm: 0 〜 2,800万円)
at
に変数のリストを渡すことで指定する 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
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):
仮説検証の結論(仮説 1)
実質的有意性:
・選挙費用が大きくなるにつれて、候補者の当選確率は大きくなる
・選挙費用が小さいとき、予測当選確率は緩やかに上昇
・選挙費用が中程度になると曲線の傾きが少しずつ急になる
・選挙費用が大きくなると再び傾きが小さくなる
実質的有意性:
・選挙費用の限界効果は、約1000万円の時が最大
・選挙費用がそれより多くても少なくても効果が小さい
→ 選挙費用が約1300万円の時、当選確率に対する影響力が最大
統計的有意性:
全ての選挙費用の範囲で統計的に有意
previous
) の値を平均値に固定した値である 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
) の記述統計を確認する
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
候補者の大多数の当選回数が 0 だとわかる
17回の当選回数を0から2回ずつ区切り、0, 2, 4, …,
16まで分割する
→ previous = seq(0, 16, by = 2)
候補者の選挙費用 (expm
)
の記述統計を確認する
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.009319 3.435080 5.899882 6.434585 8.692808 27.443685
- 選挙費用は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)
次に、当選回数が変わっても、選挙費用が当選確率に予測値に与える影響が統計的に有意なのかどうかをチェック
当選回数 (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総選挙")
previous = 15
までは、95%信頼区間が y = 0
の赤い点線に触れていないprevious = 15
までは、選挙費用が当選確率に予測値に与える影響が統計的に有意