【R】ロジット・プロビットでの限界効果とデルタメソッドによる標準誤差の計算

以前に限界効果の計算について書いた。

keita43a.hatenablog.com

前回の記事で触れていないのが、限界効果の標準誤差についてである。
仮に、ロジット（ロジスティック回帰）やプロビットの推定結果が統計的に有意であっても、
その結果によって計算した結果が同じとは限らない。
同じように有意性の検定を行うには標準誤差を計算する必要がある。

どうやって計算するのか？
- margins パッケージ
デルタ・メソッド
- ブートストラップ
デルタメソッドの実装

どうやって計算するのか？

margins パッケージ

上述の私の記事は2018年の5月8日に書いたのだが、その直後にmarginsの詳細な説明がアップされたらしく、以下の記事ではmarginsオブジェクトのsummary()によって標準誤差とp-valueまで計算してくれる。

library(sampleSelection)　　# データを使う
library(margins)                      # 限界効果の計算

# データのロード
data("Mroz87")

# ロジットの推定

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

# 限界効果の計算
logi_mar = margins(logit1)

# 限界効果の結果
summary(logi_mar)

結果は以下の通り。

> summary(logi_mar)
   factor     AME     SE       z      p   lower   upper
      age -0.0165 0.0023 -7.0591 0.0000 -0.0211 -0.0119
     educ  0.0411 0.0073  5.6025 0.0000  0.0267  0.0555
    exper  0.0217 0.0020 10.9056 0.0000  0.0178  0.0256
    kids5 -0.2608 0.0319 -8.1835 0.0000 -0.3233 -0.1984
  kids618  0.0105 0.0133  0.7941 0.4271 -0.0155  0.0366
 nwifeinc -0.0037 0.0015 -2.4801 0.0131 -0.0066 -0.0008

かんたんである。

しかし、多項ロジットを推定するmlogitパッケージ等、marginsが受け付けないオブジェクトについても計算したい時がある。
手で計算する方法は前の記事で解説したが、標準誤差（スタンダード・エラー）はどうやって計算するのだろうか？

デルタ・メソッド

そこで使うのがデルタメソッドという方法である。

たとえば二項従属変数（０か１か）のモデルの場合は、選択の確率 $p$ をモデル化するアプローチであった。（前の記事参照）
累積分布関数 $F(\cdot)$ を特定化して、それが説明変数 $X$ とパラメータ $\beta$ の条件付き確率になるようなモデルである。

$p=Pr(y|X)=F(X\beta)$

そして、限界効果は説明変数 $X$ のうち、ある一つの変数 $x_{k}$ が1単位増えたときの確率の変化であった。

$MPE = \frac{\partial P(y=1|X)}{\partial x_{k}} = \beta_{k} f(X\beta)$

ここで $f(\cdot)$ は、 $F(\cdot)$ の微分なので、確率密度関数になる。
このMPEを $\beta$ の関数として、 $MPE=h_{k}(\beta)$ とする。

説明変数の数が $K$ であるとすると、 $MPE=h_{k}(\beta)$ の勾配は

$\nabla h_{k}(\beta) = [ \beta_{k}x_{1}f'(X\beta),\beta_{k}x_{2}f'(X\beta),\ldots, \beta_{k}x_{k}f'(X\beta) + f(X\beta), \ldots, \beta_{k}x_{K}f'(X\beta) ]$

となる。ここで $f'(\cdot)$ は $f(\cdot)$ の微分である。
そして、デルタメソッドによって、限界効果の分散は

$[\nabla h_{k}(\beta)]\hat{V}[\nabla h_{k}(\beta)]'$

となる。 $\hat{V}$ は $\beta$ の分散共分散行列である。
標準誤差は、この値のルートを取れば得ることができる。

具体的な解説は、英語だがパパ・ウールドリッジと呼ばれるWooldridgeの以下の本の15章に書いてある。*1

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

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

なぜデルタメソッドが使えるのか、という証明などはここでは割愛する。*2

ブートストラップ

ここでは割愛するが、ブートストラップによって標準誤差を求めることも可能である。
とくにサンプルが小さい場合はこちらのほうがよいと思う。

デルタメソッドの実装

ここではデルタメソッドをRで実装してみる。
ざっくりいうと、AMPE.funという係数を引数にとる限界効果を計算する関数を作り、
その与えた引数の周りで微分して勾配を数値的に得る。
得た勾配を使って上述のデルタメソッドを手動で書いて計算するという順である。

# 係数を引数として、限界効果を計算する関数を作る。
AMPE.fun <- function(betas){
  temp = logit1  # 推定したglmのオブジェクトを受け渡す
  temp$coefficients <- betas   # 指定した係数を受け渡す
  P_logi = temp$fitted         # 確率を計算する。
  MPE = sapply(temp$coefficients, function(x) P_logi*(1-P_logi)*x) # MPEを計算する。手計算
  AMPE = colMeans(MPE)         # 各係数の平均限界効果を出す。
}

# 関数のチェック
tmp = AMPE.fun(logit1$coefficients)

# 勾配を数値計算
library(numDeriv)
grad <- jacobian(AMPE.fun, logit1$coefficients)

# 標準誤差の計算
ses = matrix(sqrt(diag(grad%*%vcov(logit1)%*%t(grad))), nrow = length(logit1$coefficients),byrow = TRUE) 

# 係数と標準誤差を合わせて表示
coefs = tidy(tmp) %>% as_tibble(cbind(nms = names(.), t(.))) %>%
  cbind(ses) %>%
  arrange(names)%>%
  mutate_at(c("x","ses"), list(~round(.,3))) # 小数点第三位で四捨五入

結果は以下の通り。
marginsによる計算と比べてみる。

> summary(logi_mar)
   factor     AME     SE       z      p   lower   upper
      age -0.0165 0.0023 -7.0591 0.0000 -0.0211 -0.0119
     educ  0.0411 0.0073  5.6025 0.0000  0.0267  0.0555
    exper  0.0217 0.0020 10.9056 0.0000  0.0178  0.0256
    kids5 -0.2608 0.0319 -8.1835 0.0000 -0.3233 -0.1984
  kids618  0.0105 0.0133  0.7941 0.4271 -0.0155  0.0366
 nwifeinc -0.0037 0.0015 -2.4801 0.0131 -0.0066 -0.0008
> coefs
        names       x    ses
1 (Intercept)  0.1518 0.1524
2         age -0.0165 0.0026
3        educ  0.0411 0.0078
4       exper  0.0217 0.0025
5       kids5 -0.2608 0.0365
6     kids618  0.0105 0.0133
7    nwifeinc -0.0037 0.0015