ill-identified diary

所属組織の見解などとは一切関係なく小難しい話しかしません

トップ > R > [発展編] 多項ロジットの話

[発展編] 多項ロジットの話

R 離散選択 STATA 計量経済学 教材
この記事は最終更新日から3年以上が経過しています

前回の多項ロジット (混合ロジット) の話に引き続いて, 多項ロジットが使えない場合の手法について書いておく. 今回言及するのは:

  • IIAの仮定を検証する方法
  • ネステッド・ロジット
    • 実行方法
  • 一般化極値 (GEV) モデル
  • 混合モデル (Mixture Model)
    • 実行方法

である. 今回はポンポン数式が出てくる. ある程度の知識, 特に大学学部レベルの数理統計や微積分の知識が必要な内容になった.

IIAの仮定に対する検定

そもそも IIA というのは, 多項ロジットモデルを利用する上で「単なる足かせでしかない」ということはなく, 利点を発揮することもある. 確かに IIA の仮定は, 時には現実の分析対象に対してあからさまに不自然な仮定に見えることもある. しかし一方で, もし IIA の仮定が成り立っているのなら, 観察できなかった選択肢を除いたまま多項ロジットで推定しても問題ないという利点がある. こういった利点があるため, 多項ロジットは忌避されてるわけでもないし, 多項ロジットは IIA の仮定を要求するから不完全な手法, ということにはならない.

では, どうやって IIA が成り立っているのか, いないのかを判断すればよいだろうか. その答えが Hausman and McFadden (1984) によるハウスマン検定である. まず大雑把な捉え方をするなら, ハウスマン検定は, 色違いのバスでの譬え話での説明と同じことをしている. つまり, 選択肢が増減した場合を比較し, まったく関係ない他の選択肢のオッズ比が不自然な変化をしていないかどうか, ということを確率的に検証している.

帰無仮説は, 「IIA が成り立つ」である. このとき, 全てのサンプルを使ってふつうに推定した多項ロジットのパラメータの推定値を \boldsymbol{b} とし, 選択肢のいずれかを適当に削除 (前回の解説参照. サンプルの一部を取り除くことになる) した上で推定したパラメータを \boldsymbol{b}_A とおく. 帰無仮説上はこの2つのパラメータは等価である. \boldsymbol{b}_A のほうが選択肢1つぶんパラメータが減っているので, \boldsymbol{b} から \boldsymbol{b}_A に対応するパラメータだけを抜き出した \boldsymbol{b}_H で比較する. 同様に, これらに対応する共分散行列を \mathbf{V}_H, \mathbf{V}_A とすると,

{\displaystyle [\boldsymbol{b}_H -\boldsymbol{b}_{A}]^\prime [\mathbf{V}_H-\mathbf{V}_A ] [\boldsymbol{b}_H - \boldsymbol{b}_A ] }

は漸近的に自由度 N - K のカイ2乗分布に従う. N はサンプルサイズ, K はパラメータの数である.

実行方法

STATA*1 には, IIA の仮定に対するハウスマン検定を一発でやってくれるコマンドは用意されてないため, フルサンプルの多項ロジットとサブサンプルの多項ロジットをそれぞれ推定しその結果を保存し, hausman コマンドで検定を行うことになる.

* フルサンプルで推定
mlogit y x1 x2 ....
est sto full
* サブサンプルで推定 (j には適当な選択肢の値)
qui mlogit y x1 x2 ... if (y != j )
hausman . full

このコードでは, サブサンプルでの多項ロジットは検定以外で使い道がないので, qui で結果を非表示にしている.

R の場合, mlogit パッケージに hmftest() が存在し, これで多項ロジットの IIA 仮定の検定ができる. 内部では STATA と同様, 2種類の多項ロジットを mlogit() 関数で推定した結果を与えているようが, 使用法にはバリエーションがあり,

  1. 2種類の mlogit オブジェクトを代入する. このとき, 2つは同一のモデルを使わなければならない. つまり, 片方を alt.subset= オプションを使って推定することになる.
  2. フルサンプルを推定した mlogit オブジェクトと, alt.subset= 与える.

というやり方があるので, 多少自由度がある.

IIAの仮定が棄却された場合の対処法1: ネステッド・ロジット

ネステッド, あるいは「入れ子」 (nested) ロジットモデルは, コンセプトは自体は単純で, 選択肢を選ぶ確率を定義してきた多項ロジットを拡張して, 選択肢の集合を選ぶ確率を定義する, というものである. つまり, j という点の確率でなく, 次の y \in \mathcal{G}_s という集合の中のいずれかである確率を求める:

初めに, 選択肢が S 個のグループに分けられるとする. S 個のグループのをそれぞれ \mathcal{G}_1,\cdots,\, \mathcal{G}_S という部分集合で表す. つまり, J 個の選択肢を S 個の集合 \mathcal{G}_1,\,\mathcal{G}_2,\,\cdots,\,\mathcal{G}_S に分割するということである. この部分集合を, 以降ネスト (nest) と呼ぶ. ただし, ネストに分割するには以下の条件を満たしていないとならない.

全てのネスト \mathcal{G}_s について, このネストに属する任意の2つの選択肢に対して IIA の仮定が成り立つ.

多項ロジットの IIA の仮定では, 選択肢全体のなかで, 適当に2つのペアを抜き出したとき, 常にこのペアのオッズ比が他の選択肢の影響を受けない, というものだったが, 入れ子ロジットでは, ネストが異なる場合は IIA でなくてもよい. 例えば, 前回のタクシーと色違いのバスでは, 「タクシー」という選択肢だけが含まれるネスト \mathcal{G}_1と, 「赤いバス」「青いバス」の2つの選択肢が含まれるネスト \mathcal{G}_2 の2つに分割する. すると, 「タクシー」と「赤いバス」は異なるネストだから, IIAの仮定は必要ない (= この2つのオッズ比が「青いバス」の影響を受けても前提条件に違反しない). そして, ネスト内では IIA の仮定が成り立っている必要があるので, 「赤いバス」「青いバス」の2つの選択肢のオッズ比は, 「タクシー」 (あるいは他の未知の選択肢) の影響を受けてはならない. いわば, ネステッド・ロジットで要求されるのは IIA の仮定ではなく, IIN, つまり Independent from Irrelevant Nest の仮定である.

図1: ネスト内でのみIIAが成り立つ

f:id:ill-identified:20140722014625p:plain

このような仮定の場合, それぞれのネスト \mathcal{G}_sの中で, ある選択肢が選ばれる確率は次のように表現できる.

{ \displaystyle \mathrm{P}(y=j | \boldsymbol{x},y\in\mathcal{G}_{s})=\exp(\rho_{s}^{-1}\boldsymbol{x}_{j}^{\prime}\boldsymbol{\beta})\bigg/\sum_{h\in\mathcal{G}_{s}}\exp(\rho_{s}\boldsymbol{x}_{h}^{\prime}\boldsymbol{\beta}) }

これだけではパラメータの推定には不十分である. ネスト \mathcal{G}_sが選ばれた条件のもとで, ある選択肢が選ばれる, という条件付き確率しか表現できていない. そこで, 次にネスト \mathcal{G}_s のうちいずれが選ばれるかの確率も考える必要がある.

{\displaystyle \mathrm{P}(y\in\mathcal{G}_{s}\big|\boldsymbol{x}) = \alpha_{s}\left\{ \sum_{j\in\mathcal{G}_{s}}\exp(\rho_{s}^{-1}\boldsymbol{x}_{j}^{\prime}\boldsymbol{\beta})\right\}^{\rho_{s}}\left[\sum_{r=1}^{S}\alpha_{r}\left\{ \sum_{i\in\mathcal{G}_{s}}\exp(\rho_{r}^{-1}\boldsymbol{x}_{i}^{\prime}\boldsymbol{\beta})\right\} ^{\rho_{s}}\right]^{-1} }

ここで, \alpha_s は標準化のためのパラメータである. この2つの条件付き確率を掛けることで条件を打ち消し(観察された情報の \boldsymbol{x} だけが条件として残る), 個別の選択肢の選択確率が得られる. よって, ネステッド・ロジットもこれまでと同様に最尤法として推定が可能になる.

実行方法

STATA でやる場合, 日本語だとNested logit モデルについてで紹介されているので参考にしてください.

R でやる場合, mlogit パッケージの mlogit 関数で, nests= にネスト構造を代入すると, 通常の多項ロジットではなく, ネステッド・ロジットで推定される. 具体的には, list()を使い, 各ネストごとにベクトルで表記する.

mlogit(...., n
    ests = list(taxi = 'taxi', bus = c('rad.bus', 'blue.bus'), ...)
)

というふうにやる.

さらにネストの階層を増やす場合は未確認...

対処法2: GEVモデル

ここまでが, 一般的な教科書や授業で必ずといってもいいほど説明される部分. しかし, これ以降の手法の発達が触れられることが少ない*2. かつては, 離散選択といえば Amemiya (1985) という風潮だった. たしかに雨宮先生はこの分野の先駆者だし, 実際この本も離散選択のことはほとんどなんでも書いてある, いや書いてあった. いかんせん30年近く前のものなので, 最新の研究が反映されていない. これ以降の話は主に Train (2009) を参考にしたものである. こちらは著者が自らのサイトで pdf を全文公開するという太っ腹なことをしているので, ありがたくダウンロードさせていただこう.

ネステッドロジットは, 選択肢を重複せずにきれいに分割するようにしていた. だが, もっと複雑な場合はどうすればいいか. たとえば, ミクロデータを用いた計量分析で常につきまとう, いわゆる, 観測されない個体ごとの異質性 (unobserved heterogeneity) の問題もある. このように, さらに複雑な多項選択を記述するためのモデルは, 一般化極値 (GEV) モデルと呼ばれる (よってネステッド・ロジットも GEV に含まれる) ものに集約されることが McFadden (1978) で示されている.

基本的には, GEV モデルも尤度関数で表現できれば, 最尤法でパラメータを求めることができる. しかし, GEV の場合, 多くは選択確率 P_{i,j} が closed form とならず, 解析的に解けないため, 確率を数値積分によって求める必要があり, 計算量が大幅に増えることになる.

例えば, 多項ロジット, ネステッド・ロジットいずれも, 各選択肢に対応する誤差項は同一の分布であることを仮定したが, 選択肢ごとに分散が異なる, いわゆる分散不均一性 (heteroskedasticity) が存在する場合を考慮したものは GEV モデルの中でも Heteroskedastic Extreme Value (HEV) モデルと呼ばれる. これは, 誤差項の分散が異なることから, \varepsilon_j の分散値が \mathrm{V}( \varepsilon_{j} ) = \theta_{j}^{2}\pi^{2}/6 として新たなパラメータ \{\theta_{j}\} を用いて表現する. HEVモデルは, これを標準の分散値 \pi^{2}/6 に正規化することで対処する. 極値分布 (ガンベル分布) の分布関数は一般に, \exp \{ - \exp [(x-\mu)/\theta ] \} と表し, 平均と分散はパラメータ\mu,\,\theta によって, \mathrm{E}X=\mu, \theta^{2}\pi^{2}/6 となる. よって, \mathrm{V}(\varepsilon_{j}/\theta_{j})=\pi^{2}/6 である. ここから, p_{i,j} は

{\displaystyle p_{n,i} = \int \left[\prod_{j\neq i} \exp\left( -\exp(u_{n,i} - u_{n,j} +\theta_iw)/\theta_j\right) \right] \exp\left(-\exp(-w) \right)\mathrm{e}^{-w}dw }

とする. ただし, w=\varepsilon_{n,i}/\theta_i である.

要するに, WLS, GLS と似たような発想だが, 解析的に解けるこれらと違い, 積分を シミュレーションで計算することが必要になる.

McFadden (1978) で示された一般形の GEV モデルでは, 各選択肢の効用 u_{j} の関数 G(u_1,\,u_{2},\,\cdots,\,u_{J}) が一定の条件を満たすものであるなら, 選択確率 p_{j} を {\displaystyle p_{j} = \frac{u_{j}G_{j}}{G} } で表すことができ, 極値分布モデルとして推定できるということを示している ( G_{j}:=\partial G/\partial u_{j} ). たとえば, G=\sum_{k}u_{k} とすれば多項ロジットになるし, G=\sum_{s}(\sum_{k}u_{k}^{1/\lambda_{k}})^{\lambda_{s}} とすればネステッドロジットとなる. で, このある条件というのは, ここでは長くなるので詳しくは言わないが, 効用関数や生産関数など, 経済学の最適化問題でよく用いられる目的関数の一般的な仮定とよく似たものになる. これは, 多項ロジットが経済学の実証ツールとしてよく使われる理由の1つかもしれない.

対処法3: Mixture Model

IIA の仮定を完全に緩和したロジットモデルは Mixture Model *3 と呼ばれるものである. mixture model は, GEVと同様に, 選択確率が積分の形で表現される. まず, パラメータ \boldsymbol{\beta} がある分布に従うとする. この分布の密度関数を f(\boldsymbol{\beta}) とおく. また, 所与の \boldsymbol{\beta} に対する選択肢 j の選択確率は, {\displaystyle L_{j}(\boldsymbol{\beta}):=\frac{\exp(v_j)}{\sum_{k}\exp(v_k)} } とする. ただし, 簡略化のため, 本題にあまり影響しない, 個人を表すインデックス i は省略している. すると, 選択確率 p_j は次のような積分で表すことができる. {\displaystyle p_j = \int L_j(\boldsymbol{\beta})f(\boldsymbol{\beta})d\boldsymbol{\beta} =\int\left[\frac{\exp(v_{j}(\boldsymbol{\beta}))}{\sum_{k=1}^{J}\exp(v_{k}(\boldsymbol{\beta}))}f(\boldsymbol{\beta})\right]d\boldsymbol{\beta} }

ここで, v_{j} は選択肢 j に対する効用であり, パラメータに対して線形, つまり, v_{j}:=\boldsymbol{x}_{j}^{\prime}\boldsymbol{\beta} とするのが一般的である. さらに, f(\boldsymbol{\beta}) が, \boldsymbol{\beta}=\boldsymbol{b} のとき 1に, それ以外の時にゼロとなるように定義すると, パラメータが固定されるから, 通常の多項ロジット(厳密には mixed logit)と等価になる. しかし, 一般には, 上記のように積分で表されるため, mixture model では多項ロジットのようにオッズ比 p_{j}/p_{i} が \boldsymbol{x}_{j},\,\boldsymbol{x}_{l} だけで決まるとは限らない. よって, mixture model では IIA の仮定は要求されない.

尤度関数は, 今まで同様, インジケータ関数と, 省略してきた i を用いて, {\displaystyle \ln L = \sum_{i}\sum_{j}1(y_{i} = j) \ln p_{i,j} } となる.

なお, f(\boldsymbol{\beta}) には正規分布, 対数正規分布, ガンマ分布のような連続分布が採用されることが多いが, 離散分布が用いられることもある. 特に離散分布を用いた mixture model は, latent class model (潜在クラスモデル) と呼ばれる.

計算方法

mixture model は, \boldsymbol{\beta} が分布すると仮定したので, 元々のパラメータに加え, この分布のパラメータが加えられる. これを \boldsymbol{\theta} とおくと, p_j は {\displaystyle p_j = \int L_{j}(\boldsymbol{\beta})f(\boldsymbol{\beta}|\boldsymbol{\theta})d\boldsymbol{\beta} } というふうに書き表せる. ここで, 積分によって \boldsymbol{\beta} は消えるから, p_j は \boldsymbol{\theta} の関数になる. よって, 尤度関数は \boldsymbol{\theta} の関数となるので, ここではまず \boldsymbol{\theta} を求めることになる.

実際には, p_{j} の積分は解析的に得られないので, シミュレーションによってこの部分を計算する. よって, 厳密には最尤法ではなく擬似最尤法 (シミュレーション最尤法) を用いることになる: 所与の密度関数 f(\boldsymbol{\beta}|\boldsymbol{\theta}) から R 個の乱数 \boldsymbol{\beta}^{(1)}, \boldsymbol{\beta}^{(2)}, \cdots , \boldsymbol{\beta}^{(R)}を生成し, この平均値を p_{i,j} の推定値 \check{p}_{i,j} として, 代用するのである.

実行方法

R では Train and Croissant (2013) によって mlogit パッケージでの mixture model (引用先では mixed logit model と表現) での推定方法が書かれている. これもまた mlogit() 関数で行う. rpar= およびR= 引数で, それぞれ各パラメータに対応する分布と繰り返し回数を指定できる. たとえば,

mlogit(model, data, rpar=c(x1 = 'n', x2 = 'n', ...), R=100)

と言うふうにすれば, 各パラメータが正規分布するモデルを, 100回のサンプリングで計算する. 分布は, 正規分布の 'n' のほか, 対数正規分布の 'l', 切断正規分布の 't', 一様分布の 'u' を指定できる.

STATA では不明. ユーザの作成した latent class model を EMアルゴリズムで計算するモジュールは存在するが, 少なくとも最初から用意されているわけではさそうである. また, STATA ではこのブログと同じように, mixed logit が多項ロジットの一般形を表すものとしてコマンドが実装されているので, 紛らわしい.

参考文献

  • Hausman, Jerry A. and Daniel L. McFadden (1984), ''A Specification Test for Multinomial Logit Model,'' Econometrica 52, 1219-1240
  • Wooldridge, Jeffrey M. (2010) ''Econometric Analysis of Cross Section and Panel Data'' 2nd Ed., The MIT Press
  • Amemiya, Takeshi (1985) ''Advanced Econometrics,'' Harvard University Press

*1:少なくとも STATA 12 あたりまでは

*2:Wooldridge (2010) でも, 何某で提案された方法があるよ, と書かれている程度

*3:Train (2009) などでは Mixed Logit と呼ばれるが, それでは前回紹介した混合ロジットと区別がつかないため, Wooldridge (2010) で使われている mixture という表現を採用した

ill-identified 読者になる

広告を非表示にする
コメントを書く