【A/Bテストなし】競合施策の効果を推定したい|因果推論と時系列解析
こんにちは、青の統計学です。
今回は、業務で行き詰まった部分をまとめてみました。何かの役に立てたら幸いです。
因果推論を使ったアプローチと時系列解析を使ったアプローチをご紹介します。
ケース|同業の施策の自社影響を知りたい
基本的に施策によるROIや増分効果や増分コストなどを見る場合、属性がほぼ同じのコントロール群と介入群を分けて、同じ\(t\)期のアウトカムの差分を見れば、平均処置効果がわかります。
ただ、今回は同業の施策効果なのでそのようなデータはないです。
経営陣は、同業のキャンペーン施策が自社のキャンペーン施策に影響を与え、それが結果的にKPIに影響するのではないかということを気にしています。
多分影響するとは思いますが、どれだけ影響するかによって、
「自社のキャンペーンを強化すべきか」
「強化すべきならどのように強化するのが最適か」
などが変わってきます。
どのチャネル経由でKPIに影響するのかも気になりますよね。
SEO経由のコンバージョンには大した影響はないが、SNS経由のコンバージョンには大きな影響がある、など。
| time | KPI | SEO経由CV | Lis経由CV | 自社施策-競合施策の金額差分 | 競合施策ダミー |
| 23/10 | 7000 | 4000 | 2000 | 0 | 0 |
| 23/11 | 5000 | 2000 | 1500 | -4000 | 1 |
| … | … | … | … | …. | … |
自社も同業も恒常的にキャンペーンを実施しており、今回は同業がキャンペーン金額を強化して自社キャンペーンが負けている状況を考えます。
キャッシュバックキャンペーンやポイント付与キャンペーンなど、世の中にはいろいろなバラマキ施策がありますね。
今回もその一つとして考えてください。
アプローチ1：SARIMAモデル
自分がまず考えたのが、時系列分析を使って自社と同業の金額差分を外生変数に加え、その回帰係数を求めると言うものです。
KPIには1年周期の季節性がありましたので、以下のようなSARIMAモデルを考えました。
$$(1-φ_1L)Y_t=δ+(1-Φ_1L^{12})ε_t+β_xDiff_{t}$$
上は一般形ではなく、例として、\(SARIMAX(1, 0, 0)(1, 0, 0)[12]\)モデルの数式を表現しています。
このモデルは、非季節性の自己回帰（AR）成分が1次、移動平均（MA）成分が0次、差分（I）成分が0次であり、季節性のAR成分が1次、MA成分が0次、季節性の差分成分が0次であることを意味します。
また、周期は12期です。季節性の自己回帰（SAR）成分が1つあることを意味します。
このモデルでは、12期前（通常1年前）のデータポイントが現在のデータポイントに直接的な影響を与えると仮定しています。
【時系列】ARモデルをわかりやすく解説|Yule-Walker法や最尤法も
もう少し各変数のご紹介をします。
\(Y_t\)はt期の目的変数、SEO経由CVとかですね。
\(L\)はラグ演算子です。後述します。
\(φ_1\)は非季節性自己回帰(AR)の1次の係数です。
\(Φ_1\)は季節性自己回帰(SAR)の12期の係数です
\(δ\)は定数項
\(ε_t\)はt期の誤差項です。
\(Diff\)は自社と同業のキャンペーン金額の差分ですね。今回はこの説明変数の特徴量ベクトルを求めたらいいじゃん。私は思いました。
ラグ演算子(lag operator)について
\(L\)、ラグ演算子（Lag Operator）は、時系列分析において非常に重要な概念で、時系列データの過去の値にアクセスするために使われます。
以下のような形で、過去の時点の値にアクセスすることができます。
$$L^kY_t=Y_{t-k}$$
上のように、t期時点のアウトカムをt-k時点のアウトカムを使って表せます。
これは「ラグ演算子をk回適用した状態」と呼びます。
ラグ演算子を使うと何が嬉しいのか
さて、たとえばARIMAモデルでは、非定常時系列を定常時系列に直すために差分を取りますが、ラグ演算子を使うと表現が簡単です。
つまり、モデルの安定性や因果性に関する理論的な特性を分析しやすくなります。
auto_arimaなどで探索されたパラメータによって形は変わりますが、以下のように表現できます。
$$(1-φ_1L-φ_2L^2-…-φ_pL^p)Y_t=δ+ε_t+移動平均項$$
モデルのパラメータがデータのラグされた値にどのように影響を及ぼしているかを直感的に理解しやすくなります。
また、モデルの係数が時系列の動的な挙動をどのように捉えているかの解釈が容易になります。
ラグ演算子は統計検定準1級で出ましたね。
CODE|python
＊ランダムに生成したデータを使うので、結果は参考にはなりません。
今回使うpmd_auto_arimaは、ARIMAの次数を、機械学習的なのりで、予測精度を最大限に高める次数を自動探索し決めようというライブラリです。
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pmdarima import auto_arima
from sklearn.metrics import mean_absolute_error, mean_squared_error
# データの読み込み（月次データを生成）
np.random.seed(0)
dates = pd.date_range('2020-01-01', periods=100, freq='M')
data = np.random.uniform(2000, 10000, 100)
df = pd.DataFrame(data, columns=['SEO_cv'], index=dates)
# 外生変数「diff」を追加（-5000から+5000の範囲のランダムな値）
df['diff'] = np.random.uniform(-5000, 5000, size=len(df))
# auto_arima
model = auto_arima(df['SEO_cv'], X=df[['diff']], seasonal=True, m=12,
trace=True, error_action='ignore', suppress_warnings=True)
print(model.summary())ここではランダムにデータを生成しています。
公式のドキュメントを参照し、外生変数として同業とのキャンペーン金額差分を「diff」としてX=diffとして引数に入れています。
今回diffの回帰係数を見ると、-0,0872でした。
同業のキャンペーン金額が上がった方が、SEO経由CVが上がるという解釈になります。
これは、ランダムなデータですので参考になりませんが、仮説としては同業のキャンペン金額の方が高いと、コンバージョンは落ちます。
ドメイン知識を使って分析結果は解釈していきましょう。
# モデルを使って予測を行う
train = df[:-10] # training data
test = df[-10:] # test data
model.fit(train['SEO_cv'])
forecast = model.predict(n_periods=len(test))
# 実際の値と予測値を比較する
plt.figure(figsize=(10, 6))
plt.plot(train.index, train['SEO_cv'], label='Train Data')
plt.plot(test.index, test['SEO_cv'], label='Test Data', color='gray')
plt.plot(test.index, forecast, label='Forecast', color='red')
plt.title('SEO_cv Forecast vs Actuals')
plt.xlabel('Date')
plt.ylabel('SEO_cv')
plt.legend()
plt.show()
# 評価
mae = mean_absolute_error(test['SEO_cv'], forecast)
mse = mean_squared_error(test['SEO_cv'], forecast)
print(f"Mean Absolute Error: {mae}")
print(f"Mean Squared Error: {mse}")では教師データとテストデータを使って結果を評価するコードを紹介します。
Mean Absolute Error: 2188.6810528495616 Mean Squared Error: 5854630.946243917
モデルの解釈について考えます。
出力を見ると、auto_arimaによる探索では、ARIMA(2,1,1)(0,0,0)[12]がベストらしいです。
このモデルを数式にしてみましょう。
$$(1-φ_1L-φ_2L^2)\nabla Y_t=δ+θ_1ε_{t-1}+ε_t$$
\(φ\)は自己回帰係数ですね。モデルが過去期間のデータポイントを参照していることを示しています。
\(θ_1\)は移動平均係数ですね。MAの次数が1です。
このモデルでは、季節性成分は考慮されていないため、(0,0,0)[12] の部分はモデルの数式には直接現れません。
季節性周期が12期ということは、データに年次の季節性がある場合に考慮する必要がありますが、このモデルでは季節性の効果は表現されませんでした。
時系列解析に関しては、以下のコンテンツをご覧ください。
【周期性を掴もう】pythonでコレログラムを書いてみましょう
【最短合格】統計検定準一級のチートシート|難易度や出題範囲について
アプローチ2：因果推論
因果推論的なアプローチについて、重要な論点は以下になります。
①SEO経由コンバージョンやリスティング経由コンバージョンのうち、同業施策の強化の影響を受けたものと受けなかったものはあるか。
②あった場合に、それら2群に平行トレンドの仮定は認められるか。
DID(difference in differnece)について
①②があった場合には、DID(差の差)による因果効果の推定が可能になります。
$$ATE=(Y_{post}^1-Y_{pre}^1)-(Y_{post}^0-Y_{pre}^0)$$
\(Y^1_{post}\) は、処置群（treatment group）の処置後（post-treatment）の平均アウトカムです。
\(Y^1_{pre}\) は、処置群の処置前（pre-treatment）の平均アウトカムです。
\(Y^0_{post}\) は、対照群（control group）の処置後の平均アウトカムです。
\(Y^0_{pre}\) は、対照群の処置前の平均アウトカムです。
仮にSEO経由のコンバージョンは同業施策の影響をうけ、リスティング経由のコンバージョンは同業施策の影響を受けない場合、処置群を前者、対象群を後者にして「同業施策強化」の平均処置効果を推定することができます。
図にすると以下のような感じです。
上の図のように、AとBのトレンドが同じ出ない限り、平均処置効果にはノイズが入ってしまいますね。
SEM(Structural Equation Modeling)について
今回のような、同業影響調査などはA/Bテストのように、介入群と対照群を事前に設定するわけでもないので、影響が全てのチャネルに及んでいる場合があります。
そのような場合には、DIDを使って解釈するのは難しいです。
(もちろん、同一チャネルで時期をずらせば、介入時期とコントロール時期で比較することはできます。ただ、トレンドが線形と言えなければ解釈を誤る危険があります。)
今回は、共分散構造分析を使った因果推論的なアプローチを考えます。
SEMの魅力は、ドメイン知識に基づき仮説モデルを構築でき、PDCAを回せる点です。
繰り返し繰り返しSEMをおこない、モデルを洗練させていくことができます。
パス解析や構造法的式モデリング（SEM）は、説明変数間の相関を積極的にモデルの一部に取り込み、分析者が持つ因果の仮説を表現して評価するために使用します。
つまり、単に変数間の関係性を許容するというよりも理論的な枠組みの中でその関係性を積極的に表現して、より深い因果関係の理解に繋げようというアプローチです。
そして、この手法には当然ドメイン知識などの質的なアプローチがとても役に立ちます。
観測変数と潜在変数
SEMのもう一つの利点は、潜在変数という観測できない変数をモデルに組み込むことができる点です。
今回の例で言えば、同業のキャンペーン金額が上昇したときに、「自社のキャンペーンの客観的な魅力度」が下がり、それが自社キャンペーンの通過率に影響し、各チャネルのコンバージョン数にマイナスの影響を与えると考えられます。
自社のキャンペーンの客観的な魅力度、のような観測できる変数の背後に隠れている変数を潜在変数と呼びます。
因子モデルでは、因子（あるいは共通因子）と呼びます。
CODE|共分散構造分析
semopyというライブラリを使います。
描画用にgraphvizもインストールします。
pip install semopybrew install graphvizimport semopy
import pandas as pd
import numpy as np
# モデルの定義
model_desc = """
SEO_CV ~ CP_diff
Lis_CV ~ CP_diff
KPI ~ SEO_CV + Lis_CV
"""
# 乱数を使ってデータセットを生成する
np.random.seed(42) # 再現性のためシードを設定
n_samples = 100
CP_diff = np.random.normal(size=n_samples)
SEO_CV = 0.5 * CP_diff + np.random.normal(size=n_samples)
Lis_CV = 0.3 * CP_diff + np.random.normal(size=n_samples)
KPI = 0.2 * SEO_CV + 0.4 * Lis_CV + np.random.normal(size=n_samples)
data = pd.DataFrame({'CP_diff': CP_diff, 'SEO_CV': SEO_CV, 'Lis_CV': Lis_CV, 'KPI': KPI})
# モデルのインスタンス化と推定
model = semopy.Model(model_desc)
result = model.fit(data)
# 結果の表示
print(result)Name of objective: MLW Optimization method: SLSQP Optimization successful. Optimization terminated successfully Objective value: 0.023 Number of iterations: 9 Params: 0.357 0.528 0.120 0.378 0.766 1.122 0.884
from semopy.inspector import inspect
from graphviz import Digraph
# 因果関係の情報を取得
inspect = model.inspect
semopy.semplot(model,"sem_graph.png")
p値と一緒にパス係数も出力されます。
上の状況だと、潜在変数はないのでパス解析が正確な呼び方です。
SEMの結果は通常、標準化されたパス係数として提供されます。
これは、因子と変数間の関係の相対的な強さを示します。標準化されたパス係数は、因子が一標準偏差変化したときに変数が何標準偏差変化するかを示しています。
構造方程式の結果の解釈について
では、パス係数はどう使いましょうか？
ここで意思決定者が気にしているポイントを再喝すると、「競合のキャンペーン施策がどのチャネルにどれくらい影響しているのか」という点です。
なので、各チャネルのコンバージョン数に向けたパス係数の解釈をする必要があります。
これは、そのまま回帰係数のように説明変数が1単位増えると~増えるのように使うことはできないので、各特徴量の標準偏差を使います。
たとえば、因子Aが1標準偏差変化する場合の変数Bの変化量は次のように計算されます
$$ΔA=1×σ_A$$
$$ΔB=λ_{AB}× \frac{σ_B}{σ_A}×ΔA$$
\(λ_{AB}\)はパス係数です。
上の式は因子Aと変数B間の関連の強さを実際の変化量で示しています。
標準化されたパス係数を用いることで、異なるスケールの変数間の関係を解釈しやすくなります。
因果推論に関しては、以下のコンテンツをご覧ください。
最後に
データサイエンティストの実業務において、「データが整形されてかつ数も多い」「分析基盤が整っている」「モデルをスポットで作るのではなく、オンラインで更新できるデータ基盤がある」などは、なかなかないです。
やはり、「月次データしかない中で解析をしてほしい」「来週までにある施策の因果効果が欲しい。ただA/Bテストはやっていません」など、無理難題が多いです。
「このデータだけだと分散がめっちゃ大きくなるな…」
「これ外部要因を全く考慮していないから因果効果と言えるのか…」
など、統計的な側面から言えば突っ込みたくなる状況です。
ただ、できない理由を並べるのは簡単ですが、やはり依頼には最大限答えたいので、知恵を絞って今あるデータを最大限に使って分析をしていきたいです。
