はじめに

つい最近、階層ロジスティック回帰というモデルの存在を初めて知りました。ロジスティック回帰にも階層モデルがあるんだと結構衝撃を受けました。知ってみると、意外と便利そうなモデルだったので、今回は階層ロジスティック回帰モデルを題材に記事を書いていこうと思います。
絶対この記事は伸びない自信がありますね...なぜなら自己満記事なのと、モデル自体がニッチな記事だからですね

従来のロジスティック回帰における課題

従来のロジスティック回帰モデルには、主に三つの重要な課題が存在していました。
最初の課題は、グループ構造の無視です。従来のモデルでは、データ点が独立であることを仮定していましたが、実際の多くのデータには階層構造が存在します。例えば、同じ学校の生徒のデータは互いに相関を持つ可能性が高く、この依存関係を無視することは推定の精度に影響を与えていました。
二つ目の課題は、変量効果の取り扱いです。グループ間での効果の違いを適切にモデル化できないため、各グループの特性を反映した予測が困難でした。例えば、学校ごとの教育効果の違いを考慮することができませんでした。
三つ目の課題は、予測の不確実性の過小評価です。グループ構造を無視することで、予測の信頼区間が実際よりも狭く推定される傾向がありました。これは特に、新しいグループに対する予測を行う際に問題となっていました。

階層ロジスティック回帰モデルの設計思想

まず、階層構造の明示的なモデル化です。これにより、データの自然な構造を反映した分析が可能となり、より現実に即した推論が可能となりました。
次に、変量効果と固定効果の統合です。グループレベルとグローバルレベルの効果を同時にモデル化することで、より柔軟な予測が可能となりました。
さらに、ベイズ的なアプローチの採用です。パラメータの不確実性を適切に考慮することで、より信頼性の高い予測区間を得ることができるようになりました。

階層ロジスティック回帰モデルの理論

モデルの基本構造

階層ロジスティック回帰モデルの基本構造は、以下の式で表現されます：

$logit(p_{ij}) = X_{ij}\beta + Z_{ij}u_j$

ここで、 $p_{ij}$ は $j$ 番目のグループの $i$ 番目の観測値が正のクラスに属する確率、 $X_{ij}$ は固定効果の説明変数、 $\beta$ は固定効果のパラメータ、 $Z_{ij}$ は変量効果の説明変数、 $u_j$ はグループ $j$ の変量効果を表します。

この式は、人間の意思決定プロセスに例えると理解しやすいでしょう。例えば、ある商品を購入するかどうかの決定において、商品の価格や品質といった普遍的な要因（固定効果）と、各個人の好みや経験（変量効果）の両方が影響を与えるような状況を表現しています。

変量効果の分布

変量効果は以下のような多変量正規分布に従うと仮定します：

$u_j \sim N(0, \Sigma)$

ここで、 $\Sigma$ は変量効果の分散共分散行列を表します。この行列は、グループ内での効果の相関構造を捉えています。

この仮定は、例えば、学校の教育効果を考える際に理解しやすくなります。各学校には固有の特徴がありますが、それらは完全にランダムではなく、ある程度の一貫性を持って分布していると考えるのが自然です。

推定方法

モデルのパラメータ推定は、以下の尤度関数の最大化によって行われます：

$L(\beta, \Sigma) = \prod_{j=1}^J \int \prod_{i=1}^{n_j} p_{ij}^{y_{ij}}(1-p_{ij})^{1-y_{ij}} f(u_j|\Sigma) du_j$

ここで、 $y_{ij}$ は応答変数の実測値、 $f(u_j|\Sigma)$ は変量効果の密度関数を表します。この複雑な積分は、通常、ラプラス近似や適応的ガウス求積法によって近似計算されます。

実践的応用

具体的な応用例として、全国チェーン店での購買行動分析を見てみましょう。顧客の購買確率は以下のようにモデル化されます：

$logit(p_{ijk}) = \beta_0 + \beta_1 Price_{ijk} + \beta_2 Promotion_{ijk} + u_{jk} + v_k$

ここで、 $p_{ijk}$ は店舗 $k$ における顧客 $j$ の商品 $i$ の購買確率、 $u_{jk}$ は店舗固有の顧客効果、 $v_k$ は店舗効果を表します。

このモデルは、例えば、あるコーヒーチェーン店での購買行動を考えると分かりやすくなります。お客様が商品を購入するかどうかは、商品の価格やプロモーションといった全店共通の要因（固定効果）に加えて、各店舗の雰囲気や立地といった店舗固有の要因（店舗効果）、さらには各店舗における常連客の購買パターン（店舗固有の顧客効果）によって影響を受けます。
実際の分析例として、ある全国チェーン店での新商品導入の効果分析を見てみましょう。データには以下の構造が含まれます：
　　
店舗レベル変数：
立地タイプ、店舗面積、開店からの経過年数
　　
顧客レベル変数：
年齢、性別、会員ステータス、過去の購買履歴
　　
商品レベル変数：
価格、プロモーション有無、商品カテゴリー
　　
このような複雑な階層構造を持つデータに対して、階層ロジスティック回帰モデルを適用することで、以下のような知見を得ることができます：新商品の受容性の地域差：都心型店舗と郊外型店舗での新商品の受け入れられ方の違いを定量化できます。顧客セグメント別の反応：年齢層や会員ステータスによって、新商品への反応がどのように異なるかを把握できます。店舗特性の影響：店舗の規模や経過年数が、新商品の売上にどのように影響するかを理解できます。これらの知見は、今後の商品開発やマーケティング戦略の立案に直接活用することができます。例えば、店舗特性に応じた品揃えの最適化や、顧客セグメント別のプロモーション施策の設計などに活用されています。

実装例

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.formula.api as smf
import cmdstanpy
import numpy as np

np.random.seed(42)

# ------------------------
# 1. パラメータの真の値を設定
# ------------------------
beta0_true = -1.0   # 切片
beta1_true = -2.0   # 価格の影響
beta2_true =  1.5   # プロモーションの影響

sigma_store_true = 0.5       # 店舗レベルのランダム効果の分散
sigma_store_cust_true = 0.8  # 店舗×顧客レベルのランダム効果の分散

# ------------------------
# 2. データ構造の設定
# ------------------------
K = 5       # 店舗数
J = 10      # 各店舗あたりの顧客数
N_item = 20 # 1顧客あたりの購入機会（商品数）を仮定

store_ids = np.arange(K)  # 店舗ID(0~K-1)
customer_ids_in_store = [np.arange(J) for _ in range(K)]

# ------------------------
# 3. ランダム効果の生成
#    v_k ~ Normal(0, sigma_store_true)
#    u_{jk} ~ Normal(0, sigma_store_cust_true)
# ------------------------
v_k = np.random.normal(0, sigma_store_true, size=K)  # 各店舗のランダム効果
u_jk = []
for k in range(K):
    u_jk_k = np.random.normal(0, sigma_store_cust_true, size=J)
    u_jk.append(u_jk_k)

# ------------------------
# 4. サンプルデータの生成
# ------------------------
data_list = []
for k in store_ids:
    for j in customer_ids_in_store[k]:
        for i in range(N_item):
            # Price: 0~10 の一様分布からランダムに
            price_ijk = np.random.uniform(0, 10)
            # Promotion: 0 or 1 を50%で
            promotion_ijk = np.random.binomial(1, 0.5)
            
            # ロジットを計算
            logit_p = (beta0_true 
                       + beta1_true * price_ijk
                       + beta2_true * promotion_ijk
                       + v_k[k]
                       + u_jk[k][j])
            
            # p_ijk をシグモイドで変換
            p_ijk = 1 / (1 + np.exp(-logit_p))
            
            # 購買(1) or 非購買(0) を発生させる
            y_ijk = np.random.binomial(1, p_ijk)
            
            data_list.append({
                "store_id": k,
                "customer_id": j,
                "price": price_ijk,
                "promotion": promotion_ijk,
                "purchase": y_ijk
            })

df = pd.DataFrame(data_list)

# 階層を無視するので、price と promotion のみを説明変数に使用
# ロジスティック回帰
model_simple = smf.glm(
    formula="purchase ~ price + promotion",
    data=df,
    family=sm.families.Binomial()  # 修正: statsmodels.api を使用してファミリーを指定
).fit()

print(model_simple.summary())

# まずデータを Stan に渡すための形に整形
# store_id, customer_id は Stan 内で 1-based index が必要なので +1 しておく
stan_data = {
    "N": len(df),
    "K": K,
    "J": J,
    "store_id": df["store_id"].values + 1,
    "customer_id": df["customer_id"].values + 1,
    "price": df["price"].values,
    "promotion": df["promotion"].values.astype(int),
    "purchase": df["purchase"].values.astype(int)
}

beta0 = model_simple.params["Intercept"]
beta1 = model_simple.params["price"]
beta2 = model_simple.params["promotion"]

# モデルコンパイル
stan_model = cmdstanpy.CmdStanModel(stan_file="hierarchical_logistic.stan")

# サンプリング実行
fit = stan_model.sample(
    data=stan_data,
    chains=4,
    parallel_chains=4, 
    seed=42,
    iter_warmup=1000,
    iter_sampling=1000
)

# 結果のサマリーを表示
print(fit.summary())

# 事後サンプルを取り出して真の値と比較
posterior = fit.draws_pd()
beta0_est = posterior["beta0"]
beta1_est = posterior["beta1"]
beta2_est = posterior["beta2"]

print("---- 通常ロジスティック回帰の推定結果 ----")
print(f"beta0 (Intercept):, {beta0}, true: {beta0_true}")
print(f"beta1 (price):, {beta1}, true: {beta1_true}")
print(f"beta2 (promotion):, {beta2}, true: {beta2_true}")
print("---- 階層ロジスティック回帰 (cmdstanpy) の推定結果 ----")
print(f"beta0 mean: {beta0_est.mean():.3f}, true: {beta0_true}")
print(f"beta1 mean: {beta1_est.mean():.3f}, true: {beta1_true}")
print(f"beta2 mean: {beta2_est.mean():.3f}, true: {beta2_true}")

v_columns = posterior.filter(like="v[", axis=1)  # `v[` で始まる列を選択
v_est = v_columns.mean().values  # 平均値を取得

# シミュレーションデータで店舗ごとの購買確率を推定
sim_prices = np.linspace(0, 10, 100)
sim_promotions = [0, 1]

plt.figure(figsize=(10, 6))
for k in range(K):
    v_k_effect = v_est[k]
    for promo in sim_promotions:
        p_sim = 1 / (1 + np.exp(-(beta0_est.mean() + beta1_est.mean() * sim_prices + beta2_est.mean() * promo + v_k_effect)))
        label = f"Store {k}, Promotion={promo}"
        plt.plot(sim_prices, p_sim, label=label)

plt.xlabel("Price")
plt.ylabel("Purchase Probability")
plt.title("Purchase Probability by Store and Promotion")
plt.legend()
plt.show()

これが階層ロジスティック回帰モデルの結果でこれが通常のロジスティック回帰モデルの結果です通常のロジスティック回帰モデルは、全店舗共通の固定効果（価格、プロモーション）のみを考慮。店舗間の差異を無視しているため、全ての店舗で同じ確率曲線となる。一方階層ロジスティック回帰モデルでは、店舗間のランダム効果を考慮しているため、店舗ごとに異なる確率曲線を出せる。
このような階層構造をモデルに取り込める部分はやはり利点だなと感じますね。実際の階層構造はより複雑なので、もっと実際のモデルは複雑になるはずですが、今回はあえてシンプルにしてみました。