モデル評価 - データサイエンス関連+αのメモ

コンジョイント分析で推定した部分効用モデルが、実際の選好・選択をどれだけよく説明・予測できるかを評価する観点は大きく2つある。

内的妥当性（internal validity）：モデルが推定に使ったデータ自体をどれだけよく説明できているか（当てはまりの良さ）
予測的妥当性（predictive validity）：モデルが推定に使っていない新しいデータ（実験計画（属性・水準の設計）で触れたホールドアウトタスク）をどれだけ正しく予測できるか

コンジョイント分析の目的は「未知の製品案・市場シナリオを予測すること」（マーケットシミュレーション）にあるため、実務では予測的妥当性の方が重視される。

McFaddenの疑似決定係数（pseudo R-squared）¶

線形回帰の $R^2$ に相当する当てはまりの指標として、選択モデルでは**McFaddenの疑似 $R^2$ **が使われる。

R^2_{\text{McFadden}} = 1 - \frac{\log L(\hat{\boldsymbol{\beta}})}{\log L(\mathbf{0})}

(1)

$\log L(\hat{\boldsymbol{\beta}})$ は推定したモデルの対数尤度、 $\log L(\mathbf{0})$ はすべての説明変数を除いた（選択肢を等確率で選ぶ）ヌルモデルの対数尤度である。線形回帰の $R^2$ と異なり $0〜1$ には収まらない場合もあるが、経験的には $0.2 〜 0.4$ 程度でも当てはまりが良いとされる（McFadden, 1974）。

ヒット率（Hit Rate）¶

ホールドアウトタスク（推定に使っていない選択課題）それぞれについて、モデルが最も選択確率が高いと予測した選択肢が、実際に回答者が選んだ選択肢と一致した割合を**ヒット率（hit rate）**と呼ぶ。

\text{Hit Rate} = \frac{1}{T} \sum_{t=1}^{T} \mathbb{1}\left[\hat{y}_t = y_t\right], \quad \hat{y}_t = \arg\max_{j \in C_t} \hat{P}(j \mid C_t)

(2)

選択セットの選択肢数が $J$ 個であれば、何の情報もなく当てずっぽうに選んだ場合の期待的中率（チャンスレベル）は $1/J$ になる。ヒット率がチャンスレベルをどれだけ上回っているかで、モデルの実質的な予測力を評価する。

平均絶対誤差（MAE）によるシェア予測の評価¶

個々のタスクの的中だけでなく、マーケットシミュレーションで行うような集計レベルのシェア予測の精度も重要である。ホールドアウトタスクにおける実際の選択シェアと、モデルが予測したシェア・オブ・プリファレンスとの**平均絶対誤差（Mean Absolute Error, MAE）**で評価する。

\text{MAE} = \frac{1}{J} \sum_{j=1}^{J} \left| \text{share}_j^{\text{actual}} - \text{share}_j^{\text{predicted}} \right|

(3)

実装例¶

選択型コンジョイント分析（CBC）と多項ロジットモデルと同じ設定でCBCデータをシミュレーションし、データを推定用（estimation）とホールドアウト用（holdout）に分割して、プールドロジットモデルの予測的妥当性を評価する。

import itertools
import numpy as np
import pandas as pd

rng = np.random.default_rng(0)

attributes = {
    "価格":    ["1,000円", "1,500円", "2,000円"],
    "容量":    ["500ml", "1000ml"],
    "ブランド": ["A社", "B社", "C社"],
}
true_beta = {
    "価格::1,000円": 1.2, "価格::1,500円": 0.0, "価格::2,000円": -1.2,
    "容量::500ml": -0.3, "容量::1000ml": 0.3,
    "ブランド::A社": 0.6, "ブランド::B社": 0.1, "ブランド::C社": -0.7,
}

all_profiles = list(itertools.product(*attributes.values()))

def profile_to_dummies(profile):
    price, volume, brand = profile
    return {
        "価格::1,000円": price == "1,000円", "価格::1,500円": price == "1,500円", "価格::2,000円": price == "2,000円",
        "容量::500ml": volume == "500ml", "容量::1000ml": volume == "1000ml",
        "ブランド::A社": brand == "A社", "ブランド::B社": brand == "B社", "ブランド::C社": brand == "C社",
    }

def utility(profile):
    d = profile_to_dummies(profile)
    return sum(true_beta[k] for k, v in d.items() if v)

N_RESPONDENTS = 500
N_ALTS = 3

records = []
for n in range(N_RESPONDENTS):
    choice_set = rng.choice(len(all_profiles), size=N_ALTS, replace=False)
    utilities = np.array([utility(all_profiles[j]) for j in choice_set])
    utilities_noisy = utilities + rng.gumbel(size=N_ALTS)
    chosen = np.argmax(utilities_noisy)
    for alt_idx, profile_idx in enumerate(choice_set):
        d = profile_to_dummies(all_profiles[profile_idx])
        records.append({
            "resp_id": n, "alt_id": alt_idx, "chosen": int(alt_idx == chosen), **d,
        })

df = pd.DataFrame(records)

# 回答者の8割を推定用、2割をホールドアウト用に分割
resp_ids = df["resp_id"].unique()
rng.shuffle(resp_ids)
n_holdout = int(len(resp_ids) * 0.2)
holdout_resp = set(resp_ids[:n_holdout])

df_train = df[~df["resp_id"].isin(holdout_resp)]
df_holdout = df[df["resp_id"].isin(holdout_resp)]

print(f"推定用: {df_train['resp_id'].nunique()}人, ホールドアウト用: {df_holdout['resp_id'].nunique()}人")

推定用: 400人, ホールドアウト用: 100人

from statsmodels.discrete.conditional_models import ConditionalLogit

feature_cols = [
    "価格::1,000円", "価格::1,500円",
    "容量::500ml",
    "ブランド::A社", "ブランド::B社",
]

model = ConditionalLogit(df_train["chosen"], df_train[feature_cols], groups=df_train["resp_id"])
result = model.fit()

# McFaddenの疑似R^2(推定用データに対する当てはまり)
loglik_model = result.llf
loglik_null = len(df_train["resp_id"].unique()) * np.log(1 / N_ALTS)
pseudo_r2 = 1 - loglik_model / loglik_null
print(f"log-likelihood(model) = {loglik_model:.2f}")
print(f"log-likelihood(null)  = {loglik_null:.2f}")
print(f"McFadden's pseudo R^2 = {pseudo_r2:.3f}")

log-likelihood(model) = -324.61
log-likelihood(null)  = -439.44
McFadden's pseudo R^2 = 0.261

beta_hat = result.params

def predict_probs(group):
    V = group[feature_cols].to_numpy(dtype=float) @ beta_hat.to_numpy()
    p = np.exp(V) / np.exp(V).sum()
    return p

hits = []
actual_shares = np.zeros(N_ALTS)
predicted_shares = np.zeros(N_ALTS)

for resp_id, group in df_holdout.groupby("resp_id"):
    group = group.reset_index(drop=True)
    p = predict_probs(group)
    predicted_choice = np.argmax(p)
    actual_choice = group.index[group["chosen"] == 1][0]

    hits.append(predicted_choice == actual_choice)
    actual_shares[actual_choice] += 1
    predicted_shares += p

n_holdout_tasks = df_holdout["resp_id"].nunique()
hit_rate = np.mean(hits)
chance_level = 1 / N_ALTS

actual_shares /= n_holdout_tasks
predicted_shares /= n_holdout_tasks
mae = np.mean(np.abs(actual_shares - predicted_shares))

print(f"ヒット率: {hit_rate:.1%} (チャンスレベル: {chance_level:.1%})")
print(f"シェア予測のMAE: {mae:.3f}")
pd.DataFrame({"actual_share": actual_shares, "predicted_share": predicted_shares})

ヒット率: 67.0% (チャンスレベル: 33.3%)
シェア予測のMAE: 0.036

ホールドアウトデータに対するヒット率がチャンスレベル（ $1/3 \approx 33\%$ ）を大きく上回っていれば、モデルが実質的な予測力を持っていると判断できる。また集計レベルのシェア予測についても、実際のシェアと推定シェアの差（MAE）が小さいほど、マーケットシミュレーションでの予測が信頼できることを意味する。

プールドロジット vs 階層ベイズの評価¶

階層ベイズモデルによる個人レベル部分効用の推定で個人レベルの部分効用 $\hat{\boldsymbol{\beta}}_n$ を推定した場合、個人ごとの $\hat{\boldsymbol{\beta}}_n$ を使ってホールドアウトタスクを予測すると、通常はプールドロジット（全員共通の $\hat{\boldsymbol{\beta}}$ ）よりもヒット率が高くなる。これは選好の異質性を捉えられているかどうかの実証的な確認にもなる。ただし個人内のホールドアウトタスク数が極端に少ない場合、個人レベルモデルは推定が不安定になり、かえってプールドモデルに劣ることもあるため、両者を比較検証することが望ましい。