Pythonによるマーケティングミックスモデリング
（MMM：Marketing Mix Modeling）超入門 その6

Ridge回帰モデルでMMM②（シンプルなAdStock考慮）

2621
Pythonによるマーケティングミックスモデリング（MMM：Marketing Mix Modeling）超入門 その6Ridge回帰モデルでMMM②（シンプルなAdStock考慮）

本当に売上に貢献している広告は、どの広告か？

売上と広告媒体等との関係性をモデリングし、どの広告媒体が売上にどれほど貢献していたのか分析することができます。

それが、マーケティングミックスモデリング（MMM：Marketing Mix Modeling）です。

ただ、多くの広告・販促は似たような時期に集中し実施することが多く、結果的にマルチコという線形回帰モデルを構築しうる上で非常に良くない現象が起こる危険があります。

その対策の1つに、正則化項付き線形回帰モデルの1つである、Ridge回帰でモデルを構築する、という方法があります。モデルの説明力や予測精度は若干悪化しますが、非常に有効な方法です。

3回に分けてお話しします。

  • ① AdStock非考慮　→　前回
  • ② シンプルなAdStock考慮　→　今回
  • ③ 色々なAdStock考慮

 

前回は、①のAdStockを考慮しないRidge回帰によるモデル構築についてお話ししました。

Pythonによるマーケティングミックスモデリング（MMM：Marketing Mix Modeling）超入門 その5Ridge回帰モデルでMMM①（AdStock非考慮）

今回は、②のシンプルなAdStockを考慮したRidge回帰によるモデル構築についてお話しします。

ちなみに、アドストック（Ad Stock）を考慮するとは、飽和効果（収穫逓減）とキャリーオーバー（Carryover）効果を考慮するということです。具体的には、Pythonのscikit-learn（sklearn）のパイプラインを、以下のように構築しました。

 

利用するデータセット（前回と同じ）

今回利用するデータセットの変数です。

  • Week：週
  • Sales：売上
  • TVCM：TV CMのコスト
  • Newspaper：新聞の折り込みチラシのコスト
  • Web：Web広告のコスト

以下からダウンロードできます。

MMM.csv
https://www.salesanalytics.co.jp/4zdt

 

アドストック（Ad Stock）

今回は、以下の関数でアドストック（Ad Stock）を表現しモデル構築したいと思います。

  • 定率減少型キャリーオーバー効果モデル　※ピークが広告などの投入時
  • 指数型飽和モデル（exp_Saturation）

 

　定率減少型キャリーオーバー効果モデル

ピークが広告などの投入時で徐々に低減する「定率減少型キャリーオーバー効果モデル」を、数式で表現すると、次のようになります。

\displaystyle x_t^* = \displaystyle \sum_{l=0}^{L-1} w_{t-l} \cdot x_{t-l} , \quad w_{t-l} = R^{l} x_tはt期の広告などの投入量で、x_ t^*はt期とそれ以前までの広告などの効果の累積（残存効果を足したもの）です。

ハイパーパラメータは、次の2つです。

  • L(length)：効果の続く期間　※当期含む
  • R(rate)：減衰率

グラフで表現すると、次のようになります。横軸は、広告などを投入した後の経過時間（当期、1期後、2期後、・・・）を表し、縦軸は、効果の大きさを表します。

 

　指数型飽和モデル

「指数型飽和モデル」を数式で表現すると、次のようになります。

x_ t^{**}=1-e^{-ax_ t^*} x_ t^*はt期とそれ以前までの広告などの効果の累積（残存効果を足したもの）で、x_ t^{**}は指数型飽和モデルで変換した後のt期の値です。

グラフで表現すると、次のようになります。横軸はx_ t^*で、縦軸はx_ t^{**}です。

 

この2つの関数に関しては、次の記事でもうちょっと詳しく説明しているので、気になる方は参考にしてください。

 

ライブラリーやデータセットなどの読み込み

先ずは、必要なライブラリーとデータセットを読み込みます。前回までと同じです。

 

　必要なライブラリーの読み込み

必要なライブラリーを読み込みます。

以下、コードです。

import numpy as np
import pandas as pd

from scipy.signal import convolve2d

from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score, TimeSeriesSplit
from sklearn import set_config

from sklearn.linear_model import Ridge

from optuna.integration import OptunaSearchCV
from optuna.distributions import UniformDistribution, IntUniformDistribution

import matplotlib.pyplot as plt
plt.style.use('ggplot') #グラフスタイル
plt.rcParams['figure.figsize'] = [12, 9] # グラフサイズ

 

　データセット（前回と同じ）の読み込み

前回と同じデータセットを読み込みます。

以下、コードです。

# データセット読み込み
url = 'https://www.salesanalytics.co.jp/4zdt'
df = pd.read_csv(url,
                 parse_dates=['Week'],
                 index_col='Week'
                )

# データ確認
print(df.info()) #変数の情報
print(df.head()) #データの一部

# 説明変数Xと目的変数yに分解
X = df.drop(columns=['Sales'])
y = df['Sales']

 

以下、実行結果です。

 

ハイパーパラメータ探索

Ridge回帰の最適な正則化パラメータの値を探索します。

先ずは、Pipeline用変換器を定義します。

以下、コードです。

#
# Pipeline用変換器
#

# 飽和モデル

## 関数の定義
def Saturation(X,a):
    return 1 - np.exp(-a*X)

## 変換器の定義
class ExponentialSaturation(BaseEstimator, TransformerMixin):
    
    # 初期化
    def __init__(self, a=1.0):
        self.a = a #ハイパーパラメータ
    
    # 学習
    def fit(self, X, y=None):
        return self
    
    # 処理（入力→出力）
    def transform(self, X):
        return Saturation(X,self.a)

# キャリーオーバー効果モデル

## 関数の定義
def Carryover(X: np.ndarray, rate, length):
    filter = (
        rate ** np.arange(length)
    ).reshape(-1, 1) 
    convolution = convolve2d(X, filter)
    if length > 1 : convolution = convolution[: -(length-1)]
    return convolution

## 変換器の定義
class ExponentialCarryover(BaseEstimator, TransformerMixin):
    # 初期化
    def __init__(self, rate=0.5, length=1):
        self.rate = rate #ハイパーパラメータ
        self.length = length #ハイパーパラメータ
        
    # 学習        
    def fit(self, X, y=None):
        return self
    
    # 処理（入力→出力）
    def transform(self, X: np.ndarray):
        return Carryover(X, self.rate, self.length)

 

次に、Pipelineを構築します。

以下、コードです。

#
# Pipeline構築
#

# 変換器（adstock）
adstock = ColumnTransformer(
    [
     ('TVCM_pipe', Pipeline([
                           ('carryover', ExponentialCarryover()),
                           ('saturation', ExponentialSaturation())
     ]), ['TVCM']),
     ('Newspaper_pipe', Pipeline([
                           ('carryover', ExponentialCarryover()),
                           ('saturation', ExponentialSaturation())
     ]), ['Newspaper']),
     ('Web_pipe', Pipeline([
                           ('carryover', ExponentialCarryover()),
                           ('saturation', ExponentialSaturation())
     ]), ['Web']),
    ],
    remainder='passthrough'
)

## 変換器（adstock）→学習器（regression）
MMM_pipe = Pipeline([
    ('adstock', adstock),
    ('regression', Ridge())
])

## パイプラインの確認
set_config(display='diagram')   
MMM_pipe

 

以下、実行結果です。

 

そして、Pipelineのハイパーパラメータを探索します。

以下、コードです。

#
# ハイパーパラメータ探索
#

# 探索するハイパーパラメータ範囲の設定
params = {
 ## Ridge
 'regression__alpha':
    UniformDistribution(0.01, 10000),
 ## TVCM
 'adstock__TVCM_pipe__carryover__rate':
    UniformDistribution(0, 1),
 'adstock__TVCM_pipe__carryover__length':
    IntUniformDistribution(1, 6),
 'adstock__TVCM_pipe__saturation__a':
    UniformDistribution(0, 0.01),
 ## Newspaper
 'adstock__Newspaper_pipe__carryover__rate':
    UniformDistribution(0, 1),
 'adstock__Newspaper_pipe__carryover__length':
    IntUniformDistribution(1, 6),
 'adstock__Newspaper_pipe__saturation__a':
    UniformDistribution(0, 0.01),
 ## Web
 'adstock__Web_pipe__carryover__rate':
    UniformDistribution(0, 1),
 'adstock__Web_pipe__carryover__length':
    IntUniformDistribution(1, 6),
 'adstock__Web_pipe__saturation__a':
    UniformDistribution(0, 0.01),
}

# ハイパーパラメータ探索の設定
optuna_search = OptunaSearchCV(
    estimator=MMM_pipe,
    param_distributions=params,
    n_trials=2000,
    cv=TimeSeriesSplit(),
    random_state=0,
)

# 探索実施
optuna_search.fit(X, y)

# 探索結果
optuna_search.best_params_

 

以下、実行結果です。

 

この値をハイパーパラメータに設定し、MMMを全データで構築します。

 

モデルの構築と予測

　最適ハイパーパラメータで学習

以下、コードです。

#
# 最適ハイパーパラメータで学習
#

# パイプラインのインスタンス
MMM_pipe_best = MMM_pipe.set_params(**optuna_search.best_params_)

# 全データで学習
MMM_pipe_best.fit(X, y)

# R2（決定係数）
MMM_pipe_best.score(X, y)

 

以下、実行結果です。R2（決定係数）の値です。

 

　予測式の確認

以下、コードです。

# 
# 予測式の確認
# 

# 線形回帰モデルの切片と回帰係数
intercept = MMM_pipe_best.named_steps['regression'].intercept_ #切片
coef = MMM_pipe_best.named_steps['regression'].coef_ #回帰係数
alpha = MMM_pipe_best.named_steps['regression'].alpha #正則化パラメータ

# 回帰係数をSeries（シリーズ）化
weights = pd.Series(
    coef,
    index=X.columns
)

# 結果出力（切片と係数）
print('Intercept:\n', intercept, sep='')
print()
print('Coefficients:\n',weights, sep='')
print()
print('alpha:\n', alpha, sep='')

 

以下、実行結果です。

 

　予測の実施

以下、コードです。

#
# 予測の実施
# 

# 目的変数y（売上）の予測
pred = pd.DataFrame(MMM_pipe_best.predict(X),
                    index=X.index,
                    columns=['y'])

# 各媒体による売上の予測

## 値がすべて0の説明変数
X_ = X.copy()
X_.iloc[:,:]=0

## Base
pred['Base'] = MMM_pipe_best.predict(X_)

## 媒体
for i in range(len(X.columns)):
    X_.iloc[:,:]=0
    X_.iloc[:,i]=X.iloc[:,i]
    pred[X.columns[i]] = MMM_pipe_best.predict(X_) - pred['Base']
    
print(pred) #確認

 

以下、実行結果です。y = Base + TVCM + Newspaper + Webです。

 

貢献度とマーケティングROIの算定

　貢献度の算定

以下、コードです。

#
# 貢献度の算定
#

# 予測値の補正
correction_factor = y.div(pred['y'], axis=0)   #補正係数
pred_adj = pred.mul(correction_factor, axis=0) #補正後の予測値

# 各媒体の貢献度だけ抽出
contribution = pred_adj[['Base', 'Web', 'Newspaper', 'TVCM']]

print(contribution.head()) #確認

# グラフ化
ax = (contribution
      .plot.area(
          ylabel='Sales',
          xlabel='Week')
     )
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles[::-1], labels[::-1])

 

以下、実行結果です。

 

貢献度の構成比の算出します。

以下、コードです。

#
# 貢献度の構成比の算出
#

# 媒体別の貢献度の合計
contribution_sum = contribution.sum(axis=0)

# 集計結果の出力
print('売上貢献度（円）:\n', 
      contribution_sum, 
      sep=''
     )
print()
print('売上貢献度（構成比）:\n', 
      contribution_sum/contribution_sum.sum(), 
      sep=''
     )

# グラフ化
contribution_sum.plot.pie(fontsize=24)

 

以下、実行結果です。

 

マーケティングROIの算定

以下、コードです。

#
# マーケティングROIの算定
#

# 各媒体のコストの合計
cost_sum = X.sum(axis=0)

# 各媒体のROIの計算
ROI = (contribution_sum.drop('Base', axis=0) - cost_sum)/cost_sum

print('ROI:\n', ROI, sep='') #確認

# グラフ化
ROI.plot.bar(fontsize=24)

 

以下、実行結果です。

 

まとめ

今回は、シンプルなdStockを考慮したRidge回帰によるモデル構築についてお話ししました。

  • AdStock非考慮
  • シンプルなAdStock考慮　→　今回
  • 色々なAdStock考慮　→　次回

次回は、色々なAdStockを考慮したRidge回帰によるモデル構築についてお話しします。

Pythonによるマーケティングミックスモデリング（MMM：Marketing Mix Modeling）超入門 その7Ridge回帰モデルでMMM③（色々なAdStock考慮）