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

アドストック（Ad Stock）を考慮した線形回帰モデル

9409
Pythonによるマーケティングミックスモデリング（MMM：Marketing Mix Modeling）超入門 その2アドストック（Ad Stock）を考慮した線形回帰モデル

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

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

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

前回は、「線形回帰モデルでMMMを作ろう！」というお話しをします。

Pythonによるマーケティングミックスモデリング（MMM：Marketing Mix Modeling）超入門 その1線形回帰モデルでMMMを作ろう！

単純な線形回帰モデルでも、それなりのモデルが作れます。

ただ、このモデルには、ある致命的な欠陥があります。アドストック（Ad Stock）を考慮していないということです。

アドストック（Ad Stock）を考慮したモデルにした方がいいでしょう。

ちなみに、アドストック（Ad Stock）を考慮するとは、残存効果を考慮するということです。ある日の広告宣伝活動が、その日だけに効果があるのではなく、次の日以降もその効果が続くということです。キャリーオーバー（Carryover）効果と表現されることもあります。

ということで今回は、「アドストック（Ad Stock）を考慮した線形回帰モデル」についてお話しします。

利用するデータは前回と同じです。

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

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

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

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

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

 

アドストック（Ad Stock）

　飽和効果＋キャリーオーバー効果＝アドストック

前回の線形モデルだと、広告などに触れたその日に、商品を購買するというモデルです。

広告などには通常、キャリーオーバー（Carryover）効果もしくはラグ効果、残存効果と呼ばれるものがあります。

要するに、効果が遅れてくる、もしくは、効果が後々まで残っている、というものです。

例えば、高価な商品ですと、広告などに触れた瞬間に購買するというよりも、数日もしくは数週間後に、慎重に検討した後に購買することでしょう。

消費財などの比較的安価な商品でも、何度も広告などに触れているうちに、購買したくなり、数日後に購買するということもあることでしょう。

以下は、100の効果が、徐々に薄れながらも残存して残っているイメージです。縦軸は効果で、横軸は時間軸です。詳細は、後で説明します。

 

前回の線形回帰モデルだと、ある広告にコストをかければかけるほど、売上は無限に増えていきます。

しかし、現実は異なります。

ある広告にコストをかければかけるほど、売上の上昇幅は鈍くなります。経済学でいうところの収穫逓減が起こります。売上は飽和し、いくらコストをかけても売上が伸びなくなります。

以下は、横軸は広告などの投入量で、縦軸は効果の大きさです。詳細は、後で説明します。

 

飽和効果（収穫逓減）やキャリーオーバー（Carryover）効果を考慮したモデルがいいでしょう。この2つを考慮した数値変換（変換器）の機能がアドストック（Ad Stock）です。

アドストック（Ad Stock）を考慮した線形回帰モデルで目的変数である売上（Sales）を予測するとは、広告などの説明変数のデータ（今回の例では、広告などのコスト）を数値変換（変換器）し、線形回帰モデルの新たなインプットデータを作り、そのインプットデータを使い線形回帰モデルで目的変数である売上（Sales）を予測することです。

元の説明変数のデータを数値変換（変換器）し、このデータを使い線形回帰モデルを構築します。この一連の処理の流れを機械学習っぽい用語で表現するとパイプラインと呼びます。

 

　パイプライン（入力→変換器→学習器→出力）

今回の場合ですと、2つの変換器と1つの学習器が必要になります。

  • 変換器
    • キャリーオーバー効果モデル
    • 飽和モデル
  • 学習器
    • 線形回帰モデル

今回利用するキャリーオーバー効果モデルと飽和モデルについて簡単に説明します。

 

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

キャリーオーバー効果を表現する数理モデルは、色々あります。

今回は最もシンプルな、広告などを打ったときが効果のピークで、徐々に効果が一定の割合で減衰していくモデルです。

数式で表現すると、次のようになります。

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

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

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

 

期間とは、どのくらいまで考慮するか、ということです。

例えば、減衰率50%で期間3の場合、以下のようになります。

  • 当期：100
  • 次期：50
  • 次々期：25
  • 次次々期：12.5
  • 次次次々期：0
  • 次次次次々期以降：0

 

例えば、隔期に広告を100打った場合、それが積み重なります。

  • 当期：100
  • 次期：50
  • 次々期：25+100=125
  • 次次々期：12.5+50=62.5
  • 次次次々期：0+25+100=125
  • 次次次次々期：0+12.5+50=62.5

 

　飽和モデル

飽和（収穫逓減）を表現する数理モデルは、色々あります。

今回は、最もシンプルな指数関数 1-exp(-ax) でモデル化します。a はハイパーパラメータです。

以下は、先程説明した図（横軸は広告などの投入量で、縦軸は効果の大きさ）です。

 

MMMの実施

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

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

以下、コードです。

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 LinearRegression

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] # グラフサイズ

#指数表記しない設定
np.set_printoptions(precision=3,suppress=True)
pd.options.display.float_format = '{:.3f}'.format

 

　データセット（前回と同じ）

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

以下、コードです。

# データセット読み込み
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']

# グラフ化
y.plot()
X.plot()

 

以下、実行結果です。

 

　アドストックを考慮しない線形回帰モデル

アドストックを考慮しない線形回帰モデルを構築し、どの程度の精度を持ったものだったのかを見てみます。

ちなみに、前回と同じものです。

線形回帰モデルのインスタンスを作ります。

以下、コードです。

# 線形回帰モデルのインスタンス
lr = LinearRegression()

 

線形回帰モデルが、どの程度の予測精度を持ったモデルになるのかを確かめるために、時系列のCV（クロスバリデーション）を実施します。今回は、デフォルトの5分割のCVです。

以下、コードです。

# クロスバリデーションで精度検証（R2）
np.mean(cross_val_score(lr, 
                        X, y, 
                        cv=TimeSeriesSplit()
                       )
       )

 

以下、実行結果です。

 

全データでモデルを構築し、予測精度（ R^2）を見てみます。

以下、コードです。

# 全データで精度検証（R2）
lr.fit(X, y)
lr.score(X, y)

 

以下、実行結果です。

 

　アドストックの関数

以下の2つのモデル（変換器）の関数を定義し利用することで、アドストックを表現します。

  • 飽和モデル
  • キャリーオーバー効果モデル

 

　　飽和モデル

飽和モデル（関数）を定義します。

指数関数 1-exp(-ax) です。a のハイパーパラメータで、形状が異なってきます。

以下、コードです。

# 飽和モデル（関数）の定義
def Saturation(X,a):
    return 1 - np.exp(-a*X)

 

この関数の実行例です。

以下、コードです。

# 飽和モデル（関数）の例
exp_dat = pd.DataFrame(range(500)) #入力データ

exp_sat = pd.DataFrame(index=exp_dat.index)
exp_sat['a=0.01'] = Saturation(exp_dat,0.01) 
exp_sat['a=0.02'] = Saturation(exp_dat,0.02) 
exp_sat['a=0.1'] = Saturation(exp_dat,0.1) 


print(exp_sat) #数値出力
exp_sat.plot() #グラフ

 

以下、実行結果です。

 

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

キャリーオーバー効果モデル（関数）を定義します。

減衰率（rate）と期間（ length ）がハイパーパラメータで、その値によって形状が異なってきます。

以下、コードです。

# キャリーオーバー効果モデル（関数）の定義
def Carryover(X: np.ndarray, rate, length):
    filter = (
        rate ** np.arange(length + 1)
    ).reshape(-1, 1) 
    convolution = convolve2d(X, filter)
    if length > 0 : convolution = convolution[: -length]
    return convolution

 

最初100でその後0であるシンプルなケースの、この関数の実行例です。

以下、コードです。

# キャリーオーバー効果モデル（関数）の例

## サンプルデータ
exp_dat = pd.DataFrame([100,0,0,0,0,0,0,0,0,0]) 

## キャリーオーバー効果
exp_co = Carryover(exp_dat, 0.5,2) 

print(exp_co) #数値出力
plt.bar(exp_dat.index,exp_co[:,0]) #グラフ
plt.title('rate=0.5, length=2')

 

以下、実行結果です。

 

先程の例よりも複雑なケースで、この関数の実行例を見てみます。

サンプルデータを作ります。

以下、コードです。

# サンプルデータ
exp_dat = pd.DataFrame([10,500,10,500])  #入力データ

print(exp_dat) #数値出力
plt.bar(exp_dat.index,exp_dat.values[:,0]) #グラフ

 

以下、実行結果です。

 

サンプルデータに対し、この関数の実行してみます。

以下、コードです。

# キャリーオーバー効果モデル（関数）の例
exp_co = Carryover(exp_dat, 0.5,2) 

print(exp_co) #数値出力
plt.bar(exp_dat.index,exp_co[:,0]) #グラフ

 

以下、実行結果です。

 

飽和モデルやキャリーオーバー効果モデルを単体で使うわけではありません。

入力データをキャリーオーバー効果モデルにインプットし、キャリーオーバー効果モデルのアウトプットが飽和モデルのインプットになります。飽和モデルのアウトプットは学習器（線形回帰モデル）のインプットになります。

入力データ→飽和モデル→キャリーオーバー効果モデル→学習器のインプット

先程の例で、「飽和モデル→キャリーオーバー効果モデル」を実行してみます。

以下、コードです。

# 「入力→キャリーオーバー効果モデル→飽和モデル→出力」の例
exp_sat_co = Saturation(Carryover(exp_dat, 0.5,2),0.01)

print(exp_sat_co) #数値出力
plt.bar(exp_dat.index,exp_sat_co[:,0]) #グラフ

 

以下、実行結果です。

 

　アドストックを考慮した線形回帰モデル

パイパーパラメータを設定します。

以下、コードです。

# TVCMのハイパーパラメータの設定
TVCM_carryover_rate = 0.5
TVCM_carryover_length = 4
TVCM_saturation_a = 0.000002

# Newspaperのハイパーパラメータの設定
Newspaper_carryover_rate = 0.5
Newspaper_carryover_length = 2
Newspaper_saturation_a = 0.000002

# Webのハイパーパラメータの設定
Web_carryover_rate = 0.5
Web_carryover_length = 0
Web_saturation_a = 0.000002

 

どのようなキャリーオーバー効果なのかを見てみます。

以下、コードです。

# キャリーオーバー効果モデルの出力例

## サンプルデータ
exp_dat = pd.DataFrame([100,0,0,0,0,0,0,0,0,0]) 

## キャリーオーバー効果

### TVCM
exp_co_TVCM= Carryover(
    exp_dat, 
    TVCM_carryover_rate,
    TVCM_carryover_length
) 

### Newspaper
exp_co_Newspaper = Carryover(
    exp_dat,   
    Newspaper_carryover_rate,
    Newspaper_carryover_length
) 

### Web
exp_co_Web = Carryover(
    exp_dat,    
    Web_carryover_rate,
    Web_carryover_length
) 

## グラフ
fig, axes = plt.subplots(nrows=1, ncols=3, sharex=False)

### TVCM
axes[0].bar(exp_dat.index,
            exp_co_TVCM[:,0]
           ) 
axes[0].set_title('TVCM')

### Newspaper
axes[1].bar(exp_dat.index,
            exp_co_Newspaper[:,0]
           ) 
axes[1].set_title('Newspaper')

### Web
axes[2].bar(exp_dat.index,
            exp_co_Web[:,0]
           ) 
axes[2].set_title('Web')

 

以下、実行結果です。

 

説明変数X（広告などのコスト）のデータに対し、キャリーオーバー効果モデルと飽和モデルを使い、学習器（線形回帰モデル）のインプットを作ります。

以下、コードです。

# TVの値の変換
X_TVCM = Saturation(Carryover(X[['TVCM']],
                            TVCM_carryover_rate,
                            TVCM_carryover_length),
                  TVCM_saturation_a)

# Newspaperの値の変換
X_Newspaper = Saturation(Carryover(X[['Newspaper']], 
                               Newspaper_carryover_rate,
                               Newspaper_carryover_length),
                     Newspaper_saturation_a)

# Webの値の変換
X_Web = Saturation(Carryover(X[['Web']], 
                                 Web_carryover_rate,
                                 Web_carryover_length),
                       Web_saturation_a)

# 変換した値の結合（DataFrame型へ）
X_trans = pd.DataFrame(np.concatenate([X_TVCM, 
                                       X_Newspaper,
                                       X_Web], 
                                      1))
X_trans.columns = ['TVCM','Newspaper','Web']

X_trans #確認

 

以下、実行結果です。

 

この説明変数Xを変換し作ったデータを使い、線形回帰モデルを作り精度検証してみます。

線形回帰モデルが、どの程度の予測精度を持ったモデルになるのかを確かめるために、時系列のCV（クロスバリデーション）を実施します。今回は、デフォルトの5分割のCVです。

以下、コードです。

# クロスバリデーションで精度検証（R2）
np.mean(cross_val_score(lr, 
                        X_trans, y, 
                        cv=TimeSeriesSplit()
                       )
       )

 

以下、実行結果です。

 

全データでモデルを構築し、予測精度（ R^2）を見てみます。

以下、コードです。

# 全データで精度検証（R2）
lr.fit(X_trans, y)
lr.score(X_trans, y)

 

以下、実行結果です。

 

前回な線形回帰モデルに対し、以下のように（ R^2）が変化しました。

  • CVの R^2の平均値：0.70→0.78
  • 全データ利用した場合の R^2：0.74→0.83

 

　アドストックを考慮した線形回帰モデル
（Optunaでハイパーパラメータチューニング）

変換器と学習器を繋いだパイプラインを作り、Optunaでハイパーパラメータ探索を実施し最適なモデルを作ります。

ハイパーパラメータは、変換器の以下の3つです。

  • キャリーオーバー効果モデル
    • 減衰率（rate）
    • 期間（ length ）
  • 飽和モデル
    • 形状パラメータa

 

　　パイプライン構築

飽和モデルとキャリーオーパー効果モデルを、パイプラインで利用できる変換器にします。

以下、コードです。

# Pipeline用変換器（飽和モデル）
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)

# Pipeline用変換器（キャリーオーバー効果モデル）
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の構築

## 説明変数の変換部分（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', LinearRegression())
])

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

 

以下、実行結果です。

 

　　パイパーパラメータ探索の実施

パイプラインを構築したので、パイパーパラメータを探索します。

以下、コードです。

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

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

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

 

探索結果を見てみます。

以下、コードです。

# 探索結果
optuna_search.best_params_

 

以下、実行結果です。

 

探索し得られたハイパーパラメータで、どのようなキャリーオーバー効果になるのかを見てみます。

以下、コードです。

# キャリーオーバー効果モデルの出力例

## サンプルデータ
exp_dat = pd.DataFrame([100,0,0,0,0,0,0,0,0,0]) 

## ハイパーパラメータ設定

### TVCM
TVCM_carryover_rate = optuna_search.best_params_['adstock__TVCM_pipe__carryover__rate']
TVCM_carryover_length = optuna_search.best_params_['adstock__TVCM_pipe__carryover__length']

### Newspaper
Newspaper_carryover_rate = optuna_search.best_params_['adstock__Newspaper_pipe__carryover__rate']
Newspaper_carryover_length = optuna_search.best_params_['adstock__Newspaper_pipe__carryover__length']

### Web
Web_carryover_rate = optuna_search.best_params_['adstock__Web_pipe__carryover__rate']
Web_carryover_length = optuna_search.best_params_['adstock__Web_pipe__carryover__length']

## キャリーオーバー効果

### TVCM
exp_co_TVCM= Carryover(
    exp_dat, 
    TVCM_carryover_rate,
    TVCM_carryover_length
) 

### Newspaper
exp_co_Newspaper = Carryover(
    exp_dat,   
    Newspaper_carryover_rate,
    Newspaper_carryover_length
) 

### Web
exp_co_Web = Carryover(
    exp_dat,    
    Web_carryover_rate,
    Web_carryover_length
) 

## グラフ
fig, axes = plt.subplots(nrows=1, ncols=3, sharex=False)

### TVCM
axes[0].bar(exp_dat.index,
            exp_co_TVCM[:,0]
           ) 
axes[0].set_title('TVCM(y:logarithmic scale)')
axes[0].set_yscale('log')
axes[0].set_ylim(0, 100)

### Newspaper
axes[1].bar(exp_dat.index,
            exp_co_Newspaper[:,0]
           ) 
axes[1].set_title('Newspaper')

### Web
axes[2].bar(exp_dat.index,
            exp_co_Web[:,0]
           ) 
axes[2].set_title('Web')

 

以下、実行結果です。TVCMだけ対数変換し表示しています。

 

　　最適なハイパーパラメータでパイプラインを学習

パイプラインのインスタンスを作ります。

以下、コードです。

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

 

線形回帰モデルが、どの程度の予測精度を持ったモデルになるのかを確かめるために、時系列のCV（クロスバリデーション）を実施します。今回は、デフォルトの5分割のCVです。

以下、コードです。

# クロスバリデーションで精度検証（R2）
np.mean(cross_val_score(MMM_pipe_best, 
                        X, y, 
                        cv=TimeSeriesSplit()
                       )
       )

 

以下、実行結果です。

 

全データでモデルを構築し、予測精度（ R^2）を見てみます。

以下、コードです。

# 全データで精度検証（R2）
MMM_pipe_best.fit(X, y)
MMM_pipe_best.score(X, y)

 

以下、実行結果です。

 

前回の線形回帰モデルに比べ、以下のように（ R^2）が変化しました。

  • CVの R^2の平均値：0.70→0.82
  • 全データ利用した場合の R^2：0.74→0.84

どのような線形回帰モデルなのか（切片と回帰係数）を見てみます。

以下、コードです。

# 線形回帰モデルの切片と回帰係数
intercept = MMM_pipe_best.named_steps['regression'].intercept_ #切片
coef = MMM_pipe_best.named_steps['regression'].coef_ #回帰係数

# 回帰係数をデータフレーム化
weights = pd.Series(
    coef,
    index=X.columns
)

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

 

以下、実行結果です。

 

　　売上貢献度の算出

学習し構築したパイプラインの変換器（adstock）を抽出します。

以下、コードです。

# Piplineの変換器（adstock）を抽出
adstock = MMM_pipe_best.named_steps['adstock']

 

説明変数を変換し、学習データのインプットデータを作ります。

以下、コードです。

# 説明変数Xの変換
X_trans = pd.DataFrame(adstock.transform(X),
                       index=X.
                       index,columns=X.columns)

 

先程求めた線形回帰式を使い、売上貢献度（補正前）を計算します。

以下、コードです。

# 貢献度（補正前）
unadj_contribution = X_trans.mul(weights) #Xと係数を乗算
unadj_contribution = unadj_contribution.assign(Base=intercept) #切片の追加

unadj_contribution.head() #確認

 

以下、実行結果です。

 

週ごとに各媒体の売上貢献度を合計すると、売上の予測値になります。

以下、コードです。

# 貢献度の合計（yの予測値）
y_pred = unadj_contribution.sum(axis=1)

y_pred.head() #確認

 

以下、実行結果です。

 

元の売上の実測値を見てみます。

以下、コードです。

y.head() #確認

 

以下、実行結果です。

 

予測値と実測値が乖離していることが分かります。この乖離をなくすために、補正係数（correction factor）を計算し、売上貢献度を補正します。

補正係数を計算します。

以下、コードです。

# 補正係数
correction_factor = y.div(y_pred, axis=0)

correction_factor.head() #確認

 

以下、実行結果です。

 

この補正係数を使い、売上貢献度を補正します。

以下、コードです。

# 貢献度（補正後）
adj_contribution = (unadj_contribution
                    .mul(correction_factor, axis=0)
                   ) 

# 順番の変更
adj_contribution = adj_contribution[['Base', 'Web', 'Newspaper', 'TVCM']]

#確認
adj_contribution.head()

 

以下、実行結果です。

 

週×媒体別の売上貢献度が求められたので、積み上げグラフを作成し、どのような状況になっているのかを確認してみます。

以下、コードです。

# グラフ化
ax = (adj_contribution
      .plot.area(
          figsize=(16, 10),
          linewidth=0,
          ylabel='Sales',
          xlabel='Date')
     )

handles, labels = ax.get_legend_handles_labels()
ax.legend(handles[::-1], labels[::-1])

 

以下、実行結果です。

 

媒体別に全ての週の売上貢献度を合計し、媒体別の売上貢献度（円と構成比％）と、そのグラフを作り、何がどれほど売上に貢献したのかを見てみます。

以下、コードです。

# 媒体別の貢献度の合計
contribution_sum = adj_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を計算します。

先ず、媒体別のコストを集計します。

以下、コードです。

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

cost_sum #確認

 

以下、実行結果です。

 

先程求めた売上貢献度を使い、媒体別のROIを計算しグラフ化します。

以下、コードです。

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

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

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

 

以下、実行結果です。

 

ROIは、値が大きいほど良く、最低限プラスの値である必要があります。少なくとも0以上である必要があります。

このモデルは、アドストック（Ad Stock）を考慮していますが、キャリーオーバー効果モデルの作り方から分かりますが、最初に効果のピークが来て徐々に効果が減衰するモデルです。

現実はそうではなく、効果のピークが最初に来ることもあれば、数日後や数週間後にピークが来ることもあります。そのようなモデルの方がいいでしょう。

今回は最適なハイパーパラメータを1つ探索するということをしましたが、別のハイパーパラメータの組み合わせでも同じ程度のレベルの予測精度になる可能性もあります。

ハイパーパラメータを1つとするのではなく、ある分布の代表値として見なす考え方もあります。そのためには、ベイズモデルでMMMを構築する必要があります。

 

次回

今回は、「アドストック（Ad Stock）を考慮した線形回帰モデル」というお話しをしました。

次回は、「ちょっと複雑なアドストック（Ad Stock）を考慮した線形回帰モデル」というお話しをします。

Pythonによるマーケティングミックスモデリング（MMM：Marketing Mix Modeling）超入門 その3ちょっと複雑なアドストック（Ad Stock）を考慮した線形回帰モデル