はじめに
千葉大学・株式会社Nospareの川久保です．前回の記事で，パネルデータ分析の基本的なモデルと推定法について紹介しました．今回は，各種の検定法と，Rによるデータ解析例をご紹介します．
特定化の検定
パネルデータ分析の3つの基本的なモデル
- Pooled OLS
- 変量効果モデル
- 固定効果モデル
のどれを採用するか，モデル特定化のための検定法を紹介します．
Pooled OLS vs. 固定効果モデル
固定効果モデルは，
$$
y_{it} = x_{it}^\top \beta + b_i + \varepsilon_{it}, \quad (i=1,\dots,N; \ t=1,\dots,T),
$$
といったモデルですが，固定効果がすべての$i$で等しいとき，つまり$b_1 = b_2 = \dots = b_N$の場合は，定数項を入れたPooled OLSで十分です．そこで，
\begin{align}
& H_0: b_1=b_2=\dots=b_N \\
& H_1: \mathrm{not} \ H_0
\end{align}
を検定したいですが，固定効果モデルに$N-1$個の制約を入れると，定数項を入れたPooled OLSになるため，以下の検定統計量を用いたF検定が利用できます．
F = {(\mathrm{RSS}_\mathrm{POLS} - \mathrm{RSS}_\mathrm{FE}) \ / \ (N-1) \over \mathrm{RSS}_\mathrm{FE} \ / \ (NT - N - K)}
ただし，$\mathrm{RSS}_\mathrm{POLS}$は，定数項を入れたPooled OLSの残差2乗和，$\mathrm{RSS}_\mathrm{FE}$は固定効果モデルの残差2乗和，$K$は$x_{it}$の次元（説明変数の個数）です．
固定効果モデル vs. 変量効果モデル
固定効果モデルと変量効果モデルでは，変量効果モデルの方がパラメータ数が少なく，効率的な推定ができます．一方で，$E[b_i \mid X_i] \not= 0$のときには，変量効果モデルの推定では，回帰係数を一致推定できなくなってしまうため，固定効果モデルの推定を行うのでした．そこで，
\begin{align}
& H_0: E[b_i \mid X_i] = 0 \\
& H_1: E[b_i \mid X_i] \not= 0
\end{align}
の検定を考えます．帰無仮説$H_0$のもとでは，固定効果モデルの推定量$\hat{\beta}_\mathrm{FE}$と変量効果モデルの推定量$\hat{\beta}_\mathrm{RE}$の両方が一致推定量となるため，両者は同じ値に確率収束します．この関係性を使って，Hausman検定と呼ばれる以下の検定統計量を用いた検定を利用します．
H = (\hat{\beta}_\mathrm{FE} - \hat{\beta}_\mathrm{RE})^\top [ \hat{\mathrm{Var}}(\hat{\beta}_\mathrm{FE}) - \hat{\mathrm{Var}}(\hat{\beta}_\mathrm{RE}) ]^{-1} (\hat{\beta}_\mathrm{FE} - \hat{\beta}_\mathrm{RE})
検定統計量$H$は帰無仮説$H_0$のもとで，自由度$K$のカイ二乗分布に分布収束します．帰無仮説$H_0$が棄却されれば，変量効果モデルの推定量は一致推定できないと考えられるため，固定効果モデルを採用します．
Rによるデータ解析例
plmというパッケージを用いたRによるデータ解析例を見ていきます．パッケージ内のGrunfeldという米国の企業パネルデータを読み込みます．
library(plm)
data(Grunfeld)
head(Grunfeld)でデータの冒頭を見てみます．
> head(Grunfeld)
firm year inv value capital
1 1 1935 317.6 3078.5 2.8
2 1 1936 391.8 4661.7 52.6
3 1 1937 410.6 5387.1 156.9
4 1 1938 257.7 2792.2 209.2
5 1 1939 330.8 4313.2 203.4
6 1 1940 461.2 4643.9 207.2
変数firmは企業のindexで$i=1,\dots,10$（すなわち$N = 10$），yearは西暦で1935年から1954年まで（すなわち$T = 20$），invは総投資，valueは企業価値，capitalは有形固定資産を表します．総投資を企業価値と有形固定資産に回帰するモデルを考え，Pooled OLS，固定効果モデル，変量効果モデルの3つのモデルを当てはめます．
res_p <- plm(inv ~ value + capital, data = Grunfeld, model = "pooling")
res_w <- plm(inv ~ value + capital, data = Grunfeld, model = "within")
res_r <- plm(inv ~ value + capital, data = Grunfeld, model = "random")
plm()がパネルデータ分析を実行する関数で，引数modelで当てはめるモデルを指定します．model = "pooling"でPooled OLS，model = "within"で固定効果モデル，model = "random"で変量効果モデルをそれぞれ推定することができます．
Pooled OLSの推定結果
まずはPooled OLSの結果を見てみます．
> summary(res_p)
Pooling Model
Call:
plm(formula = inv ~ value + capital, data = Grunfeld, model = "pooling")
Balanced Panel: n = 10, T = 20, N = 200
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-291.6757 -30.0137 5.3033 34.8293 369.4464
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) -42.7143694 9.5116760 -4.4907 1.207e-05 ***
value 0.1155622 0.0058357 19.8026 < 2.2e-16 ***
capital 0.2306785 0.0254758 9.0548 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Total Sum of Squares: 9359900
Residual Sum of Squares: 1755900
R-Squared: 0.81241
Adj. R-Squared: 0.8105
F-statistic: 426.576 on 2 and 197 DF, p-value: < 2.22e-16
企業価値valueと有形固定資産capitalの回帰係数はともに正で有意です．なおPooled OLSはRの標準関数を使ってlm(inv ~ value + capital, data = Grunfeld)としても同じ結果を得ることができます．
固定効果モデルの推定結果
次に固定効果モデルの結果を見ます．
> summary(res_w)
Oneway (individual) effect Within Model
Call:
plm(formula = inv ~ value + capital, data = Grunfeld, model = "within")
Balanced Panel: n = 10, T = 20, N = 200
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-184.00857 -17.64316 0.56337 19.19222 250.70974
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
value 0.110124 0.011857 9.2879 < 2.2e-16 ***
capital 0.310065 0.017355 17.8666 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Total Sum of Squares: 2244400
Residual Sum of Squares: 523480
R-Squared: 0.76676
Adj. R-Squared: 0.75311
F-statistic: 309.014 on 2 and 188 DF, p-value: < 2.22e-16
係数の推定値は若干異なっていますが，やはりともに正で有意です．固定効果の推定値を見るには，fixef()関数を用います．
> fixef(res_w)
1 2 3 4 5 6
-70.2967 101.9058 -235.5718 -27.8093 -114.6168 -23.1613
7 8 9 10
-66.5535 -57.5457 -87.2223 -6.5678
10の企業の固定効果の推定値をそれぞれ表示しています．
F検定の結果
次に，Pooled OLSで十分なのか，それとも個別効果を考慮する必要があるのか，pFtest()という関数を用いてF検定を行ってみます．
> pFtest(res_w, res_p)
F test for individual effects
data: inv ~ value + capital
F = 49.177, df1 = 9, df2 = 188, p-value < 2.2e-16
alternative hypothesis: significant effects
df1 = 9は分子の$N-1$を表し，df2 = 188は分母の$NT - N - K$（ここで$K=2$）を表しています．P値は限りなく0に近く，Pooled OLSは棄却されます．
変量効果モデルの推定結果
変量効果モデルの推定結果を見てみます．
> summary(res_r)
Oneway (individual) effect Random Effect Model
(Swamy-Arora's transformation)
Call:
plm(formula = inv ~ value + capital, data = Grunfeld, model = "random")
Balanced Panel: n = 10, T = 20, N = 200
Effects:
var std.dev share
idiosyncratic 2784.46 52.77 0.282
individual 7089.80 84.20 0.718
theta: 0.8612
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-177.6063 -19.7350 4.6851 19.5105 252.8743
Coefficients:
Estimate Std. Error z-value Pr(>|z|)
(Intercept) -57.834415 28.898935 -2.0013 0.04536 *
value 0.109781 0.010493 10.4627 < 2e-16 ***
capital 0.308113 0.017180 17.9339 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Total Sum of Squares: 2381400
Residual Sum of Squares: 548900
R-Squared: 0.7695
Adj. R-Squared: 0.76716
Chisq: 657.674 on 2 DF, p-value: < 2.22e-16
これも同じく回帰係数の推定値は正に有意です．
Hausman検定の結果
個別効果を考慮する際，変量効果で十分なのか，それとも固定効果モデルの推定を用いるべきなのか，phtest()関数を用いてHausman検定を実行してみます．
> phtest(res_w, res_r)
Hausman Test
data: inv ~ value + capital
chisq = 2.3304, df = 2, p-value = 0.3119
alternative hypothesis: one model is inconsistent
P値が0.3119と帰無仮説は棄却されず，固定効果モデルの推定を用いるべきとは言い切れません．
補足
plm()関数の引数modelには，model = "fd"というものも用意されていますが，これは1階の階差（first-difference）をとって固定効果モデルを推定するというものです．これを実行してみます．
> res_f <- plm(inv ~ value + capital, data = Grunfeld, model = "fd")
> summary(res_f)
Oneway (individual) effect First-Difference Model
Call:
plm(formula = inv ~ value + capital, data = Grunfeld, model = "fd")
Balanced Panel: n = 10, T = 20, N = 200
Observations used in estimation: 190
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-200.889558 -13.889063 0.016677 9.504223 195.634938
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) -1.8188902 3.5655931 -0.5101 0.6106
value 0.0897625 0.0083636 10.7325 < 2.2e-16 ***
capital 0.2917667 0.0537516 5.4281 1.752e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Total Sum of Squares: 584410
Residual Sum of Squares: 345460
R-Squared: 0.40888
Adj. R-Squared: 0.40256
F-statistic: 64.6736 on 2 and 187 DF, p-value: < 2.22e-16
回帰係数の推定値に，(Intercept)が含まれています．元の固定効果モデルに定数項が含まれていても，階差をとれば消えてしまうはずですが，これはどういうことでしょうか．実は，ここでは想定しているモデルが異なり，
$$
y_{it} = x_{it}^\top \beta + \delta \cdot t + b_i + \varepsilon_{it}, \quad (i=1,\dots,N; \ t = 1,\dots,T)
$$
を考えています．時刻$t$に関する線形トレンドの存在を仮定しており，その傾きのパラメータが$\delta$です．このモデルの1階の階差をとると，
$$
\Delta y_{it} = \Delta x_{it}^\top \beta + \delta + \varepsilon_{it}, \quad (i=1,\dots,N; \ t = 2,\dots,T)
$$
となり元のモデルのトレンド$\delta$が，階差モデルの定数項として出てくるわけです．このような時間トレンドを仮定しないのであれば，
plm(inv ~ value + capital - 1, data = Grunfeld, model = "fd")
としてplm()内のformulaで右辺に-1をつければ，階差モデルにおいて定数項なしで推定することができます．しかしもう1点注意しなければならない点があり，この推定値はなおmodel = "within"の推定値と一致しません．前回の記事で，「階差をとってGLS」「Within-Group推定」「最小二乗ダミー変数推定」は全て同じ結果になると説明したことと矛盾します．これはなぜかというと，どうやらplm()のmodel = "fd"推定は，「階差をとってGLS」ではなく「階差をとってOLS」推定しているようです．階差をとったデータセットに対して，lm()を当てはめると同じ結果になることから確認できます．
おわりに
株式会社Nospareには，統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospare までお問い合わせください．
