はじめに

千葉大学・株式会社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 までお問い合わせください．