因果推論-2-基本実装 ・ 学部2年 ・ 中島智哉 

目的 実践的な因果推論の知識・スキルを身につける。 今回は基本的な因果推論手法を実装する。 . 

参考書籍 因果推論: 基礎から機械学習・時 系列解析・因果探索を用いた意思 決定のアプローチ 著者：金本 拓 出版社：オーム社 ISBN：978-4-274-23123-0 発行：2024年 

分析の全体像 通常は以下のような意思決定フローを辿る。 . 

因果効果を推定するための手法の選定にあたっては、課題の前提条件によって通常以下のようなフローを辿る。 

二重にロバストな推定法まで 

回帰分析(構造方程式モデリング, SEM) 変数間の関係を線形とみなせる場合、重回帰分析により共変量を揃えることで因果効果を推定することができる。 構造方程式モデリングは因果関係を統計的に検証するツール。全ては理解していないが、今回は書籍通り実装してみる。 . 問い 炭酸飲料の売上に与える因子は何か？また、その因子がどの程度影響を与えるか？ . . 

データ収集／確認 今回はダミーデータを用意する . !pip install japanize-matplotlib !pip install semopy==2.3.11 import numpy as np import pandas as pd import matplotlib.pyplot as plt import japanize_matplotlib import seaborn as sns import semopy as sem np.random.seed(0) N=200 TV_CM = np.random.normal(50, 10, N) Temperature_C = np.random.uniform(10, 35, N) Noise = np.random.normal(0, 5, N) Soda_Sales = 10 + 0.5 * TV_CM + 2.0 * Temperature_C + Noise 

data = pd.DataFrame({ 'TV_CM': TV_CM, 'Temperature_C': Temperature_C, 'Soda_Sales': Soda_Sales }) sns.pairplot(data) plt.suptitle('変数間の関係性', y=1.02) plt.show() 

DAG の確認 炭酸飲料の売上 気温 TVCM T X Y . 識別仮定の確認 ここではデータの可視化や背景の理解から識別仮定を確認できたとする。 . 

因果効果の推定 可視化の結果から、線形性が仮定でき、多重共線性もないと判断した。 Name of objective: MLW Optimization method: SLSQP Optimization successful. Optimization terminated successfully Objective value: 0.000 Number of iterations: 12 Params: 0.513 2.043 26.145 . model_desc = """ Soda_Sales ~ TV_CM + Temperature_C """ mod = sem.Model(model_desc) res = mod.fit(data) print(res) 

結果の評価 Value DoF 3.000000e+00 DoF Baseline 5.000000e+00 chi2 1.050434e-04 chi2 p-value 9.999997e-01 chi2 Baseline 4.891441e+02 CFI 1.006196e+00 GFI 9.999998e-01 AGFI 9.999996e-01 NFI 9.999998e-01 TLI 1.010327e+00 RMSEA 0.000000e+00 AIC 5.999999e+00 BIC 1.589495e+01 LogLik 5.252169e-07 より当てはまりの良いモデルと判断した。閾値の設定には本来はドメイン知識が必要。 stats = sem.calc_stats(mod) print(stats.T) RMSEA < 0.05, GFI > 0.9 

lval op rval Estimate Std. Err z-value p-value 0 Soda_Sales ~ TV_CM 0.513118 0.035408 14.491412 0.0 1 Soda_Sales ~ Temperature_C 2.043452 0.047345 43.160776 0.0 2 Soda_Sales ~~ Soda_Sales 26.145319 2.614532 10.000000 0.0 どちらの変数についても より統計的に有意と判断した。 同様に閾値の設定には本来はドメイン知識が必要。 TVCM の指標が 1 単位増えると、炭酸飲料の売上が約 51 万円増加することが読み取れる 気温が 1℃ 上がると、炭酸飲料の売上が約 204 万円増加することが読み取れる ins = mod.inspect() print(ins) p − value < 0.05 

回帰分析(構造方程式モデリング, SEM)の課題 共変量が多い場合や交互作用項を考慮する場合には、多重共線性が発生しやすくなる 共変量が膨大な場合、任意の共変量で条件づけた際にその条件で正値性を満たしていることの確認が困難 未観測の交絡因子が存在する場合は正しい因果効果の推定ができない 

傾向スコア分析 共変量を調整して比較対象同士をより近いものにする。 → RCT を行った結果に近づけるため 傾向スコア は、未観測の交絡因子が存在しない場合において、交絡因子 で条件づけた処置確率（処置群に属する確率） として定義される。 たとえば、新薬の臨床試験において交絡因子 が与えられたときに新薬が与えられる（ ）確率が傾向スコア に 該当する。 このとき、傾向スコアは を目的変数とした回帰モデルを用いて計算される。傾向スコアを使用することで、交絡因子を 1 次元の 指標に変換できる。 . π C 0 < π(C) < 1 C T = 1 π(C) π(C) = P(T = 1∣C) T 

理想的には交絡因子が特定されており、その交絡因子のみを調整すれば十分だが、実態としては交絡因子の特定は難しい。 そのため、共変量を用いて傾向スコアを算出する。 . 問い . データ収集／確認 禁煙が体重に与える影響はどの程度か？ !pip install causallib from causallib.datasets import load_nhefs from causallib.estimation import IPW, PropensityMatching, StratifiedStandardization from causallib.evaluation import evaluate import matplotlib.pyplot as plt import numpy as np from sklearn.linear_model import LogisticRegression data = load_nhefs() 

age ：年齢（1971 年時点） race ：人種　 0：白人、1：黒人 sex ：性別　 0：男性、1：女性 smokeintensity ：1 日の喫煙本数（1971 年時点） smokeyrs ：喫煙歴、単位：年 wt71 ：体重、単位：kg（1971 年時点） 【data.X】 active ：日常生活でどの程度活発に活動するか 　　 0：非常に活発、1：適度に活発、2：活発でない exercise ：休みの日の運動量（1971 年時点） 　　 0：多くの運動、1：適度な運動、2：ほとんど運動しない education ：学歴（1971 年時点） 　　 1：8 年生以下、2：高校中退、3：高校卒業、4：大学中退、5：大 学卒業以上 qsmk ：処置変数、1971 年から 1982 年の間、禁煙をしたか 　　 1：はい、0：いいえ【data.a】 wt82_71 ：結果変数、体重の変動（1982 年の体重 − 1971 年の体重）、単位：kg【data.a】 

DAG の確認 体重 喫煙 年齢など X T Y . 識別仮定の確認 この段階では SUTVA のみ確認する。識別性仮定は傾向スコアで共変量を調整した後に検討する。 ここでは SUTVA を満たすと仮定する。 . 

注意点 ドメイン知識から重要と思われる共変量や結果変数と相関する共変量を調整対象に選ぶ。 もともとバランスが取れていた共変量に傾向スコアを導入すると却って群間のバランスが悪くなることがあるため、モデルに含め る共変量は慎重に選ぶ必要がある。 傾向スコアモデルが線形モデルである場合には多重共線性が発生し得る。GBDT の場合はそれほど問題にならない。 . 因果効果の推定 傾向スコアの推定後、これを用いて共変量を調整し因果効果を推定する。 傾向スコアの活用方法には、マッチング 、層別解析、重み付け法などがある。 . 

マッチング 処置群と対照群において傾向スコア(≒ 共変量)の値が近い者同士を組み合わせること。 以下は最近傍法で個体間でマッチングさせた例である。 Average outcome of Untreated: 1.726 Average outcome of Treated: 4.683 Average Treatment Effect: 2.956 learner = LogisticRegression( solver="liblinear", class_weight="balanced" ) pm = PropensityMatching(learner=learner) pm.fit(data.X, data.a, data.y) outcomes = pm.estimate_population_outcome(data.X, data.a) effect = pm.estimate_effect(outcomes[1],outcomes[0]) print(f"Average outcome of Untreated: {outcomes[0]:.3f}") print(f"Average outcome of Treated: {outcomes[1]:.3f}") print(f'Average Treatment Effect: {effect["diff"]:.3f}') 

出力結果より、 喫煙を続けた場合（出力 1 行目）の平均的な体重増加の推定値は 1.726 禁煙した場合（2 行目）の平均的な体重増加の推定値は 4.683 傾向スコアマッチングにより推定された ATE（3 行目）は 2.956 となり、11 年間禁煙した場合には、喫煙を続けた場合と比較して、11 年後の体重は 2.956kg 増加することが読み取れる。 . 層別解析 共変量や傾向スコアの値を層に分け、各層での結果変数の平均値と各層での割合を掛け合わせる。 各層の内部が均質であることを仮定することで、層内の ATE を算出できる。 ただし、正値性の問題が存在する。 今回は取り扱わない。 

重み付け法 観測データの各個体に適切な重みを割り当てることにより、処置群と対照群の分布が似ている疑似母集団を形成する。 重み付けは傾向スコアを用いて行う。特に逆確率重み付け(IPW)では、処置群には傾向スコアの逆数、対照群には(1-傾向ス コア)の逆数で重み付けをする。 以下は IPW を用いた例である。 learner = LogisticRegression( solver="liblinear", class_weight="balanced" ) ipw = IPW(learner=learner) ipw.fit(data.X, data.a) outcomes = ipw.estimate_population_outcome(data.X, data.a, data.y) effect = ipw.estimate_effect(outcomes[1], outcomes[0]) print(f"Average outcome of Untreated: {outcomes[0]:.3f}") print(f"Average outcome of Treated: {outcomes[1]:.3f}") print(f'Average Treatment Effect: {effect["diff"]:.3f}') 

Average outcome of Untreated: 1.563 Average outcome of Treated: 4.975 Average Treatment Effect: 3.412 出力結果より、 喫煙を続けた場合（出力 1 行目）の平均的な体重増加の推定値は 1.563 禁煙した場合（2 行目）の平均的な体重増加の推定値は 4.975 傾向スコアマッチングにより推定された ATE（3 行目）は 3.412 となり、11 年間禁煙した場合には、喫煙を続けた場合と比較して、11 年後の体重は 3.412kg 増加することが読み取れる。 

結果の評価 傾向スコアを用いて共変量を揃えることができているかを確認することで、条件付き独立性と正値性を満たしているかを評価す る。 ここでは以下の代表的な 3 つの方法を実装する。 絶対標準化平均差(ASMD) 分散比 図示 . 

絶対標準化平均差(ASMD) 共変量のバランスを評価するために絶対標準化平均差(ASMD)という指標を用いる。 : 処置群（介入群）における共変量の平均値 : 対照群（非介入群）における共変量の平均値 : 処置群と対照群を統合した標準偏差（Pooled Standard Deviation） : 処置群のサンプルサイズ : 対照群のサンプルサイズ : 処置群における共変量の分散 : 対照群における共変量の分散 ASMD = ​ ​ ​ σ ​ pool ​ − ​ X ˉ 1 X ˉ 0 σ ​ = pool 2 ​ N ​ + N ​ − 2 1 0 (N ​ − 1)σ ​ + (N ​ − 1)σ ​ 1 1 2 0 0 2 ​ X ˉ 1 ​ X ˉ 0 σ ​ pool N ​ 1 N ​ 0 σ ​ 1 2 σ ​ 0 2 

群間の共変量の平均値の差が小さいほど分布の重なりが大きく、ASMD は小さくなる。 両群を統合した分散が大きいほど分布の重なりが大きく、ASMD は小さくなる。 ASMD が 0.1 未満の場合にバランスが取れていると判断することが推奨されている。 重み付け後のすべての共変量の ASMD の値が 0.1 未満であることがわかる。 results = evaluate(ipw, data.X, data.a, data.y) fig, ax = plt.subplots(1, 1, figsize=(6, 6)) results.plot_covariate_balance(kind="love", ax=ax, thresh=0.1) 

分散比 群間の共変量の分散比を指標とする方法も提案されており、分散費が 2 未満の場合は一般的に許容可能なバランスと見做 される。 今回は実装しない。 図示 傾向スコアの広い範囲で正値性を満たしていることを確認する。 以下のように確認できた。 

傾向スコア分析の課題 未観測の交絡因子が存在する場合は正しい因果効果の推定ができない 重み付け法を用いる場合、極端に大きな重みを持つ個体には何らかの対処をする必要がある 

二重にロバストな推定法(DR) 回帰分析による因果効果の推定と IPW による推定を組み合わせた方法。 理論上、回帰分析か IPW のどちらかが正しければ他方のモデルが間違っていても正しい因果効果を推定することができる。 原理はあまり理解できなかった。実装結果を示す。 . 問い 前節(傾向スコア分析)と同様、禁煙が体重に与える影響 . . 

データ収集／確認 . DAG の確認 割愛 . !pip install causallib from sklearn.linear_model import LogisticRegression, LinearRegression from causallib.datasets import load_nhefs from causallib.estimation import IPW, Standardization, StratifiedStandardization from causallib.estimation import AIPW, PropensityFeatureStandardization, WeightedStandardization from causallib.evaluation import evaluate data = load_nhefs() . 

識別仮定の確認 回帰分析と傾向スコア分析の両方で識別仮定を満たす必要がある。 ここでは識別仮定を満たすとする。 . 因果効果の推定 qsmk 0 1.761707 1 5.205193 . ipw =IPW(LogisticRegression(solver="liblinear"), clip_min=0.05, clip_max=0.95) std = StratifiedStandardization(LinearRegression()) dr = AIPW(std, ipw) dr.fit(data.X, data.a, data.y) pop_outcome = dr.estimate_population_outcome(data.X, data.a, data.y) pop_outcome 

diff 3.443486 . 結果の評価 11 年間禁煙した場合には、喫煙を続けた場合と比較して、11 年後の体重は 約 3.4kg 増加することが読み取れる。 . 二重にロバストな推定法の課題 傾向スコアモデルが正確でも回帰分析が不正確である場合に大きなバイアスを引き起こす場合があるらしい 両方のモデルが少しでも誤っているとバイアスも分散も大きくなる可能性があるらしい 各モデルと同様、未観測の交絡因子があると正しく推定できない effect = dr.estimate_effect(pop_outcome[1], pop_outcome[0]) effect . 

自然実験(準実験) 

自然実験(準実験)の全体像 自然実験とは、RCT が実施されていない状況において、RCT が行われたかのような事象の発生を利用して因果効果を推定す る方法。 回帰分析や傾向スコア分析では未観測の交絡因子の存在を否定できない事が多く、条件付き独立性を満たしているとは言い 切れない場合が多い。 自然実験によるアプローチでは、分析手法別に特定の仮定を満たすことで RCT が行われたかのような状況を作れる。 今回は以下の手法を取り扱う。 回帰不連続デザイン 操作変数法 差分の差分法 合成コントロール法 . 

回帰不連続デザイン(RDD) 処置の有無を決める閾値(カットオフ値)付近で RCT に近い状況を作り出し、因果効果を推定する手法。 カットオフ値の右側極限と左側極限の割当はほぼランダムと考えることができる(例：合格最低点 点の生徒たちと 点の 生徒たちの間にほぼ差はないと考えられる) RDD は他の殆どの因果推論手法よりも仮定が少なく、RCT に最も近い手法であると言われている。 . 問い 顧客獲得戦略がユーザーのサービス利用時間に与える影響はどの程度か？ . ±0 −1 . 

データ収集／確認 !pip install rdrobust !pip install rddensity import numpy as np import pandas as pd import matplotlib.pyplot as plt import rdrobust import rddensity np.random.seed(0) n = 1000 X = np.random.uniform(-1, 1, size=n) threshold = 0 treatment = (X > threshold).astype(int) Y = 3 * X + 2 * treatment + np.random.normal(size=n) df = pd.DataFrame({ 'Y':Y, 'X':X, 'treatment': treatment }) 

X : 健康スコア Y : サービス利用時間 カットオフ値 : 0 健康スコアが 0 以上の人にだけ初月無料のプロモーションを行ったと想定したダミーデータ plt.scatter(df['X'], df['Y']) plt.xlabel('X') plt.ylabel('Y') plt.show() 

DAG の確認 サービス利用時間 プロモーション 健康スコア X T Y 識別仮定の確認 正値性以外の識別仮定は満たしていると見做せる事が多いが、厳密に正値性を満たすことは困難。 RDD の場合は連続性の仮定、つまり「処置がなかったのならば結果変数はカットオフ値の近傍で連続であったはずである」という ことを確認できたならばこれを正値性に換えることができる。次でその方法をまとめる。 . 

方法 1. カットオフ値付近の個体の分布の確認 横軸にカットオフの基準に使う変数をとって個体数をプロットする。 plt.hist(df[df['X'] < threshold]['X'], bins=10, alpha=0.5, label='BelowThreshold') plt.hist(df[df['X'] >= threshold]['X'], bins=10, alpha=0.5, label='AboveThreshold') plt.axvline(x=threshold, color='r', linestyle='--', label='Threshold') plt.xlabel('X') plt.ylabel('Number of Observations') plt.title('Histgram of x') plt.legend() plt.show() 

上記の出力結果は問題ないが、以下の右側のヒストグラムのようにカットオフ値近傍の観測値の数が急激に変化するような場 合の RDD の信憑性は低い傾向にある。 

方法 2. カットオフ値付近の連続性の検定 カットオフ値近傍において連続であることを帰無仮説、連続でないことを対立仮説し、帰無仮説が棄却されなければ連続である とみなす。 なぜ帰無仮説と対立仮説を逆に設定しないのかはわからなかった。 Manipulation testing using local polynomial density estimation Number of obs: 1000 Model: unrestricted Kernel: triangular BW method: estimated VCE: jackknife c = 0 Left of c Right of c Number of obs: 517 483 Eff. number of obs: 184 153 Order est. (p): 2 2 Order bias. (q): 3 3 BW est. 0.3446 0.3336 rddensity.rddensity(X, c=0) 

Method: T P > |T| Robust 0.963 0.3356 P-values of binomial tests (H0: p = [0.5] ). Window Length/2 < c >= c P>|T| 0.9616 498 468 0.3508 今回は設定した有意水準を上回り、連続性の仮定を満たすとみなせたことにする。 モヤモヤが残るが、「一般的に帰無仮説検定の非対称性により、帰無仮説を棄却できない場合、帰無仮説が正しいことを示さ れた訳では無いが、帰無仮説が間違っているという証拠もないわけであるから、便宜的に「回帰不連続デザインの適用に問題は ない」とみなす」(高橋(2022))らしい 

因果効果の推定 バンド幅の選択 バンド幅とは、因果効果の推定の対象範囲を決める幅。 バンド幅はトレードオフの関係があり、バンド幅が小さすぎる場合はカットオフ値近傍のデータを使用するためバイアスは小さくなる 一方で、使用可能なデータが減少するため分散は大きくなる。一方で、バンド幅が大きすぎる場合は推定量の分散は小さくなり やすい一方で、バイアスが大きくなる。 そのため、バンド幅を選択する際の指標としてバイアスと分散の両方を考慮できる MSE が使用されることが多く、その場合は MSE が最小化されるバンド幅が選択される。 MSE を算出するための推定値の算出には局所多項式回帰、その中でも LOWESS が用いられる。同時にカーネル密度推定を 使用することで、カットオフ値近傍に重み付けをした回帰が可能。カーネル関数には三角形関数が用いられることが多い。 

因果効果の推定 Call: rdrobust Number of Observations: 1000 Polynomial Order Est. (p): 1 Polynomial Order Bias (q): 2 Kernel: Triangular Bandwidth Selection: mserd Var-Cov Estimator: NN Left Right Number of Observations 517 483 Number of Unique Obs. 517 483 Number of Effective Obs. 112 97 Bandwidth Estimation 0.21 0.21 Bandwidth Bias 0.362 0.362 rho (h/b) 0.58 0.58 Method Coef. S.E. t-stat P>|t| 95% CI Conventional 1.194 0.304 3.922 8.767e-05 [0.597, 1.791] Robust - - 3.091 1.994e-03 [0.403, 1.799] rdrobust.rdrobust(Y, X, c=0) 

デフォルトでカーネルの種類として三角形関数、バンド幅の選択方法として MSE が設定されている。結果、バンド幅として 0.21 が選択されている。 因果効果は Coef の部分で 1.194(p 値：8.767e-05)として推定された。つまり、プロモーションを受けた新規顧客のサービス利 用時間はプロモーションを受けていない新規顧客と比較して 1 日あたり焼く 1.2 時間長いという分析結果が得られた。 

結果の評価 結果のプロット . rdrobust.rdplot(Y, X, c=0) 

プラセボテスト ~(省略)~ Method Coef. S.E. t-stat P>|t| 95% CI Conventional 0.349 0.392 0.889 3.741e-01 [-0.42, 1.118] Robust - - 0.883 3.773e-01 [-0.506, 1.337] 因果効果は 0.349 として推定されたが、p 値がよく使われる有意水準よりも大きいため、カットオフ値に 0.5 を設定した場合に は因果効果がほぼないことがわかった。 . 回帰不連続デザインの課題 カットオフ値から離れた場所での効果についての一般化は基本的に困難である カットオフ値近傍でのデータが十分に存在しない場合には適さない カットオフ値が都合よく操作された推定結果は信憑性を欠く可能性がある rdrobust.rdrobust(Y, X, c=0.5) 

操作変数法(IV 法) から への因果効果を推定したいとき、以下の条件を満たす操作変数 を導入することで、未観測の交絡因子を取り除くこ とができる。 と は関連する ( ) が を通して以外に に影響を与えない (除外制約) と は共通の原因 を持たない 最初の条件については、実際には関連するだけではなく、 と強い相関を持つ を用意する必要がある。 評価の方法としては、 と の回帰において が無関係であるという帰無仮説のもとで F 検定統計量が 10 以上であれば良い とされている。 . T Y Z Z T Cov(Z, T) =  0 Z T Y Z Y T Z Z T Z 

上記の 3 つの条件に加え、SUTVA と後述する単調性を満たす場合、 への効果(回帰係数)を への効果(回帰係数)を への効果(回帰係数)を とおくと、 となる。 Z → Y β ​ ZY Z → T β ​ ZT T → Y β ​ TY β ​ = TY ​ = Cov(T, Z) Cov(Y , Z) ​ = ​ Var(Z) Cov(T,Z) ​ Var(Z) Cov(Y ,Z) ​ β ​ ZT β ​ ZY 

単調性 個体 に対する処置割当を 、その個体が実際に処置を受けたかどうかを とすると、各個体を以下のように分類できる。 この中で、③ 反抗者のみ処置変数 が実際に処置を受ける に真逆の影響を与えるため、操作変数を導入した目的に沿わ なくなる。そのため、操作変数法では以下の単調性を条件として設定している。 この条件を設定することで、推定する処置効果は群全体で考える ATE ではなく、② 遵守者への ATE を算出することになる。こ のときの遵守者への ATE を CACE もしくは LATE と言う。 i T ​ i D ​ i T ​ i D ​ i D ​ (T ​ = i i 1) ≥ D ​ (T ​ = i i 0) 

問い . データ収集／確認 今回はダミーデータ。大規模な運動キャンペーンが地域別にランダムに展開されたと想定する。 運動量が健康指標に与える影響はどの程度か？ . import numpy as np import pandas as pd import statsmodels.api as sm from statsmodels.sandbox.regression.gmm import IV2SLS np.random.seed(0) n = 1000 Z = np.random.normal(size=n) u = np.random.normal(size=n) T = Z + np.random.normal(size=n) Y = T + u 

. DAG の確認 未観測の交絡因子 健康指標 運動量 キャンペーン規模/ 人口 Z T Y C df = pd.DataFrame({'Y': Y, 'T': T, 'Z': Z, 'u': u}) Y = df['Y'] T = sm.add_constant(df['T']) Z = sm.add_constant(df['Z']) 

識別仮定の確認 先述の仮定を満たすと考えられる。 . 因果効果の推定 . iv = IV2SLS(Y, T, Z).fit() 

結果の評価 IV2SLS Regression Results ============================================================================== Dep. Variable: Y R-squared: 0.652 Model: IV2SLS Adj. R-squared: 0.651 Method: Two Stage F-statistic: 902.4 Least Squares Prob (F-statistic): 9.68e-142 ~(中略)~ ============================================================================== coef std err t P>|t| [0.025 0.975] const 0.0105 0.031 0.342 0.733 -0.050 0.071 T 0.9679 0.032 30.040 0.000 0.905 1.031 ============================================================================== Omnibus: 1.172 Durbin-Watson: 2.079 Prob(Omnibus): 0.557 Jarque-Bera (JB): 1.243 Skew: 0.058 Prob(JB): 0.537 Kurtosis: 2.872 Cond. No. 1.36 ============================================================================== 地域の平均運動量 が ① 単位増えると、その地域の平均健康指標 が 0.9679 増えることがわかった。 print(iv.summary()) T Y 

操作変数法の課題 適切な操作変数の存在が稀で、実行が難しい。これは書籍にない私の意見だが、今回実装した例のように能動的に操作 変数を作りに行く必要がある場合が多いのではないかと考える。 

差分の差分法(DiD) 処置群と対照群でそれぞれ処置前後の差分を取り、さらにその結果の差分をとる。 → 結果変数のトレンドの影響を取り除くため 未観測の交絡因子を取り除くために以下の仮定を設定している 並行トレンド仮定：もし処置が行われなかったのなら処置群と対照群の時系列変化は同じであった 共通ショック仮定：測定期間中に目的変数に影響を与えるような別の原因が存在していない、もしくは存在したとしても両 群に同様に影響している . 

問い 新しいマーケティング施策の売上増加効果はどの程度か？ . データ収集／確認 . . !pip install causalpy from sklearn.linear_model import LinearRegression import causalpy as cp data = cp.load_data("did") data.head() 

group t unit post_treatment y 0 0 0.0 0 False 0.897122 1 0 1.0 0 True 1.961214 2 1 0.0 1 False 1.233525 3 1 1.0 1 True 2.752794 4 0 0.0 2 False 1.149207 group : その店舗が処置群か対照群かを表す t : 時点を示す変数 unit : 個々の店舗を識別するためのラベル post_treatment : 介入以前か以後かを表す y : 売上(百万円)（目的変数） 

DAG の確認 販売額 施策 T Y . 識別仮定の確認 識別仮定に加え並行トレンド仮定と共通ショック仮定が満たされていることを確認する。 ここではデータの可視化や背景の理解からこれらの仮定を確認できたとする。 . 

因果効果の推定 result = cp.DifferenceInDifferences( data, formula="y ~ 1 + group*post_treatment", time_variable_name="t", group_variable_name="group", treated=1, untreated=0, model=LinearRegression(), ) fig, ax = result.plot(); 

結果の評価 因果効果が 0.5 百万円と表示されているため、新しいマーケティング施策を行った結果 1 店舗あたりの売上を 50 万円増加さ せる効果があったと言える。 . 差分の差分法の課題 並行トレンド仮定と共通ショック仮定を満たすことを示すことが困難 . その他 そのままでは並行トレンド仮定を満たさない場合でも、傾向スコアマッチングなどの手法と組み合わせれば満たせることもあり、 それなら因果効果を測定できるはず… . 

合成コントロール法(SC) 対照群に各個体に 0~1 の重みをつけて合成対照群を作成し、この群の処置前の結果変数を処置群の処置前の結果変数に 近似させる。 → 未観測のものも含め交絡因子の影響を取り除くことができる。並行トレンド仮定も共通ショック仮定も満たす必要がない。 . 問い (コロナ禍の)緊急事態宣言の発出が商業施設の売上に与える影響はどの程度か？ . . 

データ収集／確認 a b c d e f g counterfactual causal effect actual 0 0.793234 1.277264 -0.055407 -0.791535 1.075170 0.817384 -2.607528 0.144888 -0.0 0 1 1.841898 1.185068 -0.221424 -1.430772 1.078303 0.890110 -3.108099 0.601862 -0.0 0 2 2.867102 1.922957 -0.153303 -1.429027 1.432057 1.455499 -3.149104 1.060285 -0.0 1 3 2.816255 2.424558 0.252894 -1.260527 1.938960 2.088586 -3.563201 1.520801 -0.0 1 4 3.865208 2.358650 0.311572 -2.393438 1.977716 2.752152 -3.515991 1.983661 -0.0 1 各店舗の時系列売上で構成されている。 !pip install causalpy from sklearn.linear_model import LinearRegression import causalpy as cp df = cp.load_data("sc") treatment_time = 70 df.head() 

DAG の確認 商業施設の売上 緊急事態宣言 T Y 識別仮定の確認 前提より条件付き独立性、正値性、一致性を満たすといえる。 それぞれの商業施設が遠方にあるため相互作用なしと見做せるとする。 これらの識別仮定に加え、合成コントロール法では合成対照群が処置前の処置群を近似できているかも検証する必要がある。 これは次のインが効果の推定と同時に行う。 . 

因果効果の推定 . 結果の評価 result = cp.SyntheticControl( df, treatment_time, formula="actual ~ 0 + a + b + c + d + e + f + g", model=cp.skl_models.WeightedProportion(), ) fig, ax = result.plot(plot_predictors=True) 

この決定係数 は必要十分に大きいと見做せるとする。 累積の因果効果は最大瞬間風速的なもので-4000 万円に及ぶ。 R = 2 0.92 

合成対照群を作成する際に設定した重みをプロットする。 ==================================Pre-Post Fit================================== Formula: actual ~ 0 + a + b + c + d + e + f + g Model coefficients: a 0.38 b 0.17 c 0.44 d 0 e 0 f 0 g 1e-16 主に を用いていることがわかる。 result.summary() a, b, c 

最後にプラセボテストをしようとした。対照群の 1 つを処置群と置き、もとの処置群を含めて合成対照群を作り、無いはずの因 果効果が存在しないか確かめたかった。 しかし、処置前の時系列変化をうまく近似することができず、決定係数は大きいものでも しかなかった。合成コントロ ール法を適用できる場合はあまり多くなさそうである。 R = 2 0.65 result = cp.SyntheticControl( df, treatment_time, formula="f ~ 0 + actual + a + b + c + d + e + g", model=cp.skl_models.WeightedProportion(), ) fig, ax = result.plot(plot_predictors=True) 

合成コントロール法の課題 共変量が施策の前後で大きく変化する場合には適切な推定ができない 個体の数が多いほど、時系列の長さが短いほど合成対照群の振る舞いが不安定になる