盆暗の学習記録

データサイエンスを中心として,日々学んだことの備忘録としていく予定です。初心者であり独学なので内容には誤りが含まれる可能性が大いにあります。

[R]Rでパネルデータ分析:固定効果モデル

固定効果モデルの基礎的な理論と,Rでの実行方法の基本について学んだのでメモ。

理論

パネルデータ

  • パネルデータ(panel data):複数の主体を複数の期間にわたって観察したデータ
    • 同じ人を追跡調査したデータとか
    • 国×年度の公的統計データとか
  • \(k\)個の説明変数\(X\)と被説明変数\(Y\)の主体\(i\)時間\(t\)のパネルデータ:

\[ (X_{1it}, X_{2it},…,X_{kit}, Y_{it}), \hspace{1em} i = 1,…,n \text{ and } t = 1,…,T \]

  • balanced panel:各主体と各時間の観測値がすべて揃っているパネルデータ
  • unbalanced panel:欠損のあるパネルデータ

固定効果モデル

固定効果モデルはパネルデータの分析で使われる分析手法で,固有効果(個体差とか)のうち推定の邪魔になるもの(固定効果)を除去して分析を行うものです。

固有効果,固定効果,変量効果

  • 固有効果\(F_i\):時間によって変わらない,主体に固有の要素
    • 例えば\(Y_{it} = \beta_0 + \beta_1 X_{it} + u_{it}\)というモデルを作るとき,固有効果がある場合は誤差項\(u_{it}\)に含まれ,\(u_{it} = F_i + \varepsilon_{it}\)と分解できる。
    • 「固有効果」という語は一般的でないかも。山本(2015)しか使っていないかも。
  • 固定効果(fixed effect):説明変数\(X_{it}\)と独立でない(相関がある)固有効果\(F_i\)
    • 問題:固定効果があるとき,誤差項と説明変数が独立でなくなり(\(E(u_{it}|X_{it})\neq 0\)),そのままOLSすると欠落変数バイアスにより誤った推定をしてしまう
      • \(F_i\)の効果も\(X_{it}\)の効果も混ぜ合わせて推定してしまう(奥井2015)
      • 一致性(consistency, サンプルサイズを大きくすれば,推定値がある値に収束する性質)が得られなくなる(山本2015)
    • 対処:固定効果を除去するようなモデルにする(固定効果モデルにする) → 今回の本題
  • 変量効果(random effect):説明変数\(X_{it}\)と独立の(相関がない)固有効果\(F_i\)
    • 対処:\(\varepsilon_{it}\)のほうに\(X_{it}\)との相関がなければ,誤差項\(u_{it}\)と説明変数\(X_{it}\)は独立(無相関)になり,OLS推定が可能。ただし,パネルデータの時系列的要素の関係から誤差項の自己相関は生じるため,ロバスト標準誤差(頑健標準誤差)を使う必要がある
      • 奥井(2015)やAngrist and Pischke(2008)「OLS+ロバスト標準誤差でおk」(GLSはOLSより強い仮定を要するため)
      • 山本(2015)「一般化最小二乗法(GLS)でおk」

one-way固定効果モデル

主体の固定効果モデル(within model)の推定方法

\[ Y_{it} = \beta_1 X_{it} +\alpha_i + u_{it} \]

主体の固定効果(entity fixed effects)\(\alpha_i\)の除去については,以下のような推定方法があります。

1. 差分モデル(first difference model)のOLS推定(\(T=2\)のときのみ)

\[ (Y_{i,t+1} - Y_{it}) = (\beta_0 - \beta_0) + \beta_1 (X_{i,t+1}-X_{it}) + (F_i - F_i) + (\varepsilon_{i,t+1} - \varepsilon_{it}) \]

記号を置き換えると \[ \Delta Y_{it} = \beta_1 \Delta X_{it}+ \Delta \varepsilon_{it} \]

  • 推定方法:
    1. 説明変数,被説明変数それぞれ\(t+1\)期から\(t\)期を引く
    2. 上の式をOLS推定する
2. “\(n-1\)個のダミー説明変数”を用いたOLS推定
  • 最小二乗ダミー変数推定(Least Squares Dummy Variables (LSDV) 推定)とも呼ばれるもの

\[ Y_{it} = \beta_0 + \beta_1 X_{it} + \gamma_2 D2_i + \cdots + \gamma_n Dn_i + u_{it} \\ \text{where } D2_i = \begin{cases} 1 & \text{for } i = 2\\ 0 & \text{otherwise} \end{cases} \text{, etc.} \]

  • 推定方法:
    1. ダミー変数\(D2_i, \cdots, Dn_i\)を作成する
    2. 上の式をOLS推定する
    3. 検定はロバスト標準誤差で行う
      • 均一分散が仮定できることは稀だろうから不均一分散対策としてロバスト標準誤差を使う
3. ”主体の平均除去(Entity-demeaned)”を用いたOLS推定

\[ \begin{align} \tilde{Y}_{it} &= \beta_1 \tilde{X}_{it} + \tilde{u}_{it} \\ \text{where } \tilde{Y}_{it} &= Y_{it} - \bar{Y}_i, \hspace{1em} \bar{Y}_i = \frac{1}{T} \sum^T_{t=1} Y_{it}\\ \tilde{X}_{it} &= X_{it} - \bar{X}_i, \hspace{1em} \bar{X}_i = \frac{1}{T} \sum^T_{t=1} X_{it}\\ \tilde{u}_{it} &= u_{it}- \bar{u}_i, \hspace{1em} \bar{u}_i = \frac{1}{T} \sum_{t=1}^T u_{it} \end{align} \]

  • 推定方法:
    1. 説明変数・被説明変数について,変数から期間平均を引く
    2. 上の式をOLS推定する
    3. 検定はクラスタロバスト標準誤差を使う
  • \(n-1\)個のダミー変数による推定と同じ推定量が得られる
  • ダミー変数を使う方法だと計算量が多くなりがちでスケーラビリティがないため,統計ソフトでは通常entity-demeanedによる推定が行われる

時間の固定効果モデル(between)

\[ Y_{it} = \beta_1 X_{it} +\alpha_t + u_{it} \]

時間の固定効果(time fixed effects)\(\alpha_t\)を除きたい場合も主体の場合とほぼ同様です。

  1. \(T-1\)個のダミー説明変数”を用いたOLS推定
  2. ”time-demeaned”を用いたOLS推定

two-way固定効果モデル

two-way固定効果モデルあるいは主体と時間の固定効果モデル(entity and time fixed effects model) \[ Y_{it} =\beta_1 X_{it} +\alpha_i + \alpha_t + u_{it} \] によって主体の固定効果\(\alpha_i\)と時間の固定効果\(\alpha_t\)の両方を除去したい場合,それぞれの推定方法の組み合わせになります。

  1. 差分モデルのOLS推定(T=2)
  2. ”entity demeaning”と”\(T-1\)個のダミー変数”を用いたOLS推定
  3. ”time demeaning”と“\(n-1\)個のダミー変数”を用いたOLS推定
  4. ”entity & time demeaning”を用いたOLS推定
    • 説明変数と被説明変数について,主体と時間両方の平均を引いてOLS推定

固定効果モデルの仮定

説明変数\(X\)が1つ,主体の固定効果のみのモデルの場合について表記すると \[ Y_{it} = β_1 X_{it} + α_i + u_{it}, i = 1,…,n, t = 1,…, T \]

  • 仮定1: \(E(u_{it}|X_{i1},…,X_{iT},α_i) = 0\)
  • 仮定2: \((X_{i1},…,X_{iT},u_{i1},…,u_{iT}), i =1,…,n\) は独立かつ同一の同時分布に従う
    • 主体が母集団からランダム抽出されているなら成り立つ,が,サンプリングの時間が違うのに同一母集団のはずがない ⇒ 時系列方向で相関していると考えるのが普通
    • 多くのパネルデータにおいて\(u_{it}\)は自己相関(系列相関)している可能性が高い。 ⇒クラスタロバスト標準誤差(Clustered standard error)
  • 仮定3: 大きな外れ値は存在しない:\((X_{it},u_{it})\)はゼロでない有限の4次のモーメントを持つ
  • 仮定4: 完全な多重共線性は存在しない(Xが複数ある場合)
  • 仮定5:各主体についての誤差項は,説明変数の下で,自己相関なし:\(\text{cov}(u_{it},u_{is}|X_{i1},…,X_{iT}, \alpha_i) =0, \ t \neq s\)
    • この仮定が満たされない場合は頑健標準誤差を使う

クラスタロバスト標準誤差

  • クラスタロバスト標準誤差(clustered robust standard error, clustered standard error)
    • 不均一分散と自己相関を考慮した標準誤差(heteroscedasticity and autocorrelation-consistent standard error, HAC standard error)の一種
  • その誤差項は,ある同じクラスター(同じグループ)の中では相関があり,異なるクラスターの誤差項には無相関となる
  • 留意点:
    • クラスター数が少ない(20~30個にも満たない)場合は標準誤差を過小に推定してしまうバイアスがあるので調整するらしい(太田 2013
      • 「これ使えばOK」というような簡単に調整できる手法はまだ無いようだが,「クラスターが少なくてもクラスタロバスト標準誤差が相当程度に機能する」という研究もあるとか(Angrist and Pischke 同上)。
    • two-wayクラスタロバスト手法では推定された分散が負になる可能性がある

Rによる実践

library(plm)
data("Gasoline", package="plm")
head(Gasoline)
> head(Gasoline)
  country year lgaspcar  lincomep      lrpmg  lcarpcap
1 AUSTRIA 1960 4.173244 -6.474277 -0.3345476 -9.766840
2 AUSTRIA 1961 4.100989 -6.426006 -0.3513276 -9.608622
3 AUSTRIA 1962 4.073177 -6.407308 -0.3795177 -9.457257
4 AUSTRIA 1963 4.059509 -6.370679 -0.4142514 -9.343155
5 AUSTRIA 1964 4.037689 -6.322247 -0.4453354 -9.237739
6 AUSTRIA 1965 4.033983 -6.294668 -0.4970607 -9.123903

{plm}のGasolineデータセットを使います。

変数の意味は以下の通り:

  • country:18ヶ国の国名
  • year:年(1960~1978)
  • lgaspcar:自動車1台あたりのガソリン消費量の自然対数
  • lincomep:一人あたり実質所得の自然対数
  • lrpmg:実質ガソリン価格の自然対数
  • lcarpcap:一人あたりの自動車の量(ストック)の自然対数

ガソリン価格の変化が1台あたりのガソリン消費量に与えた影響(エコカー化に与えた影響)を考える,という状況を想定してみます。

lm()による推定

ダミー変数推定(LSDV推定)

factor型にして{makedummies}でダミー変数を作るか,変数をfactor型にしてそのままlmに入れます

library(tidyverse)
library(makedummies)
# ダミー変数を作成
Gasoline_dummies <- Gasoline %>% mutate(year = as.factor(year)) %>% select(country, year) %>% makedummies()
# 作成したダミー変数を結合
data <- Gasoline %>% select(-country, -year) %>% bind_cols(Gasoline_dummies)
# 推定
lm_LSDV <- lm(lgaspcar ~ ., data = data)
summary(lm_LSDV)

f:id:nigimitama:20181105015459p:plain

(コード文中の有意水準「***」がはてなmarkdownの太字指定とみなされてレイアウトが崩れるため,summaryの返り値を画像で貼り付けています)

 \beta_1 = -0.19285と出ました。

{lmtest}で HC_1ロバスト標準誤差を計算してみます。

library(lmtest)
coeftest(lm_LSDV, vcov = vcovHC(lm_LSDV, type = "HC1"))

f:id:nigimitama:20181105015525p:plain

平均除去(demean)による推定

# 主体固定効果の除去
entity_demeaned <- Gasoline %>% group_by(country) %>% 
  mutate(lgaspcar = lgaspcar - mean(lgaspcar),
         lincomep = lincomep - mean(lincomep),
         lrpmg = lrpmg - mean(lrpmg),
         lcarpcap = lcarpcap - mean(lcarpcap)) %>% ungroup()
# 時間固定効果の除去
twoways_demeaned <- entity_demeaned %>% group_by(year) %>% 
  mutate(lgaspcar = lgaspcar - mean(lgaspcar),
         lincomep = lincomep - mean(lincomep),
         lrpmg = lrpmg - mean(lrpmg),
         lcarpcap = lcarpcap - mean(lcarpcap)) %>% ungroup()
# 推定
lm_demeaned <- lm(lgaspcar ~ -1 + lrpmg + lincomep + lcarpcap, data = twoways_demeaned)

# クラスターロバスト標準誤差
library(plm)
coeftest(lm_demeaned, 
         vcov = plm::vcovHC(lm_demeaned, method = "arellano",
                            type = "HC1", cluster = c("group", "time")))

f:id:nigimitama:20181105015541p:plain

 \beta_1 = -0.19285となりました。理論通り,ダミー変数を使っての推定と同じ結果です。

plm()による推定

まずpdata.frame()でデータをplm用の形式にします。

重要な引数は以下の2つ:

  • index:個体,時間の列名を文字列ベクトルで指定する
  • drop.index:pdata.frame変換時にindex列を除くか否か
library(plm)
Gasoline_pdf <- pdata.frame(Gasoline, index=c("country","year"), drop.index=TRUE)
head(Gasoline_pdf)
> Gasoline_pdf <- pdata.frame(Gasoline, index=c("country","year"), drop.index=TRUE)
> head(Gasoline_pdf)
             lgaspcar  lincomep      lrpmg  lcarpcap
AUSTRIA-1960 4.173244 -6.474277 -0.3345476 -9.766840
AUSTRIA-1961 4.100989 -6.426006 -0.3513276 -9.608622
AUSTRIA-1962 4.073177 -6.407308 -0.3795177 -9.457257
AUSTRIA-1963 4.059509 -6.370679 -0.4142514 -9.343155
AUSTRIA-1964 4.037689 -6.322247 -0.4453354 -9.237739
AUSTRIA-1965 4.033983 -6.294668 -0.4970607 -9.123903

そしてplm()で推定します。

重要な引数は次の2つです:

  • method:モデルを指定します
    • fd:差分モデル(first-difference model)
    • within:主体固定効果モデル
    • between:時間固定効果モデル
  • effect:推定する固定効果を指定します
    • individual:主体固定効果
    • time:時間固定効果
    • twoways:主体・時間の両方
plm_result <- plm(lgaspcar ~ lrpmg + lincomep + lcarpcap, data = Gasoline_pdf,
                  method = "within", effect = "twoways")
summary(plm_result)

f:id:nigimitama:20181105015705p:plain

coeftest(plm_result,
         vcov = plm::vcovHC(plm_result, method = "arellano",
                            type = "HC1", cluster = c("group", "time")))

f:id:nigimitama:20181105015724p:plain

 \beta_1 = -0.192850と出ました。lmの推定量と一致しました。

しかし,平均除去によるOLS推定を行ったときの結果と標準誤差の値が違いますね。なんでだろう…ここの解明は今後の課題です。

まとめ

  • 固定効果モデルは主体(個体)や時間(年度とか)に固有の効果が推定の邪魔になるときに,その固有の効果(固定効果)を除去する手法です
  • ダミー変数を作って固有の効果を吸収させたりとかします
  • Rで実行するときは{plm}パッケージを使うと楽です。

参考文献