R パッケージ
一覧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 人分サンプルから母平均(試験全体の平均値)を推定する
この母集団の特徴
最小二乗法 (OLS) | 得られたサンプルの外れ度合いを最小にする母平均は? |
最尤法 (MLE) | 得れたサンプルを生み出す蓋然性 (likelihood) を最大にする母平均は? |
クラス全体の試験の平均値(\(μ\)) がもし 5 点なら
5 個のサンプルを同時に取得する蓋然性(μ = 5
の場合)
クラス全体の試験の平均値(\(μ\)) がもし 6 点なら
- 平均 6、標準偏差 1 の正規分布の確率密度関数を使えば、母集団からサンプルを取った時に、それが「3 点である蓋然性」「4 点である蓋然性」「5 点である蓋然性」「6 点である蓋然性」「7 点である蓋然性」を計算できる
5 個のサンプルを同時に取得する蓋然性(μ = 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
における推定値 = サンプルデータを生成する蓋然性 (= 尤度)最尤推定法 (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{無制約モデルの対数尤度}{ 制約モデルの対数尤度}\]
疑似決定係数 (pseudo \(R^2\)) の限界
解決策
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
で分析するとどうなるかを示すhr96-17.csv
)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" "wl" "nocand"
[9] "seito" "j_name" "name" "term"
[13] "gender" "age" "exp" "status"
[17] "vote" "voteshare" "eligible" "turnout"
[21] "castvotes" "seshu_dummy" "jiban_seshu" "nojiban_seshu"
df1
と名前を付ける<- df %>%
df1 filter(year == 2012) %>%
select(wl, seshu_dummy, exp, eligible, status)
wl = 1
は小選挙区当選者、wl = 2
は復活当選者、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.10793 14.86928 17.18456 23.65808 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.18 | 13.50 | 0.001 | 6.11 | 23.66 | 94.59 |
変数名 | 詳細 |
wl | 選挙の当落: 1 = 当選(小選挙区+復活当選)、0 = 落選 |
seshu_dummy | 世襲候補者ダミー: 1 = 世襲、0 = 非世襲(地盤世襲 or 非世襲) |
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) |
OLS モデル(重回帰): (lm 関数
)
<- lm(wl ~ seshu_dummy + inc + exppv,
ols data = df1)
Logit モデル(ロジスティック回帰) (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.53*** |
(0.02) | (0.14) | |
Observations | 1,280 | 1,280 |
R2 | 0.29 | |
Adjusted R2 | 0.29 | |
Log Likelihood | -606.24 | |
Akaike Inf. Crit. | 1,220.48 | |
Residual Std. Error | 0.40 (df = 1276) | |
F Statistic | 177.20*** (df = 3; 1276) | |
Note: | p<0.1; p<0.05; p<0.01 |
推定値の比較
OLS モデルの推定結果
\[選挙での当落 = 0.04 + 0.36世襲 - 0.01現職 + 0.02選挙費用\]
Intercept
の係数: 0.04
が意味すること seshu_dummy = 0, inc = 0, exppv = 0
の候補者が当選する値は 0.04seshu_dummy = 0, inc = 0, exppv = 0
の候補者が当選する値は 4%Logit モデルの推定結果
\[選挙での当落 = -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.0738792
logit
におけるseshu_dummy = 0, inc = 0, exppv = 0
の候補者の当選確率は 0.075 (= 7.5%)
OLS
における seshu_dummy = 0, inc = 0, exppv = 0
の時の切片は 0.04 (= 4%)
当選回数と当落のプロット (OLS
)
<- 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)
当選回数と当選確率のプロット (Logit
)
<- 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.08436903 0.1070096 0.06172844
2 3.942531663 0.11674865 0.1421067 0.09139057
3 7.883799080 0.15939185 0.1865980 0.13218571
4 11.825066496 0.21383980 0.2420541 0.18562552
5 15.766333913 0.28067680 0.3101311 0.25122252
6 19.707601329 0.35886807 0.3918741 0.32586205
7 23.648868746 0.44535566 0.4852269 0.40548441
8 27.590136162 0.53528441 0.5835896 0.48697924
9 31.531403579 0.62297600 0.6783308 0.56762119
10 35.472670995 0.70329242 0.7621501 0.64443472
11 39.413938412 0.77274076 0.8309186 0.71456290
12 43.355205829 0.82986612 0.8838368 0.77589546
13 47.296473245 0.87495571 0.9224791 0.82743236
14 51.237740662 0.90940008 0.9495449 0.86925528
15 55.179008078 0.93506069 0.9678910 0.90223037
16 59.120275495 0.95382250 0.9800065 0.92763846
17 63.061542911 0.96735303 0.9878365 0.94686957
18 67.002810328 0.97701450 0.9928006 0.96122840
19 70.944077745 0.98386447 0.9958894 0.97183953
20 74.885345161 0.98869668 0.9977728 0.97962060
- 二つのグラフを並べて表示させる
+ LOGIT OLS
seshu_dummy
) と現職 (inc
) の値を平均値に固定した時に「有権者一人あたりに費やした選挙費用」とその「候補者の当選確率」の関係を示している predict(logit, type = "response",
newdata = data_frame(seshu_dummy = 0.09, inc = 0.31, exppv = 25))
1
0.4747446
predict(logit, type = "response",
newdata = data_frame(seshu_dummy = 0.09, inc = 0.31, exppv = 50))
1
0.8991363
seshu_dummy = 1
、inc = 1
) が、有権者一人あたりの選挙費用として 25円使うと、その候補者の当選確率は約 84 %だと予測できるpredict(logit, type = "response",
newdata = data_frame(seshu_dummy = 1, inc = 1, exppv = 25))
1
0.8394015
fit
) 比較OLS
モデル と Logit
モデルの AIC
を比較AIC(ols)
[1] 1262.83
AIC(logit)
[1] 1220.48
Logit
モデルの AIC
の値 (1223
) の方が OLS
モデルの AIC
の値 (1266
)より小さいLogit
モデルの方がよいモデルまとめ - 応答変数が 0, 1 の二値変数を OLS
モデルと Logit
モデルで分析
→ ほぼ同様の結果
- この分析では AIC
が小さい Logit
モデルの方がより良いモデルといえる
→ 質的な応答変数の場合には非線形回帰 (MLE
) を用いた方がより良い