• このセクションで使う R パッケージ一覧
library(tidyverse)
library(stargazer)
library(patchwork)
  • Plobit ModelLogit Model のような非線形モデルでは、OLS 推定は使えない
    → 最尤法 (MLE: maximum likelihood esitmation) を使って推定する

1. 最小二乗法と最尤法の違い

  • 最尤法は「さいゆうほう」と読む
  • 最も尤(もっと)もらしい方法という意味
最小二乗法 (OLS) 不偏性・効率的 (BLUE)
最尤法(MLE) 一致性・効率的・サンプル数が十分に大きければ正規分布
  • 次のような事例を考える

  • 試験の答案を返された

  • 10点満点中、自分の試験点数は「4 点」

  • これは良い点なのか、悪い点なのか?

  • 調べる方法:回りの友人達の点数を聞いてみる

  • 友人 4 人の点数を確かめた (6, 7, 3, 5点)

  • 自分の点数を含めて 5 つのサンプル(サンプルの平均値は 5)
    ← 3 + 4 + 5 + 6 + 7 = 25 ÷ 5 = 5

  • 5 人分サンプルから母平均(試験全体の平均値)を推定する  

  • この母集団の特徴

  1. 正規分布
  2. 標準偏差 = 1
  3. 母平均 (\(μ\)) は不明
  • OLS とMLE の基本的な考え方:
最小二乗法 (OLS) 得られたサンプルの外れ度合いを最小にする母平均は?
最尤法 (MLE) 得れたサンプルを生み出す蓋然性 (likelihood) を最大にする母平均は?

クラス全体の試験の平均値(\(μ\)) がもし 5 点なら

  • 平均 5、標準偏差 1 の正規分布の確率密度関数を使えば、母集団からサンプルを取った時に、それが「3 点である蓋然性」「4 点である蓋然性」「5 点である蓋然性」「6 点である蓋然性」「7 点である蓋然性」を計算できる

5 個のサンプルを同時に取得する蓋然性(μ = 5 の場合)

  • この 5 個のサンプルが無為作為に取られたのであれば、この 5 個のサンプルを同時に取得する蓋然性は次の式で求められる
(3である蓋然性) x (4である蓋然性) x (5である蓋然性) x (6である蓋然性) x (7である蓋然性)

クラス全体の試験の平均値(\(μ\)) がもし 6 点なら
- 平均 6、標準偏差 1 の正規分布の確率密度関数を使えば、母集団からサンプルを取った時に、それが「3 点である蓋然性」「4 点である蓋然性」「5 点である蓋然性」「6 点である蓋然性」「7 点である蓋然性」を計算できる


5 個のサンプルを同時に取得する蓋然性(μ = 6 の場合)

  • この 5 個のサンプルが無為作為に取られたのであれば、この 5 個のサンプルを同時に取得する蓋然性は次の式で求められる

(3である蓋然性)x (4である蓋然性) x (5である蓋然性) x (6である蓋然性) x (7である蓋然性)

クラス全体の試験の平均値(\(μ\)) が 4, 5, 6 点の場合の比較

  • これらの平均値(\(μ\))の値の中で、どれが尤(もっと)もありえそう (likely) な値か?
    → 蓋然性がもっとも大きくなる \(μ\)

  • ここでは「サンプルデータを生み出す蓋然性を最大化するような平均値(\(μ\))の値」は「5 個のサンプルの平均値」

  • これを「母集団の μ の推定値」として使う・・・最尤法 (MLE) の考え方

尤度(ゆうど)

  • 母集団からランダムに採取したサンプルが得られる蓋然性のこと
  • 尤度は私たちが推定したい母集団のパラメータ(例えば母平均 \(μ\) など)によって決まる
  • 最尤法 (MLE) は、尤度を最大化するような値をパラメータの推定値として使う

尤度関数 (likelihood function)

  • 尤度をパラメータを用いた式で表したもの
  • 実際には、計算のしやすくするため、尤度関数の対数をとって対数尤度関数を使う

最尤法 (MLE) による推定方法 ・尤度関数を求める
・尤度関数を最大化するようなパラメータの値を解く

2. 最尤法による具体的な計算方法

パラメータの計算方法 詳細
解析的 (analytical) 方法 この方法で最尤法はなかなか解けない
数値的 (numerical)方法 Newton-Raphson法、BHHH法、DFP法、BFGS法などを使って数値的 (numerical) に解く

2.1 数値演算的方法によるパラメータの計算方法

  • 最尤法による推定は、解析的 (analytical) な推定は困難
    → 数値的 (numerical) に解く
    → 単純計算が特異なコンピュータに解いてもらう

2.2 当てはまりの良さの指標

  • 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{無制約モデルの対数尤度}{ 制約モデルの対数尤度}\]

  • 無制約モデルが尤度を改善しない → 擬似決定係数は 0 に近づく
  • 無制約モデルが尤度を改善した → 擬似決定係数は 1 に近づく

疑似決定係数 (pseudo \(R^2\)) の限界

  • 決定係数 (\(R^2\)) と違って、疑似決定係数 (pseudo \(R^2\)) は「目的変数の変化のうちのどれほどを当該モデルが説明しているか」という指標として解釈できない
  • 複数の擬似決定係数 (pseudo \(R^2\)) を比較して「こちらの方がいいモデル」とはいえない
  • 調整済み決定係数 (Adjusted \(R^2\)) と違って、擬似決定係数 (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\))に相当
  • AICMLE におけるモデルの評価基準として最も広く使われている

3. 総選挙データを使った比較

3.1 OLSMLE の前提条件

  • 量的な応答変数の場合・・・線形回帰 (OLS) を用いる
  • 質的な応答変数の場合・・・非線形回帰 (MLE) を用いる
  • OLS では推定できない非線形モデルでも MLE は推定できる
  • OLS で前提としている仮定を MLE では明示的にモデル中に取り込める
  • ここでは本来 MLE で分析すべきモデルを OLS で分析するとどうなるかを示す

3.2 データの準備 (hr96-17.csv)

  • 衆議院議員総選挙の得票データ hr96-17.csv をダウンロード
    RProject folder内の data フォルダに入れる
  • dataフォルダ内から read_csv で読み取り df というデータフレーム名をつける
df <- read_csv("data/hr96-17.csv", na = ".")
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"
  • 2012年の総選挙と必要な 5 つの変数を抜き出し df1 と名前を付ける
df1 <- df %>% 
  filter(year == 2012) %>% 
  select(wl, seshu_dummy, exp, eligible, status)
  • wl = 1 は小選挙区当選者、wl = 2 は復活当選者、wl = 0 は落選者
  • これを次のように変更:
  • wl = 1 が当選者(小選挙区と比例復活), wl = 0 が落選者
df1 <- mutate(df1, wl = as.numeric(wl == 1 | wl == 2)) 
  • expeligible を使って、有権者一人あたりに使う選挙費用 (exppv) を作る
df1 <- mutate(df1, exppv = exp / eligible) # eligible は小選挙区ごとの有権者数
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 が現職
df1 <- mutate(df1, inc = as.numeric(status == 1 )) 
df1 <- df1 %>% 
  select(wl, seshu_dummy, inc, exppv)
  • これら変数の記述統計は次のとおり
  • R-Markdownhtml 表示する際にはチャンクオプションで {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 有権者一人あたりに費やす選挙費用(円)

3.3 OLS モデルMLE モデルの比較

  • ここでは2012年の総選挙における候補者の「選挙での当落 (wl)」が応答変数
  • wl は 0 が落選、1 が当選の 2 値変数(ダミー変数)

  • これに線形回帰を(無理矢理)当てはめ、\(\mathrm{\hat{Y} > 0.5}\) なら当選、その他の場合は落選と予測する

  • 2 値の変数 (wl) を OLS による回帰に使うことは道理に叶っている

  • 次のようなモデルを考える

 \(\mathrm{選挙の当落 = β_0 + β_1世襲+ β_2現職 + β_3選挙費用 +ε }\)
  • 線形回帰で得られる \(\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 関数)

ols <- lm(wl ~ seshu_dummy + inc + exppv, 
          data = df1)

Logit モデル(ロジスティック回帰) (glm 関数)

logit <- glm(wl ~ seshu_dummy + inc + exppv, data = df1,
               family = binomial(link = "logit")) # 誤差を二項分布に指定

3.4 分析結果の表示

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.04
    → これは「確率」とは言えない
    それでも、この予測値は順序付けのために使うことができる
    「確率」の大雑把な大雑把な推定値として使うことが可能
    → seshu_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%)    

3.5 グラフによる比較

当選回数と当落のプロット (OLS)

OLS <- ggplot(df1, aes(x = exppv, y = wl)) +
  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)

LOGIT <- margins::cplot(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

- 二つのグラフを並べて表示させる

OLS + LOGIT

3.6 分析結果の解釈

  • 上の右側のグラフは、世襲 (seshu_dummy) と現職 (inc) の値を平均値に固定した時に「有権者一人あたりに費やした選挙費用」とその「候補者の当選確率」の関係を示している 
  • 例えば、この平均的な候補者が有権者一人あたりに選挙費用として 25円使うと、その候補者の当選確率は約 47% だと予想できる 
  • R を使うと、次のように計算できる
predict(logit, type = "response",
        newdata = data_frame(seshu_dummy = 0.09, inc = 0.31, exppv = 25))
        1 
0.4747446 
  • 候補者一人あたりの選挙資金を50円使うと、その候補者の当選確率は約 90% まで上がる
predict(logit, type = "response",
        newdata = data_frame(seshu_dummy = 0.09, inc = 0.31, exppv = 50))
        1 
0.8991363 
  • 例えば、世襲の現職候補者 (seshu_dummy = 1inc = 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) を用いた方がより良い

参考文献