ロジット・プロビットの限界効果とRでの計算

バイナリー従属変数モデル
限界効果の定義
Rでの実装: ロジット（ロジスティック回帰）
Rでの実装: プロビット
- Marginsパッケージの利用
参考文献

回帰分析で、従属変数が二項(バイナリ)の離散値(0か1)を取る場合は、まず線形確率モデルが思いつく。
しかし、線形なので、予測値が0から1の間に収まらないという問題や、確率に与える影響が一定、不均一分散と言った問題もある。

一方で、ロジット（ロジスティック回帰）・プロビットは予測値が0-1に収まるが、非線形モデルであるために、推定されたパラメータそのものの解釈は難しい。
そこで、説明変数の1単位の変化が、従属変数の確率の変化にどの程度影響をおよぼすか、という限界効果を計算する必要がある。

バイナリー従属変数モデル

バイナリー従属変数は、以下のように、従属変数がyが0か1という二つの値を取るとして考える。この時、yが1である確率をpとする。

$y = \begin{cases} 1 & \mbox{with probability } p \\ 0 & \mbox{with probability } 1-p \end{cases}$

この確率 $p$ をパラメタライズして、モデル化するのがバイナリ変数に対するアプローチである。
具体的には、累積分布関数 $F(\cdot)$ を特定化して、それが説明変数Xとパラメータ $\beta$ の条件付き確率になるようなモデルである。

$p_i = Pr(y_i | {\bf x}) = F(x_i'\beta)$

この累積分布関数 $F(\cdot)$ が標準正規分布のものならプロビット、ロジスティック分布のものならロジットになる。

ロジット（ロジスティック回帰）
$F(x_i'\beta) = \frac{\exp(x_i' \beta)}{1+\exp(x_i' \beta)}$

プロビット
$F(x_i'\beta) = \int_{-\infty}^{x_i' \beta}\phi (z)dz$

正規分布の累積分布は明に解けないので、確率密度関数 $\phi(z)$ の積分になっている。

ちなみに、線形確率モデルならば累積分布関数はなく、 $p_i = x_i'\beta$ である。

限界効果の定義

限界確率効果 (Marginal Probability Effects, MPE)

限界効果は上述したとおり、説明変数の変化の条件付き確率に対する影響の測定である。
説明変数のベクトル $x_i$ のうちのj番目の変数 $x_{ij}$ の影響を見たいとすると、限界確率効果 (MPE) は以下の式で表される。

$\frac{\partial Pr(y_i | {\bf x})}{\partial x_{ij}} = \frac{\partial F(x_i'\beta)}{\partial x_i'\beta} \frac{\partial x_i'\beta}{\partial x_{ij}}=F'(x_i'\beta)\beta_{j}$ (1)

平均における限界確率効果 (MPE at mean, MEM)

限界確率効果はそれぞれの観察値 i について計算されるので、例えばN=1000のデータであれば1000のMPEが計算される。
それでは要約としての意味をなさず解釈できない。ではどうするか？
ひとつめは平均におけるMPEを計算する。MEMとかMPEMとか呼ばれるものである。

つまり、データの方を要約してしまい、それでMPEを計算するのである。

$\frac{\partial Pr(y | {\bf x})}{\partial \overline{x}_{j}} = \frac{\partial F(\overline{x}'\beta)}{\partial \overline{x}'\beta} \frac{\partial \overline{x}'\beta}{\partial \overline{x}_{j}}=F'(\overline{x}'\beta)\beta_{j}$

ここで $x_{i}$ というベクトルだったものが、 $\overline{x}$ に変わっている。 $\overline{x} = \{\overline{x}_1, \overline{x}_2, \ldots, \overline{x}_J \}$ というデータの平均のベクトルである。

目的に応じて、平均の代わりに最大値・最頻値・最小値、任意の値など、特定のデータの値における限界効果を計算することができる。

平均限界確率効果 (Average Marginal Probability Effects)

限界効果は、係数だけではなく変数にも依存する。そのため、平均効果を計算するのが一般的だ。
データ自体を平均するMEMとは対象的に、まず各観測値に対してMPEを計算し、その平均を取るのがAME (AMPE)である。

一般に「限界効果を計算する」と言った時はこちらが用いられることが多い。

$AMPE_{j} = \frac{1}{n}\sum^n_{i=1} F'(x_i'\beta)\beta_{j}$

AMPEの計算は他には以下の二つの方法がある。

説明変数の平均をとって、計算する。

上記のように各観測値のMPEを計算してから平均を取るのではなく、各変数の平均値を取ってから、MPEをひとつだけ計算する。

$AMPE_j = F'({\bar x_i}'\beta)\beta_j$

従属変数の平均の逆関数を線形部分の予測値として使う。

MPEの平均ではなく、従属変数の平均 ${\bar y}$ を使う。 $F(x_i\beta) = {\bar y}$ なので、 $x_i\beta = F^{-1}({\bar y})$ となる。これを使って計算すると

$AMPE_j = F'(F^{-1}({\bar y}))\beta_j$

となる。
ここでは一つ目の方法(各観測値のMPEの平均)で計算する。

ちなみにこの2つの方法はWooldridge本で詳しく説明されている。
本の中では2つ目の方法を平均における限界効果(Partial effect at the average, PEA)、1つ目の方法を平均限界効果(Average Partial Effect, APE)と呼び分けている。

Econometric Analysis of Cross Section and Panel Data, second edition (The MIT Press)

作者:Wooldridge, Jeffrey M.
発売日: 2010/10/01
メディア: ハードカバー

Rでの実装: ロジット（ロジスティック回帰）

今回は、サンプルデータを利用する。
sampleSelectionパッケージに入っているMroz87というデータを使う。
このデータは既婚女性の労働参加のデータで、今回は労働参加するかどうか(lfp)にどのような要素が影響するか推定する。
含まれる説明変数は、女性本人以外の家計への収入(nwifeinc), 教育年数(educ), 労働の経験年数(exper),女性の年齢(age), 5歳以下の子供の人数(kids5), 6歳から18歳の子供の人数(kids618)である。

# library
library(sampleSelection)

# Data load
data("Mroz87")

# Logit Estimation

logit1 = glm(lfp ~ nwifeinc + educ + exper + exper^2+age+kids5+kids618,
            family = binomial(link = "logit"), data = Mroz87)

MPEは(1)式で表されるが、ロジットなら

$\frac{\partial Pr(y_i | {\bf x})}{\partial x_{ij}} = \frac{\exp(x_i' \beta)}{1+\exp(x_i' \beta)}(1-\frac{\exp(x_i' \beta)}{1+\exp(x_i' \beta)})\beta_j$ (2)

となる。実際には、手で計算しなくても、各観測値の推定係数による予測値 $\hat{F}(x_i'\beta) = \frac{\exp(x_i' {\hat \beta})}{1+\exp(x_i' {\hat \beta})}$ はfittedで取得できる。

# get fitted probability: F(X'b)

P_logi = logit1$fitted

まず、推定された係数のベクトルを取得する。

coef = logit1$coefficients

そして、(2)式の通り計算する。sapplyを使って、係数ごとに処理を行っている。

MPE = sapply(coef, function (x) P_logi*(1-P_logi)*x)

このMPEは各観測値のMPEを計算したベクトルである。
平均限界確率効果(AMPE)を計算するにはシンプルに平均を取る。

$AMPE_j = \frac{1}{n} \sum_{i=1}^{n} F'(x_i'\beta)\beta_j$

# colMeansは列ごとの平均を計算
AMPE = colMeans(MPE)

結果は以下の通り。

 (Intercept)     nwifeinc         educ        exper          age        kids5      kids618 
 0.151837991 -0.003663443  0.041130581  0.021699206 -0.016506190 -0.260833293  0.010541658

つまり、例えばageならば、年齢が1歳上がると労働参加の確率が1.65%下がるという結果である。

Rでの実装: プロビット

プロビットでも計算可能である。
プロビットなら、累積分布関数は確率密度関数の積分だった。MPEは累積分布関数の微分と係数の積だが、積分の微分は元の関数なので、シンプルに確率密度関数とパラメータの積になる。

$\frac{\partial Pr(y_i | {\bf x})}{\partial x_{ij}} =\phi(x_i'\beta)\beta_j$

probit1 = glm(lfp ~ nwifeinc + educ + exper + exper^2+age+kids5+kids618,
             family = binomial(link = "probit"), data = Mroz87)

予測確率をfittedで取得した後、その確率に対応する分位数(Quantile)を累積分布関数の逆関数( $\Phi^{-1}(z)$ )で計算し、その分位数を用いて確率密度関数( $\phi(z)$ )から計算する。累積分布関数の逆関数はRのqnorm,確率密度関数はRのdnorm関数を使う。

# get fitted probability: F(X'b)
P_prob = probit1$fitted

# coefficients estimated
coef.p = probit1$coefficients

# MPE of Probit: 
MPE.p = dnorm(qnorm(P_prob,0,1),0,1)*coef.p["age"]

# Average MPE
AMPE.p = mean(MPE.p)

結果は

AMPE.p
[1] -0.01696685

年齢が1歳上がると労働参加の確率が約1.69%下がるという結果である。
ロジットより少しだけ大きく効果が出ている。

MPEやAMPEは、係数の推定値と、その予測値から計算するので、標準誤差が推定されない。
それはBootstrapなどで推定することが可能なようである。
時間があれば別の記事でまとめたいと思う。
https://www.r-bloggers.com/probitlogit-marginal-effects-in-r/

Marginsパッケージの利用

ここまで手で計算しておいてなんだが、marginsというmarginal effectを計算してくれるパッケージが存在する！
開発者のgithubページ

とても簡単で、プロビットやロジットの推定結果が入ったオブジェクトを関数にいれるだけである。

library("margins")

# Logit
logi_mar = margins(logit1)

#Probit
probi_mar = margins(probit1)

結果

> logi_mar
Average marginal effects
glm(formula = lfp ~ nwifeinc + educ + exper + exper^2 + age +     kids5 + kids618, family = binomial(link = "logit"), data = Mroz87)

  nwifeinc    educ  exper      age   kids5 kids618
 -0.003663 0.04113 0.0217 -0.01651 -0.2608 0.01054

> probi_mar
Average marginal effects
glm(formula = lfp ~ nwifeinc + educ + exper + exper^2 + age +     kids5 + kids618, family = binomial(link = "probit"), data = Mroz87)

  nwifeinc    educ   exper      age  kids5 kids618
 -0.003532 0.04083 0.02144 -0.01697 -0.267 0.01055