library(tidyverse)
library(stargazer)
library(patchwork)
Plobit Model
や Logit Model
のような非線形モデルでは、OLS
推定は使えない
→ 最尤法 (MLE: maximum likelihood esitmation
) を使って推定する
推定法 | 詳細 |
---|---|
最小二乗法 (OLS) | 不偏性・効率的 (BLUE) |
最尤法 (MLE) | 一致性・効率的・サンプル数が十分に大きければ正規分布 |
次のような事例を考える
試験の答案を返された
10点満点中、自分の試験点数は「4 点」
これは良い点なのか、悪い点なのか?
調べる方法:回りの友人達の点数を聞いてみる
友人 4 人の点数を確かめた (6, 7, 3, 5点)
自分の点数を含めて 5 つのサンプル(サンプルの平均値は 5)
← 3 + 4 + 5 + 6 + 7 = 25 ÷ 5 = 5
5 人分サンプルから母平均(試験全体の平均値)を推定する
この母集団の特徴
正規分布
標準偏差 = 1
母平均 (μ
) は不明
得られたサンプルの外れ度合いを最小にする母平均は?
得れたサンプルを生み出す蓋然性 (likelihood) を最大にする母平均は?
μ
) がもし 5 点ならμ = 5
の場合)\[(3である蓋然性) x (4である蓋然性) x (5である蓋然性) x (6である蓋然性) x (7である蓋然性) \]
μ
) がもし 6 点なら
- 平均 6、標準偏差 1 の正規分布の確率密度関数を使えば、母集団からサンプルを取った時に、それが「3 点である蓋然性」「4 点である蓋然性」「5 点である蓋然性」「6 点である蓋然性」「7 点である蓋然性」を計算できる
μ = 6
の場合)\[(3である蓋然性)x (4である蓋然性) x (5である蓋然性) x (6である蓋然性) x (7である蓋然性) \]
μ
) が 4, 5, 6 点の場合の比較これらの平均値(μ
)の値の中で、どれが尤(もっと)もありえそう (likely
) な値か?
ここでは「サンプルデータを生み出す蓋然性を最大化するような平均値(μ
)の値」は「5 個のサンプルの平均値」
これを「母集団の μ の推定値」として使う・・・最尤法 (MLE
) の考え方
母集団からランダムに採取したサンプルが得られる蓋然性のこと
尤度は私たちが推定したい母集団のパラメータ(例えば母平均 μ
など)によって決まる
最尤法 (MLE) は、尤度を最大化するような値をパラメータの推定値として使う
likelihood function
)尤度をパラメータを用いた式で表したもの
実際には、計算のしやすくするため、尤度関数の対数をとって対数尤度関数を使う
最尤法 (MLE) による推定方法 ・尤度関数を求める
・尤度関数を最大化するようなパラメータの値を解く
パラメータの計算方法 | 詳細 |
---|---|
解析的 (analytical) 方法 | この方法で最尤法はなかなか解けない |
数値的 (numerical)方法 | Newton-Raphson法、BHHH法、DFP法、BFGS法などを使って数値的 (numerical) に解く |
最尤法による推定は、解析的 (analytical) な推定は困難
→数値的 (numerical) に解く
→単純計算が特異なコンピュータに解いてもらう
最尤推定法 (MLE) | 最小二乗法 (OLS) |
---|---|
尤度比検定 (likelihood ratio test: LR test)・・・問題あり | F 検定 |
疑似決定係数 (pseudo \(R^2\))・・・問題あり | |
OLS における F 検定の解説 - F 検定では複数の帰無仮説を同時に検定する
- 複数の帰無仮説を同時に検定する理由:
→ 1 つずつ検定すると有意でないが、全体で検定する有意なことがありうるから
→ 2 つのモデルを作る
(1) 無制限モデル (unrestricted model
) — β1
と β2
を含むモデル
(2) 制限モデル (restricted model
) — β1
と β2
を含まないモデル
→ それぞれの \(R^2\) を計算
→ 「β1
と β2
を含むことによって、 \(R^2\) が十分に増えたか?」を検定
MLE における尤度比検定 (LR test) の解説と問題点
→ 2 つのモデルを作る
(1) 無制限モデル (unrestricted model
) — β1
と β2
を含むモデル
(2) 制限モデル (restricted model
) — β1
と β2
を含まないモデル
→ それぞれの \(R^2\) を計算
→ 「β1
と β2
を含むことによって、 尤度が十分に増えたか?」を検定
「モデル全体(説明変数の組み合わせが)無意味かどうか」というおおざっぱな評価しかできない
MLE における疑似決定係数 (pseudo\(R^2\) )の解説と問題点 - OLS における決定係数 (\(R^2\)) に類似した概念として考案
\[pseudoR^2 =1 - \frac{無制約モデルの対数尤度}{ 制約モデルの対数尤度}\]
AID: Akaike information criterion
) を使う\[AIC = -2×対数尤度+2×パラメータ数\]
- AIC が最小のモデルが「最も良いモデル」とされる
- 多くの説明変数を使っても(= パラメータ数が増えても)、それに応じた分だけ尤度が増加しなければ、そのモデルは「良いモデル」とは言えない
最尤推定法 (MLE) | 最小二乗法 (OLS) |
---|---|
尤度 | 決定係数 (\(R^2\)) |
AIC | 調整済み決定係数 (Adjusted \(R^2\)) |
MLE
における尤度は、OLS
における決定係数 (\(R^2\)) に相当MLE
における AIC
は、OLS
における調整済み決定係数 (Adjusted \(R^2\))に相当AIC
は MLE
におけるモデルの評価基準として最も広く使われているOLS
と MLE
の前提条件OLS
) を用いるMLE
) を用いるOLS
では推定できない非線形モデルでも MLE
は推定できるOLS
で前提としている仮定を MLE
では明示的にモデル中に取り込めるMLE
で分析すべきモデルを OLS
で分析するとどうなるかを示すRProject folder
内の data
フォルダに入れるdata
フォルダ内から read_csv
で読み取り df
というデータフレーム名をつける<- read_csv("data/hr96-17.csv", na = ".") df
names(df)
[1] "year" "pref" "ku" "kun"
[5] "mag" "rank" "nocand" "seito"
[9] "name" "j_name" "previous" "gender"
[13] "age" "wl" "exp" "status"
[17] "vote" "voteshare" "eligible" "turnout"
[21] "castvotes" "seshu_dummy" "jiban_dummy" "nojiban_dummy"
[25] "jiban_seshu" "nojiban_seshu"
df1
と名前を付ける<- df %>%
df1 filter(year == 2012) %>%
select(wl, seshu_dummy, exp, eligible, status)
wl = 1
は小選挙区当選者、wl = 1
は復活当選者、wl = 0
は落選者wl = 1
が当選者(小選挙区と比例復活), wl = 0
が落選者<- mutate(df1, wl = as.numeric(wl == 1 | wl == 2)) df1
exp
と eligible
を使って、有権者一人あたりに使う選挙費用 (exppv
) を作る<- mutate(df1, exppv = exp / eligible) # eligible は小選挙区ごとの有権者数 df1
summary(df1$exppv)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00126 6.08015 14.86205 17.16598 23.62160 94.59168 14
status = 0
は新人、status = 1
は現職、 status = 2
は元職status = 0
がその他、status = 1
が現職<- mutate(df1, inc = as.numeric(status == 1 )) df1
<- df1 %>%
df1 select(wl, seshu_dummy, inc, exppv)
R-Markdown
で html
表示する際にはチャンクオプションで {r, results = “asis”}
と指定する。stargazer(as.data.frame(df1),
type ="html",
digits = 2)
Statistic | N | Mean | St. Dev. | Min | Pctl(25) | Pctl(75) | Max |
wl | 1,294 | 0.33 | 0.47 | 0 | 0 | 1 | 1 |
seshu_dummy | 1,294 | 0.09 | 0.29 | 0 | 0 | 0 | 1 |
inc | 1,294 | 0.31 | 0.46 | 0 | 0 | 1 | 1 |
exppv | 1,280 | 17.17 | 13.51 | 0.001 | 6.08 | 23.62 | 94.59 |
変数名 | 詳細 |
---|---|
wl | 選挙の当落: 1 = 当選(小選挙区+復活当選)、0 = 落選 |
seshu_dummy | 世襲候補者ダミー: 1 = 世襲、0 = 非世襲 |
inc | 現職ステータス: 1 現職、0 = 非現職 |
exppv | 有権者一人あたりに費やす選挙費用(円) |
OLS モデル
と MLE モデル
の比較wl
)」が応答変数これに線形回帰を(無理矢理)当てはめ、\(\mathrm{\hat{Y} > 0.5}\) なら当選、その他の場合は落選と予測する
2 値の変数 (wl
) を OLS
による回帰に使うことは道理に叶っている
次のようなモデルを考える
線形回帰で得られる \(\mathrm{β_1,β_2, β_3}\) が \(\mathrm{Pr(当落|X)}\) の推定値
線形回帰 (OLS
) を行うと、推定値が区間[0, 1]
から外れてしまうことがある
→ これでは「確率」とは言えない
それでも、この予測値は順序付けのために使うことができる
「確率」の大雑把な大雑把な推定値として使うことが可能
OLS
では誤差項 (ε
) が正規分布しているという仮定がある
→OLS
では次の仮定を置いていることになる
\[「選挙の当落」は 「β_0 + β_1世襲+ β_2現職 + β_3選挙費用」を平均に正規分布している\]
wl
) の分布は次のとおりhist(df1$wl)
応答変数 (wl
) は 0, 1 の値をとる「二値変数」
この応答変数 (wl
) は上記の OLS
の仮定を満たさない
→ これはOLS
ではなく MLE
で分析すべき
「良いモデル」という観点から、MLE
で分析すべきモデルを OLS
で分析してみる
MLE
推定では「一般化線形モデル」(GLM: generalized linear model
) の一つ である Logit モデル(ロジスティック回帰)を使う
モデル名 | 応答変数 | 説明変数 |
---|---|---|
OLS | 選挙の当落 (wl) | 世襲 (seshu_dummy), 現職 (inc), 選挙費 (exppv) |
MLE (Logit) | 選挙の当落 (wl) | 世襲 (seshu_dummy), 現職 (inc), 選挙費 (exppv) |
lm 関数
)<- lm(wl ~ seshu_dummy + inc + exppv,
ols data = df1)
glm 関数
)<- glm(wl ~ seshu_dummy + inc + exppv, data = df1,
logit family = binomial(link = "logit")) # 誤差を二項分布に指定
stargazer(ols, logit,
type = "html",
digits = 2)
Dependent variable: | ||
wl | ||
OLS | logistic | |
(1) | (2) | |
seshu_dummy | 0.36*** | 2.04*** |
(0.04) | (0.29) | |
inc | -0.01 | -0.14 |
(0.03) | (0.16) | |
exppv | 0.02*** | 0.09*** |
(0.001) | (0.01) | |
Constant | 0.04** | -2.51*** |
(0.02) | (0.14) | |
Observations | 1,280 | 1,280 |
R2 | 0.29 | |
Adjusted R2 | 0.29 | |
Log Likelihood | -607.95 | |
Akaike Inf. Crit. | 1,223.90 | |
Residual Std. Error | 0.40 (df = 1276) | |
F Statistic | 175.52*** (df = 3; 1276) | |
Note: | p<0.1; p<0.05; p<0.01 |
\[選挙での当落 = 0.04 + 0.36世襲 - 0.01現職 + 0.02選挙費用\]
・回帰式の Intercept
の係数: 0.04
が意味すること
・seshu_dummy = 0, inc = 0, exppv = 0
の候補者が当選する値は 0.04
→ これでは「確率」とは言えない
- それでも、この予測値は順序付けのために使うことができる
- 「確率」の大雑把な大雑把な推定値として使うことが可能
→ seshu_dummy = 0, inc = 0, exppv = 0
の候補者が当選する値は 4%
\[選挙での当落 = -2.51 + 2.04世襲 - 0.14現職 + 0.09選挙費用\]
logit (=Log Odds)
Odds
と当選確率を計算する ・回帰式の Intercept
の係数: -2.51
が意味すること
・seshu_dummy = 0, inc = 0, exppv = 0
の候補者の Odds
は
exp(-2.51)
[1] 0.08126824
Odds
の値が 1 であれば、当選する確率は 50 %Odds
の値が 0.08 なので当選確率は 50% より遙かに小さいとわかるpredict(logit, type = "response",
newdata = data_frame(seshu_dummy = 0, inc = 0, exppv = 0))
1
0.07511063
・logit
におけるseshu_dummy = 0, inc = 0, exppv = 0
の候補者の当選確率は 0.075 (= 7.5%)
・参考:OLS
における seshu_dummy = 0, inc = 0, exppv = 0
の時の切片は 0.04 (= 4%)
<- ggplot(df1, aes(x = exppv, y = wl)) +
OLS geom_smooth(method = "lm", se = FALSE) +
labs(x = "選挙費用(円)", y = "選挙での当落", title = "回帰直線 (OLS)") +
geom_hline(yintercept = c(0, 1), color = "gray") +
geom_jitter(width = 0.05, height = 0.05) +
ylim(0,1) +
theme_minimal(base_family = "HiraKakuProN-W3")
print(OLS)
<- margins::cplot(logit,
LOGIT x = "exppv",
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 = "選挙費用と予測当選確率 (logit)") +
theme_minimal(base_family = "HiraKakuProN-W3")
xvals yvals upper lower
1 0.001264246 0.08590607 0.1087962 0.06301598
2 3.942531663 0.11839288 0.1439283 0.09285749
3 7.883799080 0.16100126 0.1883026 0.13369997
4 11.825066496 0.21520085 0.2434413 0.18696041
5 15.766333913 0.28152390 0.3109640 0.25208383
6 19.707601329 0.35893854 0.3919028 0.32597425
7 23.648868746 0.44447169 0.4842674 0.40467600
8 27.590136162 0.53342700 0.5816503 0.48520375
9 31.531403579 0.62030634 0.6756541 0.56495858
10 35.472670995 0.70010327 0.7591228 0.64108370
11 39.413938412 0.76936526 0.8279283 0.71080224
12 43.355205829 0.82659308 0.8811645 0.77202164
13 47.296473245 0.87198363 0.9202666 0.82370070
14 51.237740662 0.90683205 0.9478182 0.86584588
15 55.179008078 0.93292381 0.9666051 0.89924250
16 59.120275495 0.95209456 0.9790850 0.92510414
17 63.061542911 0.96598597 0.9871978 0.94477417
18 67.002810328 0.97595096 0.9923717 0.95953024
19 70.944077745 0.98304777 0.9956107 0.97048488
20 74.885345161 0.98807591 0.9975984 0.97855345
+ LOGIT OLS
seshu_dummy
) と現職 (inc
) の値を平均値に固定した時に「有権者一人あたりに費やした選挙費用」とその「候補者の当選確率」の関係を示しているpredict(logit, type = "response",
newdata = data_frame(seshu_dummy = 0.09, inc = 0.31, exppv = 25))
1
0.4735126
predict(logit, type = "response",
newdata = data_frame(seshu_dummy = 0.09, inc = 0.31, exppv = 50))
1
0.8964211
seshu_dummy = 1
、inc = 1
) が、有権者一人あたりの選挙費用として 25円使うと、その候補者の当選確率は約 84 %だと予測できるpredict(logit, type = "response",
newdata = data_frame(seshu_dummy = 1, inc = 1, exppv = 25))
1
0.8393427
fit
) 比較OLS
モデル と Logit
モデルの AIC
を比較AIC(ols)
[1] 1266.412
AIC(logit)
[1] 1223.901
Logit
モデルの AIC
の値 (1223
) の方が OLS
モデルの AIC
の値 (1266
)より小さいLogit
モデルの方がよいモデルまとめ - 応答変数が 0, 1 の二値変数を OLS
モデルと Logit
モデルで分析
→ ほぼ同様の結果
- この分析では AIC
が小さい Logit
モデルの方がより良いモデルといえる
→ 質的な応答変数の場合には非線形回帰 (MLE
) を用いた方がより良い
References
G.James他(落海他訳)『Rによる統計的学習入門』朝倉書店、2020年
浅野正彦, 矢内勇生.『Rによる計量政治学』オーム社、2018年
森田 果.『実証分析入門』日本評論社、2014年
Kosuke Imai, Quantitative Social Science: An Introduction, Princeton University Press, 2017