多項プロビットモデルの最尤推定をRで行う方法と数理、気を付けるべき点

今回は多項プロビットモデルです。少々高度な内容になります。

多項プロビットモデル
GHKアルゴリズム
Rによる多項プロビット回帰
多項ロジットとの違いは？(IIA特性)
- IIA仮定が妥当か否かの検定(ハウスマン・マクファーデン検定)
離散選択におけるその他の選択肢

多項プロビットモデル

具体例

皆さんは離散選択モデルというモデルをご存知でしょうか。

主に人間の意思決定をモデル化する際に利用されるモデルです。
例えば、コンビニに行って、何を購入するか決めるとします。
あくまで例ですが仮に「何か甘いものが食べたい」という漠然とした目的をもってコンビニに入った人だと考えましょう。

ここでこの人は普通「甘いものが売られている売り場」へ向かうはずです。
そして、パッケージ、値段、ラインナップ等々を吟味したうえでたくさんある選択肢の中から一つの商品を選び、レジへと持っていったとします。

この選択の過程を今回は多項プロビットモデルというモデルを使って解析します。

数式

さっきの話を簡単にまとめると要は

パッケージ、値段等の情報を適当な重みを付けて比較し、「最もその人にとって良いと思える商品を１つ手に取った」という状況ですね。

つまり、その人にとっての何らかの選択基準が存在し、その選択基準に基づいて最も良いと思える商品を手に取ったと考えられます。

この時の個人 $i$ の商品 $j$ に対する選択基準 $Y_{i,j}^{*}$ を以後効用と呼ぶことにします。
この効用とは経済学の概念で、例えばAさんが商品1を買うか商品2を買うか迷っているとした場合、Aさんはその商品に関する情報を使って効用 $Y_{A,1}^{*},Y_{A,2}^{*}$ を頭の中で作り上げます。

そしてもし、 $Y_{A,1}^{*}＞Y_{A,2}^{*}$ なら商品1を購入し、 $Y_{A,1}^{*}＞Y_{A,2}^{*}$ ならば商品2を購入するといった状況を考えるわけです。

ここで、「いやいや、僕はそんなよく分からない効用なんてものを考えて行動していないし、そんなモデル意味が無いよ」という批判は一旦置いておきます。*1

この効用 $Y_{i,j}^{*}$ は、誤差項 $\epsilon_{i,j}$ とその人が持ちうる情報 $I_{i,j}$ によって次のような線形関係があると考えます。

$Y_{i,j}^{*} = V_{i,j} + \epsilon_{i,j}$
ただし、 $\epsilon_{i,j}$ は期待値 $0$ 、分散共分散行列 $\Sigma$ の多変量正規分布に従うと仮定します。

更に、 $V_{i,j} = \beta_{0}^{(i)} + \beta_{1}^{(i)} X_{(i,j),1} +\beta_{2}^{(i)} X_{(i,j),2} +,\cdots ,+\beta_{k}^{(i)} X_{(i,j),k}= X^{(i,j)}\beta^{(i)}$
と仮定します。ここで $X$ とは値段や商品配置、パッケージ等の購買に関する判断基準を指します。
$\beta^{(i)}$ は個人 $i$ が暗黙的もしくは無意識的に設けているパラメータベクトルで、 $X^{(i,j)}$ は行列です。

当然添え字のつけ方で、各人共通のパラメータを設けたりすることも可能になります。
よって、 $V_{i,j}$ を条件付けた際の $Y_{i}^{*}$ の分布 $P(Y_{i}^{*} | V_{i,j},j=1,2,\cdots , k)$ は次のようになることが分かりますね。

$P(Y_{i}^{*} | V_{i,j},j=1,2,\cdots , k) = \frac{1}{(2 \pi)^\frac{k}{2}\sqrt{ |\Sigma| }}exp(-\frac{1}{2}(Y_{i}^{*}-V_{i})' \Sigma^{-1} (Y_{i}^{*}-V_{i}))$

ただし、 $Y_{i}^{*},V_{i}$ はそれぞれ効用と判断基準を並べた $j$ 次列ベクトルとします。

与えられているデータ

ここで問題になってくるのは

効用 $Y_{i}^{*}$ の値なんて観測できなくない？？ってことです。
$i$ さんにとってのチ〇ルチョコの効用は10,チョコパイの効用は20だから、 $i$ さんはチョコパイを購入したんだな～なんてことはわからなくて、実データでは、その辺の話を全部すっ飛ばして

$i$ さんは「チョコパイを購入した」という情報だけ得られます。

もし、効用まで観測できれば $P(Y_{i}^{*} | V_{i,j},j=1,2,\cdots , k)$ に対して最尤法を適用しておしまいですが、観測出来ていないので簡単にはいきません。そこで、少し工夫する形で別の尤度を考えます。

今回の $i$ さんは商品 $J$ を購入したとしましょう。

購買情報を踏まえた尤度

$i$ さんが商品 $J$ を購入したということは、次の状況が生じているはずです。

$Y_{i,J}^{*} ≧ Y_{i,j}^{*}$ for $j = 1,2,...k$
こう書き換えても問題ありません。
$Y_{i,J}^{*} = max\{Y_{i,j}^{*} |j=1,2,...k \}$

更に、 $i$ さんの購買情報 $\omega_{i}$ を次のように考えます。
$\omega_{i} = (\omega_{i,1},\cdots,\omega_{i,k})$
ここで、 $\omega_{i,J} = 1$ 、 $\omega_{i,t} = 0$
ただし、 $t≠J$

つまり、購入した商品には $1$ 、購入しなかった商品には $0$ を与えたダミー変数ベクトルを用意するわけです。

そうすると、購買情報 $\omega_{i}$ における分布 $Prob(\omega_{i} | V_{i,j},j=1,2,\cdots , k)$ は次のようになるはずですね。

$Prob(\omega_{i} | V_{i,j},j=1,2,\cdots , k)$
$= \int_{\{Y_{i,1}^{*}＜Y_{i,J}^{*}\}},\cdots \int_{\{Y_{i,J}^{*}=-\infty\}}^{\{Y_{i,J}^{*}=+\infty\}},\int_{\{Y_{i,J+1}^{*}＜Y_{i,J}^{*}\}}\cdots, \int_{\{Y_{i,k}^{*}＜Y_{i,J}^{*}\}} P(Y_{i}^{*} | V_{i,j},j=1,2,\cdots , k) dY_{i,1}^{*},\cdots,dY_{i,k}^{*}$

何かえらいこっちゃなってますが、要は多変量正規分布の密度を切断している分布になります。

まあ、当然ですがこの $Prob(\omega_{i} | V_{i,j},j=1,2,\cdots , k)$ はオープンフォームになっているため、普通の方法で最尤推定量を求めることは出来ません。

そのため、プロビットモデルは長く日の目を拝むことなく、より計算が簡単なロジットモデルが発展しました。

しかし、今日ではコンピュータの計算速度の向上により、このプロビットモデルの最尤推定値を探索的に求めることが可能になっています。

GHKアルゴリズム

最尤推定値を探索的に求める方法は様々ありますが、例えばフィッシャーのスコア法では対数尤度の二回微分に関する期待値を利用します。
しかし、肝心の尤度がゴツすぎるため、微分したものの期待値なんてもってのほか。求めることが出来ません。
そこで、一旦乱数を取り出してやり、期待値の近似値を求めます。その近似値を利用して、探索解を求めるという手順を取ります。

その際、乱数の発生させ方に一工夫必要で、その取り出し方をGHKアルゴリズムと呼びます。

もう積分記号なんか見たくないとは思いますが、再度先ほどの尤度を書いてみます。

ここで、
$Y_{i,J}^{*} ≧ Y_{i,j}^{*}$ for $j = 1,2,...k$
$Y_{i,j}^{*} = V_{i,j} + \epsilon_{i,j}$

のため、
$V_{i,J} + \epsilon_{i,J} ≧ V_{i,j} + \epsilon_{i,j}$
更に、
$V_{i,J} -V_{i,j} ≧ \epsilon_{i,j} - \epsilon_{i,J}$

ここで、 $\epsilon_{i,j} - \epsilon_{i,J} =E_{i,j^*}$ for $j^* = 1,2,...k-1$
$E[E_{i,j^*}] =0$ , $V[E_{i,j^*}] =\Sigma^*$ であるとします。

多変量正規分布には加法性があるので、 $E_{i,j^*}$ の組も多変量正規分布に従います。
つまり $E_{i}$ は多変量正規分布に従います。
さてここで、分散共分散行列 $\Sigma^*$ に対してコレスキー分解を施し、
$\Sigma^{*} = \Gamma \Gamma^{'}$ と置くと、

$E_{i} =\Gamma \xi$
ただし、 $\xi$ は各要素がi.i.dに標準正規分布に従う確率変数ベクトル。

$\Gamma$ は上三角行列なので次のように書くことが出来ますね。

$E_{i,1} =\Gamma_{1,1} \xi_1$
$E_{i,2} =\Gamma_{1,2} \xi_1+ \Gamma_{2,2} \xi_2$
$E_{i,2} =\Gamma_{1,3} \xi_1+ \Gamma_{2,3} \xi_2+ \Gamma_{3,3} \xi_3$
$\vdots$

つまり、長くなりましたが何が言いたいかというと、多変量正規分布に従う確率変数は、期待値ベクトル(今回は0)と分散共分散行列さえわかっていれば、普通の標準正規分布から得られるわけです。

こうして、対象の多変量正規分布からサンプリングを行い、更に必要な期待値計算をモンテカルロ近似によって計算し、それを用いて適当な最適化アルゴリズムを使うというのがGHKアルゴリズムの概要です。

このアルゴリズム、上の説明を見ても分かる通り
ある一つの選択肢を基準に他の

Rによる多項プロビット回帰

「GHKアルゴリズム、、なんて大変なんだ」と思った方もご安心。Rではmlogitパッケージを使って簡単に推定することが出来ます。

使用するデータ

data("TravelMode", package = "AER")

> head(TravelMode)
  individual  mode choice wait vcost travel gcost income size
1          1   air     no   69    59    100    70     35    1
2          1 train     no   34    31    372    71     35    1
3          1   bus     no   35    25    417    70     35    1
4          1   car    yes    0    10    180    30     35    1
5          2   air     no   64    58     68    68     30    2
6          2 train     no   44    31    354    84     30    2

AERパッケージのTravelModeデータです。
mlogitパッケージのExampleで使っていたので、そのまま使うことにします。
https://cran.r-project.org/web/packages/mlogit/mlogit.pdf

ちなみに各変数は
individual : 個人iの識別番号(i=1,...200)
mode : 乗り物の種類(選択肢)
choice : その乗り物を選んだか否か
wait : 乗り場における待ち時間(当然自家用車なら0)
vcost : 乗り物を使うのにかかるコスト
travel : かかる時間
income : 個人の収入
http://127.0.0.1:23970/library/AER/html/TravelMode.html

今回のデータはchoice毎に行があり、individualが重複する形の縦に長いデータとなっています。
このようなデータフレームをlong型と呼びます。ちなみに1行で個人データをひとまとめにしているデータをwide型と呼び、例えば次のようなものがあります。

> head(Fishing)
     mode price.beach price.pier price.boat price.charter catch.beach catch.pier catch.boat catch.charter   income
1 charter     157.930    157.930    157.930       182.930      0.0678     0.0503     0.2601        0.5391 7083.332
2 charter      15.114     15.114     10.534        34.534      0.1049     0.0451     0.1574        0.4671 1250.000
3    boat     161.874    161.874     24.334        59.334      0.5333     0.4522     0.2413        1.0266 3750.000
4    pier      15.134     15.134     55.930        84.930      0.0678     0.0789     0.1643        0.5391 2083.333
5    boat     106.930    106.930     41.514        71.014      0.0678     0.0503     0.1082        0.3240 4583.332
6 charter     192.474    192.474     28.934        63.934      0.5333     0.4522     0.1665        0.3975 4583.332

mlogit.data関数

mlogitを使用する前に、解析で使うデータフレームをmlogitに適用しやすい形に変換しておくことが無難です。
その際の変換に用いるのがmlogit.data関数です。

ただし、別にmlogitを使用する際に関数内で適切な指定をしてやれば変換していないデータフレームに実行をかけることも出来ます。

TM <- mlogit.data(TravelMode, choice = "choice", shape = "long",
alt.levels = c("air", "train", "bus", "car"),id.var="individual",)

choiceにはlong型の場合は選択肢を選んだか否か(yes or no ; 1 or 0)に関する変数を指定します。wide型の場合は選んだ選択肢の入った変数を指定します。
shapeは"long" か "wide" か使用するデータの型を選択します。
alt.levelsで選択肢を改めて指定しておく必要があります。long型の場合は必須ですが、wide型の場合は必要ありません。
id.varでlong型の場合に個人識別をしておく必要があります。

mlogit関数

mlogit関数で実行します。デフォルトはただのlogitなのでprobitで行う場合はprobit=TRUEとします。

pr <- mlogit(choice ~ wait + vcost | income |travel , TM,
probit = TRUE)

formulaを|で区切っていますが、これはそれぞれ次のような意味になっています。
個人iの選択jに関する効用 $U_{i,j}$ について
$U_{i,j} = \alpha_{(0,j)} + \beta_w Wait_{(i,j)} +\beta_v Vcost_{(i,j)} + \gamma_{j} Income_{i} + \omega_{j} Travel_{(i,j)}$

添え字に注目してください！

実行＆結果

> summary(pr)

Call:
mlogit(formula = choice ~ wait + vcost | income | travel, data = TM, 
    probit = TRUE)

Frequencies of alternatives:
    air   train     bus     car 
0.27619 0.30000 0.14286 0.28095 

bfgs method
22 iterations, 0h:0m:11s 
g'(-H)^-1g = 6.11E-07 
gradient close to zero 

Coefficients :
                     Estimate  Std. Error z-value  Pr(>|z|)    
train:(intercept)  0.11400922  0.50036804  0.2279 0.8197623    
bus:(intercept)   -0.25750600  0.46140096 -0.5581 0.5767789    
car:(intercept)   -1.70532857  0.67252768 -2.5357 0.0112223 *  
wait              -0.02591688  0.00596662 -4.3436 1.401e-05 ***
vcost             -0.00571812  0.00332269 -1.7209 0.0852629 .  
train:income      -0.02256556  0.00677919 -3.3287 0.0008727 ***
bus:income        -0.01125839  0.00862617 -1.3051 0.1918437    
car:income        -0.00637486  0.00788242 -0.8087 0.4186627    
air:travel        -0.01473926  0.00360134 -4.0927 4.264e-05 ***
train:travel      -0.00291735  0.00079889 -3.6517 0.0002605 ***
bus:travel        -0.00313288  0.00089814 -3.4882 0.0004863 ***
car:travel        -0.00296051  0.00089036 -3.3251 0.0008840 ***
train.bus          0.94852133  0.38768124  2.4467 0.0144190 *  
train.car          1.06025359  0.37869969  2.7997 0.0051147 ** 
bus.bus            0.38464898  0.16055087  2.3958 0.0165838 *  
bus.car            0.27405472  0.37816033  0.7247 0.4686330    
car.car            0.57842389  0.24356722  2.3748 0.0175584 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -172.56
McFadden R^2:  0.39187 
Likelihood ratio test : chisq = 222.39 (p.value = < 2.22e-16)

地味に"gradient close to zero" が重要です。これが"last step couldn't find higher value"となっている場合は問題有です。
coefficientsの最後の方に並んでいる(選択肢).(選択肢)はそれぞれの分散共分散です。一個分散が欠けているのはスケールを揃えるのに用いているためだと思われます。
McFadden R^2は一応経験的な適合度指標として使われています。あまり紋切型で考えるのもイマイチですが、0.5くらいが目安なので今回は少し小さめ。

各個人の選択確率予測値はprobabilitiesに入っています。

> head(pr$probabilities)
           air       train         bus       car
[1,] 0.1463981 0.230907011 0.069925608 0.6030067
[2,] 0.3550429 0.194959842 0.012583377 0.4531929
[3,] 0.4270799 0.025918716 0.052778939 0.4559133
[4,] 0.4102070 0.007009033 0.004554783 0.5633459
[5,] 0.2017710 0.479339802 0.114605173 0.1839117
[6,] 0.1113763 0.478200196 0.122358769 0.3097929

多項ロジットとの違いは？(IIA特性)

何故よくわからんGHKアルゴリズムを使ってまで(計算に時間をかけてまで)probitを使うのだ。多項ロジットで良いではないかと思う人もいると思うので違いを説明しておくと

多項ロジットの場合は、誤差項の共分散が無いので、各選択肢の代替関係が同等であるという仮定が入っています。
例えば車か飛行機かを選ぶにあたって、その選択肢の中に追加でバスが入ってきたとします。この場合、例えば陸路で安価に行けた方が嬉しいなら車の選択肢を選んでいた人がバスを選ぶようになる確率よりも飛行機の選択肢からバスに移る確率の方が多いように感じますし、旅行先が渋滞まみれのような場所では車で目的地に向かうのは不便なので、飛行機よりも車を選んでいた人の方がバスに流れていくはずです。このような変化に多項ロジットは対応出来ません。

これをIIA特性（無関係な選択肢からの独立）と呼びます。

これを解決するためにはネストロジットという手法もありますが、多項プロビットでも解決することが出来ます。

IIA仮定が妥当か否かの検定(ハウスマン・マクファーデン検定)

IIA仮定が妥当であれば多項ロジットを使う方が計算の点からも良く、またロジットは選択肢が実際には観測されていない場合でも機能するため、多項ロジットだって無下には出来ません。

そこで仮説検定を行います。
Hausman and McFadden (1984)
EconPapers: Specification Tests for the Multinomial Logit Model
で提案されたハウスマン・マクファーデン検定では

帰無仮説: IIA特性を持つ
対立仮説 : IIA特性は持たない

という検定を行います。検定においてはフルモデルでの推定結果と、選択肢を一個取り除いた場合の推定結果を比較します。
もしIIA特性があれば、直接関係のない選択肢が消えてもオッズ比は(誤差の範囲で)自然な値になるはずです。

mlogitパッケージではこの検定をhmftest関数として与えています。

> x <- mlogit(choice ~ wait + vcost | income |travel , TM,reflevel="car")
> g <- mlogit(choice ~ wait + vcost | income |travel , TM,reflevel="car",alt.subset=c("car","bus","train"))
> hmftest(x,g)

        Hausman-McFadden test

data:  TM
chisq = 27.572, df = 9, p-value = 0.001124
alternative hypothesis: IIA is rejected

> g2 <- mlogit(choice ~ wait + vcost | income |travel , TM,reflevel="car",alt.subset=c("car","bus","air"))
> g3 <- mlogit(choice ~ wait + vcost | income |travel , TM,reflevel="car",alt.subset=c("car","air","train"))
> hmftest(x,g2)

        Hausman-McFadden test

data:  TM
chisq = 178.24, df = 9, p-value < 2.2e-16
alternative hypothesis: IIA is rejected

> hmftest(x,g3)

        Hausman-McFadden test

data:  TM
chisq = 41.746, df = 9, p-value = 3.656e-06
alternative hypothesis: IIA is rejected

選択肢を除いたmlogitとフルモデルのmlogitで比較を行っています。
今回はどの場合でも棄却(p-value <0.05) されたため、IIA特性は妥当とは言えないということです。
よって、今回は多項プロビットモデルを利用しています(当然他のIIA仮定を緩めたロジットを利用しても良い)

離散選択におけるその他の選択肢

ここは、これから解説記事を書く中で少しずつリンクを貼ってこうと思います。
とりあえずはざっと羅列。

①多項プロビット
今回説明したもの

②多項ロジット回帰
IIA特性があるが、解釈のしやすい離散選択モデル
クローズドフォームなので計算が楽

②ネステッドロジット
IIAを緩和すべく、選択肢をネストしたロジット
選択の順序関係は考慮していないことに注意

③ネステッドロジットのIIAを更に緩めたモデル
paired combinatorial logit モデル
ordered generalized extreme value モデル
cross-nested logit モデル

④generalized nested logit
上のネストロジットを一般化したモデル
IIAはかなり緩和されるが、
パラメータの置き方や解釈の仕方が難解で計算量も多い

⑤C-logit モデルとpath size logit モデル
確定項の方で相関を表現したモデル

②～⑤までを総称してGEVモデルと呼ぶ。
GEVと多項プロビットの大きな違いは、確率項の大きさが異なっていても多項プロビットは利用できることで、GEVでは確率項の大きさが同一であることを仮定している。

⑥離散連続選択モデル
離散選択を行った後に連続値の選択を行う場合には離散連続選択モデルが利用される。これについてはmultiple discrete-continuous extreme value モデルとそれを一般化した generalized multiple discrete-continuous extreme value モデルがある。

このまとめについては下pdf参照
http://www.trans.civil.nagoya-u.ac.jp/~yamamoto/papers/TE2012.pdf

また、このpdfにはありませんが
⑦多変量プロビットモデル

というものもあります。
これについてもいつか記事にします。

*1:この記事を全部読めば、その手の批判は的外れであるということが分かるかと思います。ただ、この時点でそういう批判が出てしまう気持ちはよくわかります。