7.3 プロビットモデルの推定と解釈
7.3.1 分析の実行と結果のまとめ
プロビットモデルを構築できたら、次はモデル内のパラメータに関する推定値（\(\hat{\beta}_0, \hat{\beta}_1,...,\hat{\beta}_k\)）を求め、選択確率に関する含意を得る。プロビットモデルのパラメータ推定においては、最小二乗法ではなく、最尤（Maximum Likelihood: ML）法を用いるが、詳細については補足（7.4）を参照してほしい。
プロビットモデルを推定する場合のデータセットには、説明変数に関する列（\(x_{1i},.., x_{ki}\)）と、それに対応する個人 \(i\) の選択結果 \(y_i\)（選択していれば1、選択していなければ 0）が記録されている。choice_data.xlsx は、二つの製品に対する消費者の選択結果を捉えた人工データセットである。データには製品1の価格 \(p1\)、製品2の価格 \(p2\)（千円）、製品1のクーポン広告を受け取ったかのダミー変数 \(a1\)、製品2のクーポン広告ダミー変数 \(a2\) と、製品1の選択結果 \(y1\)、製品2の選択結果 \(y2\) が含まれている。
## # A tibble: 6 × 6
## y1 y2 p1 p2 a1 a2
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 0 10.7 16.8 0 1
## 2 0 1 11.6 9.15 1 0
## 3 1 0 8.97 9.31 1 0
## 4 1 0 4.62 8.03 0 0
## 5 1 0 7.81 19.2 1 1
## 6 1 0 7.49 8.84 0 0このchoice_data.xlsxを用いて製品1の選択に関するプロビットモデルを推定してみる。Rではglm()関数において、family = binomial(link = probit) という引数を指定することで実行が可能である14。以下では分析の実行と結果の出力を行う。
##
## Call:
## glm(formula = y1 ~ p1 + a1, family = binomial(link = probit),
## data = choice_df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.04117 0.22350 9.133 < 2e-16 ***
## p1 -0.24720 0.02213 -11.171 < 2e-16 ***
## a1 0.62763 0.08656 7.251 4.13e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1384.5 on 999 degrees of freedom
## Residual deviance: 1205.5 on 997 degrees of freedom
## AIC: 1211.5
##
## Number of Fisher Scoring iterations: 3ここまで、我々は線形確率、プロビット、ロジットという3種類のモデルでの推定方法を概観してきた。以下では、これらの異なる推定方法を用いた分析結果を併記する。
library(modelsummary)
ols1 <- lm(y1 ~ p1 + a1,
data = choice_df)
logit1 <- glm(y1 ~ p1 + a1,
family = binomial(link = logit),
data = choice_df)
models <- list()
models[["Linear probability model"]] <- ols1
models[["Probit model"]] <- probit1
models[["Logit model"]] <- logit1
modelsummary(models,
title = "モデル比較",
notes = "Values in [ ] show robust standard errors",
stars = TRUE,
statistic = "std.error",
vcov = "robust",
gof_map = "nobs") #適合度指標において、サンプルサイズ（"nobs"）のみ表示するという指示| Linear probability model | Probit model | Logit model | |
|---|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |||
| Values in [ ] show robust standard errors | |||
| (Intercept) | 1.196*** | 2.041*** | 3.332*** |
| (0.065) | (0.212) | (0.361) | |
| p1 | -0.085*** | -0.247*** | -0.404*** |
| (0.006) | (0.021) | (0.036) | |
| a1 | 0.221*** | 0.628*** | 1.035*** |
| (0.029) | (0.087) | (0.145) | |
| Num.Obs. | 1000 | 1000 | 1000 |
分析の結果、推定値の違いはあるものの、どのアプローチでも p1 は被説明変数に対して負に、a1 は正に有意な影響を与えることが示されている。しかしながら、詳しくは次項にて説明するが、係数の解釈については注意が必要である。前章で紹介した modelsummary パッケージを使えば、複数の分析アプローチを用いた結果をまとめ、併記することも簡単である。結果の頑健性チェックや、様々な比較検討のために、このような表が用いられることも多い。
7.3.2 限界効果の計算
プロビットモデルの分析には成功したが、非線形モデルで推定された係数の解釈には注意が必要である。特に、プロビットモデルによって推定されたある変数の係数値は、「他の変数の影響をコントロールした場合にその変数が選択確率に与える影響」を意味しない。そのため、OLSの分析結果のように、係数の推定値のみを見て変数の影響の程度を議論することはできず、そのような解釈を行うためには特定の変数が持つ「限界効果」を追加的に分析する必要がある。限界効果の必要性や、平均的な限界効果、個別限界効果の平均値については次節で説明を加えている。
限界効果の計算には、mfx::probitmfx() を用いる15ため、以下のようにパッケージをインストールしてほしい。
mfx パッケージを用いた限界効果の計算では、線形モデルやプロビットモデルの推定同様、以下のようにモデルを指定し、分析を行う。平均的な限界効果の計算方法については、atmean という引数で指定する。atmean = FALSE とすることで、個別限界効果の平均値を計算できる。個別限界効果の平均値とは、データ内の個人の個別限界効果を計算し、その平均値求める方法である。一方で、平均値における限界効果はatmean = TRUE と指示することで計算できる。平均値における限界効果は、モデル内の説明変数の値について、それぞれ平均値を計算し、そのうえで限界効果を求める方法である。
library(mfx)
#Average Marginal Effects (限界効果の平均)
probit1_ame <- probitmfx(y1 ~ p1 + a1,
data = choice_df, atmean = FALSE)
probit1_ame## Call:
## probitmfx(formula = y1 ~ p1 + a1, data = choice_df, atmean = FALSE)
##
## Marginal Effects:
## dF/dx Std. Err. z P>|z|
## p1 -0.0847786 0.0060787 -13.9469 < 2.2e-16 ***
## a1 0.2188866 0.0289816 7.5526 4.266e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## dF/dx is for discrete change for the following variables:
##
## [1] "a1"#Marginal Effects at Mean（平均値における限界効果）
probit1_mem <- probitmfx(y1 ~ p1 + a1,
data = choice_df, atmean = TRUE)
probit1_mem## Call:
## probitmfx(formula = y1 ~ p1 + a1, data = choice_df, atmean = TRUE)
##
## Marginal Effects:
## dF/dx Std. Err. z P>|z|
## p1 -0.0984235 0.0088074 -11.1750 < 2.2e-16 ***
## a1 0.2446609 0.0324366 7.5427 4.602e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## dF/dx is for discrete change for the following variables:
##
## [1] "a1"分析結果を比較すると、多少計算結果は異なるものの、どちらの計算方法においても、価格は製品の選択確率に負に、クーポンは正に有意な影響を与えることが伺える。また影響の程度として、限界効果の平均値に着目すると、価格が1単位（千円）上昇すると、選択確率が約\(8.5\%\)下がること、クーポンを受け取った人は受け取っていない人と比べて約\(21.9\%\)選択確率が高いことが伺える。なお、詳細については省略するが限界効果の標準誤差はデルタ法と呼ばれる方法で計算されている（Fernihough, 2019）。
また、mfx による結果も以下のように整理して出力する事ができる。下表の(1) は限界効果の平均（Average Mean Effect）、(2) は平均的個人における限界効果（Marginal Effect at Mean）をそれぞれあらわしている。
marginal <- c(
"p1 marginal" = "p1",
"a1 marginal" = "a1"
)
marg_model <- list(
"(1) Probit_AME" <- probit1_ame,
"(2) Probit_MEM" <- probit1_mem
)
marg_sum <- modelsummary(
# models to summarize side-by-side
marg_model,
# S.E. in parentheses
statistic = "std.error",
# rename and select the coefficients
coef_map = marginal,
# significance stars
stars = TRUE,
# term and component columns are combined
shape = term:component ~ model,
add_rows = data.frame("Marginal effect type", "Average Marginal Effect", "Marginal Effect at Mean"),
title = "限界効果サマリー",
# omit all goodness-of-fit statisitcs except # of observations
gof_map = "nobs")
marg_sum| (1) | (2) | |
|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
| p1 | -0.085*** | -0.098*** |
| (0.006) | (0.009) | |
| a1 | 0.219*** | 0.245*** |
| (0.029) | (0.032) | |
| Num.Obs. | 1000 | 1000 |
| Marginal effect type | Average Marginal Effect | Marginal Effect at Mean |
最後に、プロビットモデルのモデル評価指標について紹介する。プロビットモデルを最尤法で推定した場合、OLSにおける決定係数 \(R^2\) を用いてモデルの当てはまりの良さを用いることはできない。そこで、通常、疑似決定係数（Pseudo \(R^2\)）を用いる事が多い16。疑似決定係数は0から1の間の値を取り、自身の立てたモデルの当てはまりが良いほど1に近づくという特徴を持っている。Rにおいては、DescTools::PseudoR2() を用いて計算可能である。install.packages("DescTools") でインストールしてほしい。
## McFadden
## 0.12933457.3.3 限界効果がなぜ必要か？*{#whymarginal}
ここでは、プロビットモデルによって得た係数の推定値を直接的に解釈できない理由と限界効果についての詳細を説明する。パラメータの推定結果と、\(x_{1i},...,x_{ki}\) の値を得たときに \(y_i\) が 1 を取る確率は以下の様に示される。
\[ \hat{P}(y_i=1|x_{1i},...,x_{ki})=\Phi(\hat{\beta}_0+\hat{\beta}_1x_{1i}+...+\hat{\beta}_kx_{ki}) \]
ここで、説明変数 \(x_{1i}\) の限界効果は、\(x_{1i}\) が変化したときに反応確率がどのように変化するのかを表す。これは偏微分という計算方法を用いて以下の様に示すことができる。なお、以下で示す内容は \(x_{1i}\) 以外の任意の説明変数を用いても成り立つ。
\[ \frac{\partial \hat{P}(y_i=1|x_{1i},...,x_{ki})}{\partial x_{1i}}=\phi(\hat{\beta}_0+\hat{\beta}_1x_{1i}+...+\hat{\beta}_kx_{ki})\hat{\beta}_1 \]
ただし、\(\phi\) は標準正規分布の確率密度関数であり、累積分布関数を微分することで得られる。このように、OLSのような線形モデルと異なり、\(\hat{\beta}_1\) が直接的に限界効果を示しているわけではないことが伺える17。
また同式より、ある変数の変化が反応確率へ与える限界効果は、個人が持つ全ての説明変数（\(x_{1i},...,x_{ki}\)）の値によって変化することもわかる。そのため、限界効果の報告においては「平均的な」限界効果を報告することが一般的である。ただし、「何の平均を取るか」という点において、2種類の計算方法が存在する。第1に、個人の個別限界効果を計算し、その平均値求めるという以下のような方法である。
\[\begin{equation} \frac{1}{n}\sum^n_{i=1}\Big[\phi(\hat{\beta}_0+\hat{\beta}_1x_{1i}+...+\hat{\beta}_kx_{ki})\hat{\beta}_1 \Big] \tag{7.2} \end{equation}\]
この計算方法は、因果推論における平均処置効果の議論に対応しているという好ましい性質を持っている（西山ほか, 2019）。
第2の方法は、説明変数の値について、それぞれ平均値を計算し、限界効果を求めるという以下のような方法である。
\[\begin{equation} \phi(\hat{\beta_0}+\hat{\beta}_1\bar{x}_1+...+\hat{\beta}_k\bar{x}_k)\hat{\beta}_1 \tag{7.3} \end{equation}\]
これは、それぞれの説明変数について平均値を取るような「平均的個人」における限界効果として解釈可能な計算方法である。
着目する説明変数が上のデータ分析結果における a1 のように離散変数（ダミー変数）の場合、その限界効果についてこれまでのような偏微分の議論は使えない。そのため、ダミー変数が 0 の場合の反応確率と、1 の場合の反応確率の差を取る形で限界効果を捉える。例えば、クーポンを受け取っていない場合（\(a_1=0\)）と受け取った場合（\(a_1=1\)）の反応確率の差は以下の様に示すことができる。
\[\begin{equation} \Phi(\hat{\beta}_0+\hat{\beta}_1p_{1i}+\hat{\beta}_2) - \Phi(\hat{\beta}_0+\hat{\beta}_1p_{1i}) \tag{7.4} \end{equation}\]
この限界効果も、他の変数の値（例えば、\(p_{1i}\)）によって変化するため、連続変数の場合と同様に平均値を報告することが一般的である。
ロジットモデルの場合は、
family = binomial(link = logit)で計算可能である。↩︎ロジットモデルの場合は、
mfx::logitmfx()で計算可能である。↩︎この指標は、最尤推定で用いる対数尤度関数 (7.5) に実際の観測値と最尤推定値を代入して求まる対数尤度の和（\(L_{full}\)）と、定数項だけを含むプロビットモデルを推定したときの対数尤度の和（\(L_0\)）を用いて、\(\Big(1-(L_{full}/L_0) \Big)\) と定義される↩︎
\(\hat{y}_i=\hat{\beta}_0+\hat{\beta_1}x_{1i}+...+\hat{\beta}_kx_{ki}\) という線形モデルの限界効果は、以下のようになる。 \[\frac{\partial \hat{y}_i}{\partial x_{1i}}=\hat{\beta}_1 \]↩︎