• このセクションで使う R パッケージ一覧
library(tidyverse)
library(stargazer)
library(patchwork)
  • Plobit Model や Logit 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\))に相当
  • AIC は MLE におけるモデルの評価基準として最も広く使われている

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

3.1 OLS と MLE の前提条件

  • 量的な応答変数の場合・・・線形回帰 (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] "wl"            "rank"          "nocand"        "seito"        
 [9] "j_name"        "gender"        "seshu_dummy"   "jiban_seshu"  
[13] "name"          "previous"      "...15"         "age"          
[17] "exp"           "status"        "vote"          "voteshare"    
[21] "eligible"      "turnout"       "castvotes"     "...24"        
[25] "...25"         "nojiban_seshu" "mag"
  • 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))
  • exp と eligible を使って、有権者一人あたりに使う選挙費用 (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-Markdown で html 表示する際にはチャンクオプションで {r, results = “asis”} と指定する。
stargazer(as.data.frame(df1), 
          type ="html",
          digits = 2)
Statistic N Mean St. Dev. Min Max
wl 1,294 0.33 0.47 0 1
seshu_dummy 1,294 0.10 0.30 0 1
inc 1,294 0.31 0.46 0 1
exppv 1,280 17.18 13.50 0.001 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.35*** 1.94***
(0.04) (0.28)
inc -0.01 -0.12
(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 -607.20
Akaike Inf. Crit. 1,222.40
Residual Std. Error 0.40 (df = 1276)
F Statistic 176.75*** (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.07404419
  • 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")

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

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.470656
  • 候補者一人あたりの選挙資金を50円使うと、その候補者の当選確率は約 90% まで上がる
predict(logit, type = "response",
        newdata = data_frame(seshu_dummy = 0.09, inc = 0.31, exppv = 50))
        1 
0.8961336
  • 例えば、世襲の現職候補者 (seshu_dummy = 1、inc = 1) が、有権者一人あたりの選挙費用として 25円使うと、その候補者の当選確率は約 84 %だと予測できる
predict(logit, type = "response",
        newdata = data_frame(seshu_dummy = 1, inc = 1, exppv = 25))
        1 
0.8269272
  • モデルのあてはまり (fit) 比較
    OLS モデル と Logit モデルの AIC を比較
AIC(ols)
[1] 1263.786
AIC(logit)
[1] 1222.399
  • Logit モデルの AIC の値 (1223) の方が OLS モデルの AIC の値 (1266)より小さい
    => この分析では Logit モデルの方がよいモデル

まとめ - 応答変数が 0, 1 の二値変数を OLS モデルと Logit モデルで分析
→　ほぼ同様の結果
- この分析では AIC が小さい Logit モデルの方がより良いモデルといえる
→　質的な応答変数の場合には非線形回帰 (MLE) を用いた方がより良い

参考文献