本記事について
固定効果モデルはパネルデータ分析の主要な手法です。しかし、二値変数やカウントデータに対しての適用法について日本語で書かれた解説が少なかったため、この記事を書きました。
- 一般化線形モデル（ロジスティック回帰・ポアソン回帰）の紹介
- それぞれに対する固定効果モデルの推定方法
- Pythonを用いた実装
を紹介します。少しでも参考になれば嬉しいです。
一般化線形モデル
目的変数が連続変数ではなく，二値変数だったりカウント変数だったりする場合に使用するいくつかの線形回帰モデル群です。今回は二値変数である場合のロジスティック回帰、カウントデータである場合のポアソン回帰について紹介します。
ロジスティック回帰
目的変数が二値(1か0)の場合に扱う手法です。目的変数が1を取る確率を $p$ とした場合に以下のようにモデリングします。
\log{\left(\frac{p}{1-p}\right)}= \bf{\beta' x}
ポアソン回帰
目的変数がカウントデータの場合に扱う手法です。
Y_i \sim Po(\lambda_i) (i=1,...,n)
このように目的変数がポアソン分布に従っていると仮定します。そして以下のようにモデリングします。
\log{\lambda_i}= \bf{\beta' x_i}
固定効果モデルの推定方法
ロジスティック回帰
通常の線形回帰では各グループごとのダミー変数を加えてOLSすることで、固定効果モデルを推定できます。しかし、ロジスティック回帰の場合、単純にダミー変数を含んで最尤推定すると、Incidental Parameter Problem(付随パラメーター問題)によってバイアスが生まれます。特に、時間Tが個人の数Nに対して少ない場合に顕著です。
したがって、条件付き最尤推定を行う必要があります。以下ではその実装を行います。
また、条件付き最尤法の推定式については、本記事の最後に紹介しています。
実装
Pythonを用いた実装を紹介します。
まずデータセットについてですが、以下のデータセット（5行のみ表示）を人工的に作成して使います。個体500×時間30の15,000がデータ数です。
| individual_id | time | X1 | X2 | X3 | Y | |
|---|---|---|---|---|---|---|
| 0 | 0 | 0 | -1.21781 | 1.550996 | -2.61253 | 1 |
| 1 | 0 | 1 | 4.084213 | 4.851235 | -3.935622 | 0 |
| 2 | 0 | 2 | 1.681075 | 3.314259 | -2.548931 | 0 |
| 3 | 0 | 3 | 2.23045 | 2.870518 | -3.957018 | 0 |
| 4 | 0 | 4 | 3.242873 | 1.777047 | -1.110802 | 0 |
statsmodelsのConditionalLogitを使います。インスタンスの作成時、引数にグループを示すデータを入れることで推定できます。
from statsmodels.discrete.conditional_models import ConditionalLogit
# 説明変数と目的変数を設定
X = data[['X1', 'X2', 'X3']]
Y = data['Y']
group = data['individual_id']
# 固定効果モデルの設定
model = ConditionalLogit(Y, X, groups=group)
# モデルのフィット
results = model.fit()
結果は以下のように、表示されます
Conditional Logit Model Regression Results
==============================================================================
Dep. Variable: Y No. Observations: 14970
Model: ConditionalLogit No. groups: 499
Log-Likelihood: -6796.7 Min group size: 30
Method: BFGS Max group size: 30
Date: Wed, 06 Dec 2023 Mean group size: 30.0
Time: 06:46:28
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
X1 0.5019 0.012 42.880 0.000 0.479 0.525
X2 -0.3015 0.014 -22.003 0.000 -0.328 -0.275
X3 0.2111 0.020 10.573 0.000 0.172 0.250
==============================================================================
単純にダミー変数を加え、最尤推定を行った場合と比較します。今回は人工データを使っているため、真の値がわかります。以下の表に真値・グループを示すダミー変数を加えたロジスティック回帰・条件付きロジスティック回帰の推定値と標準誤差を表示します。X1の推定値が、条件付きロジスティック回帰にすると、真値に近づいていることがわかります。
| 変数名 | 推定値の種類 | 条件付きロジスティック回帰 | ダミー変数を加えたのみ | 真値 |
|---|---|---|---|---|
| X1 | params | 0.501876 | 0.521594 | 0.5 |
| 〃 | se | 0.011704 | 0.011975 | - |
| X2 | params | -0.301523 | -0.313386 | -0.3 |
| 〃 | se | 0.013703 | 0.013985 | - |
| X3 | params | 0.211083 | 0.219301 | 0.2 |
| 〃 | se | 0.019964 | 0.020357 | - |
| const | params | - | - | -1.2 |
| 〃 | se | - | - | - |
ポアソン回帰
ポアソン回帰もロジスティック回帰と同様に条件付き最尤法による推定が可能です。
ただ、ロジスティック回帰と違いポアソン回帰においては、グループを示すダミー変数を加えた最尤推定値と同一の結果を示します。
以下実装です。statsmodelsのConditionalPoissonを使用します。
from statsmodels.discrete.conditional_models import ConditionalPoisson
# 説明変数と目的変数を設定
X = data[['X1', 'X2', 'X3']]
Y = data['Y']
group = data['individual_id']
# 固定効果モデルの設定
model = ConditionalPoisson(Y, X, groups=group)
# モデルのフィット
results = model.fit()
# 結果の表示
print(results.summary())
二つの手法を比べてみましょう。ロジスティック回帰の場合と同様に、表にまとめました。推定値が同じことが確認できます。
| 条件付きポアソン回帰 | ダミー変数を加えたのみ | 真値 | ||
|---|---|---|---|---|
| X1 | params | 0.50159 | 0.50159 | 0.5 |
| se | 0.00480 | 0.00480 | - | |
| X2 | params | -0.28997 | -0.28997 | -0.3 |
| se | 0.00741 | 0.00741 | - | |
| X3 | params | 0.19951 | 0.19951 | 0.2 |
| se | 0.01030 | 0.01030 | - | |
| const | params | - | - | -1.2 |
| se | - | - | - |
おわりに
本記事では二値変数・カウントデータに対する固定効果モデルの解説と実装を行いました。実際の実証分析では、推定値の解釈の問題から線形モデルを使用することも多く、一般化線形モデルは頑健性のチェックに使われるようです。
また、本記事では計量経済学の基礎知識や数学的な部分にはあまり触れていません。その点に興味がある方は以下の参考文献ををお読みください。また、条件付き最尤法の推定式については、本記事の最後に紹介しています。
誤っている個所がありましたら、ご指摘いただければ幸いです。読んでいただきありがとうございました。
参考文献
- パネルデータ分析・超入門1
- 計量経済学, 西山慶彦, 新谷元嗣, 川口大司, 奥井亮著, 有斐閣, 2019
- 固定効果モデル, Paul D. Allison著 ; 池田裕 [ほか] 訳（計量分析one point）共立出版, 2022
- stasmodels公式ドキュメント
【補足】 条件付き最尤推定について
最後に補足として、条件付き最尤推定の推定式を説明します。
まず、グループ i 内に2個体でそれぞれ従属変数が0,1を取った場合の条件付き確率は以下のようになります。
\begin{align}
P(Y_{i1}=1,Y_{i2}=0|X_{i1},X_{i2},Y_{i1}+Y_{i1}=1)
&=\frac{P(Y_{i1}=1|X_{i1})P(Y_{i2}=0|X_{i2})}{P(Y_{i1}=1|X_{i,1})P(Y_{i2}=0|X_{i,2})+P(Y_{i1}=0|X_{i,1})P(Y_{i2}=1|X_{i,2})}\\
&=\frac{\frac{exp(\beta' X_{i1}+\alpha_i)}{1+exp(\beta'X_{i1}+\alpha_i)}
\times
\frac{1}{1+exp(\beta' X_{i2}+\alpha_i)}}
{\frac{exp(\beta' X_{i1}+\alpha_i)}{1+exp(\beta' X_{i1}+\alpha_i)}
\times
\frac{1}{1+exp(\beta ' X_{i2}+\alpha_i)}
+ \frac{exp(\beta' X_{i2}+\alpha_i)}{1+exp(\beta' X_{i2}+\alpha_i)}
\times
\frac{1}{1+exp(\beta' X_{i1}+\alpha_i)}}\\
&=
\frac{exp(\beta' X_{i1})}{exp(\beta' X_{i1})+exp(\beta' X_{i2})}
\end{align}
これを個体数kに拡張します。全体でmこの個体がある中、k個の個体が1をとった場合を考えます。この時の条件付き確率は以下のようにあらわせます。ただし、$C_k^m$ はm個の個体からk個を選択したすべての組み合わせです。
P(Y_{ij}=1 (j\leq k),Y_{ij}=0(k< j\leq m)|X_{i1},X_{i2},\sum_{j=1}^m Y_{ij}=k)
=
\frac{exp(\sum_{j=1}^k \beta' X_{ij})}{\sum_{J\in C^m_k} exp(\sum _{j\in J}\beta' X_{ij})}
各グループの条件付き確率が算出できたため、すべてのグループで積を取ることで尤度が出ます。
これを最大化するようにパラメーターを推定します。