データを活用したマーケティング戦略は、ビジネスの成功に不可欠です。その中心に位置するのが、マーケティングミックスモデリング（MMM）です。

マーケティングミックスモデリング（MMM）は、過去のデータを分析する「振り返り分析」と未来のトレンドを予測する「近未来分析」の両方で有効に活用される強力なツールです。

前回は、「振り返り分析」のための「線形回帰系モデルでのMMM構築」というお話しを、Pythonコードを交えお話ししました。

PythonによるMMM（マーケティングミックスモデリング）とビジネス活用－ 振返り分析編（その2）：線形回帰系モデルでのMMM構築 －

 

MMMの構築と後処理（貢献度やマーケティングROIの計算など）までの流れがざっくり理解できたのではないかと思います。

ただし、前回紹介したモデルには、重大な欠陥があります。

それはアドストック（Ad Stock）が考慮されていないことです。

 

ここでは、キャリーオーバー効果や飽和効果のことを、アドストック（Ad Stock）と言っています。例えば……

• キャリーオーバー効果とは「その日の広告の効果が次の日以降も継続していること」
• 飽和効果とは「その広告の投下量の増加が売上などにもたらす効果が、ある限界に達すると、それ以上の増加がほとんどないこと」

通常のMMMは、アドストック（Ad Stock）が考慮されます。ただし、アドストック（Ad Stock）にも色々なアドストック（Ad Stock）があります。

 

今回は、シンプルなアドストックを組み込んだMMMを取り上げます。

ちなみに、広告販促系のデータは、お互いの相関関係が強い傾向（同じ時期に集中し広告投下など）があり、多重共線性問題が起こる可能性が低くはありません。そのため今回は、推定器としてRidge回帰を利用します。

利用するデータセット

前回同様、以下のデータセットを前提に説明します。

• 目的変数：売上金額（Sales）
• 説明変数：TVCM、Newspaper、Webのコスト

今回利用するデータセットは、以下からダウンロードできます。

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

 

目的変数は、売上金額である必要はありません。例えば、売上点数でも受注件数、問い合せ件数、予約件数でも構いません。要は、広告販促でより良い変化をさせたい何かです。

説明変数は、コストでなくても構いません。例えば、媒体の認知や接触、露出、広告のクリック、インプレッション、配信、開封などでも構いません。ただ、コスパという概念を持ち出すとき、少なくともコストデータは必要になります。

 

MMMパイプラインとアドストック

　パイプラインの構造

MMMにアドストック（Ad Stock）を組み込むには、パイプラインを構築するといいでしょう。このようなパイプラインを、ここではMMMパイプラインと表現することにします。

このMMMパイプラインは、少なくとも2つの変換器と1つの推定器が必要になります。

• 変換器
  • キャリーオーバー効果関数（Carryover）
  • 飽和関数（Saturation）
• 推定器
  • 線形回帰系のモデル（Linear Regression）：今回はRidge回帰を利用

 

各広告販促の投下量のデータごとに（要は、各変数ごとに）、その広告販促に応じた変換（キャリーオーバー効果関数→飽和関数）を施し、その変換結果を用い推定器である線形系のモデルを構築し、売上を表現するモデルを作り上げます。

前回との大きな違いは、変換（キャリーオーバー効果関数→飽和関数）を施しすところにあります。

 

ここでは、パイプラインそのものに関する詳しい説明はいたしません。

パイプラインの概要や作り方に興味のある方は、以下を参考にしてください。

scikit-learnの機械学習パイプライン入門

 

では、今回利用するシンプルなキャリーオーバー効果関数と飽和関数について簡単に説明します。

 

　キャリーオーバー効果関数

広告販促などには通常、キャリーオーバー効果（Carryover）と呼ばれるものがあります。効果が後々まで残っているというものでラグ効果や残存効果とも言われます。

上のグラフは、広告販促の効果が、時間がたつにつれて徐々に薄れながらも残存して残っているイメージです。縦軸は効果で、横軸は時間軸です。詳細は、後で説明します。

キャリーオーバー効果を表現する関数には、色々あります。

今回は最もシンプルな、広告などを打ったときが効果のピークで、徐々に効果が一定の割合で減衰していくモデルを使います。このようなキャリーオーバー効果関数を、シンプルなキャリーオーバー効果関数とここでは呼びます。

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

シンプルなキャリーオーバー効果関数

\displaystyle x_t^* = \frac{\displaystyle \sum_{l=0}^{L-1} w_{t-l} \cdot x_{t-l}}{\displaystyle \sum_{l=0}^{L-1} w_{t-l}} , \quad w_{t-l} = R^l

 

x_t は t 期の広告などの投入量で、x_t^* は t 期とそれ以前までの広告などの効果の累積（残存効果を足したもの）です。

このように、キャリーオーバー効果関数は、単にキャリーオーバーの現象をモデル化するだけでなく、それを足し合わせるような関数です。

\displaystyle x_t \overset{Carryover}{\Rightarrow} x_t^{*}

 

この関数には、以下の2つのハイパーパラメータがあります。

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

 

この効果の続く期間 L(length)とは、どのくらいまで考慮するか、ということです。

例えば、広告などを投下した期のみ効果がでるなら、L=1になります。3日間効果がでるなら、L=3になります。

 

この減衰率 R(rate)とは、時間の経過とともにどの程度減るのかを表現したものです。

例えば、減衰率R=0.5で広告などを投下した期の効果が1とした場合、次の期にキャリーオーバーする効果（残存する効果）は0.5、次の次の期の効果にキャリーオーバーする効果（残存する効果）は0.25という感じになります。

 

Pythonのコードでこの関数を表現すると以下のようになります。

# キャリオーバー効果関数
def carryover_simplified(X, length, rate):
    X = np.append(np.zeros(length - 1), X)
    Ws = np.zeros(length)
    for l in range(length):
        Ws[length - 1 - l] = rate ** l
    carryover_X = []
    for i in range(length - 1, len(X)):
        X_array = X[i - length + 1: i + 1]
        Xi = sum(X_array * Ws) / sum(Ws)
        carryover_X.append(Xi)
    return np.array(carryover_X)

 

このコードを簡単に説明します。

• def carryover_simplified(X, length, rate):
  • この行は関数の定義を開始します。carryover_simplifiedという名前の関数は、三つの引数を取ります：
    • X：入力データの配列。
    • length：キャリオーバー効果を計算するための期間の長さ。
    • rate：キャリオーバーの減衰率。
• X = np.append(np.zeros(length - 1), X)
  • 入力配列Xの先頭にlength - 1個のゼロを追加します。
  • これは、計算の開始時にキャリオーバー効果がないことを示すためです。
• Ws = np.zeros(length)
  • 長さがlengthのゼロ配列Wsを作成します。
  • これは後で減衰係数を格納するために使用されます。
• for l in range(length): lengthの長さ分だけループを回します。
  • Ws[length - 1 - l] = rate ** l
    • 減衰率rateのl乗をWs配列の対応する位置に格納します。
    • これにより、時間の経過に伴う減衰効果が計算されます。
• carryover_X = []
  • 結果を格納するための空のリストを作成します。
• for i in range(length - 1, len(X)):入力配列Xをループし、キャリオーバー効果を計算します。
  • X_array = X[i - length + 1: i + 1]
    • 現在の位置からlength期間分のデータを取り出します。
  • Xi = sum(X_array * Ws) / sum(Ws)
    • 取り出したデータと減衰係数を掛け合わせ、その合計を減衰係数の合計で割ります。
    • これにより、その時点でのキャリオーバー効果が計算されます。
  • carryover_X.append(Xi)
    • 計算されたキャリオーバー効果をリストに追加します。
• return np.array(carryover_X)
  • 最終的なキャリオーバー効果の配列を返します。

 

この関数はパイプラインで利用するとき、クラスにして利用します。

以下、コードです。

# キャリオーバー効果を適用するカスタム変換器クラス
class CustomCarryOverTransformer(BaseEstimator, TransformerMixin):
    def __init__(self, carryover_params=None):
        self.carryover_params = carryover_params if carryover_params is not None else []

    def fit(self, X, y=None):
        return self

    def transform(self, X):
        if isinstance(X, pd.DataFrame):
            X = X.to_numpy()
        transformed_X = np.copy(X)
        for i, params in enumerate(self.carryover_params):
            transformed_X[:, i] = carryover_simplified(X[:, i], **params)
        return transformed_X

 

このコードを簡単に説明します。

• class CustomCarryOverTransformer(BaseEstimator, TransformerMixin):
  • この行でCustomCarryOverTransformerクラスを定義し、Scikit-LearnのBaseEstimatorとTransformerMixinクラスを継承しています。
• def __init__(self, carryover_params=None):
  • 初期化メソッドです。オプショナルな引数carryover_paramsを受け取ります。
  • この引数はキャリオーバー効果を計算する際に使用されるパラメータを含むリストです。
• self.carryover_params = carryover_params if carryover_params is not None else []
  • carryover_paramsが提供されていない場合、デフォルトで空のリストを使用します。
• def fit(self, X, y=None):
  • fitメソッドは、この変換器に特定の学習データを”フィット”させるために使用されます。
  • このクラスでは、特別な学習処理は行われないため、自身のインスタンスをそのまま返します。
• def transform(self, X):
  • transformメソッドは、入力データXに変換処理を適用します。
  • if isinstance(X, pd.DataFrame):
    • 入力がPandasのDataFrameである場合、それをNumPy配列に変換します。
  • transformed_X = np.copy(X)
    • 入力データのコピーを作成します。
    • これに変換処理を適用します。
  • for i, params in enumerate(self.carryover_params):
    • carryover_paramsの各変数に対してループを行い、それぞれのパラメータセットを使用して変換を行います。
    • transformed_X[:, i] = carryover_simplified(X[:, i], **params)
      • 先に定義されたcarryover_simplified関数を使用して、指定された列に対してキャリオーバー効果を適用します。
• return transformed_X
  • 変換されたデータを返します。

一見するとややこしそうに見えますが、実は非常に単純です。変換を実施するtransformメソッドに、先ほど定義したキャリオーバー効果関数を設定しているだけです。別のキャリーオーバー効果関数を定義した場合も、transformメソッドに設定する関数を入れ替えることで流用できます。

 

　飽和関数

広告販促の投下量を増やせば増やすほど売上などは上昇します。しかし、売上の上昇幅は鈍くなります。経済学でいうところの収穫逓減が起こります。

上のグラフは、横軸は広告などの投入量で、縦軸は効果の大きさです。要は、ある時点を境に売上は飽和し、いくら投下しても売上が伸びなくなるということです。

飽和を表現する関数には、色々あります。

今回は、最もシンプルな指数関数exp(*)で表現した飽和関数を利用します。

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

シンプルな飽和関数（指数型）

\displaystyle x_t^{**} = 1-exp(-dx_t^*)

 

x_t^* は、 キャリーオーバー効果を考慮したt 期の広告などの説明変数（特徴量）の値です。つまり、キャリーオーバー効果関数で変換した後の説明変数（特徴量）の値です。x_t^{**} は、この x_t^* に対し 飽和効果を考慮した値（変換した値）です。

\displaystyle x_t \overset{Carryover}{\Rightarrow} x_t^{*} \overset{Saturation}{\Rightarrow} x_t^{**}

このx_t^{**} を利用し、推定器（今回はRidge回帰）のインプットとして利用します。

 

この関数には、飽和関数の形状を表現する、以下の1つのハイパーパラメータがあります。

• d：形状パラメータ

 

Pythonのコードでこの関数を表現すると以下のようになります。

# 飽和関数（指数型）
def exponential_function(x, d):
    result = 1 - np.exp(-d * x)
    return result

 

このコードを簡単に説明します。

• def exponential_function(x, d):
  • 関数の定義です。ここではexponential_functionという名前で、二つの引数xとdを取ります。
  • xは関数の入力値です。
  • dは関数の形状を決定するパラメータです。通常は正の値を取ります。
• result = 1 - np.exp(-d * x)
  • NumPyの指数関数np.exp()を使用して、-d * xの指数を計算します。
  • この式は、xが増加するにつれて、結果が1に近づくが、1を超えることはないという特性を持っています。これは飽和特性を表しています。
  • dの値が大きいほど、飽和に達する速度が速くなります。
• return result
  • 計算された結果を返します。

 

この関数はパイプラインで利用するとき、クラスにして利用します。

以下、コードです。

# 飽和関数を適用するカスタム変換器クラス
class CustomSaturationTransformer(BaseEstimator, TransformerMixin):
    def __init__(self, curve_params=None):
        self.curve_params = curve_params if curve_params is not None else []
    def fit(self, X, y=None):
        return self
    def transform(self, X):
        if isinstance(X, pd.DataFrame):
            X = X.to_numpy()
        transformed_X = np.copy(X)
        for i, params in enumerate(self.curve_params):
            transformed_X[:, i] = exponential_function(X[:, i], **params)
        return transformed_X

 

このコードを簡単に説明します。

• class CustomSaturationTransformer(BaseEstimator, TransformerMixin):
  • この行ではCustomSaturationTransformerクラスを定義し、Scikit-LearnのBaseEstimatorとTransformerMixinクラスを継承しています。
• def __init__(self, curve_params=None):
  • 初期化メソッドです。オプショナルな引数curve_paramsを受け取ります。
  • この引数は飽和関数を計算する際に使用されるパラメータを含むリストです。
• self.curve_params = curve_params if curve_params is not None else []
  • curve_paramsが提供されていない場合、デフォルトで空のリストを使用します。
• def fit(self, X, y=None):
  • fitメソッドは、この変換器に特定の学習データを”フィット”させるために使用されます。
  • このクラスでは、特別な学習処理は行われないため、自身のインスタンスをそのまま返します。
• def transform(self, X):
  • transformメソッドは、入力データXに変換処理を適用します。
  • if isinstance(X, pd.DataFrame):
    • 入力がPandasのDataFrameである場合、それをNumPy配列に変換します。
  • transformed_X = np.copy(X)
    • 入力データのコピーを作成します。
    • これに変換処理を適用します。
  • for i, params in enumerate(self.curve_params):
    • curve_paramsの各要素に対してループを行い、それぞれのパラメータセットを使用して変換を行います。
  • transformed_X[:, i] = exponential_function(X[:, i], **params)
    • 事前に定義されたexponential_function関数を使用して、指定された列に対して飽和関数を適用します。
• return transformed_X
  • 変換されたデータを返します。

こちらも一見するとややこしそうに見えますが、実は非常に単純です。変換を実施するtransformメソッドに、先ほど定義した飽和関数を設定しているだけです。別の飽和関数を定義した場合も、transformメソッドに設定する関数を入れ替えることで流用できます。

 

準備

今回利用するモジュールやデータセットの読み込み、先ほど述べたアドストック系関数の定義などを実施していきます。

 

　モジュールなどの読み込み

　　共通利用するモジュール

以下、コードです。

import numpy as np
import pandas as pd
import optuna

from sklearn.linear_model import Ridge

from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import TimeSeriesSplit
from sklearn.model_selection import cross_val_score

from functools import partial

import warnings
warnings.simplefilter('ignore')

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

 

簡単に説明します。

• import numpy as np
  • NumPyは、数値計算を効率的に行うための基本的なライブラリ。
• import pandas as pd
  • Pandasは、データ構造とデータ分析ツールを提供するライブラリ。
• import optuna
  • Optunaは、機械学習モデルのハイパーパラメータを最適化するためのライブラリ。
• from sklearn.linear_model import Ridge
  • scikit-learnは機械学習のためのライブラリで、ここでは線形モデルの一つであるリッジ回帰（Ridge）をインポート。
• from sklearn.base import BaseEstimator, TransformerMixin
  • scikit-learnの基本クラスであるBaseEstimatorと、変換器のためのミックスインクラスであるTransformerMixin。
  • クラスを作るときに利用。
• from sklearn.pipeline import Pipeline
  • データ処理とモデルの構築を一連のステップとして管理するためのパイプライン機能。
• from sklearn.preprocessing import MinMaxScaler
  • MinMaxScalerはデータの正規化を行うためのツール。
  • 最小値が0、最大値が1になりように正規化。
• from sklearn.model_selection import TimeSeriesSplit
  • TimeSeriesSplitは時系列クロスバリデーションで利用。
• from sklearn.model_selection import cross_val_score
  • モデルの評価を行うためのcross_val_score関数をインポート。
• from functools import partial
  • 関数に事前に引数を設定するためのpartialをインポート。
• import warnings
  • Pythonの警告制御ツールをインポート。
  • warnings.simplefilter('ignore')
    • 警告を無視するように設定。
• import matplotlib.pyplot as plt
  • データの可視化を行うためのライブラリであるMatplotlibをインポートします。
  • plt.style.use('ggplot')
    • グラフのスタイルをggplot風に設定。
  • plt.rcParams['figure.figsize'] = [12, 7]
    • グラフのデフォルトサイズを設定。
  • plt.rcParams['font.size'] = 9
    • グラフのフォントサイズを設定。

 

　　前回作った関数群（MMM構築・評価など）

さらに今回は、前回作った以下の関数群を使います。

• MMM構築
  • train_and_evaluate_model:MMMの構築
  • plot_actual_vs_predicted:構築したMMMの予測値の時系列推移
• 後処理（結果の出力）
  • calculate_and_plot_contribution:売上貢献度の算出（時系列推移）
  • summarize_and_plot_contribution:売上貢献度構成比の算出
  • calculate_marketing_roi:マーケティングROIの算出

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

mmm_functions.py　※zipファイルを解凍してお使いください
https://www.salesanalytics.co.jp/oi17

以下、mmm_functions.pyの中に記載されているコードです。

import numpy as np
import pandas as pd

from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_absolute_percentage_error
from sklearn.metrics import r2_score

import warnings
warnings.simplefilter('ignore')

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

# MMM構築

def train_and_evaluate_model(model, X, y):
    """
    モデルを学習し、予測と評価を行う関数。
    
    :param model: 学習するモデルのインスタンス
    :param X: 特徴量のデータフレーム
    :param y: ターゲット変数のデータフレーム
    :return: 学習済みモデルmodel、予測値pred
    """
    # モデルの学習
    model.fit(X, y)

    # 予測
    pred = pd.DataFrame(
        model.predict(X),    
        index=X.index,
        columns=['y'],
    )

    # 精度指標の計算
    rmse = np.sqrt(mean_squared_error(y, pred))
    mae = mean_absolute_error(y, pred)
    mape = mean_absolute_percentage_error(y, pred)
    r2 = r2_score(y, pred)

    # 精度指標の出力
    print('RMSE:', rmse)
    print('MAE:', mae)
    print('MAPE:', mape)
    print('R2:', r2)

    return model,pred

# 実測値と予測値の時系列推移

def plot_actual_vs_predicted(index, actual_values, predicted_values, ylim):
    """
    実際の値と予測値を比較するグラフを描画する関数。

    :param index: データのインデックス
    :param actual_values: 実際の値の配列
    :param predicted_values: 予測値の配列
    :param ylim: y軸の表示範囲
    """
    fig, ax = plt.subplots()

    ax.plot(index, actual_values, label="actual")
    ax.plot(index, predicted_values, label="predicted", linestyle="dotted", lw=2)

    ax.set_ylim(ylim)

    plt.title('Time series of actual and predicted values')
    plt.legend()
    plt.show()

# 貢献度の算出

def calculate_and_plot_contribution(y, X, model, ylim):
    """
    各媒体の売上貢献度を算定し、結果をプロットする関数。

    :param y: ターゲット変数
    :param X: 特徴量のデータフレーム
    :param model: 学習済みモデル
    :param ylim: y軸の表示範囲
    :return: 各媒体の貢献度
    """
    # yの予測
    pred = pd.DataFrame(
        model.predict(X),    
        index=X.index,
        columns=['y'],
    )

    # 値がすべて0の説明変数
    X_ = X.copy()
    X_.iloc[:, :] = 0
    
    # Baseの予測
    base = model.predict(X_)
    pred['Base'] = base

    # 各媒体の予測
    for i in range(len(X.columns)):
        X_.iloc[:, :] = 0
        X_.iloc[:, i] = X.iloc[:, i]
        pred[X.columns[i]] = model.predict(X_) - base

    # 予測値の補正
    correction_factor = y.div(pred.y, axis=0)
    pred_adj = pred.mul(correction_factor, axis=0)
    contribution = pred_adj.drop(columns=['y'])

    # 結果の出力
    print(contribution, '\n')

    # エリアプロット
    ax = contribution.plot.area()
    handles, labels = ax.get_legend_handles_labels()
    ax.legend(reversed(handles), reversed(labels))
    ax.set_ylim(ylim)
    plt.title('Contributions over time')
    plt.show()

    return contribution

# 貢献度構成比の算出

def summarize_and_plot_contribution(contribution):
    """
    媒体別の売上貢献度の合計と構成比を計算し、結果を表示する関数。

    :param contribution: 各媒体の貢献度を含むデータフレーム
    :return: 売上貢献度の合計と構成比を含むデータフレーム
    """
    # 各媒体の貢献度の合計
    contribution_sum = contribution.sum(axis=0)

    # 各媒体の貢献度の構成比
    contribution_percentage = contribution_sum / contribution_sum.sum()

    # 結果を1つのDataFrameにまとめる
    contribution_results = pd.DataFrame({
        'contribution': contribution_sum,
        'ratio': contribution_percentage
    })

    # 結果の出力
    print(contribution_results, '\n')

    # グラフ化
    contribution_sum.plot.pie()
    plt.title('Contribution Composition Ratio')
    plt.show()

    return contribution_results

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

def calculate_marketing_roi(X, contribution):
    """
    各媒体のマーケティングROIを算定する関数。

    :param X: 各媒体のコストを含むデータフレーム
    :param contribution: 各媒体の売上貢献度を含むデータフレーム
    :return: 各媒体のROIを含むデータフレーム
    """
    # 各媒体のコストの合計
    cost_sum = X.sum(axis=0)

    # 各媒体のROIの計算
    ROI = (contribution[X.columns].sum(axis=0) - cost_sum)/cost_sum

    # 結果の出力
    print(ROI, '\n')   

    # グラフ
    ROI.plot.bar()
    plt.title('Marketing ROI')
    plt.show()

    return ROI

 

mmm_functions.pyを利用するときは、実行するPythonファイルやNotebookと同じフォルダに入れておいてください。

以下のコードで呼び出せます。

from mmm_functions import *

 

上手くいかないときは、mmm_functions.pyをメモ帳などで開き内容をコピーし、実行するPythonファイルやNotebookにコードを張り付け、Pythonで関数を作ってからMMM構築などを行ってください。

 

　データセットの読み込み

以下、コードです。

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

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

 

• データセット読み込み:
  • dataset = 'https://www.salesanalytics.co.jp/4zdt': ここで、データセットのURLが指定されています。
  • df = pd.read_csv(...): この行では、Pandasライブラリのread_csv関数を使用して、指定されたURLからデータセットを読み込んでいます。この関数は、CSV形式（カンマ区切りの値）のデータを読み込むのに広く使われています。
    • parse_dates=['Week']: この引数は、’Week’列を日付型のデータとして解析するように指示しています。
    • index_col='Week': この引数は、’Week’列をデータフレームのインデックス（行のラベル）として使用するように指示しています。
• 説明変数Xと目的変数yに分解:
  • X = df.drop(columns=['Sales']): ここでは、dropメソッドを使用して’Sales’列を除外し、残りの列を説明変数Xとして保持しています。dropは指定された列をデータフレームから削除するために使われます。
  • y = df['Sales']: この行では、’Sales’列を目的変数yとして抽出しています。目的変数は、モデルによって予測されるべき変数であり、このケースでは売上データと思われます。

 

目的変数y（売上）をプロットします。

以下、コードです。

y.plot()
plt.show()

 

以下、実行結果です。

 

• y.plot(): y（目的変数）を折れ線グラフで表示します。plot()はPandasのデフォルトのグラフ化するためのプロット機能です。
• plt.show(): 作成されたグラフを画面に出力します。

 

説明変数X（各メディアのコスト）をまとめてプロットします。

以下、コードです。

X.plot(subplots=True)
plt.show()

 

• X.plot(subplots=True): X（説明変数）を折れ線グラフで表示します。subplots=Trueは、Xの各列（各媒体）を別々のサブプロットとして表示することを意味します。これにより、各説明変数の動きを個別に確認できます。
• plt.show(): 作成されたグラフを画面に出力します。

 

以下、実行結果です。

 

　アドストック系の関数とクラス、そのグラフ化関数の定義

先ほどと重複しますが、アドストック系の関数である、シンプルなキャリーオーバー効果関数とそのクラス、シンプルな飽和関数とそのクラス、そしてそれらを視覚化（グラフ）する関数を作っていきます。

 

　　キャリオーバー効果関数

先ずは、シンプルなキャリーオーバー効果関数とそのクラスを作ります。

以下、コードです。

# キャリオーバー効果関数
def carryover_simplified(X, length, rate):
    X = np.append(np.zeros(length - 1), X)
    Ws = np.zeros(length)
    for l in range(length):
        Ws[length - 1 - l] = rate ** l
    carryover_X = []
    for i in range(length - 1, len(X)):
        X_array = X[i - length + 1: i + 1]
        Xi = sum(X_array * Ws) / sum(Ws)
        carryover_X.append(Xi)
    return np.array(carryover_X)

# キャリオーバー効果を適用するカスタム変換器クラス
class CustomCarryOverTransformer(BaseEstimator, TransformerMixin):
    def __init__(self, carryover_params=None):
        self.carryover_params = carryover_params if carryover_params is not None else []

    def fit(self, X, y=None):
        return self

    def transform(self, X):
        if isinstance(X, pd.DataFrame):
            X = X.to_numpy()
        transformed_X = np.copy(X)
        for i, params in enumerate(self.carryover_params):
            transformed_X[:, i] = carryover_simplified(X[:, i], **params)
        return transformed_X

 

　　飽和関数

次に、シンプルな飽和関数とそのクラスを作ります。

以下、コードです。

# 飽和関数（指数型）
def exponential_function(x, d):
    result = 1 - np.exp(-d * x)
    return result

# 飽和関数を適用するカスタム変換器クラス
class CustomSaturationTransformer(BaseEstimator, TransformerMixin):
    def __init__(self, curve_params=None):
        self.curve_params = curve_params if curve_params is not None else []
    def fit(self, X, y=None):
        return self
    def transform(self, X):
        if isinstance(X, pd.DataFrame):
            X = X.to_numpy()
        transformed_X = np.copy(X)
        for i, params in enumerate(self.curve_params):
            transformed_X[:, i] = exponential_function(X[:, i], **params)
        return transformed_X

 

　　アドストック視覚化（グラフ）関数

どのようなアドストックになっているのかを視覚的に理解するために、視覚化（グラフ）する関数を作ります。

以下、コードです。

def plot_carryover_effect(params, feature_name, fig, axes, i):
    max_length = max(10, params['length'])
    x = np.concatenate(([1], np.zeros(max_length - 1)))
    y = carryover_simplified(x, **params)
    y = y / max(y)
    axes[2*i].bar(np.arange(1, max_length + 1), y)
    axes[2*i].set_title(f'Carryover Effect for {feature_name}')
    axes[2*i].text(0, 1.1, params, ha='left',va='top')
    axes[2*i].set_xlabel('Time Lag')
    axes[2*i].set_ylabel('Effect')
    axes[2*i].set_xticks(range(len(y)))
    axes[2*i].set_ylim(0, 1.1)

def plot_saturation_curve(params, feature_name, fig, axes, i):
    x = np.linspace(-1, 3, 400)
    y = exponential_function(x, **params) 
    axes[2*i+1].plot(x, y, label=feature_name)
    axes[2*i+1].set_title(f'Saturation Curve for {feature_name}')
    axes[2*i+1].text(-1, max(y)* 1.1, params, ha='left',va='top')
    axes[2*i+1].set_xlabel('X')
    axes[2*i+1].set_ylabel('Transformation')
    axes[2*i+1].set_ylim(0, max(y) * 1.1)
    axes[2*i+1].set_xlim(-1, 3)

def plot_effects(carryover_params, curve_params, feature_names):
    fig, axes = plt.subplots(len(feature_names) * 2, 1, figsize=(12, 10*len(feature_names)))
    for i, params in enumerate(carryover_params):
        plot_carryover_effect(params, feature_names[i], fig, axes, i)
    for i, params in enumerate(curve_params):
        plot_saturation_curve(params, feature_names[i], fig, axes, i)
    plt.tight_layout()
    plt.show()

 

このコードは、キャリオーバー効果と飽和曲線をプロットするための関数plot_effectsを定義しています。広告や販促などの効果が時間の経過によってどのように変化するか（キャリオーバー効果）、また、ある入力に対する反応が一定のポイントを超えると増加が鈍化する（飽和）様子を視覚化するのに使用します。

大きく、キャリーオーバー効果のプロット部分と、飽和曲線のプロットの部分に分かれます。

このコードは、特定の特徴量に対するキャリーオーバー効果（時間の経過とともにどのように影響が持続するか）と飽和曲線（特定の入力値に対する応答の変化）を可視化するためのものです。

キャリーオーバー効果関数をグラフ化するplot_carryover_effect 関数と、飽和関数をグラフ化するplot_saturation_curve 関数、この2つの関数を統合しプロットするplot_effects 関数で構成されています。

簡単に説明します。

• plot_carryover_effect 関数
  • 目的:
    • 特定の特徴量に対するキャリーオーバー効果を可視化する。
  • パラメータ:
    • params: キャリーオーバー効果を計算するためのパラメータが含まれている辞書。
    • feature_name: 可視化する特徴量の名前。
    • fig: matplotlibのFigureオブジェクト。複数のプロットをまとめるコンテナです。
    • axes: matplotlibのaxesオブジェクトの配列。個々のプロットを描画するために使用されます。
    • i: 特徴量のインデックス。このインデックスに基づいて、どの軸にプロットするかを決定します。
  • 処理の流れ:
    • paramsからキャリーオーバー効果の長さを取得し、最大長さを決定します。
    • 長さに基づいて、時間の経過と共にどのように効果が減衰していくかを示すキャリーオーバー効果を計算します。
    • 計算された効果を正規化して、最大値が1になるように調整します。
    • 軸にバー図を描画して、キャリーオーバー効果を可視化します。
    • プロットのタイトル、軸ラベル、ティックを設定します。
• plot_saturation_curve 関数
  • 目的:
    • 特定の特徴量に対する飽和曲線を可視化する。
  • パラメータ:
    • plot_carryover_effect関数と同様に、パラメータ、特徴量名、Figureオブジェクト、axesオブジェクトの配列、インデックスが渡されます。
  • 処理の流れ:
    • xの範囲を定義して、その範囲内で飽和曲線を計算します。
    • 計算された飽和曲線をプロットします。
    • プロットのタイトル、軸ラベル、y軸とx軸の範囲を設定します。
• plot_effects 関数
  • 目的:
    • 複数の特徴量に対するキャリーオーバー効果と飽和曲線をまとめて可視化する。
  • 処理の流れ:
    • matplotlibのsubplotを使用して、各特徴量に対するプロット領域を設定します。
    • 各特徴量に対して、キャリーオーバー効果と飽和曲線のプロット関数を順番に呼び出します。
    • plt.tight_layout()を使用して、プロット間のスペースを適切に調整します。
    • 最後に、plt.show()を使用して、すべてのプロットを表示します。

このコードでは、np.concatenate、np.linspaceなどのNumPy関数や、plt.subplots、axes.plotなどのMatplotlibの機能を活用して、データの視覚化を行っています。

 

この関数を利用してみます。

先ず、キャリーオーバー効果関数と飽和関数のハイパーパラメータを設定します。このハイパーパラメータを、今作った視覚化（グラフ）する関数にインプットしグラフを描きます。特徴量名は何を設定しても問題ございません。

以下、コードです。

# キャリーオーバー効果関数のハイパーパラメータの設定
carryover_params = [{'length': 10, 'rate': 0.5,}]

# 飽和関数のハイパーパラメータの設定
curve_params = [{'d': 5}]

# 特徴量名の設定
feature_names = ['x']

# グラフで確認
plot_effects(carryover_params, curve_params, feature_names)

 

以下、実行結果です。

 

MMMパイプライン構築（ハイパーパラメータ手動調整）

では、MMMパイプラインを構築していきます。

 

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

アドストック（キャリーオーバー効果関数と飽和関数）とRidge回帰にハイパーパラメータがあります。あらかじめ設定します。

以下、コードです。

# 各特徴量に対するキャリーオーバー効果のパラメータ設定
carryover_params = [
    {'length': 5, 'rate': 0.5},  # TVCM
    {'length': 3, 'rate': 0.5},  # Newspaper
    {'length': 1, 'rate': 0.5},  # Web
]

# 各特徴量に対する飽和関数のパラメータ設定
curve_params = [
        {'d': 5},  # TVCM
        {'d': 3},  # Newspaper
        {'d': 1},  # Web
]

# 特徴量名の設定
feature_names = X.columns

# グラフで確認
plot_effects(carryover_params, curve_params, feature_names)

 

このコードは、特定のマーケティングチャネル（TVCM、新聞、ウェブ）に対するキャリオーバー効果と飽和効果のパラメータを設定し、それらの効果をグラフで可視化するためのプロセスを示しています。plot_effects関数を使用して、これらの効果がどのように見えるかを視覚化します。

• キャリオーバー効果のパラメータ設定
  • carryover_paramsは、各マーケティングチャネルに対するキャリオーバー効果のパラメータを含むリストです。
  • 例えば、{'length': 5, 'rate': 0.5}は、TVCMのキャリオーバー効果が5期間持続し、期間ごとに効果が半減することを意味します。
  • 新聞とウェブに対しても同様に、それぞれの効果の持続期間と減衰率が設定されています。
• 飽和関数のパラメータ設定
  • curve_paramsは、各マーケティングチャネルに対する飽和効果のパラメータを含むリストです。
  • ここでの{'d': 5}は、TVCMに対する反応が非常に急であり、少量の投資で早く飽和に達することを示します。
  • 新聞とウェブに対しても、それぞれ飽和に達する速度を示すパラメータdが設定されています。
• 特徴量名（説明変数名）の設定
  • feature_names = X.columnsは、データフレームXから特徴量名を取得しています。
  • Xは、読み込んだデータセットの特徴量（説明変数）です。
• グラフでの確認
  • 最後に、plot_effects関数を呼び出して、設定したキャリオーバー効果と飽和効果のパラメータに基づいて、それぞれのマーケティングチャネルの効果をグラフで可視化します。

 

以下、実行結果です。

 

　MMM構築

この設定でMMMパイプラインを構築していきます。

以下、コードです。

# パイプラインの構築
MMM_pipeline = Pipeline([
    ('scaler', MinMaxScaler()),
    ('carryover_transformer', CustomCarryOverTransformer(carryover_params=carryover_params)),
    ('saturation_transformer', CustomSaturationTransformer(curve_params=curve_params)),
    ('ridge', Ridge(alpha=1))
])

# パイプラインを表示
print(MMM_pipeline)

 

このコードは、Scikit-LearnのPipelineクラスを用い、複数の変換処理（キャリーオーバー効果関数→飽和関数）と最終的な予測モデル（Ridge回帰）を一連のステップとして組み込んだMMMパイプラインを構築するものです。

パイプラインの各ステップを簡単に説明します。

• MinMaxScaler:
  • このステップでは、MinMaxScalerを使用して特徴量のスケーリングを行います。
  • これは、すべての特徴量を0と1の間に正規化することによって、異なるスケールの特徴量間のバランスを取り、モデルの性能を向上させるのに役立ちます。
• CustomCarryOverTransformer:
  • カスタム変換器CustomCarryOverTransformerを使用して、特定の特徴量に対するキャリオーバー効果をデータに適用します。
  • このステップでは、carryover_paramsに基づいて各特徴量に対してキャリオーバー効果が計算されます。
• CustomSaturationTransformer:
  • カスタム変換器CustomSaturationTransformerを使用して、データに飽和効果を適用します。
  • curve_paramsに基づいて各特徴量に対する飽和曲線が適用されます。
• Ridge:
  • 最終的なステップとして、リッジ回帰モデル（Ridge）が使用されます。
  • リッジ回帰は線形回帰の一種で、重みの大きさにペナルティを課す（L2正則化）ことでモデルの過学習を防ぎます。
  • alpha=1は正則化の強さを制御します。

 

以下、実行結果です。

Pipeline(steps=[('scaler', MinMaxScaler()),
                ('carryover_transformer',
                 CustomCarryOverTransformer(carryover_params=[{'length': 5,
                                                               'rate': 0.5},
                                                              {'length': 3,
                                                               'rate': 0.5},
                                                              {'length': 1,
                                                               'rate': 0.5}])),
                ('saturation_transformer',
                 CustomSaturationTransformer(curve_params=[{'d': 5}, {'d': 3},
                                                           {'d': 1}])),
                ('ridge', Ridge(alpha=1))])

 

このMMMパイプラインを学習し、実測と予測の散布図を描いてみます。

以下、コードです。

# パイプラインを使って学習
MMM_pipeline.fit(X, y)

# 学習したモデルを使って予測
predictions = MMM_pipeline.predict(X)

# 予測結果の表示
plt.figure(figsize=(10, 6))
plt.scatter(y, predictions)
plt.plot([y.min(), y.max()], [y.min(), y.max()], 'k--', lw=3)
plt.xlabel('actual')
plt.ylabel('predicted')
plt.title('Actual vs Predicted Values')
plt.show()

 

このコードは、MMMパイプラインを学習（MMM_pipeline.fit(X, y)）し予測（predictions = MMM_pipeline.predict(X)）を行った後、実際の値（y）と予測値（predictions）を比較するためのグラフを描画しています。目的は、モデルの予測精度を視覚的に評価することです。

 

以下、実行結果です。

 

前回作った以下の関数群を使い、MMMモデル構築から貢献度算出、マーケティングROIの計算をしていきます。

• MMM構築
  • train_and_evaluate_model:MMMの構築
  • plot_actual_vs_predicted:構築したMMMの予測値の時系列推移
• 後処理（結果の出力）
  • calculate_and_plot_contribution:売上貢献度の算出（時系列推移）
  • summarize_and_plot_contribution:売上貢献度構成比の算出
  • calculate_marketing_roi:マーケティングROIの算出

 

繰り返しの説明になります。

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

mmm_functions.py　※zipファイルを解凍してお使いください
https://www.salesanalytics.co.jp/oi17

 

mmm_functions.pyを利用するときは、実行するPythonファイルやNotebookと同じフォルダに入れておいてください。

以下のコードで呼び出せます。

from mmm_functions import *

 

上手くいかないときは、mmm_functions.pyをメモ帳などで開き内容をコピーし、実行するPythonファイルやNotebookにコードを張り付け、Pythonで関数を作ってからMMM構築などを行ってください。

 

関数train_and_evaluate_modelを使い、MMMを構築します。

以下、コードです。

# MMMの構築
trained_model, pred = train_and_evaluate_model(MMM_pipeline, X, y)

 

以下、実行結果です。

RMSE: 231529.38400011233
MAE: 186164.1290630438
MAPE: 0.09520981758640282
R2: 0.8280341628644571

 

ちなみに、以下は前回（アドストックを考慮していないMMM）の結果です。

RMSE: 284667.75543375366
MAE: 237679.03680369648
MAPE: 0.12129533117154913
R2: 0.7400400171566186

 

アドストックを考慮するだけで、より良いMMMが構築できていることが分かります。

 

関数plot_actual_vs_predictedを使い、関数train_and_evaluate_modelで構築したMMMの予測値と実測値をプロットします。

以下、コードです。

# 実測値と予測値の時系列推移
plot_actual_vs_predicted(y.index, y.values, pred.y.values, (0, 4e6))

 

以下、実行結果です。

 

　後処理（結果の出力）

MMMを構築したら、貢献度やマーケティングROIなどを計算していきます。

先ずは、貢献度を計算します。

以下、コードです。

# 貢献度の算出
contribution = calculate_and_plot_contribution(y, X, trained_model,(0, 4e6))

 

以下、実行結果です。

                    Base          TVCM      Newspaper            Web
Week                                                                
2018-01-07  1.043943e+06  1.088057e+06       0.000000       0.000000
2018-01-14  8.949114e+05  7.168344e+05  477681.099799  506673.079075
2018-01-21  8.453460e+05  4.371990e+05  523596.193268  430058.870744
2018-01-28  8.776735e+05  2.607515e+05  542474.920489       0.000000
2018-02-04  1.033343e+06  1.649915e+05  377150.895351  579914.539206
...                  ...           ...            ...            ...
2021-11-28  6.530195e+05  0.000000e+00  180880.490131       0.000000
2021-12-05  1.012931e+06  0.000000e+00  597733.663017  454035.374317
2021-12-12  9.057491e+05  0.000000e+00  313965.594325  469785.316126
2021-12-19  9.653712e+05  0.000000e+00  634228.767563       0.000000
2021-12-26  7.618047e+05  0.000000e+00  315587.310590  383708.004515

[208 rows x 4 columns]

 

次に、貢献度の割合を計算します。

以下、コードです。

# 貢献度構成比の算出
contribution_results = summarize_and_plot_contribution(contribution)

 

以下、実行結果です。

           contribution     ratio
Base       1.901465e+08  0.433574
TVCM       1.012073e+08  0.230774
Newspaper  7.418071e+07  0.169148
Web        7.302152e+07  0.166504

 

最後に、各メディアのマーケティングROIを計算します。

以下、コードです。

# マーケティングROIの算出
ROI = calculate_marketing_roi(X, contribution)

 

以下、実行結果です。

TVCM         0.740528
Newspaper    0.613108
Web          1.330774
dtype: float64

 

MMMパイプライン構築（ハイパーパラメータ自動調整）

　ハイパーパラメータ自動調整

ハイパーパラメータを自動チューニングするため、Optunaの目的関数を定義します。

以下、コードです。

def objective(trial, X, y):

    carryover_params = []
    curve_params = []
    n_features = X.shape[1]

    for i in X.columns:
        carryover_length = trial.suggest_int(f'carryover_length_{i}', 1, 10)
        carryover_rate = trial.suggest_float(f'carryover_rate_{i}', 0, 1)
        carryover_params.append({'length': carryover_length, 'rate': carryover_rate})

        curve_param_d = trial.suggest_float(f'curve_param_d_{i}', 0, 10)
        curve_params.append({'d': curve_param_d})

    alpha = trial.suggest_float('alpha', 1e-3, 1e+3)

    pipeline = Pipeline(steps=[
        ('scaler', MinMaxScaler()),
        ('carryover', CustomCarryOverTransformer(carryover_params=carryover_params)),
        ('saturation', CustomSaturationTransformer(curve_params=curve_params)),
        ('ridge', Ridge(alpha=alpha))
    ])

    tscv = TimeSeriesSplit(n_splits=5)
    scores = cross_val_score(pipeline, X, y, cv=tscv, scoring='neg_mean_squared_error')
    rmse = np.mean([np.sqrt(-score) for score in scores])

    return rmse

 

このコードは、Optuna ライブラリを使用してMMMパイプラインのハイパーパラメータを最適化するための目的関数objectiveを定義しています。

• パラメータの提案: trial.suggest_*メソッドを使用して、キャリーオーバー効果のパラメータ（長さと減衰率）および飽和曲線のパラメータ（d）を提案します。また、リッジ回帰モデルの正則化強度alphaも提案されます。
• パイプラインの構築: MinMaxScaler、カスタムキャリーオーバー変換器、カスタム飽和変換器、およびリッジ回帰モデルを含む機械学習パイプラインを構築します。各変換器は、提案されたパラメータを使用してインスタンス化されます。
• クロスバリデーション: TimeSeriesSplitを使用して時系列データのクロスバリデーションを実行し、モデルの性能を評価します。ここでは、負の平均二乗誤差（neg_mean_squared_error）をスコアリングメトリックとして使用しています。
• RMSEの計算: スコア（負の平均二乗誤差）の平方根を取ることで、各分割に対するRMSE（平均二乗誤差の平方根）を計算し、その平均値を目的関数の戻り値として返します。

この目的関数は、Optuna の最適化プロセスにおいて、モデルの予測精度指標であるRMSEを最小化するようなパラメータセットを見つけるために使用されます。Optuna はこの関数を複数回呼び出し、異なるパラメータの組み合わせを試しながら、目的関数の戻り値（ここではRMSE）を最小化するパラメータを探索します。

このプロセスを通じて、時間に依存する効果（キャリーオーバー効果や飽和効果）を考慮した予測モデルのパフォーマンスを最適化することができます。

 

目的関数を最適化し、最適なハイパーパラメータを探索します。

以下、コードです。

# Optunaのスタディの作成と最適化の実行
study = optuna.create_study(direction='minimize')
objective_with_data = partial(objective, X=X, y=y)
study.optimize(objective_with_data, n_trials=1000)

print("Best trial:")
trial = study.best_trial

print(f"Value: {trial.value}")
print("Params: ")
for key, value in trial.params.items():
    print(f"{key}: {value}")

 

このコードは、Optunaを使用してMMMの予測精度を向上させる最適なパラメータセットを見つけるための探索を実行しています。

• study = optuna.create_study(direction='minimize')
  • Optunaによる最適化スタディを作成し、目標とする指標（この場合はRMSE: Root Mean Squared Error）を最小化するように指定します。
• objective_with_data = partial(objective, X=X, y=y)
  • ハイパーパラメータを探索するstudy.optimize()メソッドを使用する際、目的関数（この場合はobjective関数）が必要となります。しかし、objective関数は追加のデータセット引数（Xとy）を必要とするため、直接study.optimize()に渡すことはできません。
  • そこでpartialを使用して、これらのデータセット（Xとy）を固定した新しい関数objective_with_dataを生成し、この関数をstudy.optimize()に渡します。partialは関数の一部の引数を固定するために使用されます。
• study.optimize(objective_with_data, n_trials=1000)
  • Optunaはobjective_with_data関数を使って1000回の試行を行います。
  • 各試行では、objective_with_data関数内で定義されたハイパーパラメータの異なる組み合わせを試し、最小の戻り値（この場合はRMSE）を目指します。
  • Xとyは各試行において固定されており、partialによって事前にobjective関数に渡されています。
• print("Best trial:") / print(f"Value: {trial.value}")
  • 最適化プロセスの後、最良の試行（最小のRMSEを達成した試行）とその試行時のRMSEの値を出力します。
• for key, value in trial.params.items():
  • 最良の試行におけるハイパーパラメータの値を表示します。
  • これには、キャリオーバー効果と飽和効果の各パラメータ（長さ、減衰率、dの値）やリッジ回帰の正則化パラメータalphaが含まれます。

このプロセスを通じて、設定した試行回数内で、大量のパラメータ空間を効率的に探索し、最適なハイパーパラメータセットを自動で見つけ出すことが可能になります。

 

以下、実行結果です。

Best trial:
Value: 183761.0598208403
Params: 
carryover_length_TVCM: 6
carryover_rate_TVCM: 0.38785049743734706
curve_param_d_TVCM: 2.0488730767101266
carryover_length_Newspaper: 8
carryover_rate_Newspaper: 0.1805507173453931
curve_param_d_Newspaper: 1.4913826515787498
carryover_length_Web: 6
carryover_rate_Web: 0.10536168669865958
curve_param_d_Web: 1.5291557088843128
alpha: 0.039144675085206276

 

最適なハイパーパラメータを使いMMMを構築するため、ハイパーパラメータの値を抽出し、以下のリストなどを作成します。

• best_carryover_params：キャリーオーバー効果関数のハイパーパラメータのリスト
• best_curve_params：飽和関数のハイパーパラメータのリスト
• best_alpha：Ridge回帰の正則化パラメータ

以下、コードです。

best_trial = study.best_trial

best_carryover_params = []
best_curve_params = []
n_features = X.shape[1]

for i in X.columns:
    best_carryover_params.append({
        'length': best_trial.params[f'carryover_length_{i}'],
        'rate': best_trial.params[f'carryover_rate_{i}']
    })
    
    best_curve_params.append({
        'd': best_trial.params[f'curve_param_d_{i}']
    })

best_alpha = best_trial.params['alpha']

 

このコードは、Optunaを用いたハイパーパラメータ最適化のプロセスの後、最適な試行（best_trial）から最良のパラメータを取得し、それらを使用してモデルの構築に必要なパラメータのリストを作成しています。

• best_trial = study.best_trial
  • Optunaのstudyオブジェクトから最良の試行（最適化プロセスで最小の目的関数値を達成した試行）を取得します。
• best_carryover_paramsとbest_curve_paramsリストの初期化
  • 最適化されたキャリオーバー効果と飽和曲線のパラメータを格納するための空リストを作成します。
• n_features = X.shape[1]
  • 使用する特徴量の数をXデータセットの列数から取得します。
• for i in X.columns:ループ
  • 各特徴量について、最適化されたキャリオーバー効果のパラメータ（lengthとrate）と飽和曲線のパラメータ（d）を、最良の試行のパラメータから取得します。
  • これらのパラメータは、best_trial.params辞書から特定のキー名（例えば'carryover_length_{i}'）を使用して抽出されます。{i}はループ変数で、各特徴量に対応します。
  • 取得したパラメータは、対応するリスト（best_carryover_paramsまたはbest_curve_params）に辞書形式で追加されます。
• best_alpha = best_trial.params['alpha']
  • リッジ回帰モデルの最適化された正則化パラメータalphaを取得します。

これらのステップにより、最適化プロセスを通じて見つかった最良のハイパーパラメータを使用して、モデルの再構築や評価が可能になります。

 

キャリーオーバー効果関数のハイパーパラメータのリストであるbest_carryover_paramsと、飽和関数のハイパーパラメータのリストであるbest_curve_paramsから、どのようなアドストック（キャリーオーバー効果関数と飽和関数）になったのかをグラフで確認してみます。

以下、コードです。

# グラフで確認
plot_effects(best_carryover_params, best_curve_params, X.columns)

 

以下、実行結果です。

 

　MMM構築

このハイパーパラメータでMMMパイプラインを構築していきます。

以下、コードです。

# 最適なパイプラインの構築
MMM_pipeline = Pipeline(steps=[
    ('scaler', MinMaxScaler()),
    ('carryover', CustomCarryOverTransformer(carryover_params=best_carryover_params)),
    ('saturation', CustomSaturationTransformer(curve_params=best_curve_params)),
    ('ridge', Ridge(alpha=best_alpha))
])

 

このMMMパイプラインを学習し、実測と予測の散布図を描いてみます。

以下、コードです。

# パイプラインを使って学習
MMM_pipeline.fit(X, y)

# 学習したモデルを使って予測
predictions = MMM_pipeline.predict(X)

# 予測結果の表示
plt.figure(figsize=(10, 6))
plt.scatter(y, predictions)
plt.plot([y.min(), y.max()], [y.min(), y.max()], 'k--', lw=3)
plt.xlabel('actual')
plt.ylabel('predicted')
plt.title('Actual vs Predicted Values')
plt.show()

 

以下、実行結果です。

 

前回作った以下の関数群を使い、MMMモデル構築から貢献度算出、マーケティングROIの計算をしていきます。

関数train_and_evaluate_modelを使い、MMMを評価します。

以下、コードです。

# MMMの構築
trained_model, pred = train_and_evaluate_model(MMM_pipeline, X, y)

 

以下、実行結果です。

RMSE: 175244.63286004536
MAE: 141136.68030802492
MAPE: 0.07014970643405759
R2: 0.9014811355562776

 

ちなみに、以下は先ほどのチューニングしていないMMMの結果です。

RMSE: 231529.38400011233
MAE: 186164.1290630438
MAPE: 0.09520981758640282
R2: 0.8280341628644571

 

チューニングするだけで、より良いMMMが構築できていることが分かります。

 

関数plot_actual_vs_predictedを使い、関数train_and_evaluate_modelで構築したMMMの予測値と実測値をプロットします。

以下、コードです。

# 実測値と予測値の時系列推移
plot_actual_vs_predicted(y.index, y.values, pred.y.values, (0, 4e6))

 

以下、実行結果です。

 

　後処理（結果の出力）

MMMを構築したら、貢献度やマーケティングROIなどを計算していきます。

先ずは、貢献度を計算します。

以下、コードです。

# 貢献度の算出
contribution = calculate_and_plot_contribution(y, X, trained_model,(0, 4e6))

 

以下、実行結果です。

                    Base          TVCM      Newspaper            Web
Week                                                                
2018-01-07  1.021368e+06  1.110632e+06       0.000000       0.000000
2018-01-14  9.418685e+05  5.417976e+05  534606.406128  577827.555451
2018-01-21  9.185674e+05  2.337782e+05  534713.112114  549141.237223
2018-01-28  9.780812e+05  1.017795e+05  501692.773317   99346.509736
2018-02-04  1.215017e+06  5.006540e+04  146321.012372  743996.843338
...                  ...           ...            ...            ...
2021-11-28  6.969932e+05  0.000000e+00   75428.418789   61478.371957
2021-12-05  9.860874e+05  0.000000e+00  578500.519417  500112.052809
2021-12-12  9.701668e+05  0.000000e+00  140961.633285  578371.538870
2021-12-19  8.969053e+05  0.000000e+00  612052.612022   90642.095350
2021-12-26  8.379980e+05  0.000000e+00  153536.540014  469565.472844

[208 rows x 4 columns]

 

次に、貢献度の割合を計算します。

以下、コードです。

# 貢献度構成比の算出
contribution_results = summarize_and_plot_contribution(contribution)

 

以下、実行結果です。

           contribution     ratio
Base       1.919811e+08  0.437757
TVCM       9.511540e+07  0.216883
Newspaper  6.023040e+07  0.137338
Web        9.122923e+07  0.208022

 

最後に、各メディアのマーケティングROIを計算します。

以下、コードです。

# マーケティングROIの算出
ROI = calculate_marketing_roi(X, contribution)

 

以下、実行結果です。

TVCM         0.475873
Newspaper    0.333029
Web          1.924921
dtype: float64

 

まとめ

今回は、シンプルなアドストックを組み込んだMMMを取り上げました。

このようなシンプルなアドストックであっても組み込むことで、構築されるMMMの精度が高くなります。

そもそも、シンプルでないアドストックと何なのか？

例えば、広告や販促などの効果のピークが遅れて現れることがあります。その日には効果が薄くても数日後に多大なる効果が表れることがあります。

また、飽和曲線も指数型のシンプルなものではなく、S字曲線（ロジスティック曲線やゴンペルツ曲線など）を描くことままあります。S字曲線で描くとは、飽和現象だけでなく、広告などの投下量が少ないとほぼ効果が見らないものの、ある閾値を過ぎると効果が表れる、といった現象も捉えるということです。

次回は、今回のシンプルなキャリーオーバー効果関数ではなく、ピークが広告販促媒体によって変化する、やや高度なキャリーオーバー効果関数で、MMMを構築するお話しをします。

PythonによるMMM（マーケティングミックスモデリング）とビジネス活用－ 振返り分析編（その4）：ちょっと一般的なアドストック付き線形回帰系MMM －