【学習動機】
前回の続きで少し複雑なロジットモデルをRで推定するために, {mlogit}パッケージの使い方を覚えたい.
調べていたら, Train先生が{mlogit}の使い方のエクササイズを書いているのを発見した.
Kenneth Train’s exercises using the mlogit package for R
http://www.idg.pl/mirrors/CRAN/web/packages/mlogit/vignettes/Exercises.pdf
これに沿って学習していく. 今回は, 混合ロジットモデル(ミックスドロジットモデル(mixed logit model))の節を扱う.
ドキュメントの全てを扱うわけではない. また変数名やコードなど自分の理解しやすい形に変えているところもあるので注意されたい.
ググって探索していたら, 二週間前にこれがアップロードされていたことに気付いた. 英語が読めるならこっちを読んだ方が絶対にいい気がする.
加えてこれも発見した. 理論面はTrainのテキストに沿っていて, 補完的に読めてとてもよい.
{mlogit}パッケージ作った人の解説
Yves Croissant “Estimation of multinomial logit models in R : The mlogit Packages”
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.303.7401&rep=rep1&type=pdf
【学習内容】
3. 混合ロジットモデル
ライブラリ読み込み.
library(mlogit)
library(tidyverse)データの確認
選択実験データ(離散選択実験データ, 選好意識データ(Stated Preference Data), コンジョイント分析のデータ)と呼ばれるもの(だと思う). それぞれ異なる属性を持つ仮想の電力会社を4つ設定し, 実験参加者の個人がどの電力会社を選ぶかというデータ.
361人が実験に参加. 各個人は12個の想定パターンに回答. 12個全てに回答していない参加者もいる. そのため, 総回答数は4308となっている (全員が12個全てに回答していれば, 361*12=4332となる)
data("Electricity", package = "mlogit")
head(Electricity)## choice id pf1 pf2 pf3 pf4 cl1 cl2 cl3 cl4 loc1 loc2 loc3 loc4 wk1 wk2
## 1 4 1 7 9 0 0 5 1 0 5 0 1 0 0 1 0
## 2 3 1 7 9 0 0 0 5 1 5 0 0 1 0 1 1
## 3 4 1 9 7 0 0 5 1 0 0 0 0 0 1 0 1
## 4 4 1 0 9 7 0 1 1 0 5 0 0 1 0 1 0
## 5 1 1 0 9 0 7 0 1 0 5 1 0 0 0 0 1
## 6 4 1 0 9 0 7 0 0 1 5 0 0 1 0 0 0
## wk3 wk4 tod1 tod2 tod3 tod4 seas1 seas2 seas3 seas4
## 1 0 1 0 0 0 1 0 0 1 0
## 2 0 0 0 0 1 0 0 0 0 1
## 3 1 0 0 0 1 0 0 0 0 1
## 4 0 1 0 0 0 1 1 0 0 0
## 5 0 1 1 0 0 0 0 0 1 0
## 6 0 1 0 0 1 0 1 0 0 0変数
choice：個人の選択(1,2,3,4) id：個人のid
pfi：price fixed. 選択肢iの固定価格[cent/kWh] 価格は供給元や実験によって異なる、kWhあたりの表示セントでの固定価格. 価格はサプライヤーや実験によって異なる.
cli：contract length. 契約期間[年]. 他のサプライヤーに切り替えた場合, 違約金を支払うことになる. いつでも契約を解除できる契約も結ぶことができ, その場合は0と記録される. (携帯の2年縛りと同じようなものか.)
loci：local. サプライヤーが地元企業かどうか.
wki：well-known. サプライヤーが有名な企業家どうか.
todi：time of date. 午前8時から午後8時まで11cent/kWh, 午後8時から午前8時まで5cent/kWhというレートになるプラン.
seasi：夏季は10cent/kWh, 冬季は8cent/kWh, 春と秋は6cent/kWhというレートになるプラン. 価格は電気のみで, 地域の公益事業者による送配電は含まれていない.
wide形式から変換
Electr <- mlogit.data(Electricity, id = "id", choice = "choice",
varying = 3:26, shape = "wide", sep = "")
head(Electr)## choice id alt pf cl loc wk tod seas chid
## 1.1 FALSE 1 1 7 5 0 1 0 0 1
## 1.2 FALSE 1 2 9 1 1 0 0 0 1
## 1.3 FALSE 1 3 0 0 0 0 0 1 1
## 1.4 TRUE 1 4 0 5 0 1 1 0 1
## 2.1 FALSE 1 1 7 0 0 1 0 0 2
## 2.2 FALSE 1 2 9 5 0 1 0 0 2モデルのあてはめ
mlogit()で混合ロジットモデルするときの引数：
rpar : 変数名に分布を指定することで, その変数のパラメータに分布を当てはめる. “n”ならば正規分布, “l”ならば対数正規分布, “t”ならば切断正規分布, “u”ならば一様分布.
R : heterosc = TRUEの場合はガウス求積法に使用する関数の数. rparがNULLの場合は擬似乱数の数. (heteroscのデフォルトはNULL, Rのデフォルトは40)
print.level : printメッセージの詳細の指定. 0,1,2のいずれかをとる. 0ならば,最適化プロセスに関する情報は記載しない. 1ならば, 尤度の値,ステップおよび停止基準を記載する. 2ならば, パラメータと勾配ベクトルも記載する.
panel : rparがNULLではなく, かつデータが同じユニットの観測の繰り返しの場合に関係する. TRUEの場合、混合ロジットモデルはパネル手法を使用して推定される. (panelのデフォルトはNULL)
halton : rparがNULLでない場合のみ関係する引数. 擬似乱数の代わりにHalton sequence(ハルトン数列)が使用される. Halton数列は, 素数を選んで，(0,1)区間を分割することを繰り返して 分布範囲を均等に網羅し，かつ相互に相関などを持たない数列を準乱数として用いる方法.
使用する素数と削除する要素の数のリストを引数にとる. NAを指定したら, デフォルト値が使われれる.
Halton数列については東京大学工学部都市工学科都市生活学/ネットワーク行動学研究室の混合ロジットモデルに関する説明資料に依っている. http://bin.t.u-tokyo.ac.jp/kaken/pdf/mxl.pdf
詳しくはTrainのテキストのChap.9に説明がある. (https://eml.berkeley.edu/books/choice2.html) ここの研究室は東大の社基の羽藤研あたりが集まってるところで, 離散選択モデルに関する資料がたくさんアップロードされている. 土木・交通分野の個人の行動予測, 交通需要推計の研究が盛んみたい.
この研究分野の研究事例としては, 神奈川県のシーサイドラインが開通される前に需要予測を行っていたものなどがある.
http://library.jsce.or.jp/jsce/open/00039/1988/11-0455.pdf
using 100 draws, halton sequences and taking into account the panel data structure
混合モデル1
切片なしで, モデルの係数パラメータ6つに対して正規分布を当てはめるモデル. Halton数列を使用し, パネルデータ構造を考慮するケース.
推定には多少時間がかかる.
Elec.mxl <- mlogit(choice ~ pf + cl + loc + wk + tod + seas | 0, data = Electr,
rpar = c(pf="n", cl="n",loc="n",wk="n",tod="n",seas="n"),
R = 100, halton = NA, print.level = 0, panel = TRUE)
summary(Elec.mxl)##
## Call:
## mlogit(formula = choice ~ pf + cl + loc + wk + tod + seas | 0,
## data = Electr, rpar = c(pf = "n", cl = "n", loc = "n", wk = "n",
## tod = "n", seas = "n"), R = 100, halton = NA, panel = TRUE,
## print.level = 0)
##
## Frequencies of alternatives:
## 1 2 3 4
## 0.22702 0.26393 0.23816 0.27089
##
## bfgs method
## 21 iterations, 0h:0m:53s
## g'(-H)^-1g = 3.3E-07
## gradient close to zero
##
## Coefficients :
## Estimate Std. Error z-value Pr(>|z|)
## pf -0.973384 0.034324 -28.359 < 2.2e-16 ***
## cl -0.205557 0.013323 -15.428 < 2.2e-16 ***
## loc 2.075733 0.080430 25.808 < 2.2e-16 ***
## wk 1.475650 0.065168 22.644 < 2.2e-16 ***
## tod -9.052542 0.287219 -31.518 < 2.2e-16 ***
## seas -9.103772 0.289043 -31.496 < 2.2e-16 ***
## sd.pf 0.219945 0.010840 20.291 < 2.2e-16 ***
## sd.cl 0.378304 0.018489 20.461 < 2.2e-16 ***
## sd.loc 1.482980 0.081305 18.240 < 2.2e-16 ***
## sd.wk 1.000061 0.074182 13.481 < 2.2e-16 ***
## sd.tod 2.289489 0.110731 20.676 < 2.2e-16 ***
## sd.seas 1.180883 0.109007 10.833 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -3952.5
##
## random coefficients
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## pf -Inf -1.1217350 -0.9733844 -0.9733844 -0.82503376 Inf
## cl -Inf -0.4607190 -0.2055565 -0.2055565 0.04960589 Inf
## loc -Inf 1.0754783 2.0757333 2.0757333 3.07598832 Inf
## wk -Inf 0.8011189 1.4756497 1.4756497 2.15018054 Inf
## tod -Inf -10.5967791 -9.0525423 -9.0525423 -7.50830550 Inf
## seas -Inf -9.9002649 -9.1037717 -9.1037717 -8.30727842 Inf平均的な係数を持っている客が契約期間を１年伸ばす(減らす)のに追加的に欲しい(支払う)額を調べる.
coef(Elec.mxl)["cl"]/coef(Elec.mxl)["pf"]## cl
## 0.2111772契約期間の係数は平均-0.2, 標準偏差0.37の正規分布に従っていると推定されており, その場合, 係数が負となっている割合を計算すると
pnorm(-coef(Elec.mxl)["cl"]/coef(Elec.mxl)["sd.cl"])## cl
## 0.70656となり, 7割が契約期間が長くなることを嫌うと解釈できる.
同様に固定価格の係数についても見る.
pnorm(-coef(Elec.mxl)["pf"]/coef(Elec.mxl)["sd.pf"])## pf
## 0.9999952ほぼ100%, 係数が負となって整合的.
混合モデル3
よく知られた企業の方が好む人もいるだろうが, 好まない人も一定数いると考えがある場合, well-knownの係数の分布は正規分布ではなく一様分布であるというモデルを当てはめる
Elec.mxl3 <- update(Elec.mxl, rpar = c(cl="n", loc="n", wk="u", tod="n", seas="n"))
summary(Elec.mxl3)##
## Call:
## mlogit(formula = choice ~ pf + cl + loc + wk + tod + seas | 0,
## data = Electr, rpar = c(cl = "n", loc = "n", wk = "u", tod = "n",
## seas = "n"), R = 100, halton = NA, panel = TRUE, print.level = 0)
##
## Frequencies of alternatives:
## 1 2 3 4
## 0.22702 0.26393 0.23816 0.27089
##
## bfgs method
## 22 iterations, 0h:0m:55s
## g'(-H)^-1g = 3.71E-07
## gradient close to zero
##
## Coefficients :
## Estimate Std. Error z-value Pr(>|z|)
## pf -0.882232 0.032818 -26.883 < 2.2e-16 ***
## cl -0.217127 0.013776 -15.761 < 2.2e-16 ***
## loc 2.099319 0.081056 25.900 < 2.2e-16 ***
## wk 1.509410 0.065496 23.046 < 2.2e-16 ***
## tod -8.607002 0.282983 -30.415 < 2.2e-16 ***
## seas -8.602408 0.280671 -30.649 < 2.2e-16 ***
## sd.cl 0.381070 0.018150 20.996 < 2.2e-16 ***
## sd.loc 1.593850 0.087802 18.153 < 2.2e-16 ***
## sd.wk 1.786377 0.125764 14.204 < 2.2e-16 ***
## sd.tod 2.719078 0.119357 22.781 < 2.2e-16 ***
## sd.seas 1.945377 0.103686 18.762 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -3956.7
##
## random coefficients
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## cl -Inf -0.4741552 -0.2171273 -0.2171273 0.03990058 Inf
## loc -Inf 1.0242834 2.0993191 2.0993191 3.17435473 Inf
## wk -0.2769665 0.6162218 1.5094101 1.5094101 2.40259840 3.295787
## tod -Inf -10.4409924 -8.6070022 -8.6070022 -6.77301190 Inf
## seas -Inf -9.9145449 -8.6024084 -8.6024084 -7.29027193 Infwell-knownの係数の分布情報を確認する.
rpar(Elec.mxl3,"wk")## uniform distribution with parameters 1.509 (center) and 1.786 (span)plotで分布の形状も確認できる. 係数が負になるところは色が塗られて表題にその割合が書かれるようだ. 先のサマリーで確認できるが, -0.27から3.29の一様分布.
plot(rpar(Elec.mxl3,"wk"))
【学習予定】
プロビットモデルもある(https://cran.r-project.org/web/packages/mlogit/vignettes/e4mprobit.html)が, ロジットモデルに関してとりあえず{mlogit}の使い方が分かったので, いったん離散選択モデルの扱いの学習はここで止める. プロビットの学習やネスティッドロジットや混合ロジットのシミュレーションデータを作成と推定をしてみたりしたいところだが, 他に時間と労力を割かなければなのでまたいつか.
今回やっていて数値計算や最適化周りの知識の足りなさを感じたのでそのへんの補強も必要そう. このへんはTrainのテキストの後半パートが役立ちそう.
経済学の産業組織の分野のBLPアプローチは混合ロジットモデルを, 内生性を考慮して何かしているものっぽい. 産業組織は分からないので正しいか分からないがここにメモしておく.
