ロジスティック回帰の固定効果について
従属変数が0か1を取る離散変数だったり、３つ以上のカテゴリー変数だったりすると、まず思いつくのはロジスティック回帰 (Logistic Regression)、もしくはロジット離散選択モデル (Logit Discrete Choice Models)だと思う。
一方で、パネルデータを扱う場合、よく使われるのは個人の固定効果モデル(Individual fixed effect)だ。時間的に依らない個人の一貫した選択の傾向などがある場合は、その個人の傾向が考慮されないと誤差項に自己相関が含まれるし、個人の傾向が説明変数と相関していると推定された係数にバイアスが生じる。個人の傾向が説明変数と独立だと仮定すると変量効果(Random Effect)モデルでよいが、相関があるなら固定効果モデルで推定するとバイアスがない推定が出来る。（個人の効果を除いた残りの誤差が独立だと仮定して。）
これらのパネルデータ向けの推定は、Rだとplmとかlme4っていうパッケージがあって実際に推定できるんだけど、これらは線形モデルでの話。ロジットやプロビットと言った非線形のモデルで、そういうのはできるのか？と思って少し調べたことをまとめる。
まず、従属変数が二項の場合は、もともと入ってるglm()を使えばよい。
以下のコードでのデータはGreene(2004)*1のデータ生成過程を模したものであるが、そのコードは文末に掲載する。
glmfit = glm(y ~ x + d, data = dat, family = binomial("logit"))
このままだと、個人の効果が推定されない。直感的には個人のダミー変数を含める方法が思いつく。つまり
glmfit = glm(y ~ x + d + 0 + factor(id), data = dat, family = binomial("logit"))
と上のようにfactor(id)を含めるやり方である。
結果は以下の通りである。モデルの"真の値"はどちらも1である。Fixed effectを入れてないモデルだと、係数がおかしいようだ。統計的に違うかどうかはHausman testを検定するが、ここでは省略する。
| Dependent variable: | ||
| y | ||
| (1) | (2) | |
| x | 0.854*** | 1.078*** |
| (0.048) | (0.056) | |
| d | 0.883*** | 1.035*** |
| (0.081) | (0.091) | |
| Constant | -0.046 | |
| (0.052) | ||
| Observations | 4,000 | 4,000 |
| Log Likelihood | -2,241.378 | -1,923.132 |
| Akaike Inf. Crit. | 4,488.757 | 4,050.263 |
| Note: | *p<0.1; **p<0.05; ***p<0.01 | |
さてじゃあ固定効果含めたしこれでいいかというと、これがよくない。
Logitなど非線形なモデルに固定効果をそのまま含めて最尤法で推定すると、時間Tが個人の数Nに対して少ない場合、個人の効果の推定にバイアスが生じ、それが結果的に推定したい係数（ここではxとdの係数）にバイアスをもたらすという問題がある。これをIncidental Parameter Problemという。
じゃあどうするかというのに、答えたのがChamberlain (1980)*2の条件付きロジット(Conditional Logit）である。端的に言うと、固定効果が相殺されるような組み合わせの条件で条件付けた尤度関数を使ってパラメータを推定するというもの。これはRだとsurvivalパッケージのclogit()で実装されている。
library(survival) cl.fit = clogit(y ~ x + d + strata(id), data = dat1)
結果は以下のようになっている。推定値が少し1に近づいた。
| Dependent variable: | |||
| y | y | ||
| logistic | conditional | ||
| logistic | |||
| (1) | (2) | (3) | |
| x | 0.854*** | 1.078*** | 1.046*** |
| (0.048) | (0.056) | (0.055) | |
| d | 0.883*** | 1.035*** | 1.004*** |
| (0.081) | (0.091) | (0.090) | |
| Constant | -0.046 | ||
| (0.052) | |||
| Observations | 4,000 | 4,000 | 4,000 |
| R2 | 0.244 | ||
| Max. Possible R2 | 0.683 | ||
| Log Likelihood | -2,241.378 | -1,923.132 | -1,738.731 |
| Akaike Inf. Crit. | 4,488.757 | 4,050.263 | |
| Wald Test | 752.180*** (df = 2) | ||
| LR Test | 1,118.122*** (df = 2) | ||
| Score (Logrank) Test | 983.906*** (df = 2) | ||
| Note: | *p<0.1; **p<0.05; ***p<0.01 | ||
Conditional Logitでは、一致性を満たしたパラメータを推定できるが、固定効果自体は推定できない。
そこで、最近出たパッケージにbifeというものがある。
https://cran.r-project.org/web/packages/bife/vignettes/bife_introduction.html
このパッケージでは、Hahn and Newey (2004) *3の提案したバイアス補正 (Bias Correction)を実装している。
あと、新しいアルゴリズムが採用されていて速いらしい。
library(bife) bife(y ~ x + d | id, model = "logit", bias_corr = "ana")
結果は、clogitの結果と変わらない。
| Dependent variable: | ||||
| y | y | y | ||
| logistic | conditional | bife | ||
| logistic | ||||
| (1) | (2) | (3) | (4) | |
| x | 0.854*** | 1.078*** | 1.046*** | 1.047*** |
| (0.048) | (0.056) | (0.055) | (0.055) | |
| d | 0.883*** | 1.035*** | 1.004*** | 1.004*** |
| (0.081) | (0.091) | (0.090) | (0.091) | |
| Constant | -0.046 | |||
| (0.052) | ||||
| Observations | 4,000 | 4,000 | 4,000 | 4,000 |
| R2 | 0.244 | |||
| Max. Possible R2 | 0.683 | |||
| Log Likelihood | -2,241.378 | -1,923.132 | -1,738.731 | -1,923.132 |
| Akaike Inf. Crit. | 4,488.757 | 4,050.263 | 4,050.263 | |
| Wald Test | 752.180*** (df = 2) | |||
| LR Test | 1,118.122*** (df = 2) | |||
| Score (Logrank) Test | 983.906*** (df = 2) | |||
| Note: | *p<0.1; **p<0.05; ***p<0.01 | |||
ざっくりまとめると、二項ロジットをパネルデータで使う際に、ナイーブに個人ダミーを含めるだけで推定してはいけない。clogitを使って条件付き尤度関数で推定するか、バイアス補正を行わないといけない。
ちなみに、これは二項の場合だけど、多項ロジットでも、理論上は条件付きロジットは推定可能である。ただ、Rでそれを実装したパッケージは今のところないっぽい。Stataだとfemlogitというのがあるようだ。
付録：Greene (2004)のデータ生成過程 (Data Generating Process)
library(tidyverse) set.seed(3) # Number of observations, individual and time periods N = 100 TT = 40 # Data set for individual dat.ind = data.frame(id = 1:N) %>% mutate(a = rnorm(n(),0,1)) # Main Data Set (replicated DGP of Greene 2004) dat = expand.grid(id = 1:N,time = 1:TT) %>% # Independent variables mutate(x = rnorm(n(), 0, 1), h = rnorm(n(), 0, 1), u = runif(n(),0,1)) %>% mutate(e = log(u/(1-u)), d = ifelse(x + h > 0, 1, 0)) %>% left_join(dat.ind, by = "id") %>% group_by("id") %>% mutate(xbar = mean(x)) %>% ungroup() %>% mutate(alpha = sqrt(TT)*xbar + a) %>% mutate(w = alpha + x + d) %>% mutate(y = ifelse(w + e > 0, 1, 0)) # Only observed part dat1 = dat %>% select(y,x,d,id,time)
(4月9日追記)
Logitだと固定効果推定で係数を推定できるが、Probitだと推定できない(尤度関数の中で固定効果を相殺できない)ことが証明されているようだ。
以下のリンクの9ページ目に書いてあるところだが、元はChamberlain (2010)*4
http://www.jil.go.jp/institute/zassi/backnumber/2015/04/pdf/006-009.pdf
*1:Greene, William. 2004. “The Behaviour of the Maximum Likelihood Estimator of Limited Dependent Variable Models in the Presence of Fixed Effects.” Econometrics Journal 7 (1): 98–119.
*2:Chamberlain, Gary. 1980. “Analysis of Covariance with Qualitative Data.” The Review of Economic Studies 47 (1): 225–38.
*3:Hahn, Jinyong, and Whitney Newey. 2004. “Jackknife and Analytical Bias Reduction for Nonlinear Panel Models.” Econometrica 72 (4): 1295–1319.
*4:Chamberlain, G. (2010), Binary Response Models for Panel Data: Identification and Information. Econometrica, 78: 159-168.
