固定効果モデルの診断
太郎丸 博
Wednesday, June 03, 2015
モデルの診断とは、回帰分析などのモデルを推定した後に、モデルが適切かどうかチェックする作業のことである。データの分布のチェックもモデルの診断も、適切なモデルを選ぶための探索的な作業である。ただし、データの分布の探索が事前のチェックであるのに対して、モデルの診断は事後的なチェックである、という点で異なる。
モデルの診断はしばしば軽視されがちだが、正しいモデルを見つけるために、とても重要なので適切なモデルの目星がついてきたら、手を抜かずに必ずやってほしい。
1 Assumptions: OLS の大前提
OLS でモデルを推定する場合、以下のような条件が大前提として満たされていないと、正確な推定にならないことが知られている。
- \(\mathrm{E}(\epsilon_{it} | X) = 0\). 残差の期待値は X の値にかかわらず常にゼロ。
- \(\mathrm{E}(\epsilon_{it}X_{it}) =0\). 残差と独立変数の共分散がゼロ。
- \(\mathrm{Var}(\epsilon_{it} | X) = \sigma\). 残差の分散は X の値にかかわらず一定。
- \(\epsilon_{it}\) は正規分布にしたがう。
- 残差は互いに独立に分布。
特に重要なのは、1番目と2番目、5番目の条件で、これらが満たされていないと、推定値の期待値と真のパラメータの値が一致しない（一致推定量にならない）。3番目の条件が満たされないと、最小分散の推定量でなくなり、4番目の条件が満たされないと最尤推定量でなくなる。これらも満たされたほうがいいが、致命的ではない場合もしばしばあるのでケースバイケースで判断すればよい。
2番目の条件は重要だが、これをチェックする方法はない。なぜなら OLS では残差を説明変数とは相関しないように推定するので。5番目の条件も重要だが、これはダイナミック・モデルを紹介するときに詳しく考える。それゆえ、以下で検討する様々な問題の中で 1番重要なのは 1番目の条件のチェックである。
以下では、基本的なモデルの診断法を紹介していくが、まずは横断的データの回帰分析を例に、診断法や問題がある場合の対処法を論じていく。もしもその診断法や対処法がパネル・データの分析にそのまま応用できない場合は、パネル・データでの診断法や対処法を述べる。これは、最近のテキストではモデルの診断法が軽んじられる傾向があり、横断データの分析でも診断法があまり知られていないと感じているからである。また横断データのほうが単純で話が簡単だからでもある。
2 Residuals 残差の検討
残差の期待値が X の値にかかわらず常にゼロかどうかを検討する方法として、残差と独立変数の散布図や、残差とやモデルからの Y の予測値 \(\hat Y\) の関連を散布図にして見るという方法がある。
2.1 Simple Regression 単回帰分析の場合の残差の期待値の検討
2.1.1 Quadratic 二次曲線的なデータに線形モデルを当てはめた場合
なぜ残差の検討が適切なモデルの推定のために重要なのか、横断データを例に説明しよう。同じことがパネルにもあてはまる。話を簡単にするために、まずは説明変数が一つの場合を例に解説しよう。すべてのケースの残差のベクトルは、lm() で推定した場合、
residuals( lm()の結果 )
または
lm()の結果$residuals
で得られる。plm() で推定した場合もまったく同じやり方で残差のベクトルが得られる。
図1 の左上のグラフは、\(Y\) と \(X\) という変数の散布図を作り、単回帰分析をして回帰直線を引いたものである。
lm11 <- lm(Y ~ X, data=d01)
par(mfrow=c(2,2), mar=c(2.9, 2.9, 0.5, 0.1), mgp=c(1.7, 0.6, 0))
# mar=は下、左、上、右の順でグラフの周りに何行余白をあけるかの指定
# mgp=は、XとYのラベル、数値、目盛をグラフから何行離すかの指定
plot(X, Y, ylim=c(0,23)) # 少しグラフの上に数式を書く余白を作りたいので、ylim=c(Y軸の下限値, Y軸の上限値) でY軸の上限値を少し大きめに指定
abline(lm11, col="red") # col="色の名前" で線の色を指定
text(17, 21.5, expression(paste(hat(Y)," = 8.3 + 0.5X"))) # この切片と傾きはちょっと不正確
plot(d01$X, residuals(lm11), ylab="残差")
abline(h=0, col="blue", lty=2) # lty= で線の種類を指定。1～6まで指定できる
plot(predict(lm11), residuals(lm11), ylab="残差", xlab="Yの予測値")
abline(h=0, col="blue", lty=2)
図1: 架空の単回帰分析の結果と残差のプロット
回帰直線はモデルからの予測値を示しているが、実際の \(Y\) の値は \(X\)が小さいうちは回帰直線よりも下にあり、\(X\) が5～15ぐらいの間は、回帰直線よりも上、\(X\)が15以上だと回帰直線よりも下に来る。残差とは、実際の \(Y\) の値とこの直線の高さとの差である。すなわち、 \[ \epsilon_{i} = Y_{i} - \hat Y_{i} = Y_{i} - (b_0 + b_1 X_{i}) \] ただし、\(b_0\)と\(b_1\)は回帰直線の切片と傾きである。この残差 \(\epsilon\) と \(X\) の散布図が図1 の右上のグラフである。実際の値が予測値よりも下の場合、残差は0よりも小さく、実際の値が予測値よりも大きい場合、残差は0よりも大きくなるのがわかるだろう。図1の左下は、横軸に \(Y\) の予測値 \(\hat Y\) の散布図である。\(X\) に比例して \(\hat Y\) も大きくなるので、右上の散布図とほとんど同じ形になるのがわかる。
これらのグラフでは、\(X\)が5以下と15以上では残差の期待値が0より大きく、5～15の間では0より小さいと考えられる。それゆえ、このモデルは\(\mathrm{E}(\epsilon_{i} | X) = 0\)という条件を満たしていないと考えられる。そこで、\(X\)の二乗項をモデルに投入して図1と同じようなグラフを作ったのが、図2 である。
lm12 <- update(lm11, ~.+I(X^2))
par(mfrow=c(2,2), mar=c(2.9, 2.9, 0.5, 0.1), mgp=c(1.7, 0.6, 0))
plot(X, Y, data=d01, ylim=c(0,21))## Warning in plot.window(...): "data" はグラフィックスパラメータではありません## Warning in plot.xy(xy, type, ...): "data" はグラフィックスパラメータではありません## Warning in axis(side = side, at = at, labels = labels, ...): "data" はグラフィック
## スパラメータではありません## Warning in axis(side = side, at = at, labels = labels, ...): "data" はグラフィック
## スパラメータではありません## Warning in box(...): "data" はグラフィックスパラメータではありません## Warning in title(...): "data" はグラフィックスパラメータではありません# abline()は単回帰の時しか使えないので lines()で回帰曲線を描く
Y.hat <- predict(lm12, data.frame(X=1:20))
lines(1:20, Y.hat, col="red")
plot(d01$X, residuals(lm12), ylab="残差")
abline(h=0, col="blue", lty=2)
plot(predict(lm12), residuals(lm12), ylab="残差", xlab="Yの予測値")
abline(h=0, col="blue", lty=2)
図1.5: 図1のデータを二次曲線のモデルを当てはめた場合
図2の左上を見ると、モデルから予測される \(Y\) の曲線の上下に点が均等に分布しており、\(X\)の値によって、曲線の上や下に偏っていない。このような状態が OLS の推定が適切になるために最も重要な条件である。この時、図2の右上のグラフを見ると、残差は\(X\)の値にかかわらず、0 を中心に分布しており、このモデルは\(\mathrm{E}(\epsilon_{i} | X) = 0\)という条件を満たしていると思われる。左下のグラフに関しても同様の傾向が見られる。
2.1.2 Threshold 閾値がある場合に線形モデルを当てはめた場合
残差と\(\hat Y\)の散布図を検討していると、残差が波動する場合がある。これは独立変数の効果が非線形である可能性を示唆している。この場合、上の例のように二乗項や三乗項を投入すべきかもしれないが、閾値 (threshold value) が存在している可能性を考えるべきである。独立変数が一つの場合で閾値について解説しよう。図2の左側のようなデータに関して、単純な回帰分析をし、回帰直線を引き、左側はその結果得られた残差と独立変数の散布図である。
par(mfrow=c(2,2), mar=c(2.9, 2.9, 0.5, 0.1), mgp=c(1.7, 0.6, 0))
par(mfrow=c(1,2))
plot(Y ~ X, data=d2)
lm22 <- lm(Y ~ X, data=d2)
abline(lm22, col="red")
plot(d2$X, residuals(lm22))
abline(h=0, col="blue", lty=2)
図2: 閾値の例
図2の左側を見ると、\(X\) が 10 以下のときは \(Y\) は 5 前後の値をとっているが、\(X\) が 10 より大きくなると、20前後の値をとる。それゆえ、\(X\)に比例して \(Y\)が増加しているわけではなく、10 より大きいかどうかで \(Y\) の値が大きく変化しているのである。このようなモデルの特定の誤りを犯している場合、残差は X の値とともに波動する。それゆえ、このようなかたちの残差が得られたら、閾値の存在を疑うべきである。
閾値の存在が疑われる場合、独立変数が閾値よりも大きいとき 1 、閾値以下のとき 0 になるようなダミー変数を作る。この例の場合、\(X\)が10より大きいとき 1、10以下のとき 0 になるような変数を作ってそれに回帰させればよい。そうすると、以下のように決定係数も上昇するのがわかる。
lm23 <- lm(Y ~ X>10, data=d2) # X>10 は Xが10より大きいとき TRUE、10以下の時 FALSE を返すが、lm() 中では TRUE=1、FALSE=0のダミー変数に変換される
library(memisc)## Loading required package: lattice
## Loading required package: MASS
##
## Attaching package: 'memisc'
##
## The following objects are masked from 'package:stats':
##
## contr.sum, contr.treatment, contrasts
##
## The following objects are masked from 'package:base':
##
## as.array, trimwsmtable(lm22, lm23)##
## Calls:
## lm22: lm(formula = Y ~ X, data = d2)
## lm23: lm(formula = Y ~ X > 10, data = d2)
##
## ===================================
## lm22 lm23
## -----------------------------------
## (Intercept) 0.598 4.834***
## (1.017) (0.375)
## X 1.140***
## (0.085)
## X > 10 15.465***
## (0.531)
## -----------------------------------
## R-squared 0.648 0.896
## adj. R-squared 0.644 0.895
## sigma 4.897 2.655
## F 180.131 848.204
## p 0.000 0.000
## Log-likelihood -299.749 -238.531
## Deviance 2350.244 690.839
## AIC 605.498 483.061
## BIC 613.314 490.877
## N 100 100
## ===================================2.1.3 Dummy: 連続変数のダミー変数化の危険性
上の例のように、連続変数はいつも線形の効果を持っているとは限らないので、モデルの特定には慎重であるべきだが、それに対する対処法として、初めから連続変数をダミー変数にしてしまうという方法がある。例えば、上の2つの例ではXは 1～20の間の値を取るので、X が 1～5, 6～10, 11～15, 16～20 に対応する4カテゴリからなる変数を作り、ダミー変数としてモデルに投入するのである。例えば以下のように。
d01$Xc <- cut(d01$X, c(0, 5.5, 10.5, 15.5, 21), # 4カテゴリの変数に変換
labels=c("1-5", "6-10", "11-15", "16-20"))
xtabs(~Xc + X, data=d01) # イメージ通りのカテゴリカル変数ができているか確認## X
## Xc 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 1-5 5 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 6-10 0 0 0 0 0 5 5 5 5 5 0 0 0 0 0 0 0 0 0 0
## 11-15 0 0 0 0 0 0 0 0 0 0 5 5 5 5 5 0 0 0 0 0
## 16-20 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 5 5 5 5lm13 <- lm(Y ~ Xc, data=d01)
library(memisc)
mtable(lm11, lm12, lm13, summary.stats=c("adj. R-squared", "AIC", "BIC"))##
## Calls:
## lm11: lm(formula = Y ~ X, data = d01)
## lm12: lm(formula = Y ~ X + I(X^2), data = d01)
## lm13: lm(formula = Y ~ Xc, data = d01)
##
## =============================================
## lm11 lm12 lm13
## ---------------------------------------------
## (Intercept) 8.456*** 0.393 7.126***
## (0.775) (0.677) (0.559)
## X 0.402*** 2.601***
## (0.065) (0.149)
## I(X^2) -0.105***
## (0.007)
## Xc: 6-10/1-5 7.089***
## (0.791)
## Xc: 11-15/1-5 9.105***
## (0.791)
## Xc: 16-20/1-5 5.995***
## (0.791)
## ---------------------------------------------
## adj. R-squared 0.275 0.784 0.593
## AIC 551.132 430.881 495.376
## BIC 558.947 441.302 508.402
## =============================================plot(d01$X, residuals(lm13))
上の表を見ると、3つのダミー変数を使ったモデルは単純な線形モデルよりもずっと決定係数が大きいが、二次曲線のモデルに比べると、やはり当てはまりが悪い。これは、このモデルでは X が同じカテゴリにあるときは、Y の値も同じであると予測するからである。例えば、X が 1～5 の間は、Y=7.126 と予測されているが、実際には、X=1 のときと X=5 のときの Yの値はかなり異なっており、そのような違いをこのモデルはうまく捉えていないのである。残差をプロットしてみると、やはり若干のうねりが見られ、あまりよいモデルとは言えないことがわかる。
また、上の例では4カテゴリに分類したが、区切り値をどうするか、いくつのカテゴリに分類するかで結果が変わってくる場合もあり、注意が必要である。そもそも理論的に考えた場合、上のモデルでは、X が 1-5 の間は Y は変化しないのに、5 から 6 に変化したときに急激に変化し、またその後 10 までは変化せず、… といった段階的な上昇を Y が繰り返すと仮定しているのだが、そんなモデルを説得力をもって説明／正当化できるだろうか。
以上のような事柄を勘案すると、探索的な分析の段階で連続変数をカテゴリカル変数に変換して分析するのはいいと思うが、連続変数として扱ったほうがいい場合も少なくないので、モデルの当てはまりの良さを見ながら、どうするか決めていったほうがいいだろう。「とりあえずカテゴリに変換しておけば大丈夫」、といった分析のやり方は間違っている。
2.2 Multple Regression 重回帰分析の残差の期待値の検討
上の単回帰分析の例では、わざわざ残差を検討しなくても、散布図を見ただけで、\(X\)と\(Y\)の非線形な関係ははっきりしていた。しかし、独立変数が複数あり、独立変数間に相関があると、散布図を見ただけでは、モデルの特定の失敗に気づかない場合がある。
lm21 <- lm(Y ~ X1 + X2, data=d1)
summary(lm21)##
## Call:
## lm(formula = Y ~ X1 + X2, data = d1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.4152 -2.9772 0.2354 2.4689 8.8028
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.78005 0.83160 -3.343 0.001179 **
## X1 0.57594 0.16815 3.425 0.000902 ***
## X2 1.02535 0.07757 13.218 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4 on 97 degrees of freedom
## Multiple R-squared: 0.9422, Adjusted R-squared: 0.941
## F-statistic: 790.2 on 2 and 97 DF, p-value: < 2.2e-16par(mfrow=c(2,2), mar=c(2.9, 2.9, 0.5, 0.1), mgp=c(1.7, 0.6, 0))
plot(Y ~ X1, data=d1)
plot(residuals(lm21) ~ X1, data=d1)
abline(h=0, col="blue", lty=2)
plot(d1$X2, residuals(lm21))
abline(h=0, col="blue", lty=2)
plot(predict(lm21), residuals(lm21))
abline(h=0, col="blue", lty=2)
図3: 相関する2つの独立変数を使った重回帰分析の残差のプロット
\(Y\)は\(X_1\)と\(X_2\)に線形に回帰させられており、重回帰分析の結果は 2変数とも傾きが有意である。図2の右上のグラフを見ると、\(X_1\)と\(Y\)の関係は線形のように見える。しかし、\(X_1\)と残差のプロット（図2の右上のグラフ）を見ると、残差は\(X_1\)とともに波動しているのがわかる。左下の\(X_2\)と残差のプロットは、ほとんどそのような波動は見られず、\(X_2\)の値にかかわらず、残差の期待値は 0 であるように見える。最後に右下の \(Y\) の予測値と残差の散布図は、右上のグラフほどではないが、やはり波動しているのがわかる。この波動の仕方は、上で挙げた閾値がある場合と同じであることがわかるだろう。
このように独立変数が複数ある場合、独立変数と従属変数の散布図を見ているだけでは非線形関係の存在を見落とすことがある。残差の散布図を見ることで、非線形関係が見つかることもあるので、注意してほしい。上で見たように、独立変数と残差の関連をプロットするのがベストであるが、その場合独立変数の数だけ散布図を作らなければならないので、手間がかかる。それゆえ、\(\hat Y\) と残差の散布図をまずは見てみて、異常な兆候を見つけたら独立変数と残差のプロットを作るという方法もあろう。
2.3 Heteroscedacity 分散の不均一性
残差分散の不均一性も、残差のプロットを見ればわかる。よくある分散の不均一性は、\(\hat Y\) が大きいほど残差の分散が大きくなるという場合である。
par(mfrow=c(2,2), mar=c(2.9, 2.9, 0.5, 0.1), mgp=c(1.7, 0.6, 0))
par(mfrow=c(1,2))
plot(Y ~ X, data = d3)
lm24 <- lm(Y ~ X, data = d3)
abline(lm24, col = "red")
plot(residuals(lm24) ~ d3$X)
abline(h=0, col = "blue", lty=2)
図4: 残差の分散がXに比例して大きくなる場合
これも独立変数が一つなので、残差を見なくても X と Y の散布図だけで残差の不均一性は明らかだが、やはり重回帰分析の場合は残差を見ないとわからないので、残差をチェックする癖をつけたい。
2.3.1 WLS 重みづけ最小二乗法
モデルが正しく特定されているにもかかわらず、残差分散が不均一になる場合、重み付き最小二乗法 (Weighted Least Square: WLS) で推定することによって、OLS よりも誤差の小さい推定ができる。OLS が残差の二乗和 (\(\sum_i (Y_i- \hat Y_i)^2\)) を最小化するように係数を推定するのに対して、WLS では、個々のケースに関して残差の分散(\(\mathrm{Var}(\epsilon_i)\))の逆数で重みづけして、\(\sum_i \frac{(Y_i- \hat Y_i)^2}{\mathrm{Var}(\epsilon_i)}\)を最小化するように推定する。このように推定すれば、最小分散の推定量が得られることが知られている。
問題は、どうやって個々のケースの残差の分散を知るか、ということである。もしも事前の経験的／理論的な知識から残差の分散がどう変化するかわかるならば、lm() 関数で
weights=重みのベクトル
という引数を追加で指定すればよい。例えば、Xに比例して残差の分散が大きくなると考えられるならば、“weights=1/X” と指定する（重みは残差分散の逆数でなければならないから）。
以下は計算例。係数の推定値はほとんど同じであるが、標準誤差が小さくなっているのがわかるだろう（重みを付けたほうが決定係数が大きくなっているが、重みによるものなので、あまり比較しても意味が無い）。
lm25 <- update(lm24, weights=1/X) # Xに比例して残差の分散が大きくなるとわかっている場合
library(memisc)
mtable(lm24, lm25)##
## Calls:
## lm24: lm(formula = Y ~ X, data = d3)
## lm25: lm(formula = Y ~ X, data = d3, weights = 1/X)
##
## ===================================
## lm24 lm25
## -----------------------------------
## (Intercept) -0.088 -0.056
## (0.514) (0.232)
## X 0.503*** 0.500***
## (0.043) (0.030)
## -----------------------------------
## R-squared 0.583 0.734
## adj. R-squared 0.579 0.731
## sigma 2.476 0.675
## F 137.139 270.603
## p 0.000 0.000
## Log-likelihood -231.547 -207.472
## Deviance 600.788 44.699
## AIC 469.095 420.944
## BIC 476.910 428.760
## N 100 100
## ===================================しかし、残差の分散について事前の知識などない場合が多いだろう。そのような場合は、データから残差の分散を推定する。推定には、nlmeパッケージの gls() 関数を使う。一般的には以下のような書式に従う。
gls(モデル式, data=データフレーム名, weights=分散のモデル)
ほとんど lm() で OLS の推定をするときと同じであるが、weights= 以下が異なる。weights= 以下は、どのように残差の分散が変化するのか、指定してやる必要がある。もしも\(\hat Y\)が大きくなるにしたがって、残差の分散も大きくなると考えられるなら、以下のいずれかを引数に指定すればいいだろう。
| 引数 | 残差の分散を推定するためのモデル |
|---|---|
| weights=varExp() | \[ \mathrm{Var}(\epsilon_i) = \exp({2\gamma \hat Y_i})\] |
| weights=varPower() | \[ \mathrm{Var}(\epsilon_i) = |\hat Y_i| ^{2\gamma}\] |
| weights=varConstPower() | \[ \mathrm{Var}(\epsilon_i) = \gamma_1 + |\hat Y_i| ^{2\gamma_2} \] |
\(\hat Y\)ではなく、特定の変数の値によって残差分散が変化すると考えられるならば、weights = varExp(form = ~X) のように、分散モデルの引数に
form= ~ 残差分散の共変量
と指定すればよい。以下は gls()の使用例。
library(nlme)
gls1 <- gls(Y~X, data=d3, weights=varExp())
gls2 <- gls(Y~X, data=d3, weights=varPower())
gls3 <- gls(Y~X, data=d3, weights=varConstPower(form=~X))
library(stargazer)##
## Please cite as:
##
## Hlavac, Marek (2014). stargazer: LaTeX code and ASCII text for well-formatted regression and summary statistics tables.
## R package version 5.1. http://CRAN.R-project.org/package=stargazerstargazer(lm24, gls1, gls2, gls3, type="text")##
## ======================================================================
## Dependent variable:
## --------------------------------------------------
## Y
## OLS generalized
## least squares
## (1) (2) (3) (4)
## ----------------------------------------------------------------------
## X 0.503*** 0.501*** 0.494*** 0.493***
## (0.043) (0.033) (0.027) (0.028)
##
## Constant -0.088 -0.083 -0.020 -0.021
## (0.514) (0.197) (0.104) (0.112)
##
## ----------------------------------------------------------------------
## Observations 100 100 100 100
## R2 0.583
## Adjusted R2 0.579
## Log Likelihood -206.170 -201.402 -201.273
## Akaike Inf. Crit. 420.340 410.804 412.547
## Bayesian Inf. Crit. 430.680 421.144 425.472
## Residual Std. Error 2.476 (df = 98)
## F Statistic 137.139*** (df = 1; 98)
## ======================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01summary(gls3)## Generalized least squares fit by REML
## Model: Y ~ X
## Data: d3
## AIC BIC logLik
## 412.5469 425.4717 -201.2734
##
## Variance function:
## Structure: Constant plus power of variance covariate
## Formula: ~X
## Parameter estimates:
## const power
## 0.4451101 1.0456250
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.0209488 0.11199229 -0.187055 0.852
## X 0.4931443 0.02780285 17.737185 0.000
##
## Correlation:
## (Intr)
## X -0.652
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.42123270 -0.63629349 0.04744662 0.57934626 2.69480531
##
## Residual standard error: 0.1788628
## Degrees of freedom: 100 total; 98 residual係数の推定値にはそれほど大きな違いはないが、OLSの場合よりも標準誤差が小さくなっているのがわかるだろう。
2.3.2 GLS for Panel: パネルデータで GLS
パネルデータの場合、R で上のような分散の不均一性を扱うのは、ちょっと面倒である。plm 関数は重みを指定できない（weights= という引数を指定してもエラーは出ないが、指定しない場合とまったく同じ結果が返ってくる)。plm パッケージには pggls() という GLS で分散の重み付けをやってくれる関数があるが、これは個体間の分散の不均一性しか扱えない。それゆえ、lm() や gls() を使う必要があるが、これらは固体内偏差を使うと、残差の自由度を過大に見積もってしまうので、標準誤差や p値を過小に推定してしまう。そこで、階差をとって gls() で推定するのが、R を使う場合の最も現実的な分散の不均一性に対する対処方法といえよう。以下は gls() を使った階差の回帰分析の例。
head(d4) # 架空の残差分散が不均一なパネル・データ## id t x y
## 1 1 1 20.7 506.5
## 2 1 2 14.0 326.1
## 3 1 3 20.0 448.1
## 4 1 4 22.0 477.8
## 5 1 5 17.8 511.1
## 6 2 1 21.6 337.6plot(y ~ x, data=d4)
図5: 架空の残差分散が不均一なパネル・データ
library(plm)## Loading required package: Formulapgls1 <- pggls(y~x+as.numeric(t), data=pdata.frame(d4)) # 比較のために推定してみた
summary(pgls1)## Within model
## Random effects model
## NA
## NA
##
## Call:
## pggls(formula = y ~ x + as.numeric(t), data = pdata.frame(d4))
##
## Balanced Panel: n=100, T=5, N=500
##
## Residuals
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -354.6000 -54.7100 -0.9014 0.0000 56.5100 223.3000
##
## Coefficients
## Estimate Std. Error z-value Pr(>|z|)
## x 17.97211 2.47997 7.2469 4.264e-13 ***
## as.numeric(t) -0.32411 2.75556 -0.1176 0.9064
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Total Sum of Squares: 5735000
## Residual Sum of Squares: 3835800
## Multiple R-squared: 0.33116# 独立変数と従属変数の階差を作る
x.diff <- as.numeric(diff(d4$x)) # gls()でなぜかエラーが出るので as.numeric() で
y.diff <- as.numeric(diff(d4$y)) # 余分な属性を落としてただの数値ベクトルに変換
head(cbind(x.diff, y.diff))## x.diff y.diff
## [1,] -6.7 -180.4
## [2,] 6.0 122.0
## [3,] 2.0 29.7
## [4,] -4.2 33.3
## [5,] 3.8 -173.5
## [6,] -3.7 -74.9library(nlme)
gls4 <- gls(y.diff ~ -1 + x.diff, weights=varExp() ) # モデル式の "-1" は切片をゼロに固定するという指定。固体内偏差の回帰分析の切片は必ずゼロなので。
summary(gls4)## Generalized least squares fit by REML
## Model: y.diff ~ -1 + x.diff
## Data: NULL
## AIC BIC logLik
## 6355.558 6368.19 -3174.779
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~fitted(.)
## Parameter estimates:
## expon
## -0.0002145488
##
## Coefficients:
## Value Std.Error t-value p-value
## x.diff 19.42466 2.224421 8.732457 0
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.786125131 -0.651027273 0.004692146 0.680604919 3.213129268
##
## Residual standard error: 140.865
## Degrees of freedom: 499 total; 498 residual2.4 Normality 残差の正規性
残差が正規分布しているかどうかを視覚的にチェックしたい場合、残差のヒストグラムを作るという方法がまず考えられる。
library(plm)
p.d4 <- pdata.frame(d4)
pl1 <- plm(y ~ x, data=p.d4, model="within")
hist(residuals(pl1))
これで大雑把な分布はつかめるが、さらに詳しく検討したい場合、分位点・分位点プロット（あるいはQQプロットとも呼ばれる）という方法もある。これは qqnorm() という関数でできる。分位点・分位点プロットは以下のようなステップにしたがってグラフを作っている。
- 変数 X の個々の値が、その変数の何パーセンタイルに対応するのか計算する。
- 計算したパーセンタイルに対応する標準正規分布の値 Y を計算する。
- Y と X の散布図を描く。
もしも X が正規分布していれば、分位点・分位点プロットはきれいな右肩上がりの直線になる。下の図6は、正規分布、ポアソン分布、対数正規分布のヒストグラムとQQプロットをしめしたものである。
par(mfrow=c(3,2), mar=c(2.9, 2.9, 1, 0.1), mgp=c(1.7, 0.6, 0), cex.main=1)
x <- rnorm(1000) # 正規分布する乱数を1000作る
hist(x)
qqnorm(x)
qqline(x, col="red")
x <- rpois(1000, 1.5) # 平均1.5のポアソン分布する乱数を1000作る
barplot(table(x), xlab="x", ylab="Frequency")
qqnorm(x)
qqline(x, col="red")
x <- rlnorm(1000, 3) # 対数変換した時の平均が3の対数正規分布の乱数を1000作る
hist(x)
qqnorm(x)
qqline(x, col="red")
図6: 3種類の分布のヒストグラムとQQプロット（上段：正規分布、中段：ポアソン分布、下段：対数正規分布）
resid1 <- residuals(pl1)
qqnorm(resid1)
qqline(resid1)
図7: plm()の残差のQQプロット
図7 でプロットした残差も人工的に正規分布するように作ってあるので、この程度直線に近ければ十分な正規性があるといえる。おおむね正規分布に近い分布をしていれば、実害はないと思われるので、残差の正規性を検定する必要はあまりないと思う。
しかし、残差や変数の正規性を検定したい場合には、R では shapiro.test() と ks.test() がパッケージを特に追加しなくても使える。shapiro.test() は Shapiro-Wilk の正規性の検定のための関数で、サンプルサイズが3～5000のときしか使えないという点に注意が必要である。ks.test() は Kolmogorov-Smirnov の2変数の分布の同一性の検定をするための関数（これはサンプルサイズに制限なし）で、一変数の正規性の検定にも使える。
# Shapiro-Wilk の正規性の検定 3 < n < 5000 のときしか使えない関数
shapiro.test(x)##
## Shapiro-Wilk normality test
##
## data: x
## W = 0.6105, p-value < 2.2e-16# Kolmogorov-Smirnov 分布の同一性の検定。正規性の検定をするためには pnorm, mean=mean(ベクトル名), sd=sd(ベクトル名)の指定が必要
ks.test(resid1, pnorm, mean=mean(resid1), sd=sd(resid1))##
## One-sample Kolmogorov-Smirnov test
##
## data: resid1
## D = 0.024774, p-value = 0.9188
## alternative hypothesis: two-sided2.4.1 Nonnormal Residuals 残差が正規分布しない場合
残差が正規分布しない場合、OLS の結果は最尤推定値にならない。しかし、最小分散の一致推定量であるということにはかわりはないので、あまり神経質になる必要はない。しかし、尤度を使うようなモデルの適合度指標（例えば AIC や BIC）にはあまり意味がないかもしれない。
また、従属変数を対数変換したり、残差に関して別の分布を仮定するモデル（順序ロジットやポアソン回帰など）やトービットのようなモデルを使えば対応できるかもしれない。これに関しては、別途論じる。
3 Outlier 外れ値
外れ値 (outlier) とは、極端に大きな値や小さな値のことで、外れ値の存在は回帰分析の結果に大きな影響を与えることがある。もしも外れ値をもったケースが非常に特殊な存在で、他のサンプルと一緒に分析するのは不適切だと判断できるならば、分析から除外することも考えられる。
3.1 Outlier of Residuals: 残差の外れ値
回帰分析においては、残差の外れ値と独立変数の外れ値が問題になりうる。残差が正規分布しているならば、残差のZ得点が 1.96 以上の絶対値をとることは、全サンプルの 5% 、2.58 以上は 1%、3.29 以上は 0.1% 、3.89 以上は 0.01%程度になるはずである。例えば1000ケースのデータで回帰分析して、外れ値のZ得点が 4 になるようなケースがあったとすれば、これは明らかに異常であり、残差が正規分布していないのか、このケースだけが例外的（つまり外れ値）なのか、である。 残差のプロットをするときに、残差をZ得点に変換してからプロットすれば、どれが外れ値かわかるし、残差の分散や残差の期待値が X の値にかかわらずゼロになるかもわかるので好都合である。
## id t x y
## 1 : 5 1:100 Min. :12.00 Min. : 20.2
## 2 : 5 2:100 1st Qu.:18.50 1st Qu.:333.8
## 3 : 5 3:100 Median :20.20 Median :401.4
## 4 : 5 4:100 Mean :20.04 Mean :409.6
## 5 : 5 5:100 3rd Qu.:21.43 3rd Qu.:478.0
## 6 : 5 Max. :25.90 Max. :744.3
## (Other):470par(mfrow=c(1,2))
plot(predict(pl1), scale(as.numeric(resid1)),
ylab="残差の Z 得点") #scale()はベクトルをZ得点に変換。as.numeric() ベクトルを単純な数値に
plot(predict(pl2), scale(as.numeric(residuals(pl2)
)
),
ylab="残差の Z 得点")
図8: 残差の Z得点のプロット
図8 の左側の例は分散の均一性には問題があるが、残差の期待値には問題がみられず、残差のZ得点は おおむね -3 ～ 3 のあいだで分布しており、\(N=500\) であることを考えると特に不自然な外れ値は存在しない。一方、右側の図は -4 以下と 4 位上の値があり、外れ値と言える。
残差の外れ値は、従属変数を対数変換することで外れ値ではなくなることもある。また、サンプルサイズが大きければ外れ値を除外してもしなくてもパラメータの推定値にほとんど変化はない場合も多い。また、本当に偶然、おおきな残差が生じるという可能性も否定しきれない。外れ値を見つけた場合は、ケースバイケースで柔軟な対応したい。
3.2 Cook Distance: 独立変数の外れ値
回帰分析では、独立変数の分布はどのようなものであっても、OLS が BLUE (Best Linear Unbiased Estimator) であることには影響がない。それゆえ、独立変数に外れ値があることに対して、あまり神経質になる必要はないように思える。特にサンプル・サイズが大きい場合、外れ値がパラメータの推定に及ぼす影響は、経験則から言うと、小さいことが多い。
しかし、国や都道府県を単位とした分析の場合 \(N\) は小さく、例外的な国や地域（あるいは特定の時点）が分析全体に影響をおよぼすことはある。その場合、やはり外れ値を除外すべきか検討する必要があるかもしれない。独立変数が一つならば、外れ値を特定するのはそれほど難しいことではない。極端に大きな（あるいは小さな）値をとっている値があればそれが外れ値であり、そのケースを除外して分析したら、どの程度回帰係数が変化するか、見てみればいいのである。例えば以下のように。
head(d5)## X Y
## 1 1 3.202208
## 2 2 4.602651
## 3 3 4.933687
## 4 4 4.704860
## 5 5 4.186925
## 6 6 2.766296par(mfrow=c(1, 2))
plot(Y ~ X, data=d5)
lm4.1 <- lm(Y ~ X, data=d5)
abline(lm4.1, col="red")
lm4.2 <- update(lm4.1, subset= X < 90) # subset = 論理式 は論理式がTRUEになるケースだけを使って推定
abline(lm4.2, col="blue", lty=2)
text(70, predict(lm4.1, data.frame(X=70)) + 1.5,
"すべてのサンプルから\n 計算した回帰直線")
text(70, predict(lm4.2, data.frame(X=70)) - 1, "外れ値を除いた時の\n 回帰直線")
library(memisc)
mtable(lm4.1, lm4.2, summary.stats=c("N", "adj. R-squared", "AIC"))##
## Calls:
## lm4.1: lm(formula = Y ~ X, data = d5)
## lm4.2: lm(formula = Y ~ X, data = d5, subset = X < 90)
##
## ===================================
## lm4.1 lm4.2
## -----------------------------------
## (Intercept) 3.626*** 3.900***
## (0.137) (0.193)
## X 0.034*** 0.007
## (0.009) (0.016)
## -----------------------------------
## N 100 99
## adj. R-squared 0.121 -0.009
## AIC 274.630 268.954
## ===================================plot(d5$X, residuals(lm4.1))
\(X=100\) のケースをサンプルから除外すると、かなり劇的に係数が変化するので、このような場合検討が必要であろう。
よく、どの程度大きな（あるいは小さな）値をとったら、「外れ値」になるのか、という質問をされるが、特に定義はない。むしろ上の例のように、そのケースを除外したとき、どの程度、パラメータの推定値が変化するのか、に注目すべきであろう。そのような指標としてクック距離 (Cook’s Distance) を紹介しておこう。\(i\)番目のケースのクック距離 (\(D_i\)) は、推定したパラメータの数を \(k\)、\(i\)番目のケースをサンプルから除外して同じモデルを推定した時の\(j\)番目のケースの\(Y_j\)の予測値を\(Y_{j(-i)}\)とすると、 \[ D_i =\frac{\sum_{j=1}^N(\hat Y_j - \hat Y_{j(-i)})}{k * \mathrm{var}(\epsilon_i)} \qquad(j \neq i) \] と定義される。分子は\(i\)番目のケースを除外したときモデルからの予測値の変化を\(i\)以外のすべてのケースに関して足し合わせ、それをモデルの残差の分散と推定パラメータの数で割ったものである。クック距離がどの程度大きければ外れ値として注目すべきなのかは、いくつか説があるようであるが、Cook 自身は 1 以上という基準を示しているようである。
R の場合、
cooks.distance(lm() または glm() の結果)
でクック距離をすべてのケースに関して計算してくれる。上の例で、\(X=100\)の事例に関して Cook 距離を計算すると、以下のようになる。
Di <- cooks.distance(lm4.1)
Di[d5$X > 90] # Xが90より多いケースのクック距離## 40
## 4.817462Di[Di > 1] # クック距離が 1 より大きい事例だけ取り出す## 40
## 4.817462summary(Di) # クック距離の分布を確認## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000000 0.000845 0.003243 0.054980 0.008737 4.817000plm() の結果に関してクック距離を計算する関数は見当たらないので作ってみた。
# plm() class のオブジェクトからクック距離を計算
p.cooks.disnance <- function(model.plm){ # m は plm の結果
X <- model.matrix(model.plm) # 独立変数（とその交互作用などすべての項）のデータ行列
h <- X %*% solve(t(X) %*% X) %*% t(X) # hat 行列
h.ii <- diag(h) # テコ比
resid <- model.plm$residuals # 残差
v <- sum(resid^2) / model.plm$df.residual # 残差分散
k <- length(coef(model.plm)) # 推定したパラメータの数
h2 <- h.ii / (1 - h.ii)^2
r2 <- resid^2/ (k * v)
return(r2 * h2)
}
plot(p.cooks.disnance(pl2))
4 Multicolinearity 多重共線性
多重共線性は、日本の社会学者が唯一こだわるモデルの診断であり、VIF や独立変数のあいだの相関の強さに気を使う社会学者は多い。確かに多重共線性は標準誤差を大きくするので、注意が必要である。極端に大きな標準誤差がでている場合には、気をつけたほうがよい。
4.1 correlation 独立変数間の相関係数をみるだけではダメ
独立変数の相関しか気にしない社会学者はもう少数派だと思うが、多重共線性を気にするならば独立変数の相関だけ見ていても正確なことはわからない（が、独立変数の間の相関関係は探索的な分析としては役に立つのでやっておきたい）。以下のような架空のデータの例を考えよう。
# 架空の多重共線性のデータを作成
x1 <- 1:100
x2 <- -x1 + rnorm(100, sd=50)
x3 <- x1 + 0.5*x2 + rnorm(100)
y <- 2*x1 + 0.4*x2 + 0.3*x3 + rnorm(100, sd=20)
round(
cor(
cbind(x1, x2, x3)
),
2)## x1 x2 x3
## x1 1.00 -0.45 0.57
## x2 -0.45 1.00 0.48
## x3 0.57 0.48 1.00lm31 <- lm(y ~ x1 + x2 + x3)
summary(lm31)##
## Call:
## lm(formula = y ~ x1 + x2 + x3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -47.324 -15.270 2.315 14.854 44.089
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.7775 4.0579 0.931 0.354
## x1 -0.9235 1.8958 -0.487 0.627
## x2 -1.0609 0.9459 -1.122 0.265
## x3 3.1178 1.8976 1.643 0.104
##
## Residual standard error: 20.05 on 96 degrees of freedom
## Multiple R-squared: 0.8921, Adjusted R-squared: 0.8887
## F-statistic: 264.5 on 3 and 96 DF, p-value: < 2.2e-16library(car) # car パッケージの vif 関数を使う。##
## Attaching package: 'car'
##
## The following object is masked from 'package:memisc':
##
## recodevif(lm31)## x1 x2 x3
## 744.6258 650.6110 765.9645x1, x2, x3 の相関はけっこうあると言えるが、相関係数が 0.7 以下なら多重共線性は問題ない、と書いてある教科書もあるので、特に問題のある値とはいえない。しかし、回帰分析してみると、係数が異常な値をとっている。決定係数が 0.8 以上あるのに、変数はほとんど有意でなく、標準化係数を計算すると、絶対値が 1 以上になる。VIF を計算すると 500 以上の値をとっており、典型的な多重共線性を示している。
4.2 Definition: VIF の定義
多重共線性の診断には VIF (Variance Inflation Factor) がよく用いられる。VIF の計算は次のようなステップで行われる。
- 独立変数の中から、VIFを計算する独立変数を一つ選ぶ。これを \(X_k\) と呼んでおく。
- \(X_k\) を従属変数、残りの独立変数を独立変数とした重回帰分析を行う。
- 上の結果得られた決定係数を \(R^2_k\) とすると、 \[ \mathrm{VIF} = \frac{1}{1- R^2_k} \] で、VIF は定義される。つまり、\(R^2_k=0\) のとき VIF\(=0\) であり、\(R^2_k\) が 1 に近づくにつれて VIF は無限大に発散していく。個々の独立変数間の相関があまり大きくなくても 多重共線性が起きうるのは、多重共線性で本質的な問題は、ゼロ次の相関ではなく、重相関だからである。それゆえ VIF をチェックした方がいい。
4.3 VIF の値はいくつ以上なら問題なのか？
よく「VIF の値はいくつ以下なら問題ですか？」といった質問をされるが、この質問に対する「正解」は存在しない。VIFが大きくても、サンプルサイズが大きければ、標準誤差は小さくなるからである。上の例のように VIFが1000に近い極端な例でも、サンプルサイズが10万ぐらいあればかなり正確な推定値が得られる。例えば以下のように。
x1 <- rep(1:100, 1000)
x2 <- -x1 + rnorm(100000, sd=50)
x3 <- x1 + 0.5*x2 + rnorm(100000)
y <- 2*x1 + 0.4*x2 + 0.3*x3 + rnorm(100000, sd=20)
round(
cor(
cbind(x1, x2, x3)
),
2)## x1 x2 x3
## x1 1.0 -0.5 0.5
## x2 -0.5 1.0 0.5
## x3 0.5 0.5 1.0lm32 <- lm(y ~ x1 + x2 + x3)
summary(lm32)##
## Call:
## lm(formula = y ~ x1 + x2 + x3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -91.809 -13.519 0.004 13.496 92.923
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.19161 0.12764 -1.501 0.133
## x1 2.02625 0.06345 31.935 < 2e-16 ***
## x2 0.41175 0.03173 12.977 < 2e-16 ***
## x3 0.28019 0.06341 4.419 9.93e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.03 on 99996 degrees of freedom
## Multiple R-squared: 0.8923, Adjusted R-squared: 0.8923
## F-statistic: 2.763e+05 on 3 and 99996 DF, p-value: < 2.2e-16vif(lm32)## x1 x2 x3
## 836.1186 835.7863 835.0484最初の例でもこの例でも、x1, x2, x3 の真の傾きはそれぞれ 2, 0.4, 0.3 なのだが、サンプルを100から10万に増やせば、かなり正確な推定ができていることがわかるだろう。
私の場合、数百から数千のサンプルを扱うことが多いが、その場合、VIFが10ぐらいあっても特に問題はないと感じている。もちろん、VIFが大きくなるほど標準誤差は大きくなっていくので、VIF は小さいほど正確な推定になる（そして有意な結果が出やすくなる）のであるが、VIFが特定の値以上だから回帰分析の結果を信じられない、という考え方はあまりに短絡的といえよう。VIFの大きさは標準誤差に反映されるので、係数の標準誤差さえちゃんと見ていれば大きな失敗はふつうしないものである（それゆえ、標準誤差を論文に記載しないことのほうがよっぽど問題である）。むしろ致命的なのは、モデルの特定を誤ることであり、まじめに適切なモデルを探索・診断していくことの重要性を強調しておきたい。たまに VIF が大きくなるという理由だけで、交互作用効果をモデルに投入することを嫌がる人を見かけるが、本末転倒であろう。
4.4 VIF for plm()
car パッケージの vif() は plm クラスのオブジェクトでは正しい計算をしないようである。そこで、自分で VIF を計算する関数 p.vif() を作ってみた。
# plm() class のオブジェクトから VIF を計算
p.vif <- function(model.plm){ # m は plm の結果
X <- model.matrix(model.plm) # 独立変数（とその交互作用などすべての項）のデータ行列
k <- ncol(X) # パラメータの数
r2 <- numeric(ncol(X)) # 決定係数を格納する
names(r2) <- colnames(X)
for(i in 1:k){
r2[i] <- summary(lm(X[,i]~., data=as.data.frame(X[,-i])) )$r.squared
}
return(1/(1-r2))
}
pl3 <- plm(Y ~X1 + X2 + X3, data= p.d6)
round(cor(p.d6[, c("X1", "X2", "X3")]), 2)## X1 X2 X3
## X1 1.00 -0.46 0.83
## X2 -0.46 1.00 -0.44
## X3 0.83 -0.44 1.00p.vif(pl3) # plm クラス・オブジェクトを引数にする## X1 X2 X3
## 2.651930 1.027633 2.623106