はじめに

前回の記事では，交通手段選択を扱う多項ロジットモデルを推定しました．
データには全国旅客純流動調査のオープンデータ，推定には，PythonのBiogemeライブラリを用いました．

今回の記事では，推定にRのmlogitパッケージを用いて，同じモデルを推定してみようと思います．
推定用のデータは，前回の記事と同じで，githubで公開しています．

参考サイト

以下のサイトを参考にしました．

mlogitによる多項ロジットモデルの推定

パッケージのインポート

まずは，使用するパッケージのインストールとインポートを行います．

estimate_MNL_mlogit.R

install.packages("mlogit")

library(mlogit)
library(dfidx)
library(readr)
library(dplyr)

データの読み込み・加工

準備したデータ(詳細は前回記事参照)を読み込んで以下の加工を行います:

着目する交通手段だけに限定したデータにする．
文字を含む列を削除する．
NaNを含む行は削除する．
num列に記載されている，変数の組み合わせが同一であるサンプル数の分だけ，行を複製する(非集計データにする)．
mlogitでデータを読み込む際に，LOS変数は選択肢識別子が末尾に書いてある必要があるので，列名を変更します．

estimate_MNL_mlogit.R

# データの読み込み
df <- readr::read_csv("../data/data.csv", locale = locale(encoding = "Shift_JIS"))

# 着目する交通手段に限定（鉄道・バス・乗用車等）
df <- df[(df$mode_name == "幹線バス" | df$mode_name == "鉄道" | df$mode_name == "乗用車等"),]

# 不要な列を削除（文字列・使用しないLOS変数）
df <- df[, !(colnames(df) %in% c('O_name', 'D_name', 'purpose_name', 'sex_name', 'ship_time', 'ship_cost', 'air_time', 'air_cost'))]

# 欠損値を含む行を削除
df <- na.omit(df)

# トリップ数に応じてデータを展開（各行を'num'回繰り返す）
df <- df[rep(seq_len(nrow(df)), df$num), , drop = FALSE]

# 末尾に選択肢識別子がくるように列名を変更
df <- 
  df |> 
  rename(
    time_car = car_time, 
    time_bus = bus_time, 
    time_rail = rail_time, 
    cost_car = car_cost, 
    cost_bus = bus_cost, 
    cost_rail = rail_cost
  )

変数のスケーリング

スケーリングした変数を新たな変数として定義します．

estimate_MNL_mlogit.R

df$stime_car <- df$time_car / 60
df$stime_bus <- df$time_bus / 60
df$stime_rail <- df$time_rail / 60
df$scost_car <- df$cost_car / 10000
df$scost_bus <- df$cost_bus / 10000
df$scost_rail <- df$cost_rail / 10000

mlogitデータの生成

mlogitは列名の後半に選択肢名を付けておくと，LOS変数¹を自動で識別してくれます．
そこで，選択肢名を定義し，それを各サンプルの選択結果と対応付けます．
各選択肢の選択肢番号(5: 自動車，4: 幹線バス，2: 鉄道)と識別用の選択肢名を紐づけるcode2altリストを定義します．
code2altリストを介して，choice_alt列に各サンプルの選択結果の選択肢名を格納します．
mlogitのデータ作成時にこの列を参照することで，選択肢名を認識させます．

estimate_MNL_mlogit.R

code2alt <- c(`5`="car", `4`="bus", `2`="rail")
df$choice_alt <- unname(code2alt[as.character(df$mode_code)])
str(df)

mlogit用のデータを作成します．
mlogit.dataの

第1引数には，変換元のデータを指定します．
shapeには，データの形式の種類("wide"or"long")を指定します．それぞれのデータ形式については，参考サイトをご覧ください．
choiceには，先程作成したchoice_alt列を指定して，選択結果を選択肢名で識別させます．
varyingには，LOS変数が格納されている列番号を指定します．
sepには，LOS変数の変数名で，選択肢名とそれ以外を分けるのに用いている記号(今回は_)を指定します．これにより，例えば，stime_carという変数は，選択肢名carのstime変数であると機械的に識別してくれます．

estimate_MNL_mlogit.R

df_mlogit <- mlogit.data(
  df, 
  shape = "wide", 
  choice ="choice_alt", 
  varying = c(9:20), 
  sep="_"
)

str(df_mlogit)

MNLの推定

多項ロジットモデルを推定して，結果を表示します．.
mlogitの，

formulaで効用関数を記述します．被説明変数であるchoice_altを~の左側に書き，スケーリングした変数stimeとscostを続けます．|の右側に1を書くことで，定数項を追加できます．
dataには推定に用いるmlogit用のデータを指定します．
reflevelには，定数項のうち，基準にする選択肢を記述します．前回の記事と同じく，busを基準にします．

estimate_MNL_mlogit.R

mode_MNL <- mlogit(
  formula  = choice_alt ~ stime + scost | 1,  # | 1 で ASCs を入れる
  data     = df_mlogit,
  reflevel = "bus"
)
summary(mode_MNL)

尤度比の計算

ヌル尤度を計算して，尤度比指標を取得します．
まずは，推定したモデルの最終尤度を取得して，llに格納します．
続いて，num_sampleにサンプル数を格納し，ヌル尤度(すべてのパラメータが0であるときの尤度)を計算して，ll0に格納します．
定義に従い，尤度比指標を計算します．

estimate_MNL_mlogit.R

ll <- summary(mode_MNL)$logLik

num_samples <- nrow(df)
ll0 <- num_samples * log(1 / 3)
as.numeric( 1 - (ll / ll0) )

的中率の計算

的中率を計算します．
まずは，パラメータの推定結果をもとに，各サンプルの各選択肢の選択確率を計算します．

estimate_MNL_mlogit.R

cf <- coef(mode_MNL)
asc_car  <- unname(cf["(Intercept):car"])
asc_rail <- unname(cf["(Intercept):rail"])
b_time   <- unname(cf["stime"])
b_cost   <- unname(cf["scost"])

df$pred_V_car  <- asc_car  + b_time * df$stime_car  + b_cost * df$scost_car
df$pred_V_bus  <-            b_time * df$stime_bus  + b_cost * df$scost_bus 
df$pred_V_rail <- asc_rail + b_time * df$stime_rail + b_cost * df$scost_rail

deno <- exp(df$pred_V_car) + exp(df$pred_V_bus) + exp(df$pred_V_rail)

df$pred_P_car <- exp(df$pred_V_car) / deno
df$pred_P_bus <- exp(df$pred_V_bus) / deno
df$pred_P_rail <- exp(df$pred_V_rail) / deno

続いて，各選択肢名と予測される選択確率を結合した行列を作成します．

estimate_MNL_mlogit.R

pred_P <- cbind(
  bus  = df$pred_P_bus,
  car  = df$pred_P_car,
  rail = df$pred_P_rail
)

作成した行列の各行で最大の選択確率を取った選択肢の選択肢名(列名)を取得し，各サンプルの予測される選択肢として，pred_altに格納します．

estimate_MNL_mlogit.R

df$pred_alt <- colnames(pred_P)[max.col(pred_P, ties.method = "first")]

的中サンプル数と不的中サンプル数の和が全サンプル数と一致することを確認して，的中率を計算します．

estimate_MNL_mlogit.R

stopifnot(
  sum(df$pred_alt == df$choice_alt, na.rm = TRUE) +
  sum(df$pred_alt != df$choice_alt, na.rm = TRUE) ==
  num_samples
)

accuracy <- sum(df$pred_alt == df$choice_alt, na.rm = TRUE) / num_samples * 100
accuracy

おわりに

今回は，Rのパッケージのmlogitを用いて，多項ロジットモデルの推定を行いました．
以上のコードで推定された結果は，前回の記事の推定結果と同じになっていると思います．
多項ロジットモデルの推定コードの比較として参考にしていただければ幸いです．

データの出典

出典: 「第6回（2015年度）全国幹線旅客純流動調査」（国土交通省）を加工して作成

選択肢毎に定義される，選択肢間で値が異なる変数のことです． ↩