Rの離散選択モデルパッケージapolloを使ってみる①MNL推定
はじめに
離散選択モデルをRで推定する際にmlogitやbinが公開しているコードがよく使われていると思う。私の指導教員はapolloというパッケージを勧めてくれたのだが、これが全て英語かつ機能が豊富で、簡単なMNLを推定するのにも一苦労したので、自分のメモ代わりもかねて誰かの手助けとなればと思い、使い方の説明を残しておこうと思う。
apolloとは
Rのパッケージの1つで離散選択モデルのさまざまなモデルを推定することができます。マニュアルとサンプルコードが公開されているので、英語が得意な人はそっちを見るほうが早いと思います。
www.apollochoicemodelling.com
実践編
今回はマニュアルでも紹介されているapollo_example_3のコードを説明していこうと思います。example_1で最も基本的なMNL推定を行い、その推定結果を途中で利用しています。
起動と初期化
とくにここはいじる必要はないです。起動して諸々初期化しましょう。
# ################################################################# # #### LOAD LIBRARY AND DEFINE CORE SETTINGS #### # ################################################################# # ### Clear memory rm(list = ls()) ### Load Apollo library library(apollo) ### Initialise code apollo_initialise()
コア項目の設定
モデルの名前、説明、個人のIDがどこに入ってるかを指定します。そのほかにも色々設定できるみたいですが使ったことないのでそんなに重要ではないでしょう（多分）。
### Set core controls apollo_control = list( modelName ="Apollo_example_3", modelDescr ="MNL model with socio-demographics on mode choice SP data", indivID ="ID" )
データ読み込み
データベースを読み込みます。このサンプルコードではSPのデータだけサブセットして平均年収を計算してますがそういうのが必要なければread.csvだけで十分。
# ################################################################# # #### LOAD DATA AND APPLY ANY TRANSFORMATIONS #### # ################################################################# # database = read.csv("apollo_modeChoiceData.csv",header=TRUE) ### Use only SP data database = subset(database,database$SP==1) ### Create new variable with average income database$mean_income = mean(database$income)
モデル定義
まずはパラメータの定義から。apollo_betaでパラメータの名前と初期値を設定します。ここでは全て0になってますがそれ以外の値でも大丈夫です。
そしてapollo_fixedで初期値に固定するパラメータを指定します。0で固定するならパラメータ自体を削除しても大丈夫みたいです。
apollo_readBetaで事前に推定した結果を初期値として使うこともできます。example_1のサンプルコードで推定した結果を参照しているので同じディレクトリにexample_1の結果をファイル出力しておきましょう。初期値を既存の結果から参照する必要は必ずしもないですが、計算量が多くなって事前にパラメータを近い値に設定しておけば計算時間を短縮できるみたいです。
# ################################################################# # #### DEFINE MODEL PARAMETERS #### # ################################################################# # ### Vector of parameters, including any that are kept fixed in estimation apollo_beta=c(asc_car = 0, asc_bus = 0, asc_air = 0, asc_rail = 0, asc_bus_shift_female = 0, asc_air_shift_female = 0, asc_rail_shift_female = 0, b_tt_car = 0, b_tt_bus = 0, b_tt_air = 0, b_tt_rail = 0, b_tt_shift_business = 0, b_access = 0, b_cost = 0, b_cost_shift_business = 0, cost_income_elast = 0, b_no_frills = 0, b_wifi = 0, b_food = 0) ### Vector with names (in quotes) of parameters to be kept fixed at their starting value in apollo_beta, use apollo_beta_fixed = c() if none apollo_fixed = c("asc_car","b_no_frills") ### Read in starting values for at least some parameters from existing model output file apollo_beta = apollo_readBeta(apollo_beta, apollo_fixed, "Apollo_example_1", overwriteFixed=FALSE)
入力確認
何をやってるのかは自分でもよくわかってませんが、省略できないので実行しましょう。
# ################################################################# # #### GROUP AND VALIDATE INPUTS #### # ################################################################# # apollo_inputs = apollo_validateInputs()
特に問題なければ次のような結果が出力されます。
Several observations per individual detected based on the value of ID. Setting panelData in apollo_control set to TRUE. All checks on apollo_control completed. All checks on database completed.
効用関数の設定
一番重要（？）な効用関数の設定の部分です。
3段目の定数項と変数の計算式はこのモデルでより複雑なモデルを推定しているだけで、これをやる必要はないです。
4段目のV[ ['car'] ]で各選択肢の効用関数を定義します。アタッチしてるのでdatabase$を付ける必要はないです。選択肢の定義は次の段で決めるmnl_settingsと同じである必要があります。
5段目のmnl_settingsのalternativeで選択肢を、availで選択可能かを、choiceVarで選択結果がdatabaseのどこに格納されているかを指定します。Vはそのままで大丈夫。
7段目のapollo_panelProdはパネルデータの場合のみ使います。逆にパネルデータ以外の時にこれを実行するとエラーを吐きます。（ちなみにパネルデータを使う時は同じIDのデータは連続で並べないとうまく読み取ってくれません）
# ################################################################# # #### DEFINE MODEL AND LIKELIHOOD FUNCTION #### # ################################################################# # apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){ ### Attach inputs and detach after function exit apollo_attach(apollo_beta, apollo_inputs) on.exit(apollo_detach(apollo_beta, apollo_inputs)) ### Create list of probabilities P P = list() ### Create alternative specific constants and coefficients using interactions with socio-demographics asc_bus_value = asc_bus + asc_bus_shift_female * female asc_air_value = asc_air + asc_air_shift_female * female asc_rail_value = asc_rail + asc_rail_shift_female * female b_tt_car_value = b_tt_car + b_tt_shift_business * business b_tt_bus_value = b_tt_bus + b_tt_shift_business * business b_tt_air_value = b_tt_air + b_tt_shift_business * business b_tt_rail_value = b_tt_rail + b_tt_shift_business * business b_cost_value = ( b_cost + b_cost_shift_business * business ) * ( income / mean_income ) ^ cost_income_elast ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant V = list() V[['car']] = asc_car + b_tt_car_value * time_car + b_cost_value * cost_car V[['bus']] = asc_bus_value + b_tt_bus_value * time_bus + b_access * access_bus + b_cost_value * cost_bus V[['air']] = asc_air_value + b_tt_air_value * time_air + b_access * access_air + b_cost_value * cost_air + b_no_frills * ( service_air == 1 ) + b_wifi * ( service_air == 2 ) + b_food * ( service_air == 3 ) V[['rail']] = asc_rail_value + b_tt_rail_value * time_rail + b_access * access_rail + b_cost_value * cost_rail + b_no_frills * ( service_rail == 1 ) + b_wifi * ( service_rail == 2 ) + b_food * ( service_rail == 3 ) ### Define settings for MNL model component mnl_settings = list( alternatives = c(car=1, bus=2, air=3, rail=4), avail = list(car=av_car, bus=av_bus, air=av_air, rail=av_rail), choiceVar = choice, V = V ) ### Compute probabilities using MNL model P[["model"]] = apollo_mnl(mnl_settings, functionality) ### Take product across observation for same individual P = apollo_panelProd(P, apollo_inputs, functionality) ### Prepare and return outputs of function P = apollo_prepareProb(P, apollo_inputs, functionality) return(P) }
モデルの推定と結果の出力
あとはモデルの推定と結果を出力するだけです。apoll_estimateがモデルの推定、apollo_modelOutputでコンソールに結果を出力、apollo_modelOutputで結果をファイル出力できる。
ただし、apollo_estimamateを代入した変数の名前ではなく、apollo_controlのmodelNameで指定した名前で出力され同じ名前だと上書きされるので注意（oldとして残るが）。
# ################################################################# # #### MODEL ESTIMATION #### # ################################################################# # model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs) # ################################################################# # #### MODEL OUTPUTS #### # ################################################################# # # ----------------------------------------------------------------- # #---- FORMATTED OUTPUT (TO SCREEN) ---- # ----------------------------------------------------------------- # apollo_modelOutput(model) # ----------------------------------------------------------------- # #---- FORMATTED OUTPUT (TO FILE, using model name) ---- # ----------------------------------------------------------------- # apollo_saveOutput(model)
推定結果の出力
apollo_modelOutputを実行すると次のように結果が出力されます。
#基本情報の表示
Model name : Apollo_example_3
Model description : MNL model with socio-demographics on mode choice SP data
Model run at : 2021-08-02 18:25:18
Estimation method : bfgs
Model diagnosis : successful convergence
Number of individuals : 500
Number of rows in database : 7000
Number of modelled outcomes : 7000
Number of cores used : 1
Model without mixing
#適合度指標
LL(start) : -6359.334
LL(0) : -8196.021
LL(final) : -4830.945
Rho-square (0) : 0.4106
Adj.Rho-square (0) : 0.4085
AIC : 9695.89
BIC : 9812.4
#推定に関する情報
Estimated parameters : 17
Time taken (hh:mm:ss) : 00:00:14.06
pre-estimation : 00:00:3.36
estimation : 00:00:3.87
post-estimation : 00:00:6.83
Iterations : 27
Min abs eigenvalue of Hessian : 2.818775
#推定結果（順に推定値、標準誤差、t値、ロバスト標準誤差、ロバストt値）
Estimates:
Estimate s.e. t.rat.(0) Rob.s.e. Rob.t.rat.(0)
asc_car 0.000000 NA NA NA NA
asc_bus 0.286717 0.582949 0.4918 0.549508 0.5218
asc_air -0.903270 0.373537 -2.4182 0.361569 -2.4982
asc_rail -2.092578 0.353359 -5.9220 0.351090 -5.9602
asc_bus_shift_female 0.340222 0.132783 2.5622 0.145285 2.3418
asc_air_shift_female 0.268170 0.091506 2.9306 0.095283 2.8145
asc_rail_shift_female 0.189612 0.073759 2.5707 0.078168 2.4257
b_tt_car -0.013107 7.3440e-04 -17.8471 7.6775e-04 -17.0718
b_tt_bus -0.021266 0.001598 -13.3111 0.001518 -14.0096
b_tt_air -0.016578 0.002774 -5.9771 0.002671 -6.2074
b_tt_rail -0.007051 0.001811 -3.8930 0.001765 -3.9960
b_tt_shift_business -0.006234 6.0073e-04 -10.3773 5.9078e-04 -10.5522
b_access -0.021152 0.002865 -7.3841 0.002706 -7.8155
b_cost -0.076187 0.002097 -36.3300 0.002095 -36.3619
b_cost_shift_business 0.033379 0.002739 12.1844 0.002572 12.9783
cost_income_elast -0.613842 0.030094 -20.3976 0.030606 -20.0566
b_no_frills 0.000000 NA NA NA NA
b_wifi 1.026703 0.056142 18.2876 0.057815 17.7586
b_food 0.422063 0.055027 7.6702 0.056464 7.4749
自分でもまだよく分かってないところが多いので間違っているところがあればコメントいただければ助かります。。
