BiogemeでMNL推定～全国旅客純流動調査を用いて～

交通計画の分野において，離散選択モデルは交通手段選択行動や経路選択行動の表現など様々な用途で用いられます．
そんな離散選択モデルを組むためのPythonライブラリの1つにBiogemeがあります．
本記事では，そんなBiogemeを用いた交通手段選択行動を扱う多項ロジットモデルの推定方法をまとめます．
本記事の内容は，以下のサイトのGet started > Read the primerの資料を参考に作成しました．

データの準備

使用するデータは第6回全国旅客純流動調査で，e-Statで公開されています．

今回はモデル推定に際して，都道府県間流動表の出発地から目的地のデータを使用したいと思います．
公開ページから，2015年都道府県間（平日;代表機関;出発地-目的地）のcsvをダウンロードします．

また，OD別交通サービス水準の公開ページにおいて，OD別所要時間のexcelとOD別費用のexcelをダウンロードします．

ダウンロードしたデータを加工して，次の列を含むdata.csvを用意しました．
このデータはgithubで公開しています．

列名	説明
O_name / O_code	出発地の都道府県名とそのコード
D_name / D_code	目的地の都道府県名とそのコード
mode_name / mode_code	実際に利用した交通手段とそのコード(1: 航空，2: 鉄道，3: 幹線旅客船，4: 幹線バス，5: 自動車)
num	変数の組み合わせが同一であるサンプル数
car_time	出発地から目的地に自動車で向かう場合の所要時間[分]
rail_time	出発地から目的地に鉄道で向かう場合の所要時間[分]
bus_time	出発地から目的地に幹線バスで向かう場合の所要時間[分]
car_cost	出発地から目的地に自動車で向かう場合の所要費用[円]
rail_cost	出発地から目的地に鉄道で向かう場合の所要費用[円]
bus_cost	出発地から目的地に幹線バスで向かう場合の所要費用[円]

他の変数

配布データには，他にも以下の列を用意しました．

列名	説明
purpose_name / purpose_code	移動目的とそのコード(1: 仕事，2: 観光，3: 私用・帰省，4: その他)
sex_name / sex_code	性別とそのコード(1: Male，2: Female)
age_code	年齢コード(10: 10歳代，20: 20歳代，30: 30歳代，40: 40歳代，50: 50歳代，60: 60歳代，70: 70歳以上)
ship_time	出発地から目的地に幹線旅客船で向かう場合の所要時間[分]
air_time	出発地から目的地に航空で向かう場合の所要時間[分]
ship_cost	出発地から目的地に幹線旅客船で向かう場合の所要費用[円]
air_cost	出発地から目的地に航空で向かう場合の所要費用[円]

多項ロジットモデル推定

準備したデータを用いてBiogemeで交通手段選択モデルを推定します．
今回は，鉄道・幹線バス・自動車の3つの選択肢からの選択をモデル化します．

ライブラリのインポート

import numpy as np
import pandas as pd

import biogeme.database as db
from biogeme.expressions import Variable

import biogeme.biogeme as bio
from biogeme import models
from biogeme.expressions import Beta

データの読み込み・加工

準備したデータをpandasで読み込んで以下の加工を行います:

1. 着目する交通手段だけに限定したデータにする．
2. Biogemeでは，ヘッダー行以外で文字(string)は読み込めないので，文字を含む列を削除する．
3. Biogemeでは，NaNも読み込めないので，NaNを含む行は削除する．
4. num列に記載されている，変数の組み合わせが同一であるサンプル数の分だけ，行を複製する(非集計データにする)．
5. Biogeme形式のデータに変換して読み込む．

df = pd.read_csv("../data/data.csv", encoding='shift_jis')

# 着目する交通手段に限定（鉄道・バス・乗用車等）
df = df[(df['mode_name'] == "幹線バス") | (df['mode_name'] == "鉄道") | (df['mode_name'] == "乗用車等")].reset_index(drop=True).copy()

# 不要な列を削除（文字列・使用しないLOS変数）
df = df.drop(columns=['O_name', 'D_name', 'purpose_name', 'mode_name', 'sex_name', 'ship_time', 'ship_cost', 'air_time', 'air_cost']).copy()

# 欠損値を含む行を削除
df = df.dropna().reset_index(drop=True).copy()

# トリップ数に応じてデータを展開（各行を'num'回繰り返す）
df = df.loc[df.index.repeat(df['num'])].reset_index(drop=True).copy() 

# Biogeme形式のデータに変換
database = db.Database ('junryudou', df)

変数の定義

推定に用いる変数を定義します．
データベース内の各列をモデルの変数とするために，以下のように記述します．

mode = Variable('mode_code')  # 実際に選択した交通手段のコード
car_time = Variable('car_time')  # 自動車の所要時間
bus_time = Variable('bus_time')  # 幹線バスの所要時間
rail_time = Variable('rail_time')  # 鉄道の所要時間
car_cost = Variable('car_cost')  # 自動車の所要費用
bus_cost = Variable('bus_cost')  # 幹線バスの所要費用
rail_cost = Variable('rail_cost')  # 鉄道の所要費用

変数のスケーリングを行う場合は，つぎのように，変数を再定義します．

car_time_scaled = car_time / 60
car_cost_scaled = car_cost / 10000
rail_time_scaled = rail_time / 60
rail_cost_scaled = rail_cost / 10000
bus_time_scaled = bus_time / 60
bus_cost_scaled = bus_cost / 10000

パラメータの定義

関数Betaを用いて，パラメータを定義します．
関数Betaでは，パラメータ名，初期値，下限値，上限値，推定の要否を指定します．

# parameter = Beta('name', value, lowerBound, upperBound, status)
# status: 0=推定対象, 1=固定値

ASC_car = Beta('ASC_car', 0, None, None, 0)    # 乗用車の選択肢固有係数
ASC_rail = Beta('ASC_rail', 0, None, None, 0)  # 鉄道の選択肢固有係数
ASC_bus = Beta('ASC_bus', 0, None, None, 1)    # バスの選択肢固有定数（基準：固定）
B_time = Beta('B_time', 0, None, None, 0)      # 所要時間の選択肢共通係数
B_cost = Beta('B_cost', 0, None, None, 0)      # 費用の選択肢共通係数

効用関数の特定

以上で定義した変数とパラメータを使って，効用関数を記述します．
今回は次の式のように各選択肢の効用が各選択肢の所要時間と所要費用で表されるものとします¹．

また，パラメータは選択肢共通のものとします．

    \begin{array}
        & V_{car} = ASC_{car} + \beta_{time} T_{car} + \beta_{cost} C_{car} \\
        V_{rail} = ASC_{rail} + \beta_{time} T_{rail} + \beta_{cost} C_{rail} \\
        V_{bus} = \beta_{time} T_{bus} + \beta_{cost} C_{bus} 
    \end{array}

V_car = ASC_car + B_time * car_time_scaled + B_cost * car_cost_scaled
V_rail = ASC_rail + B_time * rail_time_scaled + B_cost * rail_cost_scaled
V_bus = ASC_bus + B_time * bus_time_scaled + B_cost * bus_cost_scaled

選択結果との対応付け・選択可能性の設定

各効用関数をデータ上の選択肢番号と対応付けます．
今回は，mode列に自動車が5，鉄道が2，幹線バスが4とされているので，各々の効用関数と対応付けます．
これには，辞書型の変数Vを用います．

# mode_codeと効用関数の対応付け
V = {5: V_car, 2: V_rail, 4: V_bus}

また，選択可能性についても各選択肢番号に対して個別に設定します．
選択可能性とは，例えば，運転免許を持っていない人は，自動車が選択肢になり得ないといった，そもそも選択肢になり得ない状況を表現するためのものです．
今回はすべてのサンプルがすべての選択肢を選択できるものとします²．

# 選択可能性(1: 選択可, 0: 選択不可)
av = {5: 1, 2: 1, 4: 1}

モデルの定義と推定

推定のためにモデルを定義します．
models.loglogit(V, av, mode)では，ロジットモデルの選択確率の対数を返します．
引数には，上記で定義した効用関数の辞書，選択可能性の辞書，選択結果の変数を渡します．

# 選択確率の対数を計算
logprob = models.loglogit(V, av, mode)

続いて，BIOGEMEオブジェクトを準備します．
biogeme公式の説明では，この変数にbiogemeという名前を設定するのは避け，たとえば，the_biogemeのような名前とすることが推奨されています．
databaseは，「データの読み込み・加工」で定義した推定用のデータを格納した変数です．

# BIOGEMEオブジェクトの準備
the_biogeme = bio.BIOGEME(database, logprob)

次のようにモデル名を定義します．
ここで定義した名前が，推定結果を記載したレポートファイルのファイル名となります．

# モデル名の定義
the_biogeme.model_name = 'mode_MNL'

各係数が0である場合の尤度はモデルの適合度評価に用いる重要な指標です．
次のような記述をすることで計算することができます．

# ヌルモデルの対数尤度を計算
the_biogeme.calculate_null_loglikelihood(av)

次のように記述することで，モデルの推定を実行します．

# モデル推定の実行
MNL_results = the_biogeme.estimate()

推定結果の表示

モデルの推定結果を表示します．
こちらのコードで表示される推定結果のほかに，ディレクトリ内に推定結果をまとめたレポートファイルが出力されます．

print(MNL_results.short_summary())

t値の計算

推定結果のレポートファイルで，推定されたパラメータのの$t$値は，ロバスト推定値のみ記載されています³．

そこで，非ロバストの$t$値も出力されるようにします．
まずは，パラメータの推定値とその分散共分散行列を取得します．

# ライブラリのインポート
from biogeme.results_processing.estimation_results import EstimationResults
from biogeme.results_processing.variance_covariance import EstimateVarianceCovariance

# MNL_results が RawEstimationResults の場合に EstimationResults でラップ
res = MNL_results if isinstance(MNL_results, EstimationResults) else EstimationResults(MNL_results)

# 推定値（β）
betas: dict[str, float] = res.get_beta_values()

# 通常（Rao–Cramer）の分散共分散行列を取得
vc = res.get_variance_covariance_matrix(
    variance_covariance_type=EstimateVarianceCovariance.RAO_CRAMER # ROBUSTを指定すると，ロバスト値を得る
)

非ロバストの$t$値と$p$値を計算し，表敬式で出力します．

# 標準誤差と t 値を計算
names = list(betas.keys())
values = np.array([betas[n] for n in names], dtype=float)
se = np.sqrt(np.diag(np.asarray(vc, dtype=float)))
t_vals = values / se

# p 値
from math import erf, sqrt
def p_from_t(t):
    Phi = 0.5 * (1.0 + erf(abs(t) / sqrt(2.0)))  # 一部erfを用いて標準正規分布の累積分布関数を設定
    return 2.0 * (1.0 - Phi)  # 両側検定

p_vals = [p_from_t(t) for t in t_vals]

# 表として出力
names = list(betas.keys())  # dict型の変数betasに対して，keysメソッドを実行
df_nonrobust = pd.DataFrame(
    {"Value": values,
     "Std err. (non-robust)": se,
     "t-test (non-robust)": t_vals,
     "p-value (non-robust)": p_vals},
    index=names
)
print(df_nonrobust)

的中率の計算

各パラメータの推定値を取得し，変数に代入します．

# パラメータの推定値の取得と変数への代入
res_ASC_car = betas['ASC_car']
res_ASC_rail = betas['ASC_rail']
res_B_time = betas['B_time']
res_B_cost = betas['B_cost']

パラメータの推定値を用いて，効用を計算します．

# 効用の計算
df['pred_V_car'] = res_ASC_car + res_B_time * df['car_time'] / 60 + res_B_cost * df['car_cost'] / 10000
df['pred_V_rail'] = res_ASC_rail + res_B_time * df['rail_time'] / 60 + res_B_cost * df['rail_cost'] / 10000
df['pred_V_bus'] = res_B_time * df['bus_time'] / 60 + res_B_cost * df['bus_cost'] / 10000

さらに，ロジットモデルの式を用いて，各選択肢の選択確率を計算します．

# 選択確率の分母
deno = np.exp(df['pred_V_car']) + np.exp(df['pred_V_rail']) + np.exp(df['pred_V_bus'])

# 各選択肢の選択確率
df['pred_P_car'] = np.exp(df['pred_V_car']) / deno
df['pred_P_rail'] = np.exp(df['pred_V_rail']) / deno
df['pred_P_bus'] = np.exp(df['pred_V_bus']) / deno

計算された各選択肢の選択確率のうち，最大であった選択肢が選ばれると予測します．

# 選択確率が最大であった選択肢に選択肢番号を対応付け
cols = ['pred_P_car', 'pred_P_rail', 'pred_P_bus']
code_map = {'pred_P_car': 5, 'pred_P_rail': 2, 'pred_P_bus': 4}
df['pred_mode'] = df[cols].idxmax(axis=1).map(code_map).astype('Int64')

実際に選択された選択肢と選択されると予測された選択肢が一致している場合には，的中，そうでない場合は，不的中であるとします．
さらにその判定結果をもとに的中率を計算します．

# 的中と不的中の割り当て
df['judge'] = np.where(df['pred_mode'].astype('Int64').eq(df['mode_code'].astype('Int64')),
                       'correct', 'wrong')
                       
# 的中となったサンプル数を数えて，的中率を計算
print( ((df["judge"] == "correct").sum() + (df["judge"] == "wrong").sum()) == len(df))
acc = (df["judge"] == "correct").sum() / ((df["judge"] == "correct").sum() + (df["judge"] == "wrong").sum()) * 100
print(f"的中率: {acc}%")

以上で推定された結果は下表です．
尤度比から一定程度の説明力を有しているモデルであると考えられ，的中率も十分に高いです．
パラメータの符号条件も自然に解釈できる結果となりました．

パラメータ	推定値	$t$値
自動車の定数項 $ASC_{car}$	2.47	472.44 **
鉄道の定数項$ASC_{rail}$	0.67	98.37 **
所要時間$\beta_{time}$	-0.93	-625.90 **
所要費用$\beta_{cost}$	-1.47	-122.73 **
サンプル数	2,216,684
尤度比$\rho^2$	0.54
的中率[%]	82.67
**: 1%有意

おわりに

今回は，PythonライブラリのBiogemeを用いて，多項ロジットモデルの推定を行いました．
多項ロジットモデルの推定以外にも用いることができるライブラリなので，幅広く活用範囲を探したいです．

データの出典

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

1. 今回の変数の組み合わせは多重共線性を生んでいるかもしれませんが，詳しくは触れません．実際に推定するときは確認してみると良いかもしれません．例えば，こちらの記事が分かりやすいと思います． ↩

2. 別途各サンプルの各選択肢の選択可能性を記入した列があれば，その列を変数として定義し，さらにその変数をここで指定することもできます． ↩

3. ロバスト推定については，例えば，こちらの記事やこちらの記事をご覧ください． ↩