階層ベイズGeoモデル・HSGP・リフトテスト統合 — PyMC-Marketing上級テクニック

はじめに

PyMC-Marketingの基本的なMMMパイプライン（Adstock → Saturation → フィッティング → ROAS → 予算最適化）は、多くのマーケティング分析課題に対応できます。しかし、実務ではそれだけでは不十分な場面が数多く存在します。

• 東京と北海道でTV広告の効果は同じなのか？
• 1年前と今でデジタル広告の効率は変わっていないか？
• A/Bテストの結果をMMMにどう取り込むか？
• 新商品は既存商品の売上を食っていないか？
• モデルの予測精度を時系列で正しく評価するには？

これらの問いに答えるのが、本記事で解説する上級テクニックです。PyMC-Marketingにはこれらの機能がすべて組み込まれており、他のMMMツール（Meridian, Robyn, LightweightMMM）にはない独自の強みとなっています。

本記事のスコープ:

1. 階層ベイズGeoモデル（multidimensional.MMM）
   → 複数地域を1つのモデルで管理、Partial Pooling

2. HSGP（ヒルベルト空間ガウス過程）
   → 時変パラメータの推定、O(n)近似

3. リフトテスト統合（add_lift_test_measurements）
   → 実験データとMMMの統合、不確実性低減

4. 顧客選択モデル（MV-ITS）
   → 新商品カニバリゼーションの定量化

5. 時系列交差検証（TimeSliceCrossValidation）
   → リーケージなしのモデル性能評価

1. 階層ベイズGeoモデル — 地域差をモデルに組み込む

なぜ階層モデルが必要なのか

複数の地域・マーケットでマーケティング活動を行う企業では、「地域ごとにMMMを構築する」か「全地域をプールして1つのMMMを構築する」かの選択に直面します。どちらにも明確な限界があります。

個別モデル（No Pooling: 地域ごとに独立）:
  メリット:
    - 地域差を完全に捉えられる
  デメリット:
    - データが少ない地域で推定が不安定（過学習リスク）
    - 地域数が多いとモデル管理が煩雑
    - 地域間の共通パターンを活用できない

統合モデル（Complete Pooling: 全地域プール）:
  メリット:
    - データ量が多く安定した推定
  デメリット:
    - 地域差を完全に無視（東京と北海道が同じパラメータ）
    - 地域固有の施策効果を捉えられない
    - Simpson's paradox（合成の誤謬）のリスク

階層モデル（Partial Pooling）:  <-- PyMC-Marketingの解法
  メリット:
    - 地域差を許容しつつ、統計的強度を共有
    - データが少ない地域は全体平均に引き寄せ（縮小推定/Shrinkage）
    - 1つのモデルで全地域を管理
    - 新規地域の追加が容易（全体の事前分布が適用される）
  考慮点:
    - 計算コストが高い（パラメータ数が増加）
    - 最低限の地域数とデータ量が必要

Partial Pooling（部分プーリング）の直感的理解

Partial Poolingは、「完全に独立」と「完全に同一」の中間に位置する考え方です。

直感的な例で理解する:

  あなたは全国チェーンのマーケティング部長です。
  東京（データ豊富）と新規出店した鳥取（データ3ヶ月分）の
  TV広告効果を推定したいとします。

  No Pooling（独立推定）:
    東京: ROAS = 2.5（信頼区間 [2.0, 3.0]）--- 安定
    鳥取: ROAS = 8.0（信頼区間 [0.5, 15.5]）--- 不安定！
    → 鳥取のデータが少なすぎて、推定がブレまくる

  Complete Pooling（完全統合）:
    全地域: ROAS = 2.8（信頼区間 [2.3, 3.3]）
    → 全地域で同じ推定。鳥取の独自性が完全に無視される

  Partial Pooling（階層モデル）:
    東京: ROAS = 2.5（信頼区間 [2.0, 3.0]）--- ほぼ独自推定
    鳥取: ROAS = 3.2（信頼区間 [1.8, 4.6]）--- 全体平均に引き寄せ
    → 東京はデータが豊富なので、独自の推定をほぼ維持
    → 鳥取はデータが少ないので、全体平均に向かって「縮小」
    → 8.0→3.2 に引き寄せられたが、東京（2.5）とは異なる推定

  これが「縮小推定（Shrinkage）」の効果:
    → データが少ない地域ほど、全体平均に強く引き寄せられる
    → データが豊富な地域は、独自のパラメータ推定を維持
    → この「引き寄せの強さ」はデータから自動的に学習される

データ準備とcheck_geo_data()

階層モデルには「地域 x 日付」のパネルデータが必要です。

import pandas as pd
import numpy as np

# パネルデータの構造（Long形式）
# date     | region   | sales | tv_spend | digital_spend | social_spend
# 2024-01-01 | Tokyo    | 50000 | 10000    | 5000          | 3000
# 2024-01-01 | Osaka    | 35000 | 8000     | 4000          | 2500
# 2024-01-01 | Hokkaido | 12000 | 3000     | 1500          | 1000
# 2024-01-08 | Tokyo    | 52000 | 11000    | 5500          | 3200

def check_geo_data(df, date_col, geo_col, channel_cols, target_col):
    """Geoモデル用パネルデータの品質チェック"""
    regions = df[geo_col].unique()
    dates = df[date_col].unique()

    print(f"=== Geoデータ品質チェック ===")
    print(f"地域数: {len(regions)}")
    print(f"日付数: {len(dates)}")
    print(f"期待レコード数: {len(regions) * len(dates)}")
    print(f"実レコード数: {len(df)}")

    # バランス確認（全地域 x 全日付の組み合わせが存在するか）
    expected = len(regions) * len(dates)
    actual = len(df)
    if actual < expected:
        missing = expected - actual
        print(f"  {missing}レコードが欠損（{missing/expected*100:.1f}%）")
    else:
        print("  パネルデータはバランス済み")

    # 地域別データ量
    print(f"\n地域別データ量:")
    for region in sorted(regions):
        n = len(df[df[geo_col] == region])
        mean_target = df[df[geo_col] == region][target_col].mean()
        print(f"  {region}: {n}行, 平均{target_col}={mean_target:,.0f}")

    # 地域別の欠損チェック
    for col in channel_cols + [target_col]:
        nulls = df[col].isnull().sum()
        if nulls > 0:
            print(f"  {col}: {nulls}件の欠損値")

check_geo_data(df_geo, "date", "region",
               ["tv_spend", "digital_spend", "social_spend"], "sales")

期待される出力:

=== Geoデータ品質チェック ===
地域数: 8
日付数: 104
期待レコード数: 832
実レコード数: 832
  パネルデータはバランス済み

地域別データ量:
  Hokkaido: 104行, 平均sales=12,500
  Kyushu: 104行, 平均sales=18,200
  Osaka: 104行, 平均sales=35,800
  Tokyo: 104行, 平均sales=52,300
  ...

Geoモデルの共通の落とし穴

Geoモデルを適用する前に、以下の条件を確認してください。

実装

from pymc_marketing.mmm.multidimensional import MMM as MultiMMM
from pymc_marketing.mmm import GeometricAdstock, LogisticSaturation

geo_mmm = MultiMMM(
    date_column="date",
    channel_columns=["tv_spend", "digital_spend", "social_spend"],
    geo_column="region",            # <- 地域列を指定
    adstock=GeometricAdstock(l_max=8),
    saturation=LogisticSaturation(),
    yearly_seasonality=2,
)

geo_mmm.fit(
    X=df_geo,
    target_column="sales",
    draws=1000,
    tune=1000,
    chains=4,
    target_accept=0.9,
)

期待される出力:

Sampling 4 chains: 100%|##########| 8000/8000 [08:30<00:00, 15.69draws/s]

解釈ガイド:

• Geoモデルは通常のMMMより3-5倍の計算時間が必要（パラメータ数が地域数に比例して増加）
• target_accept=0.9 はGeoモデルの推奨設定。divergenceが出る場合は0.93-0.95に引き上げ

地域別の効果解釈（HDIの読み方）

# 地域別のチャネル効果を可視化
geo_mmm.plot_channel_contribution_share_hdi()

# 出力例:
#
# Region: Tokyo
#   TV:      38% [30%, 46%]  <- HDIが狭い = 信頼できる推定
#   Digital: 40% [33%, 47%]
#   Social:  22% [16%, 28%]
#
# Region: Osaka
#   TV:      35% [27%, 43%]
#   Digital: 42% [34%, 50%]  <- Digital効果が東京より若干高い
#   Social:  23% [17%, 29%]
#
# Region: Hokkaido
#   TV:      45% [35%, 55%]  <- TVの効果が東京より高い
#   Digital: 35% [25%, 45%]
#   Social:  20% [12%, 28%]  <- HDIが広い = データが少なく不確実

解釈ガイド:

HDI（Highest Density Interval）の読み方:

  1. HDI幅 = 推定の不確実性
     → 幅が狭い（例: [30%, 46%]）: データが十分で信頼できる推定
     → 幅が広い（例: [12%, 28%]）: データ不足で不確実

  2. 地域間の比較:
     → HDIが重なっていれば「地域差は統計的に明確でない」
     → HDIが重なっていなければ「地域差がある可能性が高い」

  3. Shrinkageの確認:
     → データが少ない地域のHDIが全体平均付近に集中 -> Shrinkage機能中
     → データが豊富な地域のHDIが独自の位置 -> 独立推定に近い

Shrinkage推定の可視化

階層モデルの核心である縮小推定の効果を視覚的に確認します。

import arviz as az
import matplotlib.pyplot as plt

# 地域別のチャネル効果パラメータを抽出
posterior = geo_mmm.fit_result.posterior

# beta_channel の事後分布を地域別に可視化
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
channels = ["tv_spend", "digital_spend", "social_spend"]
regions = df_geo["region"].unique()

for ax, ch_idx in zip(axes, range(3)):
    for region in regions:
        # 各地域のbeta_channel事後分布
        samples = posterior["beta_channel"].sel(
            channel=channels[ch_idx], geo=region
        ).values.flatten()
        ax.hist(samples, bins=40, alpha=0.4, label=region, density=True)

    ax.set_title(f"{channels[ch_idx]}", fontsize=12)
    ax.legend()
    ax.set_xlabel("Effect Size")

plt.suptitle("Regional Channel Effects (Shrinkage Visualization)", fontsize=14)
plt.tight_layout()
plt.show()

期待される出力:

• 各チャネルについて、地域ごとの事後分布がヒストグラムで表示される
• データが少ない地域（例: 鳥取）の分布が中央に集中 -> Shrinkageが機能
• データが豊富な地域（例: 東京）の分布が独自の位置 -> 独立推定に近い

解釈ガイド:

Shrinkage可視化の読み方:

  正常なShrinkage:
    - データの少ない地域の分布が全体平均付近に集まっている
    - データの豊富な地域の分布が独自の位置にある
    - 地域間で適度な差がある

  異常なパターン:
    - 全地域が同じ位置 -> 事前分布が強すぎる（group-levelのsigmaを広げる）
    - 一部の地域が極端な位置 -> データ異常の可能性（外れ値チェック）
    - 分布が非常に広い -> データ不足（地域の統合を検討）

Geoモデルのよくある落とし穴