このセクションで使う R パッケージ一覧

library(broom)
library(car)
library(jtools)
library(normtest)
library(patchwork)
library(stargazer)
library(tidyverse)
library(zoo)

最小二乗法による回帰係数の推定

最良線形不偏推定量: BLUE

最良線形不偏推定量 (BLUE: Best Linear Unbiased Estimator) とは、次の 5 つの仮定が満たされるときの最小二乗法による回帰係数の推定量のこと

誤差項の期待値がゼロ
パラメータにおける線形性
誤差項の条件付き期待値がゼロ
多重共線性がない
誤差項の分散均一性

ここではこの BLUE に即して、社会科学において要求される最小二乗法による回帰分析の前提と回帰モデルの妥当性の診断方法を解説する

回帰分析の前提と妥当性の診断方法

回帰分析の前提	妥当性の診断方法
1. 回帰モデルは妥当なモデルである	学問上の知識に基づいて判断
2. 誤差項の期待値がゼロ
3. 加法性と線形性	Scatter plotで確認
4. 誤差が独立である	慎重に交絡変数を検討する
5. 誤差の分散が均一である	残差プロットによる診断
6. 誤差が正規分布に従う	正規QQプロットによる診断
7. 完全な多重共線性がない	分散拡大係数 (`VIF`) による診断

1. 回帰モデルは妥当なモデルである

どのような回帰モデルが「適切なのか」を判断するのはそれぞれの学問分野における知識に裏付けられて判断される
モデルの妥当性に関しては、次のセッションを参考にして下さい：
「18. 重回帰分析 1 (単回帰と重回帰）4.3.3 適切なモデルの選び方」
「26. 重回帰分析 9（コントロール変数の選び方：バックドア基準と DAG）」

2. 誤差項の期待値がゼロ　　

2.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\) の期待値はゼロと見なすことができる

2.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個生成

2.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)

2.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
R²	0.80	0.79
Adjusted R²	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\) の期待値はゼロと見なすことができる

3. 加法性と線形性

3.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\) という関係を考えることができる
線形回帰として推定可能

3.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 <- gapminder %>% 
  dplyr::select(year, lifeExp, gdpPercap)

対数変換を施した新たな変数 (log_gdp, log_lifeExp) を作り、変数を限定する

gapminder <- gapminder %>% 
  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), 
          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)

対数変数しない変数を使った散布図: `gdpPercap & lifeExp`

まず、「1 人あたり GDP (gdpPercap)」を x 軸、「平均寿命 (lifeExp)」を y 軸に指定して散布図を描いてみる

gapminder %>%
  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)

対数変数した変数を使った散布図: `log_gdpPercap & log_lifeExp`

「1 人あたり GDP: gdpPercap」の記述統計を見ると、データの範囲が 277米ドル（最低額）から 59,359 米ドル（最高額）まで非常に広範囲に及ぶ
「1 人あたり GDP」と「平均寿命」はリニアというより、対数関係のような関係
→　見やすくするため、x 軸を 対数変換 した変数 (log_gdpPercap) と「平均寿命」(log_lifeExp) との散布図を描いてみる

gapminder %>%
  ggplot() +
  geom_point(aes(x = log_gdpPercap, 
                 y = log_lifeExp)) +
  labs(x = "1 人あたりGDP (USD)の対数値", 
       y = "寿命（対数変換）") +
  theme_bw(base_family = "HiraKakuProN-W3")

変数を対数化したら、かなり綺麗な線形関係が確認できた
回帰分析をしてみる

model_4 <- lm(log_lifeExp ~ log_gdpPercap, data = gapminder)

model_1 と model_4 の分析結果を表示させる

stargazer(model_1, model_4, 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
R²	0.461	0.607
Adjusted R²	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_4:

・一人当たりGDPが 1% 増えると、寿命が 0.113% 増える
→ 一人当たりGDPが 10 % 増えると、寿命が約 1 %ポイント増える

4. 誤差が独立である

ここで、各個体 \(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. 重回帰分析（脱落変数バイアス・処置後変数バイアス）」

5. 誤差の分散が均一である

5.1 「誤差の分散が均一である」とは？

「誤差の分散が同じ」とは「全ての個体の誤差の分散が同じ」になるということ

\[全ての個体 i について、Var(ε_i) = σ^2\]

例えば、気温とビールの例だと、気温が高い夏でも寒い冬でも、ビールの出荷量のばらつきはなく一定であるということ
この事例のように分散が均一であることを均一分散性と呼ぶ
分散が均一でないことを不均一分散性と呼ぶ

`homoskerdasticity`	均一分散	分散が均一であること
`heteroskedasticity`	不均一分散	分散が均一でないこと

気温とビールの事例を図で表すと、左側のようになる
ビールの売り上げ額にかかわらず、売り上げのバラツキが一定

均一分散（左側の図）と不均一性分散（右側の図）の事例

不均一分散 (heteroskerdasticity) の事例としては「月収」と「食費」が考えられる
個人が 1 か月に使う「食費」を応答変数、その個人の「月収」を説明変数とするモデルを考える
この国で生きていくためには最低 8,000 円は毎月の食費に使う必要があるとする
月収が 1 万円の個人が食費にかける金額は 8,000 円以上 1 万円以下
これに対し、月収が 100 万円の個人が食費にかける金額は、8,000 円以上 100万円以下
つまり、月収が少ない場合、食費にかける金額を選択する余地が小さい
→　応答変数の値はあまりばらつかない
月収が増えれば選択の余地が広がり、食費を節約する人もいるし贅沢な食事を楽しむ人もいる
→　ばらつきが大きくなる
応答変数のばらつき自体は誤差のばらつきではないが、応答変数のばらつきが大きいほど説明変数で説明できない部分が大きくなる
→　誤差のばらつきも大きくなる

社会科学の事例における説明変数と誤差の分散応答変数のばらつきが大きくなるほど
→ 誤差のばらつきが大きくなる
→ 説明変数で説明できない部分が大きくなる
→　誤差に分散不均一性が生じる
→　標準誤差に影響がでて、統計的な有意性が得にくくなる
・ 誤差項の分散が不均一でも、最小二乗法による回帰係数の推定量に偏りは発生しない

5.2 不均一分散の診断

5.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.5442

`bytest()`関数の診断基準

`p-value`が 0.05以上の場合	誤差項の分散が均一である
`p-value`が 0.05以下の場合	誤差項の分散が不均一である

p-value = 0.5442 で 0.05 以上なので、「誤差項の分散が均一である」

5.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.0001846

`bytest()`関数の診断基準

`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. 不均一分散に頑健な標準誤差

もう一つの解決策として標準誤差を修正することで、不均一分散に対応する「不均一分散に頑健な標準誤差」を使う方法がある
誤差の分散が不均一の時
→「重み付け最小二乗法」や「不均一分散に頑健な標準誤差」を使うのが望ましいが、分散の不均一性は、回帰分析によって推定される説明変数と応答変数の関係には影響しないので、あまり重大な問題とは言えない

6. 誤差が正規分布に従う

このセッションでは回帰分析の前提として次の 2 つの仮定を紹介した

2. 誤差項の期待値がゼロ
6. 誤差の分散が均一である

誤差項の期待値はゼロであり、誤差の分散は均一なので、分散は定数 \(σ\) と表せる
「誤差の正規性」とは「誤差が正規分布に従う確率変数であること」
各個体の誤差が平均 0、分散 \(σ^2\) の正規分布に従うということ
すべての個体 \(i\) について次の式が成り立つということ

\[ε_i 〜 N (0, σ^2)\]

誤差項 \(ε_i\)に様々な要因が含まれていることから
→ 中心極限定理により、母集団がどのような分布であれ誤差の分布は正規分布で近似できる
→ 「誤差の正規性」仮定を正当化できる
→ 回帰直線の推定では「誤差の正規性」仮定は気にする必要はない

残差が正規分布に従っていないケース（左）と正規分布に従っているケース（右）の事例

最良線形不偏推定量 (BLUE: Best Linear Unbiased Estimator)である 5 つの仮定の中に「誤差が正規分布に従う」という仮定は含まれていない

誤差項の期待値がゼロ
パラメータにおける線形性
誤差項の条件付き期待値がゼロ
多重共線性がない
誤差項の分散均一性

誤差項が正規分布に従っていなくても、通常の最小二乗法による回帰係数の推定量は最良線形不偏推定量 (BLUE) である

ジャック・ベラの正規性検定 (`Jarque-Bera test for normality`)

Rパッケージ normtest の jp.norm.test()関数を使って検定可能

`bytest()`関数の診断基準

帰無仮説・・・「誤差項が正規分布に従っている」

`p-value`が 0.05以上の場合	誤差項が正規分布に従っている
`p-value`が 0.05以下の場合	誤差項が正規分布に従っていない

ジャック・ベラの正規性検定を使ってパットナムモデルと2009年総選挙モデルにおける誤差項の正規性を診断してみる

パットナムモデルの正規性検定

まず、パットナムモデルの残差の QQプロットを描く

putnam <- read_csv("data/putnam.csv")
fit.putnam <- lm(gov_p ~ cc, data = putnam) 
res.putnam <- resid(fit.putnam)     # 残差を取り出す

par(family = "HiraKakuProN-W3")　# 日本語文字化け防止 (Macユーザのみ)

qqnorm(res.putnam, main = "残差の正規QQプロット",
       xlab = "標準正規分布", ylab = "標準化した残差の分布")
qqline(res.putnam, col = "red")　　　　　　　　　　　　　# 正規分布に一致する場合の直線を引く

データが直線の周りに偏りなく分布しているので、正規QQプロットによる誤差が正規分布に従っているといえる
2009総選挙の残差のヒストグラムを描いてみる

hist(res.putnam)

次に bytest()関数を使って正規性を診断をする

library(normtest)

resid_putnam <- residuals(fit.putnam)
jb.norm.test(resid_putnam)


    Jarque-Bera test for normality

data:  resid_putnam
JB = 1.3857, p-value = 0.2475

p-value = 0.2475 なので帰無仮説「誤差項が正規分布に従っている」は棄却できない
→　「誤差項が正規分布でないとはいえない」
注意：ジャック・ベラの正規性検定における \(p-value\) はモンテカルロ・シミュレーションによる近似値なので、実行する度に異なる値が表示される

2009総選挙モデルの正規性検定

2009総選挙の残差の QQプロットを描く

hr09 <- read_csv("data/hr96-21.csv") %>% 
  filter(year == 2009)

fit.hr09 <- lm(voteshare ~ age, data = hr09) 
res.hr09 <- resid(fit.hr09)     # 残差を取り出す

par(family = "HiraKakuProN-W3")　# 日本語文字化け防止 (Macユーザのみ)

qqnorm(res.hr09, main = "残差の正規QQプロット",
       xlab = "標準正規分布", ylab = "標準化した残差の分布")
qqline(res.hr09, col = "red")　　　　　　　　　　　　　# 正規分布に一致する場合の直線を引く

データが直線の周りに偏りなく分布しているとはいえない
2009総選挙の残差のヒストグラムを描いてみる

hist(res.hr09)

次に bytest()関数を使って正規性を診断をする

library(normtest)

resid_hr09 <- residuals(fit.hr09)
jb.norm.test(res.hr09)


    Jarque-Bera test for normality

data:  res.hr09
JB = 73.299, p-value < 2.2e-16

p-value < 2.2e-16 なので帰無仮説「誤差項が正規分布に従っている」は棄却
→　「誤差項が正規分布に従っているとはいえない」

複数のプロット表示

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 つのプロットを表示

7. 完全な多重共線性がない

7.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\)、どちらかの説明変数だけを使う

7.2 なぜ、多重共線性が問題なのか？

重回帰分析で多重共線性は禁止されているわけではない
説明変数の間に相関があるのは当然のこと
相関係数の大きさの「程度問題」
相関係数が大きすぎることで起こること：
→　標準誤差が大きくなる
→　信頼区間の幅が大きくなる
→　帰無仮説を棄却しにくくなる
→　統計的に有意になりにくくなる
しかし、回帰係数の推定は不偏に行うことが可能
→　極端に大きな相関でない限り、分析する上で気にする必要はない
では、どの程度なら許容できるのか？

7.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}\]

7.4 共変量における多重共線性

7.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 つの説明変数のいずれか片方をモデルに含めずに分析

7.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
R²	0.035	0.743	0.793	0.794
Adjusted R²	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」は正しく推定できているので問題なし
・共変量の間に多重共線が起きていても、特に問題視する必要はない
・むしろ必要な共変量をモデルに含めないと、交絡を取り除けない
・必要と思われる共変量はモデルに含めるべき

参考文献

ロバート・パットナム（河田潤ー訳）『哲学する民主主義』NTT 出版、2001年.

森田果『実証分析入門』日本評論社、2014年.

宋財泫 (Jaehyun Song)・矢内勇生 (Yuki Yanai)「私たちのR: ベストプラクティスの探究」

高橋将宣『統計的因果推論の理論と実装』、共立出版、2022年

浅野正彦, 矢内勇生.『Rによる計量政治学』オーム社、2018年

Kosuke Imai, Quantitative Social Science: An Introduction, Princeton University Press, 2021

20. 重回帰モデルの仮定と診断

Masahiko Asano

2022-04-04

最小二乗法による回帰係数の推定

最良線形不偏推定量: BLUE

回帰分析の前提と妥当性の診断方法

1. 回帰モデルは妥当なモデルである

2. 誤差項の期待値がゼロ

2.1 証明

2.2 R を使った誤差項のシミュレーション

2.2.1 誤差項の期待値（平均値）が 0 の場合

2.2.2 誤差項の期待値（平均値）が 0 でない場合

3. 加法性と線形性

3.1 線形性とは

加法性が満たされない場合

例 1:

例 2:

3.2 変数変換の実例（対数変換）

データの準備 (gapminder)

対数変数しない変数を使った散布図: gdpPercap & lifeExp

対数変数した変数を使った散布図: log_gdpPercap & log_lifeExp

Model_1:

Model_4:

4. 誤差が独立である

誤差の独立性が満たされない一例：

誤差項の条件付き期待値がゼロということ

5. 誤差の分散が均一である

5.1 「誤差の分散が均一である」とは？

5.2 不均一分散の診断

5.2.1 パットナムのモデルでの不均一診断

bytest()関数の診断基準

5.2.2 2009年衆院選のモデルでの不均一診断

応答変数が得票率の場合

読み取り方法 1

読み取り方法 2

読み取り方法 3

bytest()関数の診断基準

誤差の分散が不均一の時の 2 つの対処法:

1. 加重最小二乗法 (WLS: Weighted least squares)

2. 不均一分散に頑健な標準誤差

6. 誤差が正規分布に従う

ジャック・ベラの正規性検定 (Jarque-Bera test for normality)

bytest()関数の診断基準

パットナムモデルの正規性検定

2009総選挙モデルの正規性検定

複数のプロット表示

7. 完全な多重共線性がない

7.1 多重共線性とは

7.2 なぜ、多重共線性が問題なのか？

7.3 多重共線性の診断方法: VIF

VIF の計算方法

7.4 共変量における多重共線性

7.4.1 イタリア地方政府に関するモデルの VIF

7.4.2 日本の2009年衆院選挙データにおける VIF

変数の記述統計

「年齢」と「ランク」の散布図を表示させる

「年齢」と「票数」の散布図を表示させる

「票数」と「得票率」の散布図を表示させる

モデルの推定

分析結果を解釈する上で大事なこと

Model_1:

Model_2 & Model_3

Model_4

Model_4 の VIF を計算する

2. 誤差項の期待値がゼロ　　

2.2.2 誤差項の期待値（平均値）が 0 でない場合　　

加法性が満たされない場合　　

データの準備 (`gapminder`)

対数変数しない変数を使った散布図: `gdpPercap & lifeExp`

対数変数した変数を使った散布図: `log_gdpPercap & log_lifeExp`

`bytest()`関数の診断基準

`bytest()`関数の診断基準

1. 加重最小二乗法 (`WLS: Weighted least squares`)

ジャック・ベラの正規性検定 (`Jarque-Bera test for normality`)

`bytest()`関数の診断基準

7.3 多重共線性の診断方法: `VIF`

7.4.1 イタリア地方政府に関するモデルの `VIF`

7.4.2 日本の2009年衆院選挙データにおける `VIF`

`Model_1`:

`Model_2` & `Model_3`

`Model_4`

`Model_4` の `VIF` を計算する