基本の固定効果モデルは、以下のような因果関係を仮定したモデルである。
\[ Y_{it} = \alpha_i + \beta X_{it} + \epsilon_{it} \qquad (1) \]
ただし、\(Y_{it}\) は個人(企業や国家など分析の単位なら個人でなくてもよい) \(i \; (i =1, \; \ldots, \; N)\) の時点 \(t \;(t=1, \; \ldots, \; T)\) における従属変数の値、\(\alpha_i\) は個人 \(i\) の切片、\(\beta\) は傾き、\(X_{it}\) は個人 \(i\) の時点 \(t\) における独立変数の値、\(\epsilon_{it}\) は誤差項である。ここでは独立変数は一つだけだが複数の場合でも考え方は同じである。切片は個人によって違っているが、\(X_{it}\)の傾き(\(\beta\))は共通である。
\(Y_{it}\)に影響をおよぼし、\(X_{it}\)と相関する変数をモデルにきちんと投入していないと、ふつうは\(\beta\)の推定値が歪む。しかし、固定効果モデルの場合、もしもその投入し損なった変数が時間とともに変化しないならば、\(\alpha_i\) の推定値が歪むだけで\(\beta\)は正しく推定できる。このような正しい推定は横断データの分析では不可能であり、パネル・データを使う大きな利点である。
固定効果モデルの推定はロング形式のデータで普通は行う。推定のやり方にはいろいろあるが、通常の最小二乗法 (Ordinary Least Square: OLS) で推定するので、lm() のようなふつうの回帰分析のための関数を使ってもできる。普通の回帰分析と違うのは、個体ID を示す変数を因子としてモデルに投入する点だけで、それ以外は横断データの回帰分析と同じである。これは個体 IDのダミー変数 の係数が \(\alpha_i\) の推定値に対応するのである。
例えば以下のように。
library(plm)
data(EmplUK)
head(EmplUK)
## firm year sector emp wage capital output
## 1 1 1977 7 5.041 13.1516 0.5894 95.7072
## 2 1 1978 7 5.600 12.3018 0.6318 97.3569
## 3 1 1979 7 5.015 12.8395 0.6771 99.6083
## 4 1 1980 7 4.715 13.8039 0.6171 100.5501
## 5 1 1981 7 4.093 14.2897 0.5076 99.5581
## 6 1 1982 7 3.166 14.8681 0.4229 98.6151
summary(EmplUK)
## firm year sector emp wage capital
## Min. : 1.0 Min. :1976 Min. :1.000 Min. : 0.104 Min. : 8.017 Min. : 0.0119
## 1st Qu.: 37.0 1st Qu.:1978 1st Qu.:3.000 1st Qu.: 1.181 1st Qu.:20.636 1st Qu.: 0.2210
## Median : 74.0 Median :1980 Median :5.000 Median : 2.287 Median :24.006 Median : 0.5180
## Mean : 73.2 Mean :1980 Mean :5.123 Mean : 7.892 Mean :23.919 Mean : 2.5074
## 3rd Qu.:110.0 3rd Qu.:1981 3rd Qu.:8.000 3rd Qu.: 7.020 3rd Qu.:27.494 3rd Qu.: 1.5010
## Max. :140.0 Max. :1984 Max. :9.000 Max. :108.562 Max. :45.232 Max. :47.1079
## output
## Min. : 86.9
## 1st Qu.: 97.1
## Median :100.6
## Mean :103.8
## 3rd Qu.:110.6
## Max. :128.4
N <- length(unique(EmplUK$firm)) # 企業数
lm1 <- lm(wage ~ factor(firm) + emp + capital + output, data=EmplUK)
summary(lm1)$coefficients[-7:-N, ] # 企業IDダミーの係数がずっと続くので途中を省略
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.09568179 1.119100158 15.276275 5.948295e-47
## factor(firm)2 6.76703783 2.467689359 2.742257 6.224627e-03
## factor(firm)3 10.76478214 1.281417140 8.400685 1.746795e-16
## factor(firm)4 4.50577113 1.380993977 3.262702 1.145672e-03
## factor(firm)5 10.91826300 2.916280309 3.743900 1.928501e-04
## factor(firm)6 9.01016984 1.175140852 7.667311 4.604913e-14
## emp -0.10941297 0.034622999 -3.160124 1.630433e-03
## capital 0.14569492 0.069578937 2.093952 3.654733e-02
## output -0.03160371 0.007941394 -3.979617 7.465870e-05
このデータは 140 社の企業のパネル・データであり、firm という変数が企業の ID を示す。これを因子として投入すれば、企業別の切片を推定したのと同じである。しかし、これは個体の数が増えると計算に非常に時間がかかるため、あまり用いられていないようである。
固定効果モデルの推定法として最もよく使われているのが、固体内偏差を使った推定法である。これは用いるすべての変数に関して固体内平均偏差を計算し、それを使って回帰分析するのである。固体内偏差については、パネル・データの探索的分析の級内相関係数と固体内偏差のプロットの節を見よ。
固体内偏差を使うとどうして固定効果モデルの推定ができるのか? それは固体内偏差の回帰分析の係数(傾き)と固定効果モデルの係数(傾き)が、必ず一致するからである。まず、個人 \(i\) の固体内平均は以下の (2) 式のように書ける。
\[ \bar Y_{i} = \alpha_i + \beta \bar X_{i} + \bar \epsilon_i \qquad (2) \]
なぜなら、個人内平均 \(\bar Y_{i}\) とは、
\[ \bar Y_{i} =\frac{Y_{i1} + Y_{i2} + \; \ldots \; + Y_{iT}}{T} \]
であるが、(1)式に個人 \(i\) に関して \(T = 1, \; 2, \ldots, \; T\) 時点の Yの値を示すと、以下のようになるが、それをたてに足して、\(T\)で割ると、(2) 式が得られる。
\[ \begin{array}{ccccc} Y_{i1} =& \alpha_i + & \beta X_{i1} + & \epsilon_{i1} & \\ Y_{i2} =& \alpha_i + & \beta X_{i2} + & \epsilon_{i2} & \\ \vdots &\vdots & \vdots & \vdots & \\ Y_{iT} =& \alpha_i + & \beta X_{iT} + & \epsilon_{iT}& \\ たてに & 合計 & すると & & \\ \sum_{t=1}^T Y_{it} =& T \alpha_i + & \beta\sum_{t=1}^T X_{it} + & \sum_{t=1}^T\epsilon_{it} & (3) \\ \end{array} \] (3) 式を \(T\)でわると、 \[ \frac{\sum_{t=1}^T Y_{it}}{T} = \alpha_i + \frac{\beta\sum_{t=1}^T X_{it}}{T} + \frac{\sum_{t=1}^T\epsilon_{it}}{T} \] となる。これから、(2) 式が得られる。
(1)式から(2)式をひくと、
\[ Y_{it} - \bar Y_{i} = \alpha_i - \alpha_i + \beta X_{it} -\beta \bar X_{i} + \epsilon_{it} - \bar \epsilon_i \] \[ Y_{it} - \bar Y_{i} = \beta (X_{it} - \bar X_{i}) + (\epsilon_{it} - \bar \epsilon_i) \qquad (4) \] となる。この (4) 式は固体内偏差を使った回帰モデルそのものであるが、その傾きは固定効果モデルのそれと同じである。それゆえ、固体内偏差を使ったモデルを推定すれば、固定効果モデルの傾きを知ることができるのである。ただし、固体内偏差を使う場合、切片はゼロになるので、個々人の切片は推定されない。しかし、個々人の切片を知りたい場合はあまりないと思われるので、この固体内偏差モデルのことを「固定効果モデル」と呼ぶこともあるようである。
EmplUK データを使って固定効果(固体内偏差)モデルを推定してみよう。固定効果モデルの推定には、plm パッケージの plm() 関数を使う。基本の書式は以下の通り。
plm(モデル式, data=pdata.frameクラスのデータ, model=“within”)
lm() で通常の回帰分析をするときとほとんど同じだが、data= の引数は pdata.frame 形式にし、model=“within” の指定も必要。計量経済学では通常、固体内偏差で回帰分析したモデルを固定効果モデルと呼び、階差を使った回帰分析とは区別しているようである。
# 自作関数を読み込む
# source("https://dl.dropboxusercontent.com/u/1144829/R/plmFunctions.R", echo=T, encoding="UTF8")
data(EmplUK)
d1 <- pdata.frame(EmplUK) # いらない変数以外はpdata.frameに変換
xtabs(~ year, data=d1)
## year
## 1976 1977 1978 1979 1980 1981 1982 1983 1984
## 80 138 140 140 140 140 140 78 35
xtabs(~xtabs(~ firm, data=d1)) # 個体数が多いのでさらに度数分布表をとる
## xtabs(~firm, data = d1)
## 7 8 9
## 103 23 14
fixed1 <- plm(wage ~ emp + capital + output, data=d1, model="within") #固体内偏差でOLS
summary(fixed1)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = wage ~ emp + capital + output, data = d1, model = "within")
##
## Unbalanced Panel: n=140, T=7-9, N=1031
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -12.600 -1.120 -0.171 0.943 17.400
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## emp -0.1094130 0.0346230 -3.1601 0.00163 **
## capital 0.1456949 0.0695789 2.0940 0.03655 *
## output -0.0316037 0.0079414 -3.9796 7.466e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 4403.4
## Residual Sum of Squares: 4233.7
## R-Squared : 0.038538
## Adj. R-Squared : 0.033193
## F-statistic: 11.8646 on 3 and 888 DF, p-value: 1.2675e-07
先ほどのlm()を使って推定した結果と同じ結果が得られているのを確認してほしい。出力も lm() とほとんど同じなので、異なる点だけ解説しよう。Call のあとに “Unbalanced Panel: n=140, T=7-9, N=1031” という記載がある。これは個体数が140で、バランスは取れておらず、個々の個体は、7時点から9時点分のデータがあるということである。“Total Sum of Squares” (TSS) は、従属変数の固体内偏差の平方和である。
sum(Within(d1$wage)^2) # Total Sum of Squares と一致しているはず
## [1] 4403.433
“Residual Sum of Squares” (RSS) は、残差の平方和である。
sum(fixed1$residuals^2) # plmの結果$residuals で残差のベクトルが得られる。
## [1] 4233.732
ちなみに \[ R^2 =1 - \frac{\mathrm{RSS}}{\mathrm{TSS}} \] である。
通常の回帰分析と同じように対数変換したり、二乗項をモデルに投入することもできる。モデルを少しだけ修正する場合、update() 関数が便利。以下は書式。
update(元のモデル, 修正する引数)
以下は使用例。
fixed2 <- update(fixed1, log(wage)~.) # update() は前に推定したモデルをちょっとだけ修正する。これは fixed1 の従属変数だけかえる
summary(fixed2)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = log(wage) ~ emp + capital + output, data = d1,
## model = "within")
##
## Unbalanced Panel: n=140, T=7-9, N=1031
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -0.73200 -0.04760 -0.00716 0.04390 0.52600
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## emp -0.00422680 0.00138711 -3.0472 0.002378 **
## capital 0.00484135 0.00278756 1.7368 0.082774 .
## output -0.00149207 0.00031816 -4.6897 3.166e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 7.1188
## Residual Sum of Squares: 6.7954
## R-Squared : 0.045429
## Adj. R-Squared : 0.039128
## F-statistic: 14.087 on 3 and 888 DF, p-value: 5.6224e-09
fixed3 <- update(fixed1, ~.+I(output^2)) # これは fixed1 に独立変数を追加
summary(fixed3)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = wage ~ emp + capital + output + I(output^2), data = d1,
## model = "within")
##
## Unbalanced Panel: n=140, T=7-9, N=1031
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -12.400 -1.140 -0.174 0.977 17.200
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## emp -0.11494963 0.03443502 -3.3382 0.0008785 ***
## capital 0.13532421 0.06919182 1.9558 0.0508043 .
## output 0.52120275 0.15611045 3.3387 0.0008768 ***
## I(output^2) -0.00261875 0.00073858 -3.5457 0.0004121 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 4403.4
## Residual Sum of Squares: 4174.6
## R-Squared : 0.051975
## Adj. R-Squared : 0.044716
## F-statistic: 12.1573 on 4 and 887 DF, p-value: 1.261e-09
固定効果モデルでは、時間とともに変化しない変数(の主効果)は扱えない。これは時間とともに変化しないと、固体内偏差をとると必ずゼロになってしまい、変数にならないからである。実際、plm() のモデル式に時定変数を投入しても無視され、係数は推定されない(というか、できない)。
一時点前の値を独立変数に使うこともできる。その場合は、lag() 関数を使う。
lag(変数名)
をモデル式の中で指定すればよい。例えば以下のように。
fixed4 <- plm(wage ~ + lag(capital) + lag(output), data=d1)
summary(fixed4)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = wage ~ +lag(capital) + lag(output), data = d1)
##
## Unbalanced Panel: n=140, T=6-8, N=891
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -13.3000 -0.9120 -0.0295 0.8720 8.8700
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## lag(capital) 0.1890520 0.0578744 3.2666 0.001138 **
## lag(output) -0.0973598 0.0072374 -13.4523 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 3069.9
## Residual Sum of Squares: 2455.4
## R-Squared : 0.20016
## Adj. R-Squared : 0.16826
## F-statistic: 93.7191 on 2 and 749 DF, p-value: < 2.22e-16
回帰分析の検定や推定では自由度が問題になる。ふつう回帰分析で自由度といえば残差の自由度(これをモデルの自由度ということもある)のことであり、係数のt検定やモデル全体の比較のために用いられる F検定でもこの自由度が用いられる。plm() の出力の最後に出てくる “DF” の直前の数値がそれである。この自由度が大きいほど誤差が小さくなり、検定結果が有意になりやすい。一般に、
自由度 = サンプル・サイズ - 推定したパラメータの数
で計算できることが多い。。例えば fixed1 と名づけた最初の固定効果モデルの場合、サンプル・サイズは1031 である。まず固体内平均を計算するが、これもパラメータなので、140 のパラメータを推定したことになる。さらに3つの変数の傾きを推定したので、合計で 143 のパラメータを推定している。それゆえ、自由度は \(1031 - 143 = 888\) となる。
上の一時点前の独立変数の値を使ったモデルでは、これまでのモデルよりも顕著に2番目の自由度が下がっている。そのため誤差も大きくなる。これはサンプル・サイズが減少しているからである。一時点前の変数を使うため、最初の時点 (\(t=1\)) の \(Y\)は予測できない。それゆえ、最初時点の\(Y\)は分析に使えないのである。それゆえ、サンプル・サイズが減少しているのである。
交互作用項を投入することも可能である。交互作用効果はモデル式の中に
変数名 : 変数名
という項を投入すればよい。交互作用効果を投入する場合、主効果も投入するのが普通なので、その場合は、
変数名 * 変数名
とすれば2つの変数の主効果も交互作用効果も投入される。以下は使用例。
fixed5 <- plm(wage ~ capital * output, data=d1) # capital + output + capital:output でも同じ
summary(fixed5)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = wage ~ capital * output, data = d1)
##
## Unbalanced Panel: n=140, T=7-9, N=1031
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -12.600 -1.160 -0.155 0.952 17.400
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## capital 0.4558454 0.2584156 1.7640 0.0780755 .
## output -0.0327935 0.0084835 -3.8656 0.0001189 ***
## capital:output -0.0035459 0.0022648 -1.5656 0.1177883
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 4403.4
## Residual Sum of Squares: 4269.6
## R-Squared : 0.030402
## Adj. R-Squared : 0.026186
## F-statistic: 9.28127 on 3 and 888 DF, p-value: 4.7799e-06
交互作用項の一方の変数は時定変数(time constant variable) でも構わない。この例では、産業は時間とともに変化しないが、産業によって資本金やアウトプットの効果が違っていても不思議はないだろう。
fixed6 <- plm(wage ~ factor(sector) : (lag(capital) + lag(output)), data=d1)
summary(fixed6)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = wage ~ factor(sector):(lag(capital) + lag(output)),
## data = d1)
##
## Unbalanced Panel: n=140, T=6-8, N=891
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -13.400 -0.754 -0.045 0.740 6.790
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## factor(sector)1:lag(capital) 0.470577 0.321944 1.4617 0.1442601
## factor(sector)2:lag(capital) -2.710023 0.636632 -4.2568 2.343e-05 ***
## factor(sector)3:lag(capital) 0.288002 0.073593 3.9135 9.945e-05 ***
## factor(sector)4:lag(capital) -0.242308 0.635901 -0.3810 0.7032792
## factor(sector)5:lag(capital) 0.139468 0.184181 0.7572 0.4491518
## factor(sector)6:lag(capital) 0.218444 0.246288 0.8869 0.3753987
## factor(sector)7:lag(capital) 0.016917 0.102872 0.1644 0.8694285
## factor(sector)8:lag(capital) 0.321451 1.390081 0.2312 0.8171878
## factor(sector)9:lag(capital) -0.333511 1.351757 -0.2467 0.8051906
## factor(sector)1:lag(output) -0.171394 0.012692 -13.5045 < 2.2e-16 ***
## factor(sector)2:lag(output) 0.057283 0.023715 2.4155 0.0159573 *
## factor(sector)3:lag(output) -0.232759 0.043812 -5.3127 1.435e-07 ***
## factor(sector)4:lag(output) -0.046020 0.012864 -3.5774 0.0003698 ***
## factor(sector)5:lag(output) -0.010024 0.059697 -0.1679 0.8667016
## factor(sector)6:lag(output) -0.238135 0.057779 -4.1214 4.195e-05 ***
## factor(sector)7:lag(output) 0.335233 0.095778 3.5001 0.0004932 ***
## factor(sector)8:lag(output) -0.109715 0.015881 -6.9087 1.063e-11 ***
## factor(sector)9:lag(output) -0.086175 0.023980 -3.5936 0.0003480 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 3069.9
## Residual Sum of Squares: 2049.7
## R-Squared : 0.33232
## Adj. R-Squared : 0.27339
## F-statistic: 20.2683 on 18 and 733 DF, p-value: < 2.22e-16
時定変数の主効果は固定効果モデルでは推定できない。それゆえ主効果はモデルに投入していない。時変変数の主効果はモデルに投入したほうがよい場合が多いが、この場合は時定変数がカテゴリカルなので、投入しなかった場合、単純に産業別に capital や output の傾きを推定してくれる。
plm() がデフォルトで計算してくれるのは調整済み決定係数ぐらいであるうえに、複数のモデルの適合度を並べて表にしてくれるような関数も見当たらなかったので、自分で作った。
aic.plm(plmモデル1, plmモデル2, …)
結果を表にまとめるだけなら stargazer パッケージの stargazer() という関数が使える。
stargazer(plmモデル1, plmモデル2, … , type=“text”)
type=“text” を省略すると latex 形式で、type=“html” にすると html 形式で出力してくれる。
round(aic.plm(fixed1, fixed3, fixed4, fixed5, fixed6), 2)
## fixed1 fixed3 fixed4 fixed5 fixed6
## adjusted R squared 0.03 0.04 0.17 0.03 0.27
## Log Likelihood -2196.57 -2189.39 -1722.22 -2200.91 -1643.38
## n of parameters 3.00 4.00 2.00 3.00 18.00
## AIC 4399.14 4386.79 3448.43 4407.82 3322.76
## BIC 4413.95 4406.54 3458.02 4422.64 3409.02
library(stargazer)
stargazer(fixed1, fixed2, fixed3, type="text")
##
## ====================================================================================
## Dependent variable:
## -----------------------------------------------------------------------
## wage log(wage) wage
## (1) (2) (3)
## ------------------------------------------------------------------------------------
## emp -0.109*** -0.004*** -0.115***
## (0.035) (0.001) (0.034)
##
## capital 0.146** 0.005* 0.135*
## (0.070) (0.003) (0.069)
##
## output -0.032*** -0.001*** 0.521***
## (0.008) (0.0003) (0.156)
##
## I(output2) -0.003***
## (0.001)
##
## ------------------------------------------------------------------------------------
## Observations 1,031 1,031 1,031
## R2 0.039 0.045 0.052
## Adjusted R2 0.033 0.039 0.045
## F Statistic 11.865*** (df = 3; 888) 14.087*** (df = 3; 888) 12.157*** (df = 4; 887)
## ====================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
固定効果モデルの係数は、固体内偏差だけでなく、階差をとることでも推定できる。(1)式を再掲すると、以下のとおりであるが、 \[ Y_{it} = \alpha_i + \beta X_{it} + \epsilon_{it} \qquad (1) \] このモデルが正しければ、これよりも一つ前の時点 \(t-1\)に関しても \[ Y_{i,t-1} = \alpha_i + \beta X_{i, t-1} + \epsilon_{i, t-1} \qquad (5) \] という関係が成り立つ。(1)式から(5)式を引くと、 \[ Y_{it} - Y_{i,t-1} = \alpha_i - \alpha_i + \beta X_{it} - \beta X_{i, t-1} + \epsilon_{it} - \epsilon_{i, t-1} \]
\[ Y_{it} - Y_{i,t-1} = \beta (X_{it} - X_{i, t-1}) + (\epsilon_{it} - \epsilon_{i, t-1}) \] となり、やはり階差をとってOLSしても固定効果モデルの係数が推定できる。plm()では下記の書式で推定できる。
plm(モデル式, data=pdata.frameクラスのデータ, model=“fd”)
例えば以下のように。
fd1 <- update(fixed1, model="fd") # update でfixed1 の model= の部分だけ変更
stargazer(fixed1, fd1, type="text")
##
## ========================================================
## Dependent variable:
## -------------------------------------------
## wage
## (1) (2)
## --------------------------------------------------------
## emp -0.109*** -0.083**
## (0.035) (0.038)
##
## capital 0.146** 0.059
## (0.070) (0.075)
##
## output -0.032*** 0.011
## (0.008) (0.013)
##
## Constant 0.095
## (0.077)
##
## --------------------------------------------------------
## Observations 1,031 891
## R2 0.039 0.006
## Adjusted R2 0.033 0.006
## F Statistic 11.865*** (df = 3; 888) 1.707 (df = 3; 887)
## ========================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
もしもモデルが正しく特定されており、パネル・データがバランスしていれば、階差を使っても固体内偏差を使っても同じ結果が推定されるが、モデルが正しくなかったり、データがバランスしていない場合、推定結果が食い違うことがある。上の例の場合、切片が推定されているが、これはデータがバランスしていないことが原因であると思われる。
階差モデルは、従属変数の変化を独立変数の変化で予測するので、その意味が理解しやすいというメリットがある。しかし、上の例のようにデータがアンバランスな場合、切片が推定されてしまったり、自由度が固体内偏差モデルよりも小さくなってしまうことがある。また、因果関係にタイムラグがあるような場合も、それを適切に特定できていないと、固体内偏差モデルよりも歪みが大きくなることがある。固定効果モデルを推定したいだけならば、固体内偏差モデルを使っておいたほうが無難なように思える。
固定効果モデルでは、時間とともに変化しない観察されない異質性はコントロールできるが、すべての擬似相関を排除できるわけではない。擬似相関の典型例の一つが、時間によって生じる擬似相関である。例えば、20人の子供に関して、小学校1年生から6年生までのあいだ繰り返し身長と学力を測定したとしよう。これは\(N=20, \; T=6\)のバランスのとれたパネル・データである。
# 人工的に作った小学校 1~6年生の身長と学力テストの結果のパネルデータ d0 の分布をチェック
par(mfrow=c(2,2), mar=c(3,3,1, 0.5), mgp=c(2, 0.7, 0))
plot(score ~ time, data=d0)
m.lines(d0$time, d0$score, d0$id)
plot(height ~ time, data=d0)
m.lines(d0$time, d0$height, d0$id)
plot(score ~ height, data=d0)
m.lines(d0$height, d0$score, d0$id)
pd0 <- pdata.frame(d0)
plot(Within(pd0$height), Within(pd0$score))
m.lines(Within(pd0$height), Within(pd0$score), pd0$id)
学力も身長も学年が上になるにしたがって、右肩上がりで上昇しているのがわかる。それゆえ、身長と学力の散布図を作ると、固体内偏差をとっても取らなくても、身長が高いほど学力も高いという相関があるのがわかる。このデータに関して、通常の固定効果モデルを推定すると、以下のように身長が学力を高めるという結果が得られる。
fixed.s1 <- plm(score ~ height, data=pd0, model="within")
summary(fixed.s1)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = score ~ height, data = pd0, model = "within")
##
## Balanced Panel: n=20, T=6, N=120
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -13.900 -2.590 0.529 3.300 11.200
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## height 1.457933 0.045748 31.869 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 30278
## Residual Sum of Squares: 2689.3
## R-Squared : 0.91118
## Adj. R-Squared : 0.75172
## F-statistic: 1015.6 on 1 and 99 DF, p-value: < 2.22e-16
しかし、常識的に考えれば身長が学力を高めているというよりも、学年が上がるほど身長が伸び、学力も高まる、という擬似相関だと考えられるだろう。このような推定の誤りは、時間をコントロールすることで回避できる。この場合、身長と学力に関して時間と線形の関係があると考えられるので、学年を連続変数としてモデルに投入してやればよい。
fixed.s2 <- update(fixed.s1, ~.+as.numeric(time)) # as.numeric() はベクトルを数値に変換する。pdata.frame では時間の変数は自動的に因子に変換されるので数値に戻している
summary(fixed.s2)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = score ~ height + as.numeric(time), data = pd0,
## model = "within")
##
## Balanced Panel: n=20, T=6, N=120
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -14.100 -2.510 0.339 2.610 10.200
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## height 0.33936 0.26588 1.2764 0.2048
## as.numeric(time) 6.89927 1.61914 4.2611 4.679e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 30278
## Residual Sum of Squares: 2268.9
## R-Squared : 0.92506
## Adj. R-Squared : 0.75547
## F-statistic: 604.883 on 2 and 98 DF, p-value: < 2.22e-16
時間は常に線形の効果を持つとは限らないので、ダミー変数として投入する場合もある。その場合、
plm(モデル式, model=“within”, effect=“twoways”)
のように effect=“twoways” という引数を指定してやると、時間のダミーをモデルに投入してくれる。
fixed.s3 <- plm(score ~ height + time, data=pd0, model="within") # 上でも下でも実質的には同じ結果が得られる
fixed.s4 <- plm(score ~ height, data=pd0, model="within", effect="twoways")
library(stargazer)
stargazer(fixed.s2, fixed.s3, fixed.s4, type="text")
##
## ===================================================================================
## Dependent variable:
## ------------------------------------------------------------------
## score
## (1) (2) (3)
## -----------------------------------------------------------------------------------
## height 0.339 0.339 0.339
## (0.266) (0.262) (0.262)
##
## as.numeric(time) 6.899***
## (1.619)
##
## time2 5.759***
## (2.139)
##
## time3 12.872***
## (3.478)
##
## time4 17.356***
## (4.972)
##
## time5 27.709***
## (6.534)
##
## time6 34.226***
## (7.966)
##
## -----------------------------------------------------------------------------------
## Observations 120 120 120
## R2 0.925 0.931 0.017
## Adjusted R2 0.755 0.729 0.014
## F Statistic 604.883*** (df = 2; 98) 209.785*** (df = 6; 94) 1.673 (df = 1; 94)
## ===================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
この例の場合、上の3つのモデルの height の係数はほぼ同じである。最後の fixed.s4 だけ顕著に決定係数が低いのは、ベースラインモデルが異なるからである。ふつうの回帰分析では、決定係数や F 値の計算をするときに、切片のみのモデルと比べてどれだけ誤差が小さいかを計算する。しかし、plm() で effect=“twoways” を指定すると、時間ダミーだけを投入したモデルと比較して、どれだけ誤差が小さいかを計算する。それゆえ、このデータの場合、時間と学力の相関が強いので、 fixed.s3 では切片のみのモデルと比較して誤差が非常に小さく、決定係数は大きい。しかし、fixed.s4 では時間ダミーのみのモデルと比較するので、身長をモデルに投入してもほとんど誤差は小さくなっていないため、決定係数は小さく F値も有意にならない。どちらも間違ってはいないが、 effect=“twoways” の計算の仕方のほうが、実質的なモデルの当てはまりを検討しているという点で好ましいと私は思う。