Rで 選択モデリング(混合ロジットモデル) を試みます。
選択モデリングのシナリオは以下の通りです。
ストーリー・シナリオ：テクノロジー・シティの通勤改革
登場人物と背景
- 舞台: テクノロジー企業がひしめく大都市「テクノロジー・シティ」。
- 住民: この街で働く多様な会社員たち。彼らは毎日の通勤で3つの選択肢を持っています。
- 車 (Car): 自由で快適だが、渋滞と高い維持費（ガソリン、駐車場）が悩み。
- 電車 (Train): 定時運行で読書などもできるが、駅までのアクセスや混雑がストレス。
- 自転車 (Bike): 健康的でエコだが、天候に左右され、坂道はきつい。
- 都市の課題: 人口急増により、朝夕の交通渋滞は限界に達し、環境汚染も深刻化しています。新市長は「スマートでグリーンな通勤都市」への変革を公約に掲げました。
都市計画チームの挑戦
市長の特命を受けた都市計画チームは、市民の通勤行動を科学的に分析し、効果的な政策を打ち出すため、選択モデリングを用いることにしました。彼らが解明したい問いは以下の3点です。
市民の価値観は？ 人々は「通勤時間」と「通勤コスト」を、どのくらいの重みで天秤にかけているのか？ 時間を1分短縮するためなら、いくらまで余分に支払うだろうか？（支払い意欲額: WTP）
価値観はみんな同じ？ 「時は金なり」と考える高給取りのエンジニアもいれば、節約を第一に考える若手社員もいるはずです。特に「時間」に対する価値観は人によって大きく異なるのではないか？ この「個人の好みの異質性（ばらつき）」を捉えるため、チームは混合ロジットモデルを選択しました。
どの政策が効く？ 分析結果をもとに、以下の政策の効果を予測したい。
- 政策A: 電車の運賃を思い切って値下げする。
- 政策B: 最新技術を導入した自転車専用ハイウェイを整備し、自転車での通勤時間を大幅に短縮する。
この分析のため、チームは100人の市民にアンケート調査（選択実験）を実施し、様々な条件（時間とコストの組み合わせ）を提示して、どの交通手段を選ぶかを回答してもらいました。
なお、統計検定における有意水準は 5%としました。
シミュレーション
ステップ1: アンケートデータの生成
混合ロジットモデルで捉えたい「個人の異質性」を持った、市民一人ひとりに異なる「時間への感度」を持った アンケートデータを生成します。
library(dplyr)
seed <- 20250617
set.seed(seed)
# --- シナリオ設定 ---
n_individuals <- 100 # 調査対象の市民の数
n_scenarios <- 10 # 各市民が回答する設問（シナリオ）の数
alternatives <- c("car", "train", "bike") # 選択肢
# --- データ生成 ---
# ここで「市民のアンケート回答」を作成します
data_list <- list()
# 市民100人分のデータを作成
for (i in 1:n_individuals) {
# 市民iさんの「時間への感度」を個人差を持たせて設定
# 平均-0.1、標準偏差0.05の正規分布から生成。値が小さいほど時間に敏感。
beta_time_i <- rnorm(1, mean = -0.1, sd = 0.05)
# コストへの感度は全員共通と仮定 (-0.01)
beta_cost <- -0.01
# 各市民が10個の設問に回答
for (s in 1:n_scenarios) {
# 選択肢（車、電車、自転車）ごとに属性値を設定
# 単位 time_values:分、cost_values:円
time_values <- c(
runif(1, 20, 70), # 車の時間
runif(1, 30, 60), # 電車の時間
runif(1, 40, 90)
) # 自転車の時間
cost_values <- c(
runif(1, 300, 800), # 車のコスト
runif(1, 200, 500), # 電車のコスト
0
) # 自転車のコストは0
for (a in 1:length(alternatives)) {
# 各選択肢の効用（魅力度）を計算
# 効用 = 時間の効用 + コストの効用 + 観測できない誤差
utility <- beta_time_i * time_values[a] + beta_cost * cost_values[a] + rlogis(1)
data_list[[length(data_list) + 1]] <- data.frame(
id = i, # 個人ID
scenario_id = paste0(i, "_", s), # 設問ID
alt = alternatives[a], # 選択肢名
time = time_values[a],
cost = cost_values[a],
utility = utility
)
}
}
}
# リストを一つのデータフレームに変換
df <- do.call(rbind, data_list)
# 各設問で、最も効用が高い選択肢を「選ばれた(chosen=1)」とする
df <- df %>%
group_by(scenario_id) %>%
mutate(chosen = ifelse(utility == max(utility), 1, 0)) %>%
ungroup()
cat("--- アンケートデータの一部を確認---\n\n")
cat("--- head(df)---\n")
head(df)
cat("\n--- tail(df)---\n")
tail(df)--- アンケートデータの一部を確認---
--- head(df)---
# A tibble: 6 × 7
id scenario_id alt time cost utility chosen
<int> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 1 1_1 car 49.2 798. -9.37 0
2 1 1_1 train 46.6 318. -3.62 0
3 1 1_1 bike 87.6 0 0.901 1
4 1 1_2 car 57.6 649. -2.54 0
5 1 1_2 train 56.5 455. -4.08 0
6 1 1_2 bike 89.6 0 4.12 1
--- tail(df)---
# A tibble: 6 × 7
id scenario_id alt time cost utility chosen
<int> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 100 100_9 car 46.7 564. -12.7 0
2 100 100_9 train 41.4 347. -8.86 0
3 100 100_9 bike 53.9 0 -7.26 1
4 100 100_10 car 52.7 541. -13.3 0
5 100 100_10 train 53.9 223. -11.7 0
6 100 100_10 bike 62.3 0 -9.58 1ステップ2: 混合ロジットモデルの推定
続いて、市民の意思決定メカニズムを解き明かすモデルを構築します。mlogitパッケージを使い、「時間」の係数に個人差（正規分布）を仮定して推定します。
library(mlogit)
# mlogitパッケージが扱える形式にデータを変換
mlogit_data <- mlogit.data(
df,
choice = "chosen",
shape = "long",
alt.var = "alt",
id.var = "id",
chid.var = "scenario_id"
)
cat("--- mlogit_data の一部を確認---\n\n")
print(mlogit_data)--- mlogit_data の一部を確認---
~~~~~~~
first 10 observations out of 3000
~~~~~~~
# A tibble: 10 × 8
idx$chid $id $alt id scenario_id alt time cost utility chosen
<chr> <int> <fct> <int> <chr> <fct> <dbl> <dbl> <dbl> <lgl>
1 1_1 1 bike 1 1_1 bike 87.6 0 0.901 TRUE
2 1_1 1 car 1 1_1 car 49.2 798. -9.37 FALSE
3 1_1 1 train 1 1_1 train 46.6 318. -3.62 FALSE
4 1_10 1 bike 1 1_10 bike 43.3 0 1.70 TRUE
5 1_10 1 car 1 1_10 car 37.5 680. -4.59 FALSE
6 1_10 1 train 1 1_10 train 40.6 367. -2.67 FALSE
7 1_2 1 bike 1 1_2 bike 89.6 0 4.12 TRUE
8 1_2 1 car 1 1_2 car 57.6 649. -2.54 FALSE
9 1_2 1 train 1 1_2 train 56.5 455. -4.08 FALSE
10 1_3 1 bike 1 1_3 bike 74.3 0 0.524 TRUE
# A tibble: 3,000 × 3
chid id alt
<chr> <int> <fct>
1 1_1 1 bike
2 1_1 1 car
3 1_1 1 train
4 1_10 1 bike
5 1_10 1 car
6 1_10 1 train
7 1_2 1 bike
8 1_2 1 car
9 1_2 1 train
10 1_3 1 bike
# ℹ 2,990 more rows混合ロジットモデルを推定します。
なお、引数は以下の通りとします。
-
formula = chosen ~ time + cost | 0[選択結果変数] ~ [選択肢に依存する変数] | [個人に依存する変数]-
chosen: 「選択結果変数（被説明変数）」になります。mlogit_dataの中で、各設問（scenario_id）において、回答者が実際に選んだ選択肢（alt）には1（TRUE）、選ばなかった選択肢には0（FALSE）が入っています。モデルは、この1が選ばれる確率を最大化するようにパラメータを推定します。 -
~ time + cost: 「説明変数」になります。ここでは、選択肢の魅力度（効用）がtime（所要時間）とcost（費用）によって説明されることを意味します。これらの変数は、選択肢ごと（車、電車、自転車）に値が異なるため、「選択肢に依存する変数」になります。 -
| 0:|の右側には、「個人に依存する変数」（例えば、個人の年収や性別など、どの選択肢を選んでも値が変わらない変数）を指定します。今回はそのような変数をモデルに含めていないため、0を指定して「該当なし」を明示しています。 - まとめると: 「
chosen（選択されたかどうか）は、各選択肢のtimeとcostによって決まる」というモデルを定義しています。
-
rpar = c(time = "n")-
rparは “random parameters”（ランダムパラメータ）の略です。 - 意味: 「特定のパラメータ（係数）が、回答者間で一定ではなく、ある確率分布に従ってばらついている（ランダムである）」とモデルに指定します。
-
c(time = "n"):-
time =: 「timeという変数の係数（時間に対する感度）」をランダムパラメータとして扱います。 -
"n": その係数が、回答者間で正規分布 (Normal distribution) に従うと仮定します。
-
- もしこの引数がなければ:
mlogit(chosen ~ time + cost | 0, ...)のようにrparを指定しない場合、このモデルは「すべての回答者の時間とコストに対する感度は同じである」と仮定する、より単純な多項ロジットモデル（Multinomial Logit Model, MNL）になります。 - 推定されるもの: この指定により、
mlogitは以下の2つの値を推定します。-
time係数の平均値: 人々の時間感度の平均的な値。 -
time係数の標準偏差 (sd): 人々の時間感度が、その平均値からどれくらいばらついているか。標準偏差が大きいほど、人々の好みが多様であることを意味します。
-
-
-
R = 100- これは、ランダムパラメータを推定するためのシミュレーションの繰り返し回数を指定する引数です。（“Replications” や “Random draws” のRです）。
- 背景: 混合ロジットモデルは、解析的に解を求めるのが非常に困難です。そのため、シミュレーションを用いてパラメータを推定します（最尤シミュレーション法）。
- 処理の流れ: モデルは、
rparで指定された正規分布からランダムに係数の値をR回抽出し（ドロー）、各個人・各選択肢の効用を計算し、尤度（そのデータが観測されるもっともらしさ）を評価します。このプロセスを繰り返して、尤度が最大となる分布の平均値と標準偏差を探し出します。 - 値の目安:
Rの値が大きいほど、推定結果は安定し精度が上がりますが、計算時間は長くなります。今回の例では100としています。
-
halton = NA- これは、シミュレーションで用いる乱数の種類を指定する引数です。
-
halton = NA(デフォルト): これは「ハルトン数列 (Halton sequence)」と呼ばれる疑似乱数を使用することを意味します。- 通常の乱数（擬似乱数）:
rnorm()などで生成される乱数は、完全にランダムに見えますが、偶然に偏りが生じることがあります（例えば、最初の10個がすべてプラスの値になるなど）。 - ハルトン数列: より均一に空間をカバーするように設計された数列です。これにより、少ないシミュレーション回数（小さい
R）でも、より安定した推定結果が得られる傾向があります。
- 通常の乱数（擬似乱数）:
-
halton = NULL: 通常の擬似乱数を使用します。
# chosen ~ time + cost: 選択(chosen)は時間(time)と費用(cost)で決まる
# rpar = c(time = "n"): timeの係数が個人間で正規分布("n")に従って変動すると仮定
model <- mlogit(
chosen ~ time + cost | 0,
data = mlogit_data,
rpar = c(time = "n"), # "n" for normal distribution
R = 100, # シミュレーションの繰り返し回数
halton = NA # 準乱数/ハルトン数列
)
cat("--- 結果の確認 ---\n")
summary(model)--- 結果の確認 ---
Call:
mlogit(formula = chosen ~ time + cost | 0, data = mlogit_data,
rpar = c(time = "n"), R = 100, halton = NA)
Frequencies of alternatives:choice
bike car train
0.641 0.097 0.262
bfgs method
9 iterations, 0h:0m:6s
g'(-H)^-1g = 1.21E-06
successive function values within tolerance limits
Coefficients :
Estimate Std. Error z-value Pr(>|z|)
time -0.07387138 0.00633553 -11.6599 < 2.2e-16 ***
cost -0.00758946 0.00045848 -16.5536 < 2.2e-16 ***
sd.time 0.04340723 0.01122239 3.8679 0.0001098 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Log-Likelihood: -690.19
random coefficients
Min. 1st Qu. Median Mean 3rd Qu. Max.
time -Inf -0.1031491 -0.07387138 -0.07387138 -0.04459365 Infステップ3: 結果の解釈と政策提言
モデルの推定結果から、都市計画チームは市長に報告する洞察を得ます。
【市長への報告内容】
- 市民の平均的な価値観:
-
costの係数（-0.00758946）とtimeの係数（-0.07387138）は、どちらもマイナス かつ p値は設定した有意水準を下回っています。つまり、市民はコストと時間がかかるほどその選択肢を嫌うことを示しています。 - 時間の価値（支払い意欲額, WTP）は
|time係数 / cost係数|で計算できます。abs(-0.07387138 / -0.00758946) = 9.73これは「通勤時間を1分短縮するためなら、市民は平均して約9.7円を追加で支払う意思がある」ことを意味します。
-
- 価値観の多様性（混合ロジットモデルの成果）:
-
sd.time（時間の係数の標準偏差）が 0.04340723 と推定され、こちらもp値は設定した有意水準を下回っています。 - これは「時間に対する価値観は、人によって異なる」という仮説を裏付ける結果です。
-
ステップ4: 政策シミュレーション
推定したモデルを使って、政策Aと政策Bの効果を予測します。
# --- 政策シミュレーション ---
# 1. 現状のシェアを予測
original_data <- mlogit_data
pred_original <- apply(predict(model, newdata = original_data), 2, mean)
# 2. 政策A: 電車の運賃を100円値下げ
policy_A_data <- original_data
policy_A_data$cost[policy_A_data$alt == "train"] <- policy_A_data$cost[policy_A_data$alt == "train"] - 100
pred_A <- apply(predict(model, newdata = policy_A_data), 2, mean)
# 3. 政策B: 自転車の所要時間を平均15分短縮
policy_B_data <- original_data
policy_B_data$time[policy_B_data$alt == "bike"] <- policy_B_data$time[policy_B_data$alt == "bike"] - 15
pred_B <- apply(predict(model, newdata = policy_B_data), 2, mean)
# 4. 結果を比較
results <- data.frame(
Status_Quo = pred_original,
Policy_A_Train_Discount = pred_A,
Policy_B_Bike_Highway = pred_B
)
print(round(results * 100, 2)) # パーセント表示 Status_Quo Policy_A_Train_Discount Policy_B_Bike_Highway
bike 63.97 54.76 79.37
car 10.54 8.48 6.24
train 25.49 36.76 14.39【最終政策提言】
シミュレーション結果から、以下の提言ができます。
- 政策A（電車の運賃を100円値下げ）: 電車のシェアが 25%から37%へと増加 し、バイクのシェアは 64% → 55% 、車のシェアは 11% → 8% と、両者ともに低下することが予測されます。
- 政策B（自転車道整備により自転車の所要時間を平均15分短縮）: 自転車のシェアが 64%から79%に増加 し、車のシェアは 11% → 6% 、電車のシェアは 25% → 14% と、両者ともに低下することが予測されます。
以上です。
