機械学習から統計に足を踏み入れた身としては、推定や予測の話になると、すぐさま機械学習の枠組みにはめてしまいがちです。何でも機械学習状態です。しかしながら、データによっては注意が必要な場合があります。特にパネルデータでは、複数の観測個体の複数の時点が集まったデータであるため、被説明変数には個体と時間の両方の影響があります。
パネルデータにおいて、被説名変数に対する説明変数の影響度合いを見ることはよくありますが、このとき個体特有の要素が加味されていないと、欠落変数の影響で適切に推定されません。そこで、パネルデータにおいては、個体特有の影響を加味して推定する必要があり、これを固定効果モデルと呼びます。
今回は、固定効果モデルの数式を使いながら理論的背景を見ていくとともに、Pythonによる実践も行います。今回のコードはこちらに置いてあります。
固定効果モデルを理解する
パネルデータとは、複数個体、複数時点のデータが集まっているデータのことです。そのため、通常通りの線形回帰モデルで推定を行うと、問題がある場合があります。それでは、何が問題なのかを見ていきます。
線形回帰での最小二乗法における推定の問題
推定の話は単純化して考えたほうが理解しやすいので、単回帰モデルで考えます。単回帰モデルは以下の式で表され、得られたデータ \(y, x\) を使ってパラメータ \(a, b\) を推定します。
パラメータ \(a, b\) を推定するためには、最小二乗法を使います。最小二乗法のやり方に関しては、他で詳細に解説されているので、ここでは割愛します。今回大事なことは、最小二乗法によって適切な推定量が必ず求まるということではないという点です。求めた推定量が適切かを条件に照らし合わせて確認できます。そして、この条件を満たした推定量のことをBLUE: Best Linear Unbiased Estimatorと呼びます。日本語で読むと、最良線形不偏推定量です。
BLUEとなるためには以下の４つの特性を保たなければいけません。
- 線形性
- 不編成
推定量の期待値が真の値に等しい - 効率性
推定量の分散が最小 - 一致性
サンプルサイズを大きくすれば、推定量がある値に収束する
上の４つの特性を保つためには、誤差項においてガウスマルコフの定理を満たす必要があります。ガウスマルコフの定理は、「均一分散」「共分散ゼロ」「説明変数と独立」の３つです。３つありますが今回の話で大事なのは、「誤差項と説明変数が独立」です。誤差項と説明変数が独立ではない、つまり相関すると、BLUEの特性である一致性が満たされなくなります。
誤差項 \(\epsilon\) は被説名変数 \(y\) の決定要因であるため、両者は相関します。そうすると、被説名変数 \(y\) に影響された誤差項 \(\epsilon\) によって、誤差項 \(\epsilon\) と相関している場合説名変数 \(x\) も影響を受けてしまいます。これは、被説名変数によって説明変数が決まるという、因果が逆になってしまっている状態です（同時決定バイアス）。被説名変数と誤差項が説明変数に影響を与えるため、説名変数のパラメータの一致性がなくなります。
パネルデータに関する推定の問題
説名変数と誤差項が相関していると、最小二乗法によって適切な推定量は求まらないことがわかりました。これをパネルデータを使って推定する場合で考えてみます。パネルデータは複数個体の複数時点を集めたデータですが、時間一定の個体に関する変数は欠落してしまっていることが多いです。以下の式の \(\mu_i\) の部分です。
この変数が欠落したまま、 \( y_{it} = \beta x_{it} + \epsilon_{it} \) として推定すると、欠落変数の影響で、誤差項 \(\epsilon_{it}\) と 説名変数 \(x_{it}\) が相関してしまい、適切な推定量が求まりません。つまり、パネルデータにおける推定の問題とは、個体を表す変数の欠落です。これを解決するモデルは、時間依存しない \(\mu_i \) いわば固定された要素を考慮したモデルということで固定効果モデルと呼ばれています。
固定効果モデルの導出
上記によって、パネルデータにおけるモデルのパラメータの推定での問題は、誤差項と説明変数が相関してしまうことだとわかりました。適切な推定量を求めるようにするため、この問題に対処します。考え方としては、「相関する部分を消してしまえばええやん」です。個体 \(i\) の要素の \( \mu_i \) さえ消してしまえば、誤差項と説明変数が相関しなくなるのでOKということです。では、どうやって \( \mu_i \) を消すかというと、両辺からそれぞれ時間平均したものを引きます。
\begin{align*}
y - \bar y &= \beta x_{it} + \mu_i + \epsilon_{it} - (\beta \bar x_i + \mu_i + \bar \epsilon_i) \\
&= \beta (x_{it} - \bar x_i) + \epsilon_{it} - \bar \epsilon_i
\end{align*}
ここで、\(\bar y = \sum_{t=1}^T \frac {y_{it}} {T} \)、 \(\bar x = \sum_{t=1}^T \frac {x_{it}} {T} \)、 \(\bar \epsilon = \sum_{t=1}^T \frac {\epsilon_{it}} {T} \) としています。上式を見てもらえば分かる通り、両辺から時間平均を引くことで、個体特有の要素を消すことができます。これは、個体特有の要素は時間依存していないからです。
以上により、固定効果モデルのやり方が分かりました。それでは、実際にPythonを使って、固定効果モデルを実践していきます。
固定効果モデルの実践
上記までで、パネルデータにおける適切な推定の仕方がわかりました。やり方は何のことはなく、変数の時間平均を引いた上で、線形回帰モデルで推定するというまさしくいシンプルイズザベストですね。早速この方法をPythonで試してみましょう。
準備
まず必要なモジュールをインポートします。statsmodelsで線形回帰モデルを構築し、使用データはstatsmodelsのデータセットに含まれている「grunfeld」です。grunfeldはアメリカの有名企業の投資額に関するデータであり、企業価値や有形固定資産も入っています。このデータは、それぞれの企業（10社）に対して20期分あり、パネルデータになっています。
from pathlib import Path from statsmodels.formula.api import ols from statsmodels.datasets import grunfeld import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt result_dir_path = Path('result') if not result_dir_path.exists(): result_dir_path.mkdir(parents=True)
データの取得
grunfeldデータを取得し、中身を見てみます。
grunfeld_data = grunfeld.load_pandas().data grunfeld_data.head()
上図を見れば分かる通り、企業名(firm)、年(year)、投資額(invest)、企業価値(value)、有形固定資産(capital)のデータになっていることがわかります。それでは、それぞれの企業の投資額を年を横軸、投資額を縦軸として見てみます。
plt.figure(figsize=(8, 5)) sns.lineplot( data=grunfeld_data, x='year', y='invest', hue='firm' )
上図を見ると、General MotorsとUS Steelは他の企業に比べて投資額高いまま推移しています。それぞれの企業特有の要素があることがわかりますが、その要素は今回のデータには含まれていません。そのため、このまま線形回帰モデルを行うと、説明変数と誤差項が相関してしまうと思われます。
通常の線形回帰モデルによる推定
上記グラフでそれぞれの企業の投資額を見たことで、時間依存しないそれぞれ企業の特有の要素が投資額に関係していることがわかりました。その効果を取り除かないといけないのですが、先に通常の線形回帰モデルを行うとどうなるのかを見てみます。
model_result = ols(
data=grunfeld_data,
formula='invest ~ +1 + capital + value '
).fit()
model_result.summary()
上図の結果を見ると、企業価値より有形固定資産の方が投資額には影響があり、有形固定資産の限界効果は0.23であることがわかります。それぞれの変数の結果は有意なので、このままの結果を採用してしまいそうですね。
固定効果モデルの実践
まずは、個体に対する時間平均を引いたデータを作ります。
- 企業名、年以外の変数を一つのカラムとして扱うために縦長のデータに変換
- それぞれの変数でgroup byを行って平均を求める
- 2. の結果を企業名、変数ごとに結合し、もとの値から平均を引く
- 元のデータの形式に戻す
target_data = grunfeld_data.loc[
:, ['firm', 'year', 'invest', 'capital', 'value']
]
target_data = target_data.set_index(['firm', 'year'])
target_data = target_data.stack()
target_data = target_data.reset_index()
target_data.columns = ['firm', 'year', 'colname', 'value']
target_data_mean = target_data.groupby(['firm', 'colname'])['value'].mean()
target_data_mean = target_data_mean.reset_index()
target_data_mean.rename(columns={'value':'value_mean'}, inplace=True)
target_data = pd.merge(target_data, target_data_mean, on=['firm', 'colname'], how='left')
target_data['value'] = target_data['value'] - target_data['value_mean']
target_data = pd.pivot_table(
data=target_data.loc[:, ['firm', 'year', 'colname', 'value']],
index=['firm', 'year'],
columns='colname'
)
target_data.columns = [col[1] for col in target_data.columns]
target_data = target_data.reset_index()
それでは、上記のデータを使って、固定化効果モデルを構築します。
model_result = ols(
data=target_data,
formula='invest ~ -1 + capital + value '
).fit()
model_result.summary()
上記の結果を見ると、有形固定資産の限界効果は0.31となっており、通常の線形回帰モデルで求めたときの限界効果は0.23でした。つまり、通常の線形回帰モデルでの推定では、0.08ほど過小に見積もられていたことがわかりました。
おまけ Rを使って固定効果モデルを実践
Rには固定効果モデル用のパッケージがあります。それを使って、Pythonで行った結果と同じになるか確かめてみます。
library(plm)
library(dplyr)
grunfeld_data <- read.csv('grunfeld_data.csv')
grunfeld_data <- grunfeld_data %>%
dplyr::mutate(firm, as.numeric(as.factor(firm))) %>%
dplyr::select(firm, year, invest, capital, value)
result<- plm(
invest~value+capital,
data=grunfeld_data,
model="within"
)
summary(result)
Pythonでの結果と同じになりました。
終わりに
今回はPythonを使って固定効果モデルを実践してみました。なぜこのテーマにしたかというと、ここ最近因果推定に関連して計量経済学を勉強しており、線形回帰周りを今まで雑に理解してきたと感じたからです。計量経済学や因果推定の理解が進んだのは以下の本のおかげです。具体例が多くわかりやすかったです。ありがとうございます。
