一般最小二乗法（OLS）では、独立で同一の分布に従う誤差を仮定していました。しかし、観測値の誤差がすべて同一の分布に従うと見做せないことは多いです。不均一分散です。
集計データでは大きく異なる人口の行政区が同時に集計されていますが、行政区ごとの平均値の標準誤差は人口によって異なります。観察できない偶発的な事象が、ある集団の間で同じように影響する一方、他の集団には影響しないこともあります。個体ごとに複数の状態を考えるとき、ある個体ごとの観察できない状態は、その個体のすべての状態に影響します。無作為に介入群と対称群を割り当てる実験でも、介入の有無でサブサンプルの分布が異なってしまうことはあります。
不均一分散は二つの意味で問題を引き起こします。一つは、有効推定量にならないことです。小サンプルにおいて点推定に信頼がおけないことになります。一つは、標準誤差にバイアスが入ることです。区間推定と統計的仮説検定に信頼がおけないことになります。大標本でパラメーターの点推定量だけが問題であればさほど気にしなくてもよいのですが、小標本であるなど不均一分散が結果に大きく影響しうる場合で、かつその影響の緩和策が知られている場合は対処する方が望ましいです。
1 一般化最小二乗法（GLS）
不均一分散を前提とした最小二乗法を一般化最小二乗法（GLS）と呼びます。
\[ y = X \beta + \Gamma \epsilon \]
\(y\)は従属変数のベクトル、\(X\)は説明変数の行列、\(\epsilon\)はOLSと同様のi.i.dの誤差項です。\(\Gamma\)が不均一分散をあらわす行列となります。\(\Gamma\)をどう仮定するかで、様々な不均一分散の入り方を表現することができます。
例えばサンプルサイズ\(4\)、説明変数の数を切片項も入れて\(3\)、1行目と2行目の観測値と、3行目と4行目の観測値がそれぞれ同一集団で、観測値は属する集団の他の観測値に影響している偶発的事象にも影響されるとすると、
\[ \begin{pmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \end{pmatrix} = \begin{pmatrix} 1 & x_{11} & x_{12} \\ 1 & x_{21} & x_{22} \\ 1 & x_{31} & x_{32} \\ 1 & x_{41} & x_{42} \end{pmatrix}\begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \end{pmatrix} + \begin{pmatrix} 1 & \gamma_{12} & 0 & 0 \\ \gamma_{22} & 1 & 0 & 0 \\ 0 & 0 & 1 & \gamma_{34} \\ 0 & 0 & \gamma_{43} & 1 \end{pmatrix}\begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \epsilon_3 \\ \epsilon_4 \end{pmatrix} \]
と書くことができます。1行目の観測値は
\[ y_1 = \beta_0 + \beta_1 x_{11} + \beta_2 x_{12} + 1\cdot\epsilon_1 + \gamma_{12}\epsilon_2 + 0\cdot\epsilon_3 + 0\cdot\epsilon_4 \]
になるわけですが、自身に生じた誤差（と言うかモデル外の要因）の影響を\(1\)だけ受け、同一集団の2行目の観測値に生じた誤差の影響を\(\gamma_{12}\)だけ受け、異なる集団の3行目と4行目に生じた誤差の影響は受けないことになっているのがわかります。
推定は\(\Gamma\)がわかっていて、かつ正則であれば簡単です。逆行列\(\Gamma^{-1}\)を両辺にかければ
\[ \Gamma^{-1} y = \Gamma^{-1} X \beta + \epsilon \]
となりますから、OLSの計算がそのまま使えます。係数の推定量は
\[ \hat{\beta} = ({}^t\! X {}^t\! \Gamma^{-1} \Gamma^{-1} X)^{-1} {}^t\! X {}^t\! \Gamma^{-1} \Gamma^{-1} y \]
となり、ここで\(\Omega = {}^t\! \Gamma^{-1} \Gamma^{-1}\)と置いて、
\[ \hat{\beta} = ({}^t\! X \Omega X)^{-1} {}^t\! X \Omega y \]
と書けます。標準誤差も同様に導出できます。一般化と称されているだけあって、線形回帰モデルに入り込む不均一分散はすべてこの形で表現することができます。
ただし、仮定した\(\Gamma\)を観測できるかは別です。一般には観測できず、そもそも不均一分散の入り方を直接観察することもできないため、不均一分散への完全な対処方法はありません。多かれ少なかれ分析者が置く仮定に依存することになります。
2 ウェイト付き最小二乗法（WLS）
もっとも単純な不均一分散は、観測値ごとの誤差の分散が異なると言うものです。この場合、
\[\begin{align} \Gamma = \begin{pmatrix} \gamma_{11} & 0 & \cdots & 0 & 0 \\ 0 & \gamma_{22} & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & \gamma_{(n-1)(n-1)} & 0 \\ 0 & 0 & \cdots & 0 & \gamma_{nn} \end{pmatrix} \end{align}\]
と書け、だいたい\(\gamma_{ii}\ne\gamma_{jj}\ (i \ne j)\)となります。両辺に乗じるのは\(\Gamma^{-1}\)ですが、この対角成分をウェイトと呼びます。
\(\Gamma\)が分からないので、このままでは推定できません。何かの仮定をおいて\(\Gamma\)を推定します。
2.1 観察可能な変数をウェイトにする
何か観察可能な変数（を単純に変換したもの）と\(\Gamma\)が比例するのであれば、それをウェイトにする事ができます。行政区ごとの平均値を回帰する場合であれば、行政区の人口\(N_i\)の平方根に比例していると考えてよい場合もあるでしょう1。その場合は、\(\gamma{ii}=1/\sqrt{N_i}\)と置き、ウェイトを求めます。
2.2 残差を推定して予測値からウェイトをつくる
観察可能な変数をどう変換すればよいのか分からない場合は、OLSの残差の大きさを推定し、予測値からウェイトをつくります。FGLSと呼ばれる手法になります。
残差の大きさの予測モデルは様々なものが考えられますが、ここではOLSの残差\(\rho\)の二乗を対数化し、説明変数で回帰します。手間がかかりますが、ウェイトが非負になることが保証されます。
\[\begin{align} \log{\rho^2} &= x_i \alpha + \eta_i\\ \hat{\rho_i} &= \sqrt{\exp(x_i \alpha)} \end{align}\]
ここでの\(x_i\)は説明変数の横ベクトルで、\(\alpha\)は係数の縦ベクトル、\(\eta\)は誤差項です。
\(\hat{\rho_i}\)を\(\gamma_{ii}\)として用います。
2.2.1 データセット
不均一分散の入ったデータセットをつくります。
set.seed(1234)
df01 <- with(new.env(), {
n <- 50
x <- runif(n, 0, 1)
z <- c(rep(1, n - 3), rep(10, 3))
y <- 1 + 2*x + 0*z + rnorm(n, sd = z)
data.frame(y = y, x = x, z = z)
})推定に用いるモデルとは、誤差の作り方が異なることに注意してください。
2.2.2 行列演算による推定
早速、推定してみましょう。
model <- y ~ x + z
mf <- model.frame(model, df01)
y <- model.response(mf)
X <- model.matrix(mf, mf)
# OLS
beta <- solve(t(X) %*% X) %*% t(X) %*% y
sigma <- sqrt(sum((y - X %*% beta)^2)/(nrow(X) - ncol(X)))
ols_vcov <- sigma^2 * solve(t(X) %*% X)
# WLS
R <- y - X %*% beta
P <- X %*% solve(t(X) %*% X) %*% t(X) %*% log(R^2)
W <- 1/sqrt(exp(c(P)))
Gamma <- diag(1/W)
Omega <- t(solve(Gamma)) %*% solve(Gamma)
w_beta <- solve(t(X) %*% Omega %*% X) %*% t(X) %*% Omega %*% y
w_sigma <- sqrt(sum((W*(y - X %*% w_beta))^2)/(nrow(X) - ncol(X)))
w_vcov <- w_sigma^2 * solve(t(X) %*% Omega %*% X)
# Comparing WLS with OLS
(CMP <- matrix(c(c(1, 2, 0), beta, sqrt(diag(ols_vcov)), w_beta, sqrt(diag(w_vcov))), ncol = 5,
dimnames = list(colnames(X), c("Coef.(DGP/True)", "Coef.(OLS)", "S.E.(OLS)", "Coef.(WLS)", "S.E.(WLS)")))) Coef.(DGP/True) Coef.(OLS) S.E.(OLS) Coef.(WLS) S.E.(WLS)
(Intercept) 1 -0.7607539 0.8307883 -0.1264285 0.6529731
x 2 3.2196465 1.4706592 2.3943684 0.6625233
z 0 0.8993327 0.1797359 0.6736195 0.5921787OLS推定量よりも、WLS推定量の方がデータ生成プロセスの係数に近いです。
2.2.3
lmを使った推定
実践ではビルトイン関数lmを使いましょう。短く記述でき、計算時間もかかりません。
r_lm <- lm(y ~ x + z, df01)
r_lgr <- lm(log(residuals(r_lm)^2) ~ x + z, df01)
r_wls <- lm(model, df01, weight = 1/exp(predict(r_lgr)))
coef(summary(r_wls)) Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.1264285 0.6529731 -0.1936198 0.8473079596
x 2.3943684 0.6625233 3.6140136 0.0007322129
z 0.6736195 0.5921787 1.1375274 0.26108385192.3 不均一分散の存在の検定
無視できないほど大きなWLSで対処できる不均一分散がありそうか検定をかけることができます。
OLSで推定して、分散を計算しておきます。
教科書的には残差と分散の比がχ二乗分布に従うことを使うBreusch-Pagan検定を用います。
g <- residuals(r_lm)^2 / sigma2
TSS <- sum((g - mean(g))^2) # mean(g)は1
r_g <- lm(g ~ x + z, df01)
RSS <- sum(residuals(r_g)^2)
(LM <- (TSS - RSS)/2)[1] 273.2341[1] 4.65542e-60実用ではどうもstudentized Breusch-Pagan testの方が標準的に使われており、こちらは残差の二乗と分散の差の二乗がχ二乗分布に従うことを使っています。
w <- residuals(r_lm)^2 - sigma2
TSS <- sum((w - mean(w))^2) # mean(w)は0
r_w <- lm(w ~ x + z, df01)
RSS <- sum(residuals(r_w)^2)
(s.LM <- length(residuals(r_lm)) * (TSS - RSS) / sum(w^2))[1] 32.86126[1] 7.315922e-083 見かけ上無関係な回帰（SUR）
WLSより手の混んだGLSに分類される手法としてSURを紹介します。
以下の連立方程式を同時推定することを考えます。
\[ \begin{eqnarray} \left\{ \begin{array}{l} y_{1} = X_{1}\beta_1 + \rho_{1} \\ y_{2} = X_{2}\beta_2 + \rho_{2} \end{array} \right. \end{eqnarray} \]
\(y_1\)と\(y_2\)は従属変数のベクトルで、\(X_{1}\)と\(X_{2}\)は説明変数の行列、\(\beta_1\)と\(\beta_2\)はそれらの係数のベクトル、\(\rho_1\)と\(\rho_2\)は誤差項のベクトルです。
\(\rho_1\)と\(\rho_2\)は同一の分布に従うとは限りません。\(X_{1}\)と\(X_{1}\)は同じ変数が入っていてもよいですし、\(\beta_1\)と\(\beta_2\)も同じ値の部分があると仮定することもできます。
二つの式はそれぞれはOLSの仮定を満たしているので、一つ一つ推定しても良い気がしますが、\(\mbox{COV}(\rho_1, \rho_2) \neq 0\)の場合は推定量を改善する余地があります。
最初は\(\mbox{COV}(\rho_1, \rho_2)\)の情報が無いので、以下をOLSで推定します。
\[ \begin{pmatrix} y_{1} \\ y_{2} \end{pmatrix} = \begin{pmatrix} X_{1} & 0 \\ 0 & X_{2} \end{pmatrix} \begin{pmatrix} \beta_{1} \\ \beta_{2} \end{pmatrix} + \begin{pmatrix} \rho_1 \\ \rho_2 \end{pmatrix} \]
今回は\(\beta_1\)と\(\beta_2\)に共通の値はないと仮定しています。
上の式をOLSで推定すると残差を計算できるので、以下の分散共分散行列\(\Sigma\)がつくれます。
\[ \Sigma = \begin{pmatrix} \mbox{COV}(\rho_1, \rho_1) & \mbox{COV}(\rho_1, \rho_2) \\ \mbox{COV}(\rho_1, \rho_2) & \mbox{COV}(\rho_2, \rho_2) \end{pmatrix} \]
\(\rho_1\)と\(\rho_2\)は期待値\(0\)で分散\(1\)の誤差項のベクトル\(\epsilon_1\)と\(\epsilon_2\)の間とは以下のような関係があります。
\[ \begin{pmatrix} \rho_1 \\ \rho_2 \end{pmatrix} = \Sigma^\frac{1}{2} \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \end{pmatrix} \]
両辺の分散共分散行列をつくると一致するからです。
\[\begin{align} \frac{1}{n} \begin{pmatrix} \rho_1 \\ \rho_2 \end{pmatrix}\begin{pmatrix} \rho_1 & \rho_2 \end{pmatrix} &= \frac{1}{n} \Sigma^\frac{1}{2} \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \end{pmatrix} \begin{pmatrix} \epsilon_1 & \epsilon_2 \end{pmatrix} {}^t\! \Sigma^\frac{1}{2} \\ &= \frac{1}{n} \Sigma^\frac{1}{2} \begin{pmatrix} n & 0 \\ 0 & n \end{pmatrix}{}^t\! \Sigma^\frac{1}{2} \\ &= \Sigma \end{align}\]
二つの式を合併して、以下のように表記を簡素化します。
\[ y = X\beta + \rho \]
\(\rho\)を\(\Sigma^{\frac{1}{2}} \otimes I_n \epsilon\)に置き換えると、GLSの形式にできます。
\[ y = X\beta + \Sigma^{\frac{1}{2}} \otimes I_n \epsilon \]
\(\otimes\)はクロネッカー積で、\(I_n\)は\(n\)行\(n\)列の恒等行列です。見慣れない演算子かも知れませんが、例えば\(n=3\)のときは、\(\Sigma^{\frac{1}{2}}\)の\(i\)行\(j\)列の成分を\(s_{ij}\)として、
\[\begin{align} \Sigma^{\frac{1}{2}} \otimes I_3 \epsilon &= \begin{pmatrix} s_{11} & s_{12} \\ s_{12} & s_{22} \end{pmatrix} \otimes \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{pmatrix} \begin{pmatrix} \epsilon_{11} \\ \epsilon_{21} \\ \epsilon_{31} \\ \epsilon_{12} \\ \epsilon_{22} \\ \epsilon_{32} \end{pmatrix}\\ &= \begin{pmatrix} s_{11} & 0 & 0 & s_{12} & 0 & 0 \\ 0 & s_{11} & 0 & 0 & s_{12} & 0 \\ 0 & 0 & s_{11} & 0 & 0 & s_{12} \\ s_{12} & 0 & 0 & s_{22} & 0 & 0 \\ 0 & s_{12} & 0 & 0 & s_{22} & 0 \\ 0 & 0 & s_{12} & 0 & 0 & s_{22} \end{pmatrix}\begin{pmatrix} \epsilon_{11} \\ \epsilon_{21} \\ \epsilon_{31} \\ \epsilon_{12} \\ \epsilon_{22} \\ \epsilon_{32} \end{pmatrix} \\ &= \begin{pmatrix} s_{11} \epsilon_{11} + s_{12} \epsilon_{12}\\ s_{11} \epsilon_{21} + s_{12} \epsilon_{22}\\ s_{11} \epsilon_{31} + s_{12} \epsilon_{32}\\ s_{12} \epsilon_{11} + s_{22} \epsilon_{12}\\ s_{12} \epsilon_{21} + s_{22} \epsilon_{22}\\ s_{12} \epsilon_{31} + s_{22} \epsilon_{32} \end{pmatrix} \end{align}\]
と展開されます。つまり、左の行列の要素に右の行列を乗じたものでブロック行列をつくって連結します。
連結前の1番目式の\(i\)行目の誤差項と2番目の式と\(i\)行目の誤差項が、それぞれ\(\epsilon_{i1}\)と\(\epsilon_{i2}\)の関数となっており、ゼロではない共分散を持つようにできたことがわかります。
\(\Omega\)は以下のようになります。テキストの定義ではクロネッカー積を乗じる前の逆行列になっていると思うので注意してください。
\[\begin{align} \Omega &= {}^t\!(\Sigma^\frac{1}{2} \otimes I_n)(\Sigma^\frac{1}{2} \otimes I_n) \\ &= (\Sigma \otimes I_n)^{-1} \\ &= \Sigma^{-1} \otimes I_n \end{align}\]
係数の分散共分散行列を計算します。頑強ではなく簡便な手法を用います。
\[\begin{align} \hat{\beta}_{SUR} - \beta &\approx ({}^t \! X \Omega X)^{-1} {}^t \! X \Omega y - ({}^t \! X \Omega X)^{-1} {}^t \! X \Omega (y - \rho) \\ &= ({}^t \! X \Omega X)^{-1} {}^t \! X \Omega \rho \\ (\hat{\beta} - \beta) {}^t\!(\hat{\beta} - \beta) &= ({}^t \! X \Omega X)^{-1} {}^t \! X \Omega \rho {}^t \! \rho {}^t\! \Omega X {}^t \! ({}^t \! X \Omega X)^{-1} \end{align}\]
方程式が同じであれば等分散を仮定しているので\(E[\rho {}^t \! \rho]\)は\(\Sigma \otimes I_n\)と見なすことができ、\(\Sigma\)と\(\Omega\)は対称行列であることに注意すると、
\[ \Sigma \otimes I_n {}^t\! \Omega = \Sigma \otimes I_n \Omega = (\Sigma \otimes I_n)(\Sigma \otimes I_n)^{-1} = I_{2n} \]
なので、
\[\begin{align} &\approx ({}^t \! X \Omega X)^{-1} {}^t \! X \Omega I_{2n} X {}^t \! ({}^t \! X \Omega X)^{-1} \\ &= ({}^t \! X \Omega X)^{-1} {}^t \! X \Omega X {}^t \! ({}^t \! X \Omega X)^{-1}\\ &= {}^t \! ({}^t \! X \Omega X)^{-1} = ({}^t \! X \Omega X)^{-1} \end{align}\]
になります。
3.1 データセット
二つの方程式を連結したデータセットを作成します。
set.seed(20181103)
df02 <- with(new.env(), {
n <- 50 # 標本サイズが小さい方が、OLSとSURの差が出る
l <- rnorm(n, sd=2) # 共分散を持つようにする
e1 <- rnorm(n, sd=1) + l
e2 <- rnorm(n, sd=1) + l
x1 <- runif(n, min=0, max=2)
x2 <- runif(n, min=-2, max=2)
y1 <- 1 + 2*x1 + e1
y2 <- 3 + 2*x2 + e2
X1 <- matrix(c(rep(1, n), x1, rep(0, 2*n)), n, 4)
X2 <- matrix(c(rep(0, 2*n), rep(1, n), x2), n, 4)
X <- rbind(X1, X2)
colnames(X) <- c("Const.1", "x1", "Const.2", "x2")
y <- c(y1, y2)
data.frame(y, X)
})3.2 行列による推定
実際のSUR推定量の計算を行います。browseVignettes("systemfit")をすると統計学的に詳しい説明も参照できます。
mf <- model.frame(y ~ 0 + Const.1 + x1 + Const.2 + x2, df02)
y <- model.response(mf)
X <- model.matrix(mf, mf)
n <- nrow(X)/2
# 1段階目のOLS推定量を計算
beta_ols <- solve(t(X)%*%X)%*%t(X)%*%y
# 誤差項を計算
r_e <- y - X %*% beta_ols
r_e1 <- r_e[1:n]
r_e2 <- r_e[(n+1):(2*n)]
# ∑を計算
df <- n - 2 # 自由度は複数の方程式の最も小さい値になる
Sigma <- matrix(c(sum(r_e1^2)/df, sum(r_e1*r_e2)/df, sum(r_e1*r_e2)/df, sum(r_e2^2)/df), 2, 2)
Omega <- solve(Sigma) %x% diag(n)
# SUR推定量を計算
beta_sur <- solve(t(X) %*% Omega %*% X) %*% t(X) %*% Omega %*% y
# 分散共分散行列を計算
(vcov_sur <- solve(t(X) %*% Omega %*% X)) # Usual Variance Matrix (Wooldridge (2002) p.161) Const.1 x1 Const.2 x2
Const.1 0.188358645 -9.421814e-02 8.287258e-02 0.0036769042
x1 -0.094218144 9.776368e-02 8.868483e-05 -0.0038152701
Const.2 0.082872582 8.868483e-05 1.024398e-01 -0.0004940602
x2 0.003676904 -3.815270e-03 -4.940602e-04 0.02125474313.3
systemfitで推定
実践においてはパッケージを使う方が無難です。
library(systemfit)
# データセットをほどく（通常は不要）
y1 <- y[1:n]
y2 <- y[(n+1):(2*n)]
X1 <- X[1:n, 1:2]
X2 <- X[(n+1):(2*n), 3:4]
# 連立方程式を書く
eq1 <- y1 ~ 0 + X1
eq2 <- y2 ~ 0 + X2
system <- list(eq1, eq2)
r_systemfit <- systemfit(system, method = "SUR") # 気づくと関数が統合されて使い方が微妙に変わっていた
summary(r_systemfit)
systemfit results
method: SUR
N DF SSR detRCov OLS-R2 McElroy-R2
system 100 96 484.321 7.34881 0.418894 0.713533
N DF SSR MSE RMSE R2 Adj R2
eq1 50 48 236.737 4.93203 2.22082 0.134333 0.116298
eq2 50 48 247.584 5.15800 2.27112 0.557864 0.548653
The covariance matrix of the residuals used for estimation
eq1 eq2
eq1 4.87787 4.14790
eq2 4.14790 5.12141
The covariance matrix of the residuals
eq1 eq2
eq1 4.93203 4.2533
eq2 4.25330 5.1580
The correlations of the residuals
eq1 eq2
eq1 1.000000 0.843282
eq2 0.843282 1.000000
SUR estimates for 'eq1' (equation 1)
Model Formula: y1 ~ 0 + X1
Estimate Std. Error t value Pr(>|t|)
X1Const.1 0.942839 0.434003 2.17242 0.034793 *
X1x1 1.993500 0.312672 6.37569 6.6711e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.220817 on 48 degrees of freedom
Number of observations: 50 Degrees of Freedom: 48
SSR: 236.73726 MSE: 4.932026 Root MSE: 2.220817
Multiple R-Squared: 0.134333 Adjusted R-Squared: 0.116298
SUR estimates for 'eq2' (equation 2)
Model Formula: y2 ~ 0 + X2
Estimate Std. Error t value Pr(>|t|)
X2Const.2 3.068725 0.320062 9.5879 9.861e-13 ***
X2x2 1.886858 0.145790 12.9423 < 2.22e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.271124 on 48 degrees of freedom
Number of observations: 50 Degrees of Freedom: 48
SSR: 247.584146 MSE: 5.158003 Root MSE: 2.271124
Multiple R-Squared: 0.557864 Adjusted R-Squared: 0.548653 eq1_X1Const.1 eq1_X1x1 eq2_X2Const.2 eq2_X2x2
eq1_X1Const.1 0.188358645 -9.421814e-02 8.287258e-02 0.0036769042
eq1_X1x1 -0.094218144 9.776368e-02 8.868483e-05 -0.0038152701
eq2_X2Const.2 0.082872582 8.868483e-05 1.024398e-01 -0.0004940602
eq2_X2x2 0.003676904 -3.815270e-03 -4.940602e-04 0.02125474314 計量経済学の変量効果モデル
計量経済学で言う変量効果モデル2もFGLSに分類されます。観測値ごとの誤差\(u_{it}\)に加えて、個体ごとの誤差\(c_{i}\)がある線形回帰モデル\(y_{it} = X\beta + c_{i} + u_{it}\)です。\(y\)は被説明変数、\(X\)は説明変数、\(\beta\)は係数、\(i\)は個体を示す添字、\(t\)は時点を示す添字になります。観察できない偶発的な事象が、ある個体の観測値で同じように影響する一方、他の個体には影響しないことを表しています。
例えば個体ごとに\(5\)期間のデータがある場合、分散共分散\(\Sigma\)は以下のように表されます。
\[ \Sigma = \begin{bmatrix} \sigma^2_u + \sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_c \\ \sigma^2_c & \sigma^2_u + \sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_c \\ \sigma^2_c & \sigma^2_c & \sigma^2_u + \sigma^2_c & \sigma^2_c & \sigma^2_c \\ \sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_u + \sigma^2_c & \sigma^2_c \\ \sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_u + \sigma^2_c \\ \end{bmatrix} \]
このとき、\(n\)を個体の数として、\(\Omega = I_n \otimes \Sigma^{-1}\)とすると、
\[ \hat{\beta}_{RE} = (^t\!X \Omega X)^{-1} {^t\!X} \Omega y \]
と推定することができます。計量モデルと言うよりはデータの並びの問題で、クロネッカー積の左右がSURと逆なことに注意してください。
4.1 データセット
モデルに沿って不均一分散が入ったデータセットをつくります。
set.seed(1303)
mkdata <- function(n_of_t, n_of_i){
# n_of_t: 時点の数, n_of_i: 個体の数
# 個体番号をつくる
id <- rep(1:n_of_i, each = n_of_t, min=-0.5, max=0.5)
# 時点変数をつくる
t <- rep(1:n_of_t, n_of_i)
# 説明変数x,zをつくる
x <- runif(n_of_i*n_of_t)
z <- runif(n_of_i*n_of_t)
# 誤差項をつくる
e_ind <- rep(rnorm(n_of_i, sd = 2), each = n_of_t)
e_ids <- rnorm(n_of_i*n_of_t)
e <- e_ind + e_ids
# 従属変数yをつくる
y <- 1 + 2*x - 3*z + e
data.frame(id = as.factor(id), t = as.factor(t), y, x, z)
}
df03 <- mkdata(5, 100)4.2 行列による推定
定番テキストのWooldridgeのGSLとRandom Effect Modelの説明に沿って計算します。ただし、\(\Omega\)の定義が異なるので注意してください。
frml <- y ~ x + z
X <- model.matrix(frml, data = df03)
mf <- model.frame(frml, data = df03)
y <- model.response(mf)
beta_ols <- solve(t(X)%*%X) %*% t(X) %*% y
residuals <- X %*% beta_ols - y
df <- nrow(X) - ncol(X)
sigma2_ols <- sum(residuals^2)/df
T <- length(levels(df03$t))
n <- length(levels(df03$id))
# total error variance
sigma2_v <- sum(residuals^2)/df
# individual error variance
sigma2_c <- with(df03, {
sum(sapply(levels(id), function(i){
r <- residuals[id == i]
cp <- outer(r, r)
sum(cp[upper.tri(cp, diag = FALSE)])
# 以下と同じ結果
# sum <- 0
# for(t in 1:(T-1)){
# for(k in (t+1):T){
# sum <- sum + r[t]*r[k]
# }
# }
# sum
}))/((n*T*(T-1)/2) - ncol(X))
})
# idiosyncratic error variance
sigma2_u <- sigma2_v - sigma2_c
sigma <- matrix(sigma2_c, T, T) + sigma2_u*diag(T)
sigma_inv <- solve(sigma)
omega <- diag(n) %x% sigma_inv
beta_gls <- solve(t(X) %*% omega %*% X) %*% t(X) %*% omega %*% y
# Usual Variance Matrix Estimator
usual_vcov <- solve(t(X) %*% omega %*% X)
# Robust Variance Matrix Estimator
rbst_vcov <- with(df03, {
dfs <- split(df03, id)
A <- matrix(0, ncol(X), ncol(X))
B <- matrix(0, ncol(X), ncol(X))
n <- length(dfs)
for(i in 1:n){
mf <- model.frame(frml, data = dfs[[i]])
X <- model.matrix(mf, mf)
y <- model.response(mf)
u <- y - X %*% beta_gls
A <- A + t(X) %*% sigma_inv %*% X
B <- B + t(X) %*% sigma_inv %*% u %*% t(u) %*% sigma_inv %*% X
}
solve(A) %*% B %*% solve(A)
})自由度がややこしいですね。
4.3
plmで推定
比較のために、plmパッケージで推定します。ぐっと手間暇が減っていますね。
library(plm)
pd01 <- pdata.frame(df03, index=c("id"))
r_plm <- plm(frml, model = "random", data = pd01)パネルデータ分析に使えるパッケージは多数あり、plmにこだわる必要はありません。一方、plmのヴィネットはパネルデータ分析を行う前に一読を推奨します。
4.4 推定結果を比較
係数の推定量と標準誤差を比較してみましょう。
(beta <- matrix(
c(1, 2, -3, beta_ols, beta_gls, coef(summary(r_plm))[, 1]),
ncol(X), dimnames = list(rownames(beta_gls), c("true", "OLS", "GLS", "plm/random")))) true OLS GLS plm/random
(Intercept) 1 0.7393101 1.145877 1.147142
x 2 2.6266132 2.086162 2.084467
z -3 -2.7875092 -3.041145 -3.041919(sigma <- matrix(
c(sqrt(sigma2_ols*diag(solve(t(X) %*% X))), sqrt(diag(usual_vcov)), sqrt(diag(rbst_vcov)), coef(summary(r_plm))[, 2]),
ncol(X), dimnames = list(rownames(beta_gls), c("OLS", "GLS/Usual", "GLS/Robust", "plm/random")))) OLS GLS/Usual GLS/Robust plm/random
(Intercept) 0.2877172 0.2514105 0.2571024 0.2508719
x 0.3630824 0.1858884 0.2043492 0.1823705
z 0.3511725 0.1754158 0.1586268 0.1720832FGLSで推定した係数は、OLSよりも真の値に近くなっており、係数も標準誤差もplmパッケージとほとんど同じです。
ほとんど同じと言うところに引っかかると思いますが、plmパッケージでは以下の連立方程式を解いて個体内変化と個体間相違を計算して、idiosyncratic error varianceとindividual error varianceにしており、分散分解の方法が教科書のものと異なっているためです。
- \(\mbox{withinの残差平方和} = \mbox{個体内変化}((T-1)*N-K+1) + \mbox{個体間相違}*0\)
- \(\mbox{betweenの残差平方和} = \mbox{個体内変化}*(N-K)/T + \mbox{個体間相違}*(N-K)\)
なお、定番パッケージの作成者は研究動向に詳しい統計学者である事が多いので、実用面ではパッケージの手法を信じた方が無難です。
5 まとめ
GLSの例として、WLS、SUR、計量経済学のランダム効果モデルをざっと説明してきましたが、一般化と言いつつ不均一分散の構造を特定してから話がはじまることが分かりました。誤差の構造と推定方法は他にも様々あり、もっと高度な手法も広く用いられています。
実用では特殊ケースに応じたパッケージに頼る事になるので、もっと高度な手法でもコーディングの手間暇はほとんどありませんが、どういう構造を想定しているかを把握することが結果の解釈では重要になります。分かりやすい手法から、構造をしっかり把握していきましょう。