🎯 Phase 1: 統計的因果推論の理論とPython実装
Week 2-5（4週間） | 因果推論の基礎を理解し、交絡を調整した効果推定をPythonで実装する
📊 Phase 1で達成すること
潜在的アウトカム、ATEを定義し、因果の言語で議論できる
交絡構造を特定し、調整すべき変数を判断できる
交絡調整のツールとして回帰分析を使いこなせる
Propensity Scoreを算出し、マッチングまたは層別化を実装できる
📚 使用教材
主教材
- 理論書: 金本『因果推論』（第1-4章）
- 因果推論の基本的な枠組み、DAG、回帰分析による推定
- 統計的因果推論の理論的基礎を体系的に学習
- 実装教材: Udemy『データサイエンスの為の因果推論入門』（Python利用）
- Python（pandas, statsmodels, scikit-learn）での実装例
- 傾向スコア、回帰による推定を中心に学習
- 実データを使った実践的なハンズオン
副読書・資料
- 概念書: 『統計的因果推論の基礎』(Pearl著) イントロダクションおよびDAG部分
- 実装資料: statsmodels公式ドキュメント（OLS, Logit, WLS）
- 補助教材: 『効果検証入門』（技術評論社）対応箇所
- ✅ Python環境（Jupyter Notebook, pandas, numpy, statsmodels）が動作すること
- ✅ 基本的な統計（平均、分散、回帰分析の係数の解釈）を理解していること
- ✅ Phase 0で「相関≠因果」の雰囲気を掴んでいること
- ✅ 統計検定準1級レベルの確率・統計の基礎知識
- 広告効果測定: 「特定の広告（施策）がCVRに与える純粋な影響」を年齢や過去購買履歴の**交絡**を除去して推定する
- コンテンツ効果: 「ブログ記事Aを読んだこと（施策）が、その後のサイト滞在時間に与える影響」を推定する
- 価格弾力性: 価格変動が売上に与える影響を、競合の動向などの交絡を調整して計測する
- サイト適用: オフィス物件の問い合わせ数に対する各種マーケティング施策の真の効果を推定
📅 Week 2: 潜在的アウトカムとDAG（因果グラフ）の基礎
Week 2 の学習目標: 因果推論の「言語」を習得する
- **潜在的アウトカム**（反実仮想）を理解し、因果効果を数学的に定義できる
- **セレクション・バイアス**の正体を言語化し、なぜ単純比較が誤るかを説明できる
- **DAG（有向非巡回グラフ）**の基本ルールを理解し、因果構造を可視化できる
- **交絡因子**、**媒介変数**、**合流点（コライダー）**をDAG上で識別できる
- **バックドア基準**を使って、調整すべき変数セットを決定できる
📖 月曜・火曜（各3時間）
- 理論（金本 第1章）: 潜在的アウトカム、Rubinの因果モデル（RCM）
- 因果効果の定義: 個別因果効果、**平均因果効果 (ATE)**、処置群の平均因果効果 (ATT)
- 問題の核心: **「反実仮想（もし〜だったら）」**の観測不能性（根本問題）
- RCTと無視可能性: RCTがなぜバイアスを打ち消すのかを数学的に理解する
- セレクション・バイアス: 処置を受ける人と受けない人が系統的に異なる場合の問題
各個体 i に対して:
- Yi(1): 処置を受けた場合のアウトカム（潜在的）
- Yi(0): 処置を受けなかった場合のアウトカム（潜在的）
- 個別因果効果: τi = Yi(1) - Yi(0)
- ATE（平均因果効果）: τ = E[Yi(1) - Yi(0)] = E[Yi(1)] - E[Yi(0)]
根本問題: Yi(1) と Yi(0) を同時に観測することは不可能
💻 Pythonによるシミュレーションデータ作成
import numpy as np
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
import seaborn as sns
np.random.seed(42)
# 1. 交絡因子 Z: ユーザーの習熟度 (正規分布)
n = 1000
Z = np.random.normal(5, 2, n)
# 2. 施策 X: 広告接触 (Zが高い人ほど広告に接触しやすい)
P_X = 1 / (1 + np.exp(-(-1 + 0.3 * Z))) # ロジスティック関数
X = np.random.binomial(1, P_X) # 0: 広告なし, 1: 広告あり
# 3. アウトカム Y: CVR (真の因果効果: 広告でCVRが+3%増加)
Y0 = 0.1 + 0.05 * Z + np.random.normal(0, 0.05, n) # 施策なしのCVR (Zに依存)
Y1 = Y0 + 0.03 # 施策ありのCVR (+0.03が真の因果効果)
Y = Y1 * X + Y0 * (1 - X) # 観測されるCVR (反実仮想)
df = pd.DataFrame({'Z': Z, 'X': X, 'Y': Y, 'Y0': Y0, 'Y1': Y1})
# 真のATEを計算（本来は観測不可能）
true_ate = (Y1 - Y0).mean()
# 単純な比較 (セレクション・バイアスあり)
mean_Y1 = df[df['X']==1]['Y'].mean()
mean_Y0 = df[df['X']==0]['Y'].mean()
naive_ate = mean_Y1 - mean_Y0
print("=" * 60)
print("因果推論の基本問題: 反実仮想の観測不能性")
print("=" * 60)
print(f"真のATE (観測不可能): {true_ate:.4f}")
print(f"単純比較 (ナイーブATE): {naive_ate:.4f}")
print(f"バイアス: {naive_ate - true_ate:.4f}")
print(f"バイアスの割合: {(naive_ate - true_ate) / true_ate * 100:.1f}%")
# 可視化: 交絡因子の分布
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.hist(df[df['X']==0]['Z'], bins=30, alpha=0.6, label='広告なし (X=0)', color='#2196f3')
plt.hist(df[df['X']==1]['Z'], bins=30, alpha=0.6, label='広告あり (X=1)', color='#9c27b0')
plt.xlabel('習熟度 (Z)')
plt.ylabel('頻度')
plt.title('交絡因子の分布の違い\n→ 広告を見る人は元々習熟度が高い', fontweight='bold')
plt.legend()
plt.grid(alpha=0.3)
plt.subplot(1, 2, 2)
plt.scatter(df[df['X']==0]['Z'], df[df['X']==0]['Y'], alpha=0.4, label='広告なし', s=20, color='#2196f3')
plt.scatter(df[df['X']==1]['Z'], df[df['X']==1]['Y'], alpha=0.4, label='広告あり', s=20, color='#9c27b0')
plt.xlabel('習熟度 (Z)')
plt.ylabel('CVR (Y)')
plt.title('習熟度とCVRの関係\n→ 習熟度が両方に影響（交絡）', fontweight='bold')
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()
上のコードでは、単純な平均比較**（ナイーブATE）**は真のATE（0.03）から大きくズレます。これは、**交絡因子Z（習熟度）**が「広告接触（X）」と「CVR（Y）」の両方に影響しているためです。このズレこそが**セレクション・バイアス**です。
- Z → X: 習熟度の高いユーザーは広告をクリックしやすい
- Z → Y: 習熟度の高いユーザーは元々CVRが高い
- 結果: 広告の効果を過大評価してしまう
解決策: 交絡因子Zを調整する手法（回帰分析、傾向スコアなど）を使う
📖 水曜（3時間）- 演習と復習
- 実装演習: 上記のシミュレーションコードを改変し、交絡の強さを変えた場合のバイアスの変化を観察
- サイトへの適用: 「物件問い合わせ」に影響する交絡因子（企業規模、業種、地域など）をリストアップ
- 課題: RCTがなぜバイアスを除去できるかを、Pythonシミュレーションで確認
📖 木曜・金曜（各3時間）
- 理論（金本 第2章/Udemy）: DAG、**交絡因子**、**媒介変数**、**合流点（コライダー）**の役割
- バックドア基準: どの変数を調整すれば交絡が除去できるかを判断するルール
- d-分離: DAG上で因果関係の有無を判定する理論的ツール
- フロントドア基準: 交絡因子が観測不能な場合の代替手法（advanced）
| 構造 | パス | 因果推論での役割 | 調整の要否 |
|---|---|---|---|
| 交絡（Fork） | X ← Z → Y | 偽の相関を生む（バイアスの原因） | 必要（バックドアを閉じる） |
| 媒介（Chain） | X → M → Y | 真の因果が伝わる経路 | 不要（調整すると因果効果を消す） |
| 合流点（Collider） | X → C ← Y | 通常は関係なし | 不要（調整すると偽の相関を生む） |
| 直接因果 | X → Y | 推定したい因果効果 | 調整対象ではない |
処置XからアウトカムYへの因果効果を推定するために、変数セットZを調整すべき条件:
- 条件1: ZはXの子孫（descendant）を含まない
- 条件2: ZはX→Yを含まないすべてのバックドアパスをブロックする
実務での使い方: DAGを描き、Xからバックドアで辿れるすべてのパスを特定 → それらを閉じる最小の変数セットを見つける
import matplotlib.pyplot as plt
from matplotlib.patches import FancyBboxPatch
# DAGの作成と可視化（サイトの広告効果分析の例）
def draw_dag(title):
G = nx.DiGraph()
# ノードの定義
nodes = {
'企業規模': (0, 2),
'業種': (0, 1),
'広告': (2, 2),
'サイト訪問': (4, 2),
'問い合わせ': (6, 2),
'季節性': (3, 0)
}
# エッジの定義（因果関係）
edges = [
('企業規模', '広告'), # 大企業ほど広告予算が大きい
('企業規模', '問い合わせ'), # 大企業ほど問い合わせしやすい（交絡）
('業種', '広告'), # 業種によって広告の選好が異なる
('業種', '問い合わせ'), # 業種によって問い合わせ率が異なる（交絡）
('広告', 'サイト訪問'), # 広告でサイト訪問が増える
('サイト訪問', '問い合わせ'), # 訪問から問い合わせへ（媒介）
('季節性', 'サイト訪問'), # 季節で訪問数が変動
('季節性', '問い合わせ') # 季節で問い合わせ数が変動（交絡）
]
G.add_nodes_from(nodes.keys())
G.add_edges_from(edges)
# 描画
plt.figure(figsize=(12, 8))
ax = plt.gca()
# ノードの色分け
node_colors = []
for node in G.nodes():
if node == '広告':
node_colors.append('#9c27b0') # 処置変数（紫）
elif node == '問い合わせ':
node_colors.append('#f44336') # アウトカム（赤）
elif node == 'サイト訪問':
node_colors.append('#ffc107') # 媒介変数（黄）
else:
node_colors.append('#4caf50') # 交絡因子（緑）
nx.draw(G, nodes,
node_color=node_colors,
node_size=3000,
with_labels=True,
font_size=11,
font_weight='bold',
font_family='Yu Gothic',
arrows=True,
arrowsize=20,
arrowstyle='->',
edge_color='#666',
width=2,
ax=ax)
plt.title(title, fontsize=14, fontweight='bold', pad=20)
# 凡例
legend_elements = [
plt.Line2D([0], [0], marker='o', color='w', markerfacecolor='#9c27b0',
markersize=10, label='処置変数'),
plt.Line2D([0], [0], marker='o', color='w', markerfacecolor='#f44336',
markersize=10, label='アウトカム'),
plt.Line2D([0], [0], marker='o', color='w', markerfacecolor='#ffc107',
markersize=10, label='媒介変数'),
plt.Line2D([0], [0], marker='o', color='w', markerfacecolor='#4caf50',
markersize=10, label='交絡因子')
]
plt.legend(handles=legend_elements, loc='upper left', fontsize=10)
plt.tight_layout()
plt.show()
draw_dag('サイトの広告効果分析DAG\n調整すべき変数: 企業規模、業種、季節性')
- ステップ1: 処置変数（X）とアウトカム（Y）を明確にする
- ステップ2: ドメイン知識で「何が何に影響するか」をリストアップ
- ステップ3: リストを矢印で繋いでDAGを描く
- ステップ4: X→Y以外のパスを特定し、どれがバックドアか確認
- ステップ5: バックドアパスを閉じる最小の変数セットを決定
注意: 媒介変数（X→M→Y）は調整しない！調整すると因果効果が消える
📖 土曜・日曜（演習と復習）
- 実務演習: あなたの会社のA/Bテストデータ（または想定データ）の因果構造をDAGで図示してみる
- 課題1: 「広告費（X）→ 売上（Y）」の間に、「景気（Z）」と「ブランド認知度（M）」を加えたDAGを作成し、調整すべき変数を特定せよ
- 課題2: サイトの「サイトリニューアル → 問い合わせ数」のDAGを作成し、少なくとも5つの変数を含めよ
- Week 2振り返り: 潜在的アウトカム、ATE、DAG、バックドア基準を自分の言葉でまとめる
✅ Week 2 完了チェックリスト
📅 Week 3: 回帰分析による因果効果の推定（Adjustment Formula）
Week 3 の学習目標: 回帰分析を「交絡調整ツール」として使う
- 回帰分析（OLS）が因果推論に適用できる条件を理解する
- **調整公式（Adjustment Formula）**を回帰式として実装する
- カテゴリ変数、連続変数での因果効果推定をマスターする
- 回帰分析の限界（線形性の仮定、外挿の問題）を理解する
- 交互作用項を使った異質性（Heterogeneity）の分析ができる
📖 月曜・火曜（各3時間）
- 理論（金本 第3章）: 条件付き独立性、回帰モデルと層別化の関係
- 回帰式の解釈: Y = β₀ + β₁X + β₂Z + ε における β₁ の意味
- 調整公式: E[Y|do(X=x)] = ∑_z E[Y|X=x, Z=z]P(Z=z) の回帰版
- 標準誤差と信頼区間: 推定値の不確実性を評価する
基本式: Y = β₀ + β₁X + β₂Z + ε
- β₁: 交絡因子Zを調整した後のXのYに対する因果効果（ATE）
- β₂: 交絡因子Zの効果（これ自体は因果効果とは限らない）
- 前提条件:
- 線形性（X、Zと の関係が線形）
- すべての交絡因子がモデルに含まれている
- 誤差項が独立同分布（i.i.d.）
💻 回帰分析による因果効果（ATE）推定のコード
import statsmodels.api as sm
# Week 2で作成した交絡を含むデータ（df）を使用
# 1. ナイーブ回帰モデル (交絡Zを調整しない)
model_naive = smf.ols('Y ~ X', data=df).fit()
print("=" * 70)
print("ナイーブ回帰 (交絡Zを無視)")
print("=" * 70)
print(f"Xの係数 (ナイーブATE): {model_naive.params['X']:.4f}")
print(f"95%信頼区間: [{model_naive.conf_int().loc['X', 0]:.4f}, {model_naive.conf_int().loc['X', 1]:.4f}]")
print(f"p値: {model_naive.pvalues['X']:.4f}")
print(f"真のATE: 0.0300")
print(f"バイアス: {model_naive.params['X'] - 0.03:.4f}")
# 2. 調整回帰モデル (交絡Zを共変量として追加)
model_adj = smf.ols('Y ~ X + Z', data=df).fit()
print("\n" + "=" * 70)
print("調整回帰 (交絡Zを調整)")
print("=" * 70)
print(f"Xの係数 (調整ATE): {model_adj.params['X']:.4f}")
print(f"95%信頼区間: [{model_adj.conf_int().loc['X', 0]:.4f}, {model_adj.conf_int().loc['X', 1]:.4f}]")
print(f"p値: {model_adj.pvalues['X']:.4f}")
print(f"真のATE: 0.0300")
print(f"推定誤差: {model_adj.params['X'] - 0.03:.4f}")
# 3. モデルの詳細サマリー
print("\n" + "=" * 70)
print("調整回帰モデルの詳細")
print("=" * 70)
print(model_adj.summary())
# 4. 可視化: 回帰直線の比較
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# ナイーブモデル
axes[0].scatter(df['X'], df['Y'], alpha=0.5, c='#9c27b0', s=30)
x_vals = np.array([0, 1])
y_vals_naive = model_naive.params['Intercept'] + model_naive.params['X'] * x_vals
axes[0].plot(x_vals, y_vals_naive, 'r-', linewidth=3, label=f'傾き={model_naive.params["X"]:.4f}')
axes[0].set_xlabel('広告接触 (X)')
axes[0].set_ylabel('CVR (Y)')
axes[0].set_title('ナイーブ回帰\n→ バイアスあり', fontweight='bold')
axes[0].legend()
axes[0].grid(alpha=0.3)
# 調整モデル（Zの影響を除去した後）
# 残差プロット: YとXからZの影響を除去
residY = smf.ols('Y ~ Z', data=df).fit().resid
residX = smf.ols('X ~ Z', data=df).fit().resid
axes[1].scatter(residX, residY, alpha=0.5, c='#2196f3', s=30)
model_resid = smf.ols('residY ~ residX', data=pd.DataFrame({'residY': residY, 'residX': residX})).fit()
x_resid = np.linspace(residX.min(), residX.max(), 100)
y_resid = model_resid.params['Intercept'] + model_resid.params['residX'] * x_resid
axes[1].plot(x_resid, y_resid, 'r-', linewidth=3, label=f'傾き={model_adj.params["X"]:.4f}')
axes[1].set_xlabel('X（Z調整後の残差）')
axes[1].set_ylabel('Y（Z調整後の残差）')
axes[1].set_title('調整回帰（Frisch-Waugh-Lovell定理）\n→ 真の因果効果に近い', fontweight='bold')
axes[1].legend()
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
- 線形性: 回帰分析はあくまで**線形モデル**。非線形な関係（例：広告費が一定以上で効果激減）は多項式項や対数変換で対応
- 変数選択: **調整する変数（Z）**を理論（DAG）に基づいて慎重に選ぶこと。たくさん入れれば良いわけではない（Colliderを調整すると逆にバイアス）
- 多重共線性: 交絡因子同士が強く相関している場合、係数が不安定になる → VIF（分散拡大係数）で確認
- 外挿の危険: データの範囲外での予測は信頼できない
- 未観測の交絡: モデルに含まれていない交絡因子があると、依然としてバイアスが残る
- 関数形の誤設定: 真の関係が非線形なのに線形モデルを使うとバイアス
- 外生性の仮定: 誤差項εがXと無相関という仮定が崩れると因果推論は失敗
- 解釈の注意: 回帰係数が「因果効果」と解釈できるのは、DAGで正しく変数を選択した場合のみ
📖 水曜・木曜（各3時間）
- 実装演習（Udemy）: ロジスティック回帰による施策効果の推定、アウトカムが二値変数の場合の対応
- 異質性（Heterogeneity）: 調整因子Zの水準によって効果（Xの係数）が異なるケース（交互作用項の導入）
- カテゴリ変数の扱い: ダミー変数化と解釈
- 診断: 残差プロット、VIF、モデルの適合度（R²、AIC）
# 「習熟度（Z）が高い人ほど、広告の効果が大きい」という仮説を検証
# 交互作用項を含むモデル
model_interaction = smf.ols('Y ~ X + Z + X:Z', data=df).fit()
print("=" * 70)
print("交互作用項を含むモデル（異質性の分析）")
print("=" * 70)
print(model_interaction.summary())
print("\n【解釈】")
print(f"Xの主効果: {model_interaction.params['X']:.4f}")
print(f"X:Zの交互作用: {model_interaction.params['X:Z']:.4f}")
print(f"→ Zが1単位増えると、Xの効果が{model_interaction.params['X:Z']:.4f}だけ変化する")
# 異質性の可視化
z_levels = [3, 5, 7] # 習熟度の低・中・高
colors = ['#e57373', '#9c27b0', '#42a5f5']
plt.figure(figsize=(10, 6))
for z, color in zip(z_levels, colors):
effect_at_z = model_interaction.params['X'] + model_interaction.params['X:Z'] * z
plt.bar(z, effect_at_z, color=color, alpha=0.7, width=0.5,
label=f'Z={z}: 効果={effect_at_z:.4f}')
plt.xlabel('習熟度 (Z)')
plt.ylabel('広告の効果（Xの係数）')
plt.title('異質性の分析: 習熟度別の広告効果\n→ 習熟度が高いほど広告が効く', fontweight='bold')
plt.legend()
plt.grid(alpha=0.3)
plt.show()
✅ Week 3 完了チェックリスト
📅 Week 4: 傾向スコア（Propensity Score）による推定
Week 4 の学習目標: 高次元の交絡を単一変数で調整する
- **傾向スコア（Propensity Score）**の定義と役割を理解する
- ロジスティック回帰で傾向スコアを算出できる
- **傾向スコア・マッチング**と**逆確率重み付け（IPW）**の概念を理解し、実装の準備をする
- **共通サポート**の概念と可視化を理解する
- バランスチェックで調整の妥当性を評価できる
📖 月曜・火曜（各3時間）
- 理論（金本 第4章/Udemy）: 傾向スコアの定義、**強い無視可能性（Unconfoundedness）**、**共通サポート（Common Support）**
- 算出方法: ロジスティック回帰 e(Z) = P(X=1|Z) によるスコア算出
- Rosenbaumの定理: なぜ傾向スコア1つで多次元の交絡を調整できるのか
- バランス診断: 標準化差分（Standardized Difference）の計算
定義: e(Z) = P(X=1|Z) = 観測された共変量Zが与えられたとき、処置を受ける確率
なぜ強力か:
- 交絡因子が10個、20個あっても、傾向スコアという1つの変数に圧縮できる
- 傾向スコアが同じ個体同士を比較すれば、すべての交絡因子が調整される（Rosenbaumの定理）
- 高次元の呪いを回避できる
前提条件:
- 強い無視可能性: すべての交絡因子が観測されている
- 共通サポート: すべての傾向スコアの値で、処置群と対照群の両方が存在する
💻 傾向スコア算出の実装
from sklearn.neighbors import NearestNeighbors
import matplotlib.pyplot as plt
import seaborn as sns
# ロジスティック回帰モデルで傾向スコアを算出
model_ps = LogisticRegression(solver='liblinear', random_state=42)
model_ps.fit(df[['Z']], df['X'])
# 施策を受ける確率 (傾向スコア) を計算
ps_scores = model_ps.predict_proba(df[['Z']])[:, 1]
df['PS'] = ps_scores
print("=" * 70)
print("傾向スコアの算出")
print("=" * 70)
print(f"傾向スコアの範囲: [{ps_scores.min():.4f}, {ps_scores.max():.4f}]")
print(f"処置群の平均PS: {df[df['X']==1]['PS'].mean():.4f}")
print(f"対照群の平均PS: {df[df['X']==0]['PS'].mean():.4f}")
# 傾向スコアの分布を可視化（共通サポートの確認）
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# プロット1: 傾向スコアの分布
axes[0, 0].hist(df[df['X']==1]['PS'], bins=30, alpha=0.6, label='処置群 (X=1)', color='#9c27b0', density=True)
axes[0, 0].hist(df[df['X']==0]['PS'], bins=30, alpha=0.6, label='対照群 (X=0)', color='#2196f3', density=True)
axes[0, 0].axvline(df[df['X']==1]['PS'].min(), color='#9c27b0', linestyle='--', linewidth=2, label='処置群の最小値')
axes[0, 0].axvline(df[df['X']==0]['PS'].max(), color='#2196f3', linestyle='--', linewidth=2, label='対照群の最大値')
axes[0, 0].set_xlabel('傾向スコア (PS)')
axes[0, 0].set_ylabel('密度')
axes[0, 0].set_title('傾向スコアの分布と共通サポート\n→ 重なり部分が共通サポート', fontweight='bold')
axes[0, 0].legend()
axes[0, 0].grid(alpha=0.3)
# プロット2: 習熟度Zと傾向スコアの関係
axes[0, 1].scatter(df[df['X']==0]['Z'], df[df['X']==0]['PS'], alpha=0.4, label='対照群', s=30, color='#2196f3')
axes[0, 1].scatter(df[df['X']==1]['Z'], df[df['X']==1]['PS'], alpha=0.4, label='処置群', s=30, color='#9c27b0')
z_range = np.linspace(df['Z'].min(), df['Z'].max(), 100)
ps_pred = model_ps.predict_proba(z_range.reshape(-1, 1))[:, 1]
axes[0, 1].plot(z_range, ps_pred, 'r-', linewidth=2, label='ロジスティック曲線')
axes[0, 1].set_xlabel('習熟度 (Z)')
axes[0, 1].set_ylabel('傾向スコア (PS)')
axes[0, 1].set_title('交絡因子と傾向スコアの関係\n→ Zが高いほどPSが高い', fontweight='bold')
axes[0, 1].legend()
axes[0, 1].grid(alpha=0.3)
# プロット3: 傾向スコアの箱ひげ図
ps_data = [df[df['X']==0]['PS'], df[df['X']==1]['PS']]
bp = axes[1, 0].boxplot(ps_data, labels=['対照群', '処置群'], patch_artist=True)
bp['boxes'][0].set_facecolor('#2196f3')
bp['boxes'][1].set_facecolor('#9c27b0')
axes[1, 0].set_ylabel('傾向スコア (PS)')
axes[1, 0].set_title('傾向スコアの分布（箱ひげ図）\n→ 処置群の方がPS高い', fontweight='bold')
axes[1, 0].grid(alpha=0.3, axis='y')
# プロット4: 共通サポートの詳細確認
common_support_min = max(df[df['X']==1]['PS'].min(), df[df['X']==0]['PS'].min())
common_support_max = min(df[df['X']==1]['PS'].max(), df[df['X']==0]['PS'].max())
in_support = (df['PS'] >= common_support_min) & (df['PS'] <= common_support_max)
axes[1, 1].hist(df[in_support & (df['X']==1)]['PS'], bins=20, alpha=0.6,
label=f'処置群（共通サポート内: {in_support[df["X"]==1].sum()}人）', color='#9c27b0')
axes[1, 1].hist(df[in_support & (df['X']==0)]['PS'], bins=20, alpha=0.6,
label=f'対照群（共通サポート内: {in_support[df["X"]==0].sum()}人）', color='#2196f3')
axes[1, 1].axvspan(common_support_min, common_support_max, alpha=0.2, color='green', label='共通サポート範囲')
axes[1, 1].set_xlabel('傾向スコア (PS)')
axes[1, 1].set_ylabel('頻度')
axes[1, 1].set_title(f'共通サポート: [{common_support_min:.3f}, {common_support_max:.3f}]\n除外: {(~in_support).sum()}人',
fontweight='bold')
axes[1, 1].legend()
axes[1, 1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
print(f"\n共通サポート範囲: [{common_support_min:.4f}, {common_support_max:.4f}]")
print(f"共通サポート内のサンプル数: {in_support.sum()}人 ({in_support.mean()*100:.1f}%)")
print(f"共通サポート外のサンプル数: {(~in_support).sum()}人 ({(~in_support).mean()*100:.1f}%)")
共通サポート（ヒストグラムが重なる範囲）がない場合、その範囲のユーザーは**比較不可能**です。例えば：
- 処置群にPS=0.9のユーザーがいるが、対照群にPS=0.9のユーザーがいない場合、そのユーザーは除外（または重み調整）する必要がある
- 共通サポート外のサンプルを強引にマッチングすると、バイアスが大きくなる
- 除外されたサンプルの特徴を必ず確認し、結果の一般化可能性を議論すること
傾向スコアによる調整が成功しているかを確認する指標
標準化差分 = (処置群の平均 - 対照群の平均) / √[(処置群の分散 + 対照群の分散) / 2]
判断基準:
- |標準化差分| < 0.1: バランス良好
- 0.1 ≤ |標準化差分| < 0.25: やや不均衡
- |標準化差分| ≥ 0.25: バランス不良（調整方法の見直しが必要）
def standardized_difference(df, var, treatment_col='X'):
"""標準化差分を計算"""
treated = df[df[treatment_col]==1][var]
control = df[df[treatment_col]==0][var]
mean_diff = treated.mean() - control.mean()
pooled_std = np.sqrt((treated.var() + control.var()) / 2)
return mean_diff / pooled_std if pooled_std > 0 else 0
# 調整前のバランス
std_diff_before = standardized_difference(df, 'Z')
print("=" * 70)
print("バランス診断（傾向スコア調整前）")
print("=" * 70)
print(f"習熟度Zの標準化差分: {std_diff_before:.4f}")
if abs(std_diff_before) < 0.1:
print("→ バランス良好")
elif abs(std_diff_before) < 0.25:
print("→ やや不均衡")
else:
print("→ バランス不良（調整が必要）")
# 傾向スコアで層別化した後のバランス確認（簡易版）
# 傾向スコアを5つの層に分割
df['PS_quintile'] = pd.qcut(df['PS'], q=5, labels=False, duplicates='drop')
print("\n" + "=" * 70)
print("傾向スコア層別化後のバランス")
print("=" * 70)
for q in sorted(df['PS_quintile'].unique()):
df_q = df[df['PS_quintile'] == q]
if len(df_q[df_q['X']==1]) > 0 and len(df_q[df_q['X']==0]) > 0:
std_diff_q = standardized_difference(df_q, 'Z')
print(f"第{q+1}層（PS: {df_q['PS'].min():.3f}-{df_q['PS'].max():.3f}）")
print(f" 標準化差分: {std_diff_q:.4f}")
print(f" 処置群: {len(df_q[df_q['X']==1])}人, 対照群: {len(df_q[df_q['X']==0])}人")
📖 水曜・木曜（各3時間）
- 理論/実装（Udemy）: **傾向スコア・マッチング**の原理（Nearest Neighbor Matchingなど）
- マッチング手法: 1対1マッチング、1対多マッチング、カリパーマッチング
- IPWの概念: **逆確率重み付け（Inverse Probability Weighting）**によるATE推定の理論
- IPW実装準備: 重み W の計算 W = X/PS + (1-X)/(1-PS) の理解
- 層別化: 傾向スコアで層を作り、各層内でATE推定
| 手法 | 原理 | 長所 | 短所 |
|---|---|---|---|
| マッチング | PSが近い対照群を見つけてペアを作る | 直感的でわかりやすい 外れ値に頑健 |
サンプルを捨てる可能性 マッチングの質に依存 |
| 層別化 | PSで層を作り、各層内で比較 | シンプル すべてのサンプルを使用 |
層の数の選択が必要 層内の不均衡が残る可能性 |
| IPW | PSの逆数で重み付け | すべてのサンプルを使用 理論的に最も効率的 |
極端なPSで重みが不安定 実装がやや複雑 |
from sklearn.neighbors import NearestNeighbors
# 共通サポート内のデータのみ使用
df_support = df[in_support].copy()
# 処置群と対照群に分割
treated = df_support[df_support['X'] == 1].copy()
control = df_support[df_support['X'] == 0].copy()
print("=" * 70)
print("傾向スコア・マッチング（1対1 Nearest Neighbor）")
print("=" * 70)
print(f"処置群: {len(treated)}人")
print(f"対照群: {len(control)}人")
# Nearest Neighbor探索（カリパー付き）
caliper = 0.2 * df_support['PS'].std() # カリパー = PSの標準偏差の0.2倍
print(f"カリパー: {caliper:.4f}")
nn = NearestNeighbors(n_neighbors=1, metric='euclidean')
nn.fit(control[['PS']].values)
distances, indices = nn.kneighbors(treated[['PS']].values)
# カリパー内のマッチのみ採用
valid_matches = distances.flatten() < caliper
matched_treated_indices = treated.index[valid_matches]
matched_control_indices = control.index[indices.flatten()[valid_matches]]
print(f"\nマッチング成功: {valid_matches.sum()}ペア")
print(f"マッチング失敗（カリパー外）: {(~valid_matches).sum()}人")
# マッチング後のデータセット
df_matched = pd.concat([
df_support.loc[matched_treated_indices],
df_support.loc[matched_control_indices]
])
# マッチング後のバランス確認
std_diff_after = standardized_difference(df_matched, 'Z')
print("\n" + "=" * 70)
print("マッチング後のバランス診断")
print("=" * 70)
print(f"調整前の標準化差分: {std_diff_before:.4f}")
print(f"調整後の標準化差分: {std_diff_after:.4f}")
print(f"改善: {abs(std_diff_before) - abs(std_diff_after):.4f}")
# マッチング後のATE推定
ate_matched = (
df_matched[df_matched['X']==1]['Y'].mean() -
df_matched[df_matched['X']==0]['Y'].mean()
)
print("\n" + "=" * 70)
print("因果効果の推定")
print("=" * 70)
print(f"真のATE: 0.0300")
print(f"マッチング後のATE: {ate_matched:.4f}")
print(f"推定誤差: {ate_matched - 0.03:.4f}")
# 可視化: マッチング前後の比較
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# マッチング前
axes[0].scatter(df_support[df_support['X']==0]['PS'],
df_support[df_support['X']==0]['Z'],
alpha=0.5, label='対照群', s=30, color='#2196f3')
axes[0].scatter(df_support[df_support['X']==1]['PS'],
df_support[df_support['X']==1]['Z'],
alpha=0.5, label='処置群', s=30, color='#9c27b0')
axes[0].set_xlabel('傾向スコア (PS)')
axes[0].set_ylabel('習熟度 (Z)')
axes[0].set_title(f'マッチング前\n標準化差分={std_diff_before:.3f}', fontweight='bold')
axes[0].legend()
axes[0].grid(alpha=0.3)
# マッチング後
axes[1].scatter(df_matched[df_matched['X']==0]['PS'],
df_matched[df_matched['X']==0]['Z'],
alpha=0.5, label='対照群（マッチング後）', s=30, color='#2196f3')
axes[1].scatter(df_matched[df_matched['X']==1]['PS'],
df_matched[df_matched['X']==1]['Z'],
alpha=0.5, label='処置群（マッチング後）', s=30, color='#9c27b0')
axes[1].set_xlabel('傾向スコア (PS)')
axes[1].set_ylabel('習熟度 (Z)')
axes[1].set_title(f'マッチング後\n標準化差分={std_diff_after:.3f}', fontweight='bold')
axes[1].legend()
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
- カリパーの設定: PSの標準偏差の0.1〜0.25倍が一般的。厳しすぎるとマッチング数が減る
- 1対多マッチング: 処置群が少ない場合、1人に対して複数の対照群をマッチングする
- 置換の有無: 同じ対照群を複数回使うか（置換あり）、1回のみか（置換なし）
- マッチング後のバランスチェック必須: 標準化差分で必ず確認
✅ Week 4 完了チェックリスト
📅 Week 5: 準実験デザインの概要とPhase 1総復習
Week 5 の学習目標: 観測データから因果を炙り出すテクニックを知る
- **差分の差分法（DiD）**の基本原理と適用例を理解する
- **回帰不連続デザイン（RDD）**の基本原理を理解する
- DiDをPythonで実装し、因果効果を推定できる
- Phase 1の全ての概念を整理し、Phase 2への学習準備を完了する
- 実務でどの手法を選ぶべきかの判断基準を持つ
📖 月曜・火曜（各3時間）
- DiD（金本 第5章/Udemy）: DiDの数式と「平行トレンドの仮定」の重要性
- WEBマーケでのDiD例: 「特定地域へのテレビ広告投入」の前後で、広告を打たなかった地域と比較する
- RDD: 「LTV上位10%に達したユーザーへの特別施策」のように、**閾値**を利用するケース
- 準実験デザイン: RCTができない場合の代替手法群
基本アイデア: 処置群の前後変化 - 対照群の前後変化 = 因果効果
DiD = [E[Y|処置群,後] - E[Y|処置群,前]] - [E[Y|対照群,後] - E[Y|対照群,前]]
なぜ強力か:
- 時間不変の交絡因子（処置群と対照群の元々の差）が自動的に除去される
- 共通トレンド（両群に共通の時間変動）も除去される
- 残るのは処置の純粋な効果のみ
重要な仮定:
- 平行トレンドの仮定: 処置がなければ、両群のトレンドは平行だったはず
- この仮定が崩れるとDiDは失敗する
💻 DiDモデルの実装
import pandas as pd
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
# DiDシミュレーションデータ作成
# サイト: 東京エリア（処置群）でサイトリニューアル、大阪エリア（対照群）は変更なし
np.random.seed(42)
N = 400
df_did = pd.DataFrame({
'Group': np.repeat([0, 1], N/2), # 0: 対照（大阪）, 1: 処置（東京）
'Time': np.tile(np.repeat([0, 1], N/4), 2), # 0: リニューアル前, 1: リニューアル後
'Y': np.zeros(N)
})
# データ生成（真のDiD効果=10件）
df_did['Y'] = (
50 # ベースライン（問い合わせ数）
+ 5 * df_did['Group'] # グループ差（東京は元々5件多い）
+ 8 * df_did['Time'] # 時間トレンド（全体で8件増加）
+ 10 * df_did['Group'] * df_did['Time'] # 真のDiD効果 (10件)
+ np.random.normal(0, 3, N) # ランダムノイズ
)
print("=" * 70)
print("DiD（差分の差分法）分析: サイトサイトリニューアル効果")
print("=" * 70)
# 各グループ・期間の平均を計算
means = df_did.groupby(['Group', 'Time'])['Y'].mean()
print("\n【各グループ・期間の平均問い合わせ数】")
print(f"大阪（対照群）・リニューアル前: {means[0, 0]:.2f}件")
print(f"大阪（対照群）・リニューアル後: {means[0, 1]:.2f}件")
print(f"東京（処置群）・リニューアル前: {means[1, 0]:.2f}件")
print(f"東京（処置群）・リニューアル後: {means[1, 1]:.2f}件")
# 第1の差分（各グループの前後比較）
diff_control = means[0, 1] - means[0, 0]
diff_treated = means[1, 1] - means[1, 0]
print(f"\n【第1の差分（前後比較）】")
print(f"大阪の変化: {diff_control:.2f}件")
print(f"東京の変化: {diff_treated:.2f}件")
# 第2の差分（DiD効果）
did_effect = diff_treated - diff_control
print(f"\n【DiD効果（差分の差分）】")
print(f"DiD = {diff_treated:.2f} - {diff_control:.2f} = {did_effect:.2f}件")
print(f"真の因果効果: 10.00件")
print(f"推定誤差: {did_effect - 10:.2f}件")
# DiD回帰モデル: Y = β0 + β1*Group + β2*Time + β3*(Group*Time) + ε
# β3がDiD効果に相当する
model_did = smf.ols('Y ~ Group * Time', data=df_did).fit()
print("\n" + "=" * 70)
print("DiD回帰モデルによる推定")
print("=" * 70)
print(model_did.summary())
print(f"\n【回帰モデルでの推定結果】")
print(f"DiD効果（Group:Timeの係数）: {model_did.params['Group:Time']:.4f}")
print(f"95%信頼区間: [{model_did.conf_int().loc['Group:Time', 0]:.4f}, {model_did.conf_int().loc['Group:Time', 1]:.4f}]")
print(f"p値: {model_did.pvalues['Group:Time']:.4f}")
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# プロット1: 時系列推移
for group, label, color in [(0, '大阪（対照群）', '#2196f3'), (1, '東京（処置群）', '#9c27b0')]:
df_group = df_did[df_did['Group'] == group]
group_means = df_group.groupby('Time')['Y'].mean()
axes[0].plot([0, 1], group_means.values, 'o-', linewidth=3, markersize=10,
label=label, color=color)
# 反実仮想（東京がリニューアルしなかった場合）
if group == 1:
counterfactual= [means[1, 0], means[1, 0] + diff_control]
axes[0].plot([0, 1], counterfactual, 'o--', linewidth=2, markersize=8,
alpha=0.5, color=color, label='東京（反実仮想）')
axes[0].axvline(x=0.5, color='red', linestyle='--', linewidth=2, alpha=0.7, label='リニューアル実施')
axes[0].set_xticks([0, 1])
axes[0].set_xticklabels(['リニューアル前', 'リニューアル後'])
axes[0].set_ylabel('問い合わせ数（件）')
axes[0].set_title('DiDの可視化\n→ 反実仮想との差がDiD効果', fontweight='bold')
axes[0].legend()
axes[0].grid(alpha=0.3)
# プロット2: DiD効果の分解
categories = ['大阪\n（対照群）', '東京\n（処置群）']
before_means = [means[0, 0], means[1, 0]]
after_means = [means[0, 1], means[1, 1]]
x = np.array([0, 1])
width = 0.35
bars1 = axes[1].bar(x - width/2, before_means, width, label='リニューアル前',
color='lightblue', alpha=0.8)
bars2 = axes[1].bar(x + width/2, after_means, width, label='リニューアル後',
color='orange', alpha=0.8)
# 変化量を矢印で表示
for i, (before, after, diff) in enumerate([(means[0,0], means[0,1], diff_control),
(means[1,0], means[1,1], diff_treated)]):
axes[1].annotate('', xy=(i, after), xytext=(i, before),
arrowprops=dict(arrowstyle='->', lw=2, color='black'))
axes[1].text(i + 0.15, (before + after)/2, f'+{diff:.1f}',
fontsize=11, fontweight='bold', color='red' if i==1 else 'blue')
axes[1].set_ylabel('問い合わせ数（件）')
axes[1].set_xticks(x)
axes[1].set_xticklabels(categories)
axes[1].set_title(f'DiD効果の分解\nDiD = {diff_treated:.1f} - {diff_control:.1f} = {did_effect:.1f}件',
fontweight='bold')
axes[1].legend()
axes[1].grid(axis='y', alpha=0.3)
plt.tight_layout()
plt.show()
DiDが機能するためには、**「もし処置がなければ、処置群と対照群のアウトカムのトレンドは平行だった」**という仮定が必要です。
検証方法:
- プラセボテスト: 処置前の複数期間でDiDを実施し、効果がゼロかを確認
- グラフ検証: 処置前の時系列グラフで、両群のトレンドが平行かを目視確認
- イベントスタディ: 処置前後の各期間での効果を推定し、処置前はゼロかを確認
サイトでの例: 東京と大阪で、リニューアル前3ヶ月間の問い合わせ数のトレンドが平行かを確認
- 予期された処置: リニューアル予告で事前に行動変化 → 処置前から差が出る
- 同時期の他の変化: 処置群だけに別のイベントが発生 → DiD効果に混入
- 構成変化: 処置群のユーザー構成が処置後に変化 → 平行トレンドが崩れる
- スピルオーバー: 対照群が処置の影響を受ける → 純粋な対照群ではない
📖 水曜（3時間）- RDDの基礎
- RDD（金本 第6章/Udemy）: 回帰不連続デザインの原理
- 基本アイデア: 閾値の前後で処置が決まる場合、閾値付近で比較すれば因果効果が推定できる
- Sharp RDD vs Fuzzy RDD: 閾値で処置が完全に決まるか、確率的に決まるか
- サイトでの応用例: 「問い合わせ数100件以上の企業に特別対応」の効果測定
基本アイデア: 何らかのスコアに閾値があり、閾値を超えると処置が決まる場合、閾値前後で比較
例（サイト）:
- 企業規模が従業員100人以上 → プレミアム対応
- 閾値付近（95〜105人）の企業を比較すれば、ほぼランダム化と同じ
重要な仮定:
- 連続性: 閾値付近で潜在的アウトカムが連続（ジャンプしない）
- 操作不可能性: 企業が閾値を超えるために従業員数を操作していない
np.random.seed(42)
n_rdd = 500
# Running Variable: 企業規模（従業員数）
company_size = np.random.uniform(50, 150, n_rdd)
# 閾値: 100人
threshold = 100
# 処置: 100人以上でプレミアム対応
treatment_rdd = (company_size >= threshold).astype(int)
# アウトカム: 成約数（閾値で10件ジャンプ = RDD効果）
outcome_rdd = (
50 + # ベースライン
0.3 * (company_size - threshold) + # 企業規模の線形効果
10 * treatment_rdd + # 真のRDD効果（10件）
np.random.normal(0, 3, n_rdd) # ノイズ
)
df_rdd = pd.DataFrame({
'CompanySize': company_size,
'Treatment': treatment_rdd,
'Outcome': outcome_rdd,
'Distance': company_size - threshold # 閾値からの距離
})
print("=" * 70)
print("RDD（回帰不連続デザイン）分析")
print("=" * 70)
# バンド幅内（閾値±20人）でのRDD推定
bandwidth = 20
df_rdd_bw = df_rdd[abs(df_rdd['Distance']) <= bandwidth].copy()
print(f"\n使用サンプル: 全体{len(df_rdd)}社 → バンド幅内{len(df_rdd_bw)}社")
print(f"閾値: 従業員{threshold}人")
print(f"バンド幅: ±{bandwidth}人（{threshold-bandwidth}〜{threshold+bandwidth}人）")
# 閾値直前と直後の平均比較
mean_below = df_rdd_bw[df_rdd_bw['Treatment']==0]['Outcome'].mean()
mean_above = df_rdd_bw[df_rdd_bw['Treatment']==1]['Outcome'].mean()
rdd_effect_naive = mean_above - mean_below
print(f"\n【単純比較】")
print(f"閾値未満（<100人）の平均成約数: {mean_below:.2f}件")
print(f"閾値以上（≥100人）の平均成約数: {mean_above:.2f}件")
print(f"差（RDD効果）: {rdd_effect_naive:.2f}件")
# 局所線形回帰によるRDD推定
model_rdd = smf.ols('Outcome ~ Treatment + Distance + Treatment:Distance',
data=df_rdd_bw).fit()
print(f"\n【局所線形回帰】")
print(f"RDD効果（Treatmentの係数）: {model_rdd.params['Treatment']:.4f}")
print(f"95%信頼区間: [{model_rdd.conf_int().loc['Treatment', 0]:.4f}, "
f"{model_rdd.conf_int().loc['Treatment', 1]:.4f}]")
print(f"p値: {model_rdd.pvalues['Treatment']:.4f}")
print(f"真の因果効果: 10.00件")
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# プロット1: RDDの散布図
axes[0].scatter(df_rdd[df_rdd['Treatment']==0]['CompanySize'],
df_rdd[df_rdd['Treatment']==0]['Outcome'],
alpha=0.4, s=30, color='#2196f3', label='処置なし（<100人）')
axes[0].scatter(df_rdd[df_rdd['Treatment']==1]['CompanySize'],
df_rdd[df_rdd['Treatment']==1]['Outcome'],
alpha=0.4, s=30, color='#9c27b0', label='処置あり（≥100人）')
# 回帰直線（バンド幅内）
x_below = np.linspace(threshold - bandwidth, threshold, 100)
x_above = np.linspace(threshold, threshold + bandwidth, 100)
y_below = (model_rdd.params['Intercept'] +
model_rdd.params['Distance'] * (x_below - threshold))
y_above = (model_rdd.params['Intercept'] +
model_rdd.params['Treatment'] +
(model_rdd.params['Distance'] + model_rdd.params['Treatment:Distance']) *
(x_above - threshold))
axes[0].plot(x_below, y_below, 'r-', linewidth=3, label='回帰直線（閾値未満）')
axes[0].plot(x_above, y_above, 'r-', linewidth=3, label='回帰直線（閾値以上）')
axes[0].axvline(x=threshold, color='green', linestyle='--', linewidth=2, label='閾値')
axes[0].axvspan(threshold - bandwidth, threshold + bandwidth, alpha=0.1, color='yellow')
axes[0].set_xlabel('企業規模（従業員数）')
axes[0].set_ylabel('成約数（件）')
axes[0].set_title('RDD: 閾値でのジャンプ\n→ ジャンプ幅=処置効果', fontweight='bold')
axes[0].legend()
axes[0].grid(alpha=0.3)
# プロット2: バンド幅内の拡大図
axes[1].scatter(df_rdd_bw[df_rdd_bw['Treatment']==0]['CompanySize'],
df_rdd_bw[df_rdd_bw['Treatment']==0]['Outcome'],
alpha=0.6, s=50, color='#2196f3', label='処置なし')
axes[1].scatter(df_rdd_bw[df_rdd_bw['Treatment']==1]['CompanySize'],
df_rdd_bw[df_rdd_bw['Treatment']==1]['Outcome'],
alpha=0.6, s=50, color='#9c27b0', label='処置あり')
axes[1].plot(x_below, y_below, 'r-', linewidth=3)
axes[1].plot(x_above, y_above, 'r-', linewidth=3)
axes[1].axvline(x=threshold, color='green', linestyle='--', linewidth=2)
# RDD効果を矢印で表示
y_at_threshold_below = model_rdd.params['Intercept']
y_at_threshold_above = model_rdd.params['Intercept'] + model_rdd.params['Treatment']
axes[1].annotate('', xy=(threshold, y_at_threshold_above),
xytext=(threshold, y_at_threshold_below),
arrowprops=dict(arrowstyle='<->', lw=3, color='red'))
axes[1].text(threshold + 2, (y_at_threshold_below + y_at_threshold_above)/2,
f'RDD効果\n{model_rdd.params["Treatment"]:.1f}件',
fontsize=11, fontweight='bold', color='red')
axes[1].set_xlabel('企業規模（従業員数）')
axes[1].set_ylabel('成約数（件）')
axes[1].set_title(f'バンド幅内の拡大図（±{bandwidth}人）\n→ 閾値付近はほぼランダム', fontweight='bold')
axes[1].legend()
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
📖 木曜・金曜（各3時間）
- 総復習: 潜在的アウトカム、DAG、回帰分析、傾向スコア、DiD、RDDの全てを整理
- 手法の選択基準: どの状況でどの手法を使うべきかのフローチャート作成
- Phase 2の予習: ベイズ統計学で因果推論をどう表現するかを軽く確認
- 実務課題: サイトの実データ（または仮想データ）で複数手法を試す
| 状況 | 推奨手法 | 理由 |
|---|---|---|
| ランダム化可能 | RCT（A/Bテスト） | 最も信頼性が高い（ゴールドスタンダード） |
| すべての交絡因子が観測済み | 回帰分析 or 傾向スコア | 交絡を調整できる。高次元なら傾向スコア |
| 処置前後データ + 対照群あり | DiD | 時間不変の交絡を自動除去。平行トレンド要確認 |
| 閾値で処置が決まる | RDD | 閾値付近はランダム化に近い。操作不可要確認 |
| 未観測の交絡あり + 操作変数あり | IV（操作変数法） | 未観測交絡に対応可能。Phase 1では詳細未学習 |
| 対照群がない | 合成コントロール法 or 中断時系列 | 複数の対照候補から合成。Phase 1では詳細未学習 |
| 手法 | 主な仮定 | 長所 | 短所 |
|---|---|---|---|
| 回帰分析 | ・すべての交絡観測済み ・線形性 |
シンプル 解釈が容易 |
未観測交絡に弱い 非線形に対応困難 |
| 傾向スコア | ・すべての交絡観測済み ・共通サポート |
高次元交絡に対応 バランスが視覚的 |
未観測交絡に弱い モデル依存 |
| DiD | ・平行トレンド ・SUTVA |
時間不変交絡を除去 直感的 |
平行トレンド検証必須 処置前後データ必要 |
| RDD | ・閾値での連続性 ・操作不可能性 |
未観測交絡に頑健 説得力がある |
閾値付近のみ推定 外的妥当性が限定的 |
📖 土曜・日曜（総合演習とまとめ）
- 総合演習: サイトの架空データで、4つの手法（回帰、PS、DiD、RDD）を実装し、結果を比較
- レポート作成: 経営層向けに「Phase 1で習得したスキル」をまとめる
- Phase 2への準備: ベイズ統計学の基礎（事前分布、事後分布）を予習
- 振り返り: つまずいたポイント、よく理解できた概念を整理
✅ Week 5 完了チェックリスト
🎓 Phase 1 完了チェックリスト総まとめ
🎉 Phase 1のまとめと次のステップ
この4週間で、あなたは「相関」と「因果」の壁を完全に越えました。単なるデータ分析者から、**因果構造を設計できる戦略家**へと進化しました。
- ✅ 因果効果を定義する**潜在的アウトカム**の思考回路
- ✅ 交絡を特定する**DAG**という最強のツール
- ✅ 交絡を調整する**回帰分析・傾向スコア**という実用的な武器
- ✅ 時間と閾値を活用する**DiD・RDD**という準実験デザイン
- ✅ 実務でどの手法を選ぶべきかの**判断力**
次のPhase 2では、**ベイズ統計学**を導入します。Phase 1で学んだ手法は「点推定」（効果がX%）でしたが、ベイズ統計学では効果を「確率分布」で推定し、**不確実性**を考慮した意思決定が可能になります。
- 課題1: 回帰分析の係数（ATE）の**信頼区間**を、ベイズ統計でより柔軟に表現する
- 課題2: 処置効果がゼロでない**確率**を算出することで、意思決定の確度を高める
- 課題3: 事前知識（ドメイン知識）を**事前分布**として組み込む
Phase 1の統計学的基礎が、Phase 2での柔軟なベイズモデリングを可能にします。
- 未観測の交絡: すべての手法が「観測された交絡因子が全て」という前提に依存
- 外的妥当性: 特にRDDやマッチングは、一部のサンプルでの推定 → 全体への一般化に注意
- 仮定の検証: 平行トレンド、共通サポート、操作不可能性などの仮定を必ず確認
- 因果推論は万能ではない: 複雑な因果構造、フィードバックループ、動的な関係には限界がある
次のステップ: Phase 6で学ぶDoWhyなどのツールで、仮定の検証と感度分析を強化