パネルデータ分析・超入門2

# はじめに
千葉大学・株式会社Nospareの[川久保](https://sites.google.com/site/ykawakubostat/)です．[前回の記事](https://qiita.com/ykawakubo/items/0dcfacaa00188b48737d)で，パネルデータ分析の基本的なモデルと推定法について紹介しました．今回は，各種の検定法と，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で十分です．そこで，

```{math}
\begin{align}
& H_0: b_1=b_2=\dots=b_N \\
& H_1: \mathrm{not} \ H_0
\end{align}
```

を検定したいですが，固定効果モデルに$N-1$個の制約を入れると，定数項を入れたPooled OLSになるため，以下の検定統計量を用いたF検定が利用できます．

```{math}
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$のときには，変量効果モデルの推定では，回帰係数を一致推定できなくなってしまうため，固定効果モデルの推定を行うのでした．そこで，

```{math}
\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検定と呼ばれる以下の検定統計量を用いた検定を利用します．

```{math}
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`という米国の企業パネルデータを読み込みます．

```R
library(plm)
data(Grunfeld)
```

`head(Grunfeld)`でデータの冒頭を見てみます．

```R
> 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つのモデルを当てはめます．

```R
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の結果を見てみます．

```R
> 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)`としても同じ結果を得ることができます．


### 固定効果モデルの推定結果
次に固定効果モデルの結果を見ます．

```R
> 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()`関数を用います．

```R
> 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検定を行ってみます．

```R
> 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は棄却されます．


### 変量効果モデルの推定結果
変量効果モデルの推定結果を見てみます．

```R
> 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検定を実行してみます．

```R
> 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）をとって固定効果モデルを推定するというものです．これを実行してみます．

```R
> 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$が，階差モデルの定数項として出てくるわけです．このような時間トレンドを仮定しないのであれば，

```R
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](https://bit.ly/38WtGZH) までお問い合わせください．
<a href="https://bit.ly/38WtGZH">
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1023506/d5c625c0-d6f4-680b-d702-3c40f11d6c39.jpeg" alt="株式会社Nospare" width=400>
</a>