16. 重回帰モデルの仮定と診断
Masahiko Asano
2023-07-04
- このセクションで使う
R パッケージ一覧
library(broom)
library(car)
library(jtools)
library(patchwork)
library(stargazer)
library(tidyverse)
library(zoo)最小二乗法による回帰係数の推定
最良線形不偏推定量: BLUE
- 最良線形不偏推定量
(
BLUE: Best Linear Unbiased Estimator) とは、次の 5 つの仮定が満たされるときの最小二乗法による回帰係数の推定量のこと
- 誤差項の期待値がゼロ
- パラメータにおける線形性
- 誤差項の条件付き期待値がゼロ
- 多重共線性がない
- 誤差項の分散均一性
- ここではこの
BLUEに即して、社会科学において要求される最小二乗法による回帰分析の前提と回帰モデルの妥当性の診断方法を解説する
回帰分析の前提と妥当性の診断方法
| 回帰分析の前提 | 妥当性の診断方法 |
|---|---|
| 0. 回帰モデルは妥当なモデルである | 学問上の知識に基づいて判断 |
| 1. 誤差項の期待値がゼロ | |
| 2. 加法性と線形性 | Scatter plotで確認 |
| 3. 誤差が独立である | 慎重に交絡変数を検討する |
| 4. 誤差の分散が均一である | 残差プロットによる診断 |
| 5. 誤差が正規分布に従う | 正規QQプロットによる診断 |
| 6. 完全な多重共線性がない | 分散拡大係数 (VIF) による診断 |
・因果推論するためには仮定 1 と仮定 2
を満たすことが特に重要
・この仮定を完全に満たすことは困難 → 傾向スコアを使った回帰分析
傾向スコアが使えない場合の因果推論の方法（準実験）
- 回帰不連続デザイン
(
RDD: Regression Discontinuity Design)
- 差の差分法 (
DID: Difference-in-differences)
- 操作変数法 (
IV: Instrumental Variables)
0. 回帰モデルは妥当なモデルである
- どのような回帰モデルが「適切なのか」を判断するのはそれぞれの学問分野における知識に裏付けられて判断される
- モデルの妥当性に関しては、次のセッションを参考にして下さい：
「18. 重回帰分析 1 (単回帰と重回帰）4.3.3 適切なモデルの選び方」
「26. 重回帰分析 9（コントロール変数の選び方：バックドア基準と DAG）」
1. 誤差項の期待値がゼロ
1.1 証明
- 誤差項は観測されない
- 「誤差項の期待値がゼロである」という仮定を満たすことは難しいことではない
- ここで次のようなモデルを考える
\[Y_i = \alpha_0 + β_1X_1i +u_i\]
\(u_i\) はこのモデルにおける誤差項
「誤差項の期待値がゼロでない」場合でも、誤差項の期待値がゼロになることを示す
つまり、誤差項の期待値 \(E(u_i) = σ（シグマ）、σ≠0\) ということ
ここでは、誤差項 \(u_i\) の期待値がゼロでない場合から始めても、誤差項 \(ε_i\)の期待値はゼロとなり仮定は満たされることを示してみる
| モデル | 解説 |
|---|---|
| \(Y_i = \alpha_0 + β_1X_1i +u_i\) | 誤差項 \(u_i\) の期待値がゼロでないと仮定 |
| \(Y_i = \alpha_0 + β_1X_1i +u_i + \sigma - \sigma\) | \(\sigma=0\)なので、計算上右辺に \(\sigma\) を足して引く |
| \(Y_i = (\alpha_0 + \sigma) + β_1X_1i + (u_i - \sigma)\) | 新たな \(Y\) 切片は \((\alpha_0 + \sigma)\)、新たな誤差項は \((u_i - \sigma) になる\) |
| \(\beta_0 = (\alpha_0 + \sigma), ε_i = (u_i - \sigma)\)と定義する | |
| \(Y_i = \beta_0 + β_1X_1i + ε_i\) | 誤差項 \(u_i\) の期待値がゼロでない場合、\(Y\) 切片の値は \(\alpha_0\) から \(\beta_0 = \alpha_0 + σ\) に変わるが、傾きは \(\beta_1\) のまま変化せず |
| \(E[ε_i]=E[u_i-\sigma]\) | 再定義した誤差項 \(ε_i = u_i - \sigma\) の期待値をとる |
| \(E[ε_i]=E[u_i] - E[\sigma]\) | 期待値の加法性 \(E[X+Y] = E[X] + E[Y]\) より |
| \(E[ε_i]=\sigma - \sigma\) | \(E[u_i] = \sigma\)であり、\(E[定数] = 定数\)だから（\(\sigma\) は定数） |
| \(E[ε_i] = 0\) | 従って、誤差項 \(ε_i\)の期待値はゼロとなり、仮定は満たされた |
ポイント
・誤差項 \(u_i\)
の期待値がゼロでない場合、\(Y\) 切片の値は \(\alpha_0\) から \(\beta_0 = \alpha_0 + σ\)
に変わるが、傾きは \(\beta_1\)
のまま変化しない
・新たに定義し直した誤差項 \(ε_i = u_i - \sigma\)
の期待値はゼロと見なすことができる
1.2 R を使った誤差項のシミュレーション
- 人工的にデータを作り、誤差項の期待値が 0 の場合と、誤差項の期待値が 0 でない場合に、それぞれ正しくモデルの傾きを推定できるかどうか確かめてみる
set.seed(20220225) # seed の設定（数字は任意）
n <- 1000 # 生成するデータ数を n = 1000 と設定
a0 <- 1 # モデルの切片 (α0)を 1 に設定
b1 <- 2 # モデルの傾き (β1) を 2 に設定
x1 <- rnorm(n) # x1 の値をランダムに 1000個生成1.2.1 誤差項の期待値（平均値）が 0 の場合
u1 <- rnorm(n, 0, 1) # 誤差項の期待値を 0 と設定
# 標準偏差を 1 と設定- 生成した誤差項 \(u1\) を可視化する
hist(u1)
- \(u1\) の記述統計を確認
summary(u1) Min. 1st Qu. Median Mean 3rd Qu. Max.
-2.93364 -0.62634 0.02522 0.03083 0.68300 3.02237 当初指定したとおり、誤差項 \(u1\) の期待値 \(E(u1)\) はほぼゼロ: 0.03083
応答変数 \(y1\) に説明変数 \(x1\) を回帰させて結果を表示させる
y1 <- a0 + b1*x1 + u1 # 応答変数 (y1) を定義する
model_1 <- lm(y1 ~ x1)1.2.2 誤差項の期待値（平均値）が 0 でない場合
u2 <- rnorm(n, 10, 1) # 誤差項の期待値を 10 と設定
# 標準偏差を 1 と設定- 生成した誤差項 \(u2\) を可視化する
hist(u2)
- \(u2\) 記述統計を確認
summary(u2) Min. 1st Qu. Median Mean 3rd Qu. Max.
6.688 9.267 9.986 9.976 10.644 12.987 当初指定したとおり、誤差項 \(u2\) の期待値 \(E(u2)\) はほぼ 10: 9.976
応答変数 \(y2\) に説明変数 \(x1\) を回帰させて結果を表示させる
y2 <- a0 + b1*x1 + u2 # 応答変数 (y2) を定義する
model_2 <- lm(y2 ~ x1)- 分析結果を表示させる
stargazer(model_1, model_2,
type = "html",
digits = 2)| Dependent variable: | ||
| y1 | y2 | |
| (1) | (2) | |
| x1 | 1.98*** | 1.96*** |
| (0.03) | (0.03) | |
| Constant | 1.03*** | 10.97*** |
| (0.03) | (0.03) | |
| Observations | 1,000 | 1,000 |
| R2 | 0.80 | 0.79 |
| Adjusted R2 | 0.80 | 0.79 |
| Residual Std. Error (df = 998) | 1.00 | 1.03 |
| F Statistic (df = 1; 998) | 3,930.37*** | 3,649.38*** |
| Note: | p<0.1; p<0.05; p<0.01 | |
傾きはほとんど同じ: \(1.98\) と \(1.96\)
切片は異なる: \(1.03\) と \(10.97\)
二つのモデルを可視化してみる
シミュレーションの結果
・「誤差項の期待値がゼロでない」モデル \(y_2\) の場合でも、影響を受けるのは \(Y\) 切片 (1.03→10.97) だけであり、傾き
\(\beta_i = 1.98\)
にはほとんど影響はない
・誤差項 \(ε_i\)
の期待値はゼロと見なすことができる
2. 加法性と線形性
2.1 線形性とは
回帰モデルにおけるパラメータ（母数）が線形という仮定
「線形回帰モデルでは、誤差を除く応答変数の決定要因が各説明変数の線形関数で表されなければならない」という仮定
この仮定は、最小二乗法による回帰係数の推定量が「不偏」であるために必要
応答変数を \(y\)、説明変数を \(x_m (m = 1, 2, ..., k )\)とすると、次の関係があるということ
\[y = β_1x_1 + β_2x_2 +...+ β_kx_k\]
- 切片 \(\beta_0\)
は説明変数とは関係ないので決定要因には含めない
- 切片を含めた \(y = β_0 + β_1x_1 + β_2x_2 + ... + β_kx_k(β_0 ≠ 0)\) は一次関数ではあるが線形関数ではない
加法性が満たされない場合
例 1:
\[y = x_1x_2x_3\]
- この式は加法性を満たさない
★ 解決策 1 → 両辺の対数をとる
\[logy=log(x_1x_2x_3)=logx_1 +logx_2 +logx_3\]
- ここで \(Y=log_y, X_j = logx_j(j = 1, 2,3)\) と新たに定義すると
\[Y = X_1 + X_2 + X_3\]
- 加法性を満たす式となり、説明変数と応答変数の関係をモデル化できる
★ 解決策 2 → \(z = x_1x_2x_3\) と定義
\[y = z\]
- \(y=z\) という線形モデルができるので線形回帰として推定可能
例 2:
- 次のような関係を考えたい
\[y=logx\]
★ 解決策 → 新たに \(t=logx\) という変数を作る
log関数で変数変換し、\(y=t\) という関係を考えることができる
- 線形回帰として推定可能
2.2 変数変換の実例（対数変換）
- ここでは対数変換を使った分析方法と結果の結果の解釈方法を紹介する
- ここでは I と IV のモデルを比較して解説する
- モデル I は説明変数と応答変数のどちらも対数変換しない場合
- モデル IV は説明変数と応答変数の両方を対数変換した場合
| モデル | 解釈 |
|---|---|
| I. \(\hat{y}_i = \hat{\beta_0} + \hat{\beta_1}x_{1i}\) | \(x_1\)が1単位変化すると、\(y\) は\(\hat{\beta_1}単位変化する\) |
| II. \(\hat{y}_i = \hat{\beta_0} + \hat{\beta_1}log(x_{1i})\) | \(x_1\)が1%変化すると、\(y\) は\(\hat{\beta_1}/100単位変化する\) |
| III. \(\widehat{log(y_i)} = \hat{\beta_0} + \hat{\beta_1}x_{1i}\) | \(x_1\)が1単位変化すると、\(y\) は100\(\hat{\beta_1}\)%変化する |
| IV. \(\widehat{log(y_i)} = \hat{\beta_0} + \hat{\beta_1}log(x_{1i})\) | \(x_1\)が1%変化すると、\(y\) は\(\hat{\beta_1}\)%変化する |
gapminderパッケージに含まれている「GDP」と「寿命」変数を使う
データの準備 (gapminder)
gapminderパッケージをダウンロードする
library(gapminder)gapminderには次の変数が含まれている
names(gapminder)[1] "country" "continent" "year" "lifeExp" "pop" "gdpPercap"- 必要な変数だけを選ぶ:
lifeExp,gdpPercap
gapminder_1 <- gapminder |>
dplyr::select(year, lifeExp, gdpPercap) - 対数変換を施した新たな変数 (
log_gdp,log_lifeExp) を作り、変数を限定する
gapminder_1 <- gapminder_1 |>
mutate(log_gdpPercap = log(gdpPercap)) |>
mutate(log_lifeExp = log(lifeExp)) |>
round(digits = 2) |> # 表示を 2 桁に指定
dplyr::filter(year == 2007) # 2007年のデータだけを残すgapminderの記述統計を示す
library(stargazer)| 変数名 | 詳細 |
|---|---|
| lifeExp | 寿命 |
| log_lifeExp | 寿命（対数変換） |
| gdpPercap | 1 人あたり GDP, 2005年の米ドル基準 |
| log_gdpPercap | 1 人あたり GDP（対数変化） |
- チャンクオプションで
{r, results = "asis"}と指定する
stargazer(as.data.frame(gapminder_1),
type ="html",
digits = 2)| Statistic | N | Mean | St. Dev. | Min | Max |
| year | 142 | 2,007.00 | 0.00 | 2,007 | 2,007 |
| lifeExp | 142 | 67.01 | 12.07 | 39.61 | 82.60 |
| gdpPercap | 142 | 11,680.07 | 12,859.94 | 277.55 | 49,357.19 |
| log_gdpPercap | 142 | 8.62 | 1.36 | 5.63 | 10.81 |
| log_lifeExp | 142 | 4.19 | 0.20 | 3.68 | 4.41 |
gapminder の中身を確認する
DT::datatable(gapminder_1)対数変数しない変数を使った散布図:
gdpPercap & lifeExp
- まず、「1 人あたり GDP (
gdpPercap)」を x 軸、「平均寿命 (lifeExp)」を y 軸に指定して散布図を描いてみる
gapminder_1 |>
ggplot() +
geom_point(aes(x = gdpPercap,
y = lifeExp)) +
labs(x = "1 人あたりGDP (USD)",
y = "寿命") +
theme_bw(base_family = "HiraKakuProN-W3")
- リニアな関係とは言えないが、回帰分析をしてみる
model_1 <- lm(lifeExp ~ gdpPercap, data = gapminder_1)- 結果は、
model_4と一緒に後に示す
対数変数した変数を使った散布図:
log_gdpPercap & log_lifeExp
- 「1 人あたり
GDP:gdpPercap」の記述統計を見ると、データの範囲が 277米ドル（最低額）から 49,359 米ドル（最高額）まで非常に広範囲に及ぶ
- 「1 人あたり
GDP」と「平均寿命」はリニアというより、対数関係のような関係
→ 見やすくするため、x 軸を対数変換した変数 (log_gdpPercap) と「平均寿命」(log_lifeExp) との散布図を描いてみる
gapminder_1 |>
ggplot() +
geom_point(aes(x = log_gdpPercap,
y = log_lifeExp)) +
labs(x = "1 人あたりGDP (USD)の対数値",
y = "寿命（対数変換）") +
theme_bw(base_family = "HiraKakuProN-W3")
- 変数を対数化したら、かなり綺麗な線形関係が確認できた
- 回帰分析をしてみる
model_2 <- lm(log_lifeExp ~ log_gdpPercap, data = gapminder_1)model_1とmodel_2の分析結果を表示させる
stargazer(model_1, model_2, type = "html")| Dependent variable: | ||
| lifeExp | log_lifeExp | |
| (1) | (2) | |
| gdpPercap | 0.001*** | |
| (0.0001) | ||
| log_gdpPercap | 0.113*** | |
| (0.008) | ||
| Constant | 59.566*** | 3.211*** |
| (1.010) | (0.067) | |
| Observations | 142 | 142 |
| R2 | 0.461 | 0.607 |
| Adjusted R2 | 0.457 | 0.605 |
| Residual Std. Error (df = 140) | 8.898 | 0.124 |
| F Statistic (df = 1; 140) | 119.556*** | 216.516*** |
| Note: | p<0.1; p<0.05; p<0.01 | |
対数変換した回帰分析 (model_4)
の解釈
Model_1:
・一人当たりGDPが 1ドル増えると、寿命が 0.001歳増える
→ 一人当たりGDPが 1000ドル増えると、寿命が 1 歳増える
Model_2:
・一人当たりGDPが 1% 増えると、寿命が 0.113%ポイント増える→ 一人当たりGDPが 10 % 増えると、寿命が約 1 %ポイント増える
3. 誤差が独立である
- ここで、各個体 \(i\) について次の関係を想定してみる
\[Y_i = \beta_0 + \beta_1x_{1i} + .... + \beta_kx_{ki} + ε_i\]
- 誤差項 \(ε_i\)：母集団の回帰直線（予測値）と個体
\(i\) の応答変数（観測値）のズレ
- 誤差項 \(ε_j\)：母集団の回帰直線（予測値）と個体 \(j\) の応答変数（観測値）のズレ
「誤差が独立」というのは「別々の個体 \(i\) と \(j\) の誤差の間に何の関係もない」という意味
上に書いた回帰モデルでは、応答変数のうち説明変数に含まれていない変数の影響は全て誤差項に含まれる
従って、誤差が独立であるというのは「説明変数に含まれるもの以外で応答変数に影響する要因が互いに独立である（＝無関係である）」ということ
「誤差が独立である」ことを式で表現すると
\[Cov (ε_i, ε_j) = 0\]
\[i ≠ j であるようなすべての (i , j) について、ε_i と ε_j の共分散が 0 になる\]
これは「系列相関がない」と同じ意味
例えば、個人のビールの消費量 \(Beer\) が個人の収入 \(Income\) によって影響を受けるという次のようなモデルを考える
\[Beer_i = β_0 - β_1Income_i + ε_i\]
「誤差が独立である」＝「2 つの個体、\(i\) と \(j\) の誤差の間に何の関係もない」
上記の回帰モデルでは、応答変数 \(Beer\) のうち説明変数 \(Income\) に含まれていない変数の影響は誤差項 \(ε_i\)に含まれる
従って、「誤差が独立である」＝「説明変数に含まれるもの以外で、応答変数に影響する要因が互いに独立である（＝無関係である）」
誤差の独立性が満たされない一例：
☆ \(i\) さんと \(j\)
さんはゼミの飲み友達で、一緒にビールを飲みに行く
→ \(i\)
さんのビール消費量が増えると、\(j\)
さんのビール消費量も増える
- これが事実だとすると \(Beer\) のうちで収入 \(Income\) では説明できない部分が \(Beer\) と正の相関を持つことになり、誤差の独立性は満たされないことになる
「平均独立」と「条件付き平均独立」
誤差項の条件付き期待値がゼロということ
・誤差項と共変量 \(X\) は独立であると仮定する・つまり \(E[ε_i|X] = E[ε_i]\)
・これが満たされているとき「誤差項 \(ε_i\) は共変量 \(X\) と平均独立 (
mean independent) であるという・また、誤差項 \(ε_i\) の期待値がゼロなので \(E[ε_i]=0\)
・共変量 \(X\) が与えられたとき、誤差項 \(ε_i\) の期待値はゼロ
・このとき、誤差項 \(ε_i\) は共変量 \(X\) と条件付き平均独立 (
conditional mean independent) であるという\[E[ε_i|X] = 0\]
＊「平均独立」の詳細に関しては「16. ランダム化比較試験:RCT」を参照
「誤差が独立」であるという仮定は最小二乗法による回帰係数の推定量が不偏であるために必要
→ 「統制すべき交絡因子が十分にモデルに含まれていなければならない」ということにつながるどのような変数を交絡因子としてモデルに含めるべきか（含めるべきではないか）という問題は次のセッションを参照：
「26. 重回帰分析 9（コントロール変数の選び方：バックドア基準と DAG）」不要な変数をモデルに含めた場合と、必要な変数をモデルに含めなかった場合にどういうことが起きるかという問題は次のセッションを参照：
「25. 重回帰分析（脱落変数バイアス・処置後変数バイアス）」
4. 誤差の分散が均一である
4.1 「誤差の分散が均一である」とは？
- 「誤差の分散が同じ」とは「全ての個体の誤差の分散が同じ」になるということ
\[全ての個体 i について、Var(ε_i) = σ^2\]
- 例えば、気温とビールの例だと、気温が高い夏でも寒い冬でも、ビールの出荷量のばらつきはなく一定であるということ
- この事例のように分散が均一であることを均一分散性と呼ぶ
- 分散が均一でないことを不均一分散性と呼ぶ
homoskerdasticity |
均一分散 | 分散が均一であること |
heteroskedasticity |
不均一分散 | 分散が均一でないこと |
- 気温とビールの事例を図で表すと、左側のようになる
- ビールの売り上げ額にかかわらず、売り上げのバラツキが一定
均一分散（左側の図）と不均一性分散（右側の図）の事例
- 不均一分散 (
heteroskerdasticity) の事例としては「月収」と「食費」が考えられる
- 個人が 1 か月に使う「食費」を応答変数、その個人の「月収」を説明変数
とするモデルを考える
- この国で生きていくためには最低 8,000 円は
毎月の食費に使う必要があるとする
- 月収が 1 万円の個人が食費にかける金額は 8,000 円以上 1
万円以下
- これに対し、月収が 100 万円の個人が食費にかける金額は、8,000 円以上
100万円以下
- つまり、月収が少
ない場合、食費にかける金額を選択する余地が小さい
→ 応答変数の値はあまりばらつかない
- 月収が増えれば選択の余地が広がり、食費を節約する人もいるし贅沢な食事を楽しむ人もいる
→ ばらつきが大きくなる
- 応答変数のばらつき自体は誤差のばらつきではないが、応答変数のばらつきが大きいほど説明変数で説明できない部分が大きくなる
→ 誤差のばらつきも大きくなる
社会科学の事例における説明変数と誤差の分散
応答変数のばらつきが大きくなるほど
→ 誤差のばらつきが大きくなる
→ 説明変数で説明できない部分が大きくなる
→ 誤差に分散不均一性が生じる
→ 標準誤差に影響がでて、統計的な有意性が得にくくなる
・ 誤差項の分散が不均一でも、最小二乗法による回帰係数の推定量に偏りは発生しない
4.2 不均一分散の診断
4.2.1 パットナムのモデルでの不均一診断
Breusch-Pagan Testによって、誤差項の分散が均一かどうか検定できる政府のパフォーマンスに関するデータ (putnam.csv) をダウンロードし、
RProjectフォルダー内のdataフォルダに保存するデータのダウンロードが終わったら、データを読み込み
putnamと名前を付ける
putnam <- read_csv("data/putnam.csv")gov_pを応答変数とした単回帰分析を行う
fit.putnam <- lm(gov_p ~ cc, data = putnam) - 誤差 = 残差 (
residuals) をプロットする
res.putnam <- resid(fit.putnam) # 残差を取り出す
pred.putnam <- fitted(fit.putnam) # 予測値を取り出す
plot(pred.putnam,
res.putnam,
main = "Residuals vs Fitted",
xlab = "Fitted Values",
ylab = "residuals")
abline(h = 0, col = "red") # 残差 = 0 に平行な赤線を引く
応答変数
Fitted Valuesの大小によって、誤差項が偏っている（= 不均一分散）ようには見えないRパッケージ
lmtestのbytest()関数を使った不均一分散の診断
library(lmtest)bptest(fit.putnam)
studentized Breusch-Pagan test
data: fit.putnam
BP = 0.36787, df = 1, p-value = 0.5442bytest()関数の診断基準
p-valueが 0.05以上の場合 |
誤差項の分散が均一である | |
p-valueが 0.05以下の場合 |
誤差項の分散が不均一である | |
p-value = 0.5442で 0.05 以上なので、「誤差項の分散が均一である」
4.2.2 2009年衆院選のモデルでの不均一診断
応答変数が得票率の場合
Breusch-Pagan Testによって、誤差項の分散が均一かどうか検定できるデータのダウンロード方法
RProjectフォルダ内に data という名称のフォルダを作成するhr96-21.csv をクリックしてデータをパソコンにダウンロード
RProjectフォルダ内に data という名称のフォルダを作成する
ダウンロードした
hr96-21.csvを手動でRProjectフォルダ内にある data フォルダに入れる選挙データの読み取り方法
次のいずれかの方法で
hr96-21.csvを読み取る
読み取り方法 1
na = "."というコマンドは「欠損値をドットで置き換える」という意味
- 欠損値を空欄のまま残すと、本来「数値
(
numeric)」型のデータが「」文字型 (character)」として認識されるなど、エラーの原因になるため、読み取る時点で事前に対処する
hr <- read_csv("data/hr96-21.csv",
na = ".") 読み取り方法 2
- 読み取った値の日本語が文字化けする場合
locale()関数を使って日本語エンコーディング (cp932) を指定する
hr <- read.csv("data/hr96-21.csv",
na = ".",
locale = locale(encoding = "cp932"))読み取り方法 3
hr <- read.csv("data/hr96-21.csv",
na = ".") - 読み取ったデータの変数を確認
names(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"- 2009年の総選挙データだけを抜き出す
hr09 <- hr |>
dplyr::filter(year == 2009) |>
dplyr::select(voteshare, age)voteshareを応答変数とした単回帰分析を行う
fit.hr09 <- lm(voteshare ~ age, data = hr09) - まず、回帰分析から残差 (
residuals) のプロットを描いてみる
fit.hr09 <- lm(voteshare ~ age, data = hr09)
res.hr09 <- resid(fit.hr09) # 残差を取り出す
pred.hr09 <- fitted(fit.hr09) # 予測値を取り出す
plot(pred.hr09,
res.hr09,
main = "Residuals vs Fitted",
xlab = "Fitted Values",
ylab = "residuals")
abline(h = 0, col = "red") # 残差 = 0 に平行な赤線を引く
Fitted Valuesが大きくなるにつれて、Risiduals（赤の平行点線）を中心に下の方がデータが広がっているので、誤差の独立性は保たれていないようにみえるRパッケージ
lmtestのbytest()関数を使った不均一分散の診断
library(lmtest)bptest(fit.hr09)
studentized Breusch-Pagan test
data: fit.hr09
BP = 13.982, df = 1, p-value = 0.0001846bytest()関数の診断基準
p-valueが 0.05以上の場合 |
誤差項の分散が均一である | |
p-valueが 0.05以下の場合 |
誤差項の分散が不均一である | |
p-value = 0.00509で 0.05 以下なので、「誤差項の分散が不均一である」
誤差の分散が不均一の時の 2 つの対処法:
| 1. 加重最小二乗法 | 重み付けたOLS | 不均一にしている変数とその関数形が分かる場合 |
| 2. 不均一分散に頑健な標準誤差 | 標準誤差を修正 | 不均一にしている変数とその関数形が分からない場合 |
1. 加重最小二乗法 (WLS: Weighted least squares)
- 各観測点の分散の逆数(これを精度
precisionと呼ぶ)に比例し た重みを使った「重み付き最小二乗法」の推定量を使う
- この推定法法の方が、最小二乗法による推定量より「効率的」
-「効率的」= 「標準誤差が小さい」という意味
2. 不均一分散に頑健な標準誤差
- もう一つの解決策として標準誤差を修正することで、不均一分散に対応する「不均一分散に頑健な標準誤差」を使う方法がある
- 誤差の分散が不均一の時
→「重み付け最小二乗法」や「不均一分散に頑健な標準誤差」を使うのが望ましいが、分散の不均一性は、回帰分析によって推定される説明変数と応答変数の関係には影響しないので、あまり重大な問題とは言えない
5. 誤差が正規分布に従う
- このセッションでは回帰分析の前提として次の 2 つの仮定を紹介した
2. 誤差項の期待値がゼロ
6. 誤差の分散が均一である
- 誤差項の期待値はゼロであり、誤差の分散は均一なので、分散は定数 \(σ\) と表せる
- 「誤差の正規性」とは「誤差が正規分布に従う確率変数であること」
- 各個体の誤差が平均 0、分散 \(σ^2\)
の正規分布に従うということ
- すべての個体 \(i\) について次の式が成り立つということ
\[ε_i 〜 N (0, σ^2)\]
- 誤差項 \(ε_i\)に様々な要因が含まれていることから
→ 中心極限定理により、母集団がどのような分布であれ誤差の分布は正規分布で近似できる
→ 「誤差の正規性」仮定を正当化できる
→ 回帰直線の推定では「誤差の正規性」仮定は気にする必要はない
残差が正規分布に従っていないケース（左）と正規分布に従っているケース（右）の事例
- 最良線形不偏推定量
(
BLUE: Best Linear Unbiased Estimator)である 5 つの仮定の中に「誤差が正規分布に従う」という仮定は含まれていない
- 誤差項の期待値がゼロ
- パラメータにおける線形性
- 誤差項の条件付き期待値がゼロ
- 多重共線性がない
- 誤差項の分散均一性
- 誤差項が正規分布に従っていなくても、通常の最小二乗法による回帰係数の推定量は最良線形不偏推定量
(
BLUE) である
複数のプロット表示
- R では
par( )関数を使って 2 行 2 列のグラフィック環境を整えてから、plot( )関数を使って、4 つのプロットを並べて表示することもできる
fit.putnam <- lm(gov_p ~ cc + econ, data = putnam) # gov_p を応答変数とした重回帰分析を行う par(mfrow = c(2,2)) # 2 行 2 列のグラフィック環境を作る plot(fit.putnam) # 4 つのプロットを表示
6. 完全な多重共線性がない
6.1 多重共線性とは
多重共線性 (
multicollinearity) の定義：
「説明変数の間の相関係数が 1 ではないものの、強い相関があること」次のようなモデルを考える
\[Y_i = \beta_0 + \beta_1X_{1i} + \beta_2X_{2i} + ε_i\]
- 応答変数が \(Y\)
- 説明変数が \(X_1\) と \(X_2\)
- 通常の回帰分析では左図ように、説明変数である \(X_1\) と \(X_2\) の間にはある程度の相関があるが、相関係数は 1 ではない
出典：高橋(2022)p.107を参考に著者が作成した図
- 多重共線性が問題になるのは、右図のように、説明変数である \(X_1\) と \(X_2\) の相関係数が 1 の時
- このような場合の解決策：
→ \(X_1\) と \(X_2\)、どちらかの説明変数だけを使う
6.2 なぜ、多重共線性が問題なのか？
- 重回帰分析で多重共線性は禁止されているわけではない
- 説明変数の間に相関があるのは当然のこと
- 相関係数の大きさの「程度問題」
- 相関係数が大きすぎることで起こること：
→ 標準誤差が大きくなる
→ 信頼区間の幅が大きくなる
→ 帰無仮説を棄却しにくくなる
→ 統計的に有意になりにくくなる
- しかし、回帰係数の推定は不偏に行うことが可能
→ 極端に大きな相関でない限り、分析する上で気にする必要はない
- では、どの程度なら許容できるのか？
6.3 多重共線性の診断方法: VIF
- 通常、私たちは多変量解析を行う
→ 説明変数のペアごとの相関ではなく、説明変数全体での相関を確かめる必要がある
- そのために、多重共線性は分散拡大要因
(
VIF: variance inflation factor) という指標を使って診断する
VIF の計算方法
- 次のような説明変数の数は \(p\) 個のモデルを考える
\[\hat{Y} = \hat{\beta_0} + \hat{\beta_1X_{1i}} + \hat{\beta_2X_{2i}} +...+ \hat{\beta_pX_{pi}}\]
- \(p\) 個の説明変数 \(X\) のうち、\(j\) 番目の値 \(X_{ji}\) を応答変数として、残りの \(X\) を右辺に持ってきた時の重回帰モデルは次のとおり
\[\hat{X_{ji}} = \hat{γ_0} + \hatγX_{1i} ... + \hatγX_ {pi} \]
- この回帰式の決定係数が \(R^2_j\)
\[VIF = \frac{1}{1-R^2_j}\]
6.4 共変量における多重共線性
6.4.1 イタリア地方政府に関するモデルの VIF
Rパッケージ carを使ってモデルのVIFを計算できる
putnam <- read_csv("data/putnam.csv")
fit.putnam <- lm(gov_p ~ cc + econ,
data = putnam) # gov_p を応答変数とした重回帰分析を行うlibrary("car")
vif(fit.putnam) cc econ
5.632172 5.632172 VIF > 10なら、多重共線性の疑いがある
VIF > 10というのは説明変数間の相関係数が0.9を超えるような場合
ccとeconの VIF は 5.6 程度なので問題ない
VIF > 10の時の対処法
→ 相関の強い 2 つの説明変数のいずれか片方をモデルに含めずに分析
6.4.2 日本の2009年衆院選挙データにおける VIF
- 2009年衆院選データにおける 4 つの変数を使って回帰分析してみる
- ここで知りたいのは候補者の年齢
ageが小選挙区でのランクrankに与える影響
- 応答変数：小選挙区でのランク
rank
- 説明変数：候補者の年齢
age
- 候補者の年齢
ageが最も興味のある重要な変数
| 変数名 | 詳細 | |
|---|---|---|
rank |
小選挙区でのランク | 応答変数 |
age |
候補者の年齢 | 説明変数 |
vote |
候補者が得た票数 | 共変量 |
voteshare |
候補者が得た得票率 (%) | 共変量 |
- 共変量として「候補者が得た票数
vote」と「候補者が得た得票率voteshare」をモデルに含める
voteもvoteshareも「交絡因子」ではない
- 票数や得票率が年齢に影響を与えるとは考えられない
- しかし、どちらの変数も「ランク
rank」に影響を与えるし、年齢とは「相関」があると思われるため便宜上モデルに含める
- 年齢と票数との相関係数は −0.19（かなり小さい）
- 票数と得票率との相関係数は 0.96（かなり大きい）
- ここで使う変数の関係をベン図で示してみる
- 総選挙データを読み取る
hr <- read_csv("data/hr96-21.csv")- 必要な変数だけに絞る
hr09 <- hr |>
dplyr::filter(year == 2009) |> # 2009年総選挙に限る
dplyr::select(age, vote, voteshare, rank) # 必要な変数を選ぶ変数の記述統計
stargazer(as.data.frame(hr09),
type = "html",
digits = 2)| Statistic | N | Mean | St. Dev. | Min | Max |
| vote | 1,139 | 61,939.59 | 53,880.34 | 177 | 201,461 |
| voteshare | 1,139 | 26.34 | 22.20 | 0.10 | 95.30 |
| rank | 1,139 | 2.50 | 1.24 | 1 | 9 |
「年齢」と「ランク」の散布図を表示させる
- ここで知りたいのはこの 2 つの変数の関係
hr09 |>
ggplot(aes(age, as.factor(rank))) +
geom_point() +
stat_smooth(method = lm, se = FALSE) +
theme_bw()
- 負の相関がありそう
- 年配の候補者ほどランクが良くなる（= 1位に近い）
「年齢」と「票数」の散布図を表示させる
hr09 |>
ggplot(aes(age, vote)) +
geom_point() +
stat_smooth(method = lm, se = FALSE) +
theme_bw()
- 年齢が増えるほど票が増える
- ほとんど相関はない: \(r = 0.19\)
hr09$age <- as.numeric(hr09$age)
hr09$vote <- as.numeric(hr09$vote)with(hr09, cor.test(age, vote))
Pearson's product-moment correlation
data: age and vote
t = 6.5822, df = 1133, p-value = 7.072e-11
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.1352367 0.2473405
sample estimates:
cor
0.1919146 「票数」と「得票率」の散布図を表示させる
hr09 |>
ggplot(aes(vote, voteshare)) +
geom_point() +
stat_smooth(method = lm, se = FALSE) +
theme_bw()
- 強い相関がある: \(r = 0.96\)
with(hr09, cor.test(vote, voteshare))
Pearson's product-moment correlation
data: vote and voteshare
t = 117.15, df = 1137, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.9562771 0.9651931
sample estimates:
cor
0.960984 モデルの推定
- ここでは次の 4 つのモデルを推定する
m1 <- lm(rank ~ age, data = hr09)
m2 <- lm(rank ~ age + vote, data = hr09)
m3 <- lm(rank ~ age + voteshare, data = hr09)
m4 <- lm(rank ~ age + vote + voteshare, data = hr09)stargazer(m1, m2, m3, m4,
type = "html")| Dependent variable: | ||||
| rank | ||||
| (1) | (2) | (3) | (4) | |
| age | -0.021*** | -0.003 | -0.001 | -0.001 |
| (0.003) | (0.002) | (0.002) | (0.002) | |
| vote | -0.00002*** | -0.00000 | ||
| (0.00000) | (0.00000) | |||
| voteshare | -0.050*** | -0.046*** | ||
| (0.001) | (0.003) | |||
| Constant | 3.538*** | 3.847*** | 3.845*** | 3.847*** |
| (0.166) | (0.086) | (0.077) | (0.077) | |
| Observations | 1,135 | 1,135 | 1,135 | 1,135 |
| R2 | 0.035 | 0.743 | 0.793 | 0.794 |
| Adjusted R2 | 0.034 | 0.743 | 0.793 | 0.793 |
| Residual Std. Error | 1.221 (df = 1133) | 0.630 (df = 1132) | 0.566 (df = 1132) | 0.566 (df = 1131) |
| F Statistic | 41.288*** (df = 1; 1133) | 1,637.589*** (df = 2; 1132) | 2,168.995*** (df = 2; 1132) | 1,448.713*** (df = 3; 1131) |
| Note: | p<0.1; p<0.05; p<0.01 | |||
分析結果を解釈する上で大事なこと
- 私たちの主要な関心は「候補者の年齢
age」が候補者の「ランク」に与える影響
- 共変数である「候補者が得た票数
vote」と「候補者が得た得票率voteshare」の偏回帰係数は気にする必要なし
- 共変量は交絡調整のためにモデルに含められるため、これらの偏回帰係数は多重共線の影響を受け、正しく推定できない
- 共変量の偏回帰係数は、主要な推定対象ではないので無視してかまわない
- 注目すべきは「候補者の年齢
age」偏回帰係数
Model_1:
- 候補者の年齢
ageが 1 歳増えると、候補者のランクが 0.02 ポイントだけ良くなる（= 1 位に近くなる）
- しかし、
Model_1には統制変数が含まれていないので、この結果は候補者の年齢ageを過大評価していると思われる
→ 統制変数（共変量）を含むモデル (Model_2&Model_3) を推定する必要あり
Model_2 & Model_3
- 「候補者の票」
voteをモデルに含めると、候補者の年齢ageの影響力は消える - 「候補者の得票率」
voteshareをモデルに含めると、候補者の年齢ageの影響力は消える
→Model_1の結果は候補者の年齢ageを課題評価していることが裏付けられた
Model_4
- 「候補者の票」
voteと「候補者の得票率」voteshareを両方同時にモデルに含めても、候補者の年齢ageの影響力は消え、しかも、Model_2, 3, 4と、候補者の年齢ageの係数が変わらない
→ 「ランクに対する候補者の年齢ageの影響力はない」という推定が正しそう
Model_4 の VIF を計算する
carパッケージを使って、モデル4のVIFを計算する
library(car)
vif(m4) age vote voteshare
1.043137 13.127862 13.187917 多重共線のポイント
・ 回帰モデルが含む 2 つの共変量 vote と
voteshare の VIF が 13 以上
・共変量の間に多重共線が起きている
→ VIF が 10 以上ではあるが、これら
2 つの変数をモデルに加えても「候補者の年齢
age」は正しく推定できているので問題なし
・共変量の間に多重共線が起きていても、特に問題視する必要はない
・むしろ必要な共変量をモデルに含めないと、交絡を取り除けない
・必要と思われる共変量はモデルに含めるべき
7. 仮定が満たされない場合の対処法
7.1 誤差が均一でない場合
解決策：
- 誤差項の不均一分散に対して頑健な標準誤差（
heteroskedasticity robust standard error）を考慮した回帰分析を行う
- 不均一分散に対して頑健な標準誤差（
heteroskedasticity robust standard error）を考慮した回帰分析を行う
7.2 自己相関の問題がある場合
解決策：
- 固定効果 (fixed effect) をモデルに入れる
- パネルデータで観察されない個人や時間の特徴が従属変数に影響を与える場合
- 各時点で同じ個人の効果が効き続ける
→ この効果を捕捉しないと、推定誤差に系列相関が生まれてしまう - 簡単なのは
lm()の中で、個人や時点などをダミー変数としてモデルに含める
事例 1：GDPで寿命を予測する (1952-2007)
gapminderパッケージに含まれている「GDP」と「寿命」変数を使う
データの準備 (gapminder)
gapminderパッケージをダウンロードする
library(gapminder)gapminderには次の変数が含まれている
names(gapminder)[1] "country" "continent" "year" "lifeExp" "pop" "gdpPercap"- 必要な変数だけを選ぶ:
lifeExp,gdpPercap
gapminder <- gapminder |>
dplyr::select(year, lifeExp, gdpPercap) - 対数変換を施した新たな変数 (
log_gdp,log_lifeExp) を作り、変数を限定する
gapminder_2 <- gapminder |>
mutate(log_gdpPercap = log(gdpPercap)) |>
mutate(log_lifeExp = log(lifeExp)) |>
round(digits = 2) # 表示を 2 桁に指定gapminderの記述統計を示す
library(stargazer)| 変数名 | 詳細 |
|---|---|
| year | 年度 |
| lifeExp | 寿命 |
| log_lifeExp | 寿命（対数変換） |
| gdpPercap | 1 人あたり GDP, 2005年の米ドル基準 |
| log_gdpPercap | 1 人あたり GDP（対数変化） |
- チャンクオプションで
{r, results = "asis"}と指定する
stargazer(as.data.frame(gapminder_2),
type ="html",
digits = 2)| Statistic | N | Mean | St. Dev. | Min | Max |
| year | 1,704 | 1,979.50 | 17.27 | 1,952 | 2,007 |
| lifeExp | 1,704 | 59.47 | 12.92 | 23.60 | 82.60 |
| gdpPercap | 1,704 | 7,215.33 | 9,857.45 | 241.17 | 113,523.10 |
| log_gdpPercap | 1,704 | 8.16 | 1.24 | 5.49 | 11.64 |
| log_lifeExp | 1,704 | 4.06 | 0.23 | 3.16 | 4.41 |
gapminder の中身を確認する
DT::datatable(gapminder_2)対数変数しない変数を使った散布図:
gdpPercap & lifeExp
- まず、「1 人あたり GDP (
gdpPercap)」を x 軸、「平均寿命 (lifeExp)」を y 軸に指定して散布図を描いてみる
gapminder_2 |>
ggplot() +
geom_point(aes(x = gdpPercap,
y = lifeExp)) +
labs(x = "1 人あたりGDP (USD)",
y = "寿命") +
theme_bw(base_family = "HiraKakuProN-W3")
- リニアな関係ではないので、対数変数したデータを使ってみる
対数変数した変数を使った散布図:
log_gdpPercap & log_lifeExp
- 「1 人あたり
GDP:gdpPercap」の記述統計を見ると、データの範囲が 241米ドル（最低額）から 113,523 米ドル（最高額）まで非常に広範囲に及ぶ
- 「1 人あたり
GDP」と「平均寿命」はリニアというより、対数関係のような関係
→ 見やすくするため、x 軸を対数変換した変数 (log_gdpPercap) と「平均寿命」(log_lifeExp) との散布図を描いてみる
gapminder_2 |>
ggplot() +
geom_point(aes(x = log_gdpPercap,
y = log_lifeExp)) +
labs(x = "1 人あたりGDP (USD)の対数値",
y = "寿命（対数変換）") +
theme_bw(base_family = "HiraKakuProN-W3")
- 変数を対数化したら、かなり綺麗な線形関係が確認できた
- 回帰分析をしてみる
model_22 <- lm(log_lifeExp ~ log_gdpPercap, data = gapminder_2)model_22の分析結果を表示させる
stargazer(model_22, type = "html")| Dependent variable: | |
| log_lifeExp | |
| log_gdpPercap | 0.147*** |
| (0.003) | |
| Constant | 2.864*** |
| (0.023) | |
| Observations | 1,704 |
| R2 | 0.613 |
| Adjusted R2 | 0.613 |
| Residual Std. Error | 0.145 (df = 1702) |
| F Statistic | 2,698.395*** (df = 1; 1702) |
| Note: | p<0.1; p<0.05; p<0.01 |
対数変換した回帰分析
(model_22) の解釈
Model_22:
・一人当たりGDPが 1% 増えると、寿命が 0.147% 増える→ 一人当たりGDPが 10 % 増えると、寿命が約 1.47 %ポイント増える
- ここで、パネルデータで観察されない年度の特徴が従属変数に影響を与えると考える
- 各時点で同じ年度の効果が効き続ける
→ この効果を捕捉しないと、推定誤差に系列相関が生まれてしまう
- モデルに
year ダミーを含める
→factor(year)をモデルに追加する
model_33 <- lm(log_lifeExp ~ log_gdpPercap + factor(year),
data = gapminder_2)stargazer(model_22, model_33,
type = "html")| Dependent variable: | ||
| log_lifeExp | ||
| (1) | (2) | |
| log_gdpPercap | 0.147*** | 0.135*** |
| (0.003) | (0.003) | |
| factor(year)1957 | 0.033** | |
| (0.016) | ||
| factor(year)1962 | 0.058*** | |
| (0.016) | ||
| factor(year)1967 | 0.079*** | |
| (0.016) | ||
| factor(year)1972 | 0.095*** | |
| (0.016) | ||
| factor(year)1977 | 0.116*** | |
| (0.016) | ||
| factor(year)1982 | 0.144*** | |
| (0.016) | ||
| factor(year)1987 | 0.171*** | |
| (0.016) | ||
| factor(year)1992 | 0.184*** | |
| (0.016) | ||
| factor(year)1997 | 0.186*** | |
| (0.016) | ||
| factor(year)2002 | 0.185*** | |
| (0.016) | ||
| factor(year)2007 | 0.185*** | |
| (0.016) | ||
| Constant | 2.864*** | 2.841*** |
| (0.023) | (0.023) | |
| Observations | 1,704 | 1,704 |
| R2 | 0.613 | 0.684 |
| Adjusted R2 | 0.613 | 0.681 |
| Residual Std. Error | 0.145 (df = 1702) | 0.131 (df = 1691) |
| F Statistic | 2,698.395*** (df = 1; 1702) | 304.605*** (df = 12; 1691) |
| Note: | p<0.1; p<0.05; p<0.01 | |
- 変数
yearの値を確認してみる
table(gapminder_2$year)
1952 1957 1962 1967 1972 1977 1982 1987 1992 1997 2002 2007
142 142 142 142 142 142 142 142 142 142 142 142 - 1952年から2007年まで 12 個ある
- 結果には
factor(year)1957からfactor(year)2007まで 11 個しかない
→ 非表示の 1952年はベースカテゴリーとして認識されている
→ 例えば 1957年の係数0.033
→ 「ベースである 1952 年と比較して」「1957年は0.033だけ増える」と解釈する
固定効果を含むモデル(model_33)
の解釈
Model_33:
・一人当たりGDPが 1% 増えると、寿命が 0.135% 増える→ 一人当たりGDPが 10 % 増えると、寿命が約 1.35 %ポイント増える
estimatr パッケージの
lm_robust()関数を使っても可能
library(estimatr)fixed_robust <- lm_robust(log_lifeExp ~ log_gdpPercap,
fixed_effects = ~ year,
se_type = "stata", # ロバストな標準誤差のタイプを選択
data = gapminder_2)summary(fixed_robust)
Call:
lm_robust(formula = log_lifeExp ~ log_gdpPercap, data = gapminder_2,
fixed_effects = ~year, se_type = "stata")
Standard error type: HC1
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
log_gdpPercap 0.1347 0.002678 50.31 0 0.1295 0.14 1691
Multiple R-squared: 0.6837 , Adjusted R-squared: 0.6815
Multiple R-squared (proj. model): 0.6065 , Adjusted R-squared (proj. model): 0.6037
F-statistic (proj. model): 2531 on 1 and 1691 DF, p-value: < 2.2e-16fixest パッケージの
feols()関数を使っても可能
library(fixest)feols <- feols(log_lifeExp ~ log_gdpPercap,
cluster = ~year,
data = gapminder_2)summary(feols)OLS estimation, Dep. Var.: log_lifeExp
Observations: 1,704
Standard-errors: Clustered (year)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.863995 0.081979 34.9359 1.2684e-12 ***
log_gdpPercap 0.146576 0.007996 18.3315 1.3578e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 0.144431 Adj. R2: 0.612989事例 2：総選挙におけるパネルデータの分析
- 1996-2021までの総選挙データを読み込む
hr <- read_csv("data/hr96-21.csv",
na = ".")hr$exp <- as.numeric(hr$exp)- 選挙費用 (
exp) をチェック
hr |>
drop_na(exp) |>
group_by(year) |>
summarise(mean_exp = mean(exp),
obs = n())# A tibble: 8 × 3
year mean_exp obs
<dbl> <dbl> <int>
1 1996 9136316. 1198
2 2000 8388889. 1164
3 2003 7935408. 1009
4 2005 8142244. 985
5 2009 6118181. 1124
6 2012 5769988. 1280
7 2014 7440127. 39
8 2017 9298783. 30- 2014年と2017年の
expは欠損値があり不完全なので使わない
hr96_12 <- hr |>
filter(year < 2014)- 必要な変数を整備する
- district = ku + kun
hr96_12 <- hr96_12 |>
mutate(district = str_c(ku, kun, sep = "_")) |>
mutate(inc = if_else(status == 1, 1, 0)) |> # = 1 if incumbent, = 0 otherwise
mutate(exppv = exp/eligible) |> # 一人あたり選挙費用 exppv
select(year, seito, district, voteshare, exppv, nocand, previous, age, inc)hr96_12# A tibble: 6,908 × 9
year seito district voteshare exppv nocand previous age inc
<dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1996 新進 aichi_1 40 28.3 7 2 47 1
2 1996 自民 aichi_1 25.7 26.9 7 2 72 0
3 1996 民主 aichi_1 20.1 26.6 7 2 53 1
4 1996 共産 aichi_1 13.3 6.28 7 0 43 0
5 1996 文化フォーラム aichi_1 0.4 NA 7 0 51 0
6 1996 国民党 aichi_1 0.3 NA 7 0 51 0
7 1996 無所 aichi_1 0.2 NA 7 0 45 0
8 1996 新進 aichi_2 32.9 38.2 8 2 51 1
9 1996 自民 aichi_2 26.4 48.8 8 0 71 0
10 1996 民主 aichi_2 25.7 33.8 8 0 30 0
# ℹ 6,898 more rows1996年から2012年までの6回の総選挙を対象
次の回帰分析を実行する
- 応答変数・・・voteshare
- 説明変数・・・exppv（有権者一人あたりに費やすy選挙費用）
- コントロール変数・・・nocand, previous, age, inc
- 固定効果・・・ year, seito
- clusters・・・district
通常の重回帰分析
model_1 <- lm(voteshare ~ exppv +
nocand +
previous +
age +
inc,
data = hr96_12)通常の重回帰分析 + 固定効果 (year)
model_2 <- lm(voteshare ~ exppv +
nocand +
previous +
age +
inc +
factor(year),
data = hr96_12)通常の重回帰分析 + 固定効果 (year, seito)
model_3 <- lm(voteshare ~ exppv +
nocand +
previous +
age +
inc +
factor(year) +
factor(seito),
data = hr96_12)通常の重回帰分析 + 固定効果 (year, seito) + clusters(distict)
model_4 <- estimatr::lm_robust(voteshare ~
exppv +
nocand +
previous +
age +
inc +
factor(year) +
factor(seito),
clusters = district,
se_type = "stata",
data = hr96_12)model_1 から model_3 までの結果
stargazer(model_1,
model_2,
model_3,
type = "html")| Dependent variable: | |||
| voteshare | |||
| (1) | (2) | (3) | |
| exppv | 0.440*** | 0.471*** | 0.159*** |
| (0.009) | (0.010) | (0.009) | |
| nocand | -3.184*** | -3.004*** | -3.341*** |
| (0.138) | (0.148) | (0.118) | |
| previous | 2.256*** | 2.229*** | 1.693*** |
| (0.089) | (0.088) | (0.069) | |
| age | -0.214*** | -0.209*** | -0.118*** |
| (0.015) | (0.015) | (0.012) | |
| inc | 9.397*** | 8.830*** | 3.664*** |
| (0.424) | (0.423) | (0.338) | |
| factor(year)2000 | 1.405*** | 2.308*** | |
| (0.487) | (0.407) | ||
| factor(year)2003 | 3.605*** | 1.754*** | |
| (0.524) | (0.432) | ||
| factor(year)2005 | 2.888*** | 0.690 | |
| (0.533) | (0.438) | ||
| factor(year)2009 | 4.673*** | 5.248*** | |
| (0.507) | (0.449) | ||
| factor(year)2012 | 4.519*** | 1.139*** | |
| (0.486) | (0.432) | ||
| factor(seito)みんな | 15.592* | ||
| (9.027) | |||
| factor(seito)世界経済共同体党 | -0.089 | ||
| (10.991) | |||
| factor(seito)保守 | 20.701** | ||
| (9.260) | |||
| factor(seito)保守新 | 18.144* | ||
| (9.380) | |||
| factor(seito)公明 | 23.214** | ||
| (9.064) | |||
| factor(seito)共産 | 4.405 | ||
| (8.975) | |||
| factor(seito)国民新党 | 14.885 | ||
| (9.187) | |||
| factor(seito)安楽死党 | 3.299 | ||
| (12.685) | |||
| factor(seito)市民新党にいがた | 2.301 | ||
| (12.691) | |||
| factor(seito)幸福 | -6.934 | ||
| (8.992) | |||
| factor(seito)当たり前 | 1.831 | ||
| (12.683) | |||
| factor(seito)改革 | -7.779 | ||
| (12.694) | |||
| factor(seito)改革クラブ | 11.558 | ||
| (10.042) | |||
| factor(seito)政事公団太平会 | 14.308 | ||
| (12.701) | |||
| factor(seito)政治団体代表 | 12.272 | ||
| (10.988) | |||
| factor(seito)新党さきがけ | 14.155 | ||
| (9.321) | |||
| factor(seito)新党大地 | 12.384 | ||
| (9.515) | |||
| factor(seito)新党尊命 | 5.904 | ||
| (12.694) | |||
| factor(seito)新党日本 | 12.772 | ||
| (9.461) | |||
| factor(seito)新社会 | 5.031 | ||
| (9.099) | |||
| factor(seito)新進 | 23.789*** | ||
| (9.001) | |||
| factor(seito)日本新進 | 0.821 | ||
| (12.690) | |||
| factor(seito)日本未来 | 7.715 | ||
| (9.011) | |||
| factor(seito)日本維新の会 | 17.787** | ||
| (9.000) | |||
| factor(seito)民主 | 23.134*** | ||
| (8.978) | |||
| factor(seito)民主改革連合 | 25.549** | ||
| (10.996) | |||
| factor(seito)沖縄社会大衆党 | 31.737** | ||
| (12.691) | |||
| factor(seito)無所 | 11.260 | ||
| (8.986) | |||
| factor(seito)無所属の会 | 22.100** | ||
| (9.465) | |||
| factor(seito)社民 | 8.903 | ||
| (8.991) | |||
| factor(seito)緑の党 | 3.035 | ||
| (12.693) | |||
| factor(seito)自民 | 27.768*** | ||
| (8.980) | |||
| factor(seito)自由 | 6.763 | ||
| (9.060) | |||
| factor(seito)自由連合 | 0.250 | ||
| (8.999) | |||
| factor(seito)諸派 | -1.572 | ||
| (10.038) | |||
| factor(seito)青年自由 | 1.611 | ||
| (12.693) | |||
| Constant | 33.816*** | 29.530*** | 21.265** |
| (0.951) | (1.077) | (9.012) | |
| Observations | 6,756 | 6,756 | 6,756 |
| R2 | 0.613 | 0.620 | 0.781 |
| Adjusted R2 | 0.613 | 0.620 | 0.779 |
| Residual Std. Error | 11.883 (df = 6750) | 11.772 (df = 6745) | 8.968 (df = 6709) |
| F Statistic | 2,136.532*** (df = 5; 6750) | 1,101.649*** (df = 10; 6745) | 519.493*** (df = 46; 6709) |
| Note: | p<0.1; p<0.05; p<0.01 | ||
model_4 の結果
tidy(model_4) term estimate std.error statistic
1 (Intercept) 21.2646680 1.18285804 17.9773627
2 exppv 0.1594163 0.01376447 11.5817259
3 nocand -3.3409380 0.19770321 -16.8987546
4 previous 1.6930050 0.10528271 16.0805616
5 age -0.1184489 0.01414018 -8.3767577
6 inc 3.6640448 0.41684652 8.7899132
7 factor(year)2000 2.3075938 0.29918886 7.7128336
8 factor(year)2003 1.7538862 0.36532928 4.8008367
9 factor(year)2005 0.6902432 0.40074709 1.7223909
10 factor(year)2009 5.2476748 0.38348112 13.6843107
11 factor(year)2012 1.1393473 0.37721627 3.0204087
12 factor(seito)みんな 15.5915074 1.33430401 11.6851237
13 factor(seito)世界経済共同体党 -0.0885252 0.38982810 -0.2270878
14 factor(seito)保守 20.7008395 3.29429748 6.2838404
15 factor(seito)保守新 18.1443088 4.48664963 4.0440664
16 factor(seito)公明 23.2144967 1.45799443 15.9222122
17 factor(seito)共産 4.4054811 0.26170479 16.8337807
18 factor(seito)国民新党 14.8850695 2.95080739 5.0444057
19 factor(seito)安楽死党 3.2986367 0.41759442 7.8991398
20 factor(seito)市民新党にいがた 2.3014537 0.41366528 5.5635650
21 factor(seito)幸福 -6.9335745 0.41762398 -16.6024337
22 factor(seito)当たり前 1.8313532 0.09450484 19.3784075
23 factor(seito)改革 -7.7794500 0.58561965 -13.2841343
24 factor(seito)改革クラブ 11.5578391 1.69574992 6.8157686
25 factor(seito)政事公団太平会 14.3075697 0.90985251 15.7251528
26 factor(seito)政治団体代表 12.2718002 2.45154317 5.0057451
27 factor(seito)新党さきがけ 14.1549995 3.11200819 4.5485097
28 factor(seito)新党大地 12.3841303 2.14869192 5.7635672
29 factor(seito)新党尊命 5.9044980 0.56907360 10.3756314
30 factor(seito)新党日本 12.7723456 6.09692855 2.0948820
31 factor(seito)新社会 5.0313204 0.68132766 7.3845827
32 factor(seito)新進 23.7890163 0.72518594 32.8040232
33 factor(seito)日本新進 0.8211593 0.37295727 2.2017515
34 factor(seito)日本未来 7.7147413 0.55827700 13.8188414
35 factor(seito)日本維新の会 17.7871031 0.71562985 24.8551721
36 factor(seito)民主 23.1343992 0.50186191 46.0971410
37 factor(seito)民主改革連合 25.5487484 0.89199772 28.6421678
38 factor(seito)沖縄社会大衆党 31.7370846 0.44320085 71.6088084
39 factor(seito)無所 11.2603994 0.92632161 12.1560366
40 factor(seito)無所属の会 22.0997488 3.78631455 5.8367440
41 factor(seito)社民 8.9025608 0.92275061 9.6478514
42 factor(seito)緑の党 3.0347489 0.53629269 5.6587550
43 factor(seito)自民 27.7675600 0.60035088 46.2522184
44 factor(seito)自由 6.7629082 1.14627237 5.8999138
45 factor(seito)自由連合 0.2495205 0.60159801 0.4147629
46 factor(seito)諸派 -1.5721258 0.74939389 -2.0978631
47 factor(seito)青年自由 1.6109601 0.47942662 3.3601807
p.value conf.low conf.high df outcome
1 9.894119e-50 18.93704214 23.59229383 304 voteshare
2 6.112057e-26 0.13233062 0.18650202 304 voteshare
3 1.235101e-45 -3.72997801 -2.95189799 304 voteshare
4 1.568682e-42 1.48582993 1.90018016 304 voteshare
5 2.032114e-15 -0.14627392 -0.09062386 304 voteshare
6 1.127372e-16 2.84377495 4.48431457 304 voteshare
7 1.772277e-13 1.71885057 2.89633712 304 voteshare
8 2.482677e-06 1.03499195 2.47278051 304 voteshare
9 8.601578e-02 -0.09834622 1.47883252 304 voteshare
10 1.530056e-33 4.49306135 6.00228821 304 voteshare
11 2.738595e-03 0.39706184 1.88163280 304 voteshare
12 2.642997e-26 12.96586644 18.21714834 304 voteshare
13 8.205081e-01 -0.85562822 0.67857782 304 voteshare
14 1.142998e-09 14.21832706 27.18335200 304 voteshare
15 6.664971e-05 9.31548800 26.97312960 304 voteshare
16 6.239590e-42 20.34545793 26.08353537 304 voteshare
17 2.179799e-45 3.89049887 4.92046328 304 voteshare
18 7.838654e-07 9.07847620 20.69166278 304 voteshare
19 5.177046e-14 2.47689520 4.12037825 304 voteshare
20 5.800948e-08 1.48744389 3.11546342 304 voteshare
21 1.646746e-44 -7.75537419 -6.11177479 304 voteshare
22 4.892965e-55 1.64538678 2.01731966 304 voteshare
23 4.560247e-32 -8.93183130 -6.62706878 304 voteshare
24 5.044809e-11 8.22094556 14.89473262 304 voteshare
25 3.472666e-41 12.51716366 16.09797581 304 voteshare
26 9.440708e-07 7.44765808 17.09594226 304 voteshare
27 7.814916e-06 8.03119560 20.27880346 304 voteshare
28 2.020606e-08 8.15593835 16.61232228 304 voteshare
29 8.617682e-22 4.78467598 7.02431994 304 voteshare
30 3.700886e-02 0.77482092 24.76987026 304 voteshare
31 1.479212e-12 3.69060514 6.37203574 304 voteshare
32 6.981790e-102 22.36199675 25.21603583 304 voteshare
33 2.843380e-02 0.08725462 1.55506388 304 voteshare
34 4.861916e-34 6.61616485 8.81331772 304 voteshare
35 3.309903e-75 16.37888797 19.19531814 304 voteshare
36 3.178609e-139 22.14683630 24.12196217 304 voteshare
37 2.434604e-88 23.79347693 27.30401985 304 voteshare
38 2.278449e-192 30.86495478 32.60921439 304 voteshare
39 5.615580e-28 9.43758543 13.08321331 304 voteshare
40 1.364171e-08 14.64904608 29.55045155 304 voteshare
41 2.176661e-19 7.08677386 10.71834772 304 voteshare
42 3.523962e-08 1.97943316 4.09006467 304 voteshare
43 1.300116e-139 26.58619061 28.94892929 304 voteshare
44 9.689372e-09 4.50727557 9.01854086 304 voteshare
45 6.786080e-01 -0.93430290 1.43334399 304 voteshare
46 3.674257e-02 -3.04678169 -0.09746987 304 voteshare
47 8.782908e-04 0.66754528 2.55437489 304 voteshare結果
通常の重回帰分析の場合 (model_1)
・一人あたり選挙費用を 1 円使うと、得票率が 0.44%ポイント増える
通常の重回帰分析 + 固定効果 (year) の場合
(model_2)
・一人あたり選挙費用を 1 円使うと、得票率が 0.47%ポイント増える
通常の重回帰分析 + 固定効果 (year, seito) の場合
(model_3)
・一人あたり選挙費用を 1 円使うと、得票率が 0.159%ポイント増える
通常の重回帰分析 + 固定効果 (year, seito) の場合
(model_4)
・一人あたり選挙費用を 1 円使うと、得票率が 0.159%ポイント増える
・model_1 と model_2 では exppv の影響力を過大評価している: 0.44 と 0.47
・model_3 と model_4 ではより適切に exppv の影響力を評価している: 0.159
7.3 外生性が満たされない場合
解決策：
- 操作変数 (
intsrumental variable) を使って分析する
- 詳細は 31. 重回帰分析 13（操作変数法：IV 基礎）を参照