Największym wyzwaniem w przejściu na bardziej złożone modele w MMM był brak narzędzi wyjaśniających wpływ poszczególnych kanałów medialnych. Chociaż społeczność zajmująca się uczeniem maszynowym szeroko korzysta z podejść do wyjaśniania modeli, takich jak SHAP, co sugerują setki artykułów i przemówień konferencyjnych, nadal bardzo trudno jest znaleźć przykłady użycia SHAP w kontekście MMM. Ten świetny artykuł łączy MMM z SHAP i wyjaśnia, jak możemy interpretować wyniki marketingu mix. Zainspirowany tym artykułem napisałem niemal ogólne rozwiązanie do modelowania miksu marketingowego, łącząc pomysły metodologii Robyn dotyczące dekompozycji trendów i sezonowości, używając estymatora Random Forest (który można łatwo zmienić na inne algorytmy) oraz optymalizując zasoby reklamowe i model -specyficzne parametry przy użyciu Optuna (struktura optymalizacji hiperparametrów). Rozwiązanie pozwala na przełączanie pomiędzy optymalizacją jednego celu, jak zwykle ma to miejsce w MMM, a optymalizacją wielu celów, jak robi to Robyn.

Nadal korzystam ze zbioru danych udostępnionego przez Robyn na licencji MIT, tak jak w moim pierwszym artykule, w celu zapewnienia spójności i porównań, a także wykonuję te same kroki przygotowania danych, stosując Prophet do rozkładania trendów, sezonowości i świąt.

Zbiór danych obejmuje 208 tygodni przychodów (od 2015–11–23 do 2019–11–11) obejmujących:

• 5 kanałów wydatków na media: tv_S, ooh_S, print_S, facebook_S, search_S
• 2 kanały medialne, które zawierają również informacje o ekspozycji (wyświetlenia, kliknięcia): facebook_I, search_clicks_P (nieużywane w tym artykule)
• Media organiczne bez wydatków: newsletter
• Zmienne kontrolne: wydarzenia, święta, sprzedaż konkurencji (competitor_sales_B)

Okno analityczne obejmuje 92 tygodnie od 2016–11–21 do 2018–08–20.

Niezależnie od algorytmu modelowania, materiał reklamowy odgrywa w MMM ważną rolę. Dlatego musimy zdecydować, z jakim rodzajem adstocku będziemy eksperymentować i jakie są minimalne i maksymalne wartości, jakie może on potencjalnie mieć dla każdego kanału medialnego (przejrzyj mój poprzedni artykuł, aby zapoznać się z różnymi funkcjami adstocku). Algorytm optymalizacji wypróbuje każdą wartość adstock z zakresu zdefiniowanych wartości, aby znaleźć najlepszą, która minimalizuje kryteria optymalizacji.

Używam geometrycznej funkcji adstock zaimplementowanej w scikit-learn w następujący sposób:

from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.utils import check_array
from sklearn.utils.validation import check_is_fittedclass AdstockGeometric(BaseEstimator, TransformerMixin):
def __init__(self, alpha=0.5):
self.alpha = alpha
def fit(self, X, y=None):
X = check_array(X)
self._check_n_features(X, reset=True)
return self
def transform(self, X: np.ndarray):
check_is_fitted(self)
X = check_array(X)
self._check_n_features(X, reset=False)
x_decayed = np.zeros_like(X)
x_decayed[0] = X[0]
for xi in range(1, len(x_decayed)):
x_decayed[xi] = X[xi] + self.alpha* x_decayed[xi - 1]
return x_decayed

Wspomniałem już, że modele liniowe nie są w stanie uchwycić nieliniowych zależności pomiędzy różnymi poziomami wydatków na reklamę a wynikami. Dlatego przed modelowaniem do kanałów medialnych zastosowano różne transformacje nieliniowe, takie jak moc, ujemna wykładnicza i Hill.

Algorytmy oparte na drzewach są w stanie wychwycić nieliniowości. Dlatego nie stosuję jawnie żadnej transformacji nieliniowej i pozwalam modelowi samodzielnie uczyć się nieliniowości.

Modelowanie składa się z kilku etapów:

Parametry Adstocka

To, jak długo reklama może oddziaływać, zależy od kanału medialnego. Ponieważ szukamy optymalnego tempa zaniku materiału reklamowego, musimy realistycznie ocenić możliwe zakresy tego parametru. Wiadomo na przykład, że reklama telewizyjna może działać długotrwale, a drukowana – krócej. Musimy więc mieć elastyczność w definiowaniu realistycznych hiperparametrów dla każdego kanału medialnego. W tym przykładzie używam dokładnych zakresów zaproponowanych przez Robyn w jej pliku demonstracyjnym.

adstock_features_params = {}
adstock_features_params["tv_S_adstock"] = (0.3, 0.8)
adstock_features_params["ooh_S_adstock"] = (0.1, 0.4)
adstock_features_params["print_S_adstock"] = (0.1, 0.4)
adstock_features_params["facebook_S_adstock"] = (0.0, 0.4)
adstock_features_params["search_S_adstock"] = (0.0, 0.3)
adstock_features_params["newsletter_adstock"] = (0.1, 0.4)

Walidacja krzyżowa szeregów czasowych

Chcemy znaleźć parametry, które dobrze generalizują na niewidoczne dane. Musimy podzielić nasze dane na zbiór uczący i testowy. Ponieważ nasze dane reprezentują wydatki i przychody występujące na osi czasu, musimy zastosować weryfikację krzyżową szeregów czasowych, tak aby zbiór szkoleniowy składał się wyłącznie ze zdarzeń, które wystąpiły przed zdarzeniami w zestawie testowym.

Algorytmy uczenia maszynowego działają najlepiej, gdy są szkolone na dużych ilościach danych. Algorytm Random Forest nie jest wyjątkiem i aby wychwycić nieliniowości i interakcje pomiędzy zmiennymi, należy go trenować na dużej liczbie danych. Jak wspomniałem wcześniej, w sumie mamy tylko 208 punktów danych i 92 punkty danych w oknie analizy. Potrzebujemy kompromisu między możliwością uogólniania a zdolnością modelu do uczenia się.

Po kilku eksperymentach zdecydowałem się zastosować 3 podziały CV, przydzielając dane z 20 tygodni (około 10%) jako zbiór testowy.

from sklearn.model_selection import TimeSeriesSplit
tscv = TimeSeriesSplit(n_splits=3, test_size = 20)

Każdy kolejny podział zestawu treningowego jest większy niż poprzedni.

tscv = TimeSeriesSplit(n_splits=3, test_size = 20)
for train_index, test_index in tscv.split(data):
print(f"train size: {len(train_index)}, test size: {len(test_index)}")#train size: 148, test size: 20
#train size: 168, test size: 20
#train size: 188, test size: 20

Optymalizacja hiperparametrów przy użyciu Optuna

Optymalizacja hiperparametrów składa się z szeregu eksperymentów lub prób. Każdą próbę można z grubsza podzielić na trzy etapy.

• Zastosuj transformację Adstock w kanałach medialnych, korzystając z zestawu parametrów Adstock

for feature in adstock_features:
adstock_param = f"{feature}_adstock"
min_, max_ = adstock_features_params[adstock_param]
adstock_alpha = trial.suggest_uniform(f"adstock_alpha_{feature}", min_, max_)
adstock_alphas[feature] = adstock_alpha
#adstock transformation
x_feature = data[feature].values.reshape(-1, 1)
temp_adstock = AdstockGeometric(alpha = adstock_alpha).fit_transform(x_feature)
data_temp[feature] = temp_adstock

#Random Forest parameters
n_estimators = trial.suggest_int("n_estimators", 5, 100)
min_samples_leaf = trial.suggest_int('min_samples_leaf', 1, 20)
min_samples_split = trial.suggest_int('min_samples_split', 2, 20)
max_depth = trial.suggest_int("max_depth", 4,7)
ccp_alpha = trial.suggest_uniform("ccp_alpha", 0, 0.3)
bootstrap = trial.suggest_categorical("bootstrap", [False, True])
criterion = trial.suggest_categorical("criterion",["squared_error"])

• Sprawdź krzyżowo i zmierz średni błąd we wszystkich zestawach testowych

for train_index, test_index in tscv.split(data_temp):
x_train = data_temp.iloc[train_index][features]
y_train =  data_temp[target].values[train_index]
x_test = data_temp.iloc[test_index][features]
y_test = data_temp[target].values[test_index]
#apply Random Forest
params = {"n_estimators": n_estimators, 
"min_samples_leaf":min_samples_leaf, 
"min_samples_split" : min_samples_split,
"max_depth" : max_depth, 
"ccp_alpha" : ccp_alpha, 
"bootstrap" : bootstrap, 
"criterion" : criterion
}
#train a model      
rf = RandomForestRegressor(random_state=0, **params)
rf.fit(x_train, y_train)
#predict test set
prediction = rf.predict(x_test)
#RMSE error metric       
rmse = mean_squared_error(y_true = y_test, y_pred = prediction, squared = False)
#collect errors for each fold
scores.append(rmse)#finally return the average of the cv error
return np.mean(scores)

Każda próba zwraca zasoby reklamowe, parametry modelu i metryki błędów jako atrybut użytkownika. Umożliwia to łatwe wyszukanie parametrów w najlepszej próbie.

trial.set_user_attr("scores", scores)
trial.set_user_attr("params", params)
trial.set_user_attr("adstock_alphas", adstock_alphas)

Główną funkcją rozpoczynającą optymalizację jest optuna_optymalizacja. Zwraca obiekt Optuna Study ze wszystkimi próbami, w tym najlepszą próbą (posiadającą minimalny średni błąd RMSE)

tscv = TimeSeriesSplit(n_splits=3, test_size = 20)adstock_features_params = {}
adstock_features_params["tv_S_adstock"] = (0.3, 0.8)
adstock_features_params["ooh_S_adstock"] = (0.1, 0.4)
adstock_features_params["print_S_adstock"] = (0.1, 0.4)
adstock_features_params["facebook_S_adstock"] = (0.0, 0.4)
adstock_features_params["search_S_adstock"] = (0.0, 0.3)
adstock_features_params["newsletter_adstock"] = (0.1, 0.4)OPTUNA_TRIALS = 2000#experiment is an optuna study object
experiment = optuna_optimize(trials = OPTUNA_TRIALS, 
data = data, 
target = target, 
features = features, 
adstock_features = media_channels + organic_channels, 
adstock_features_params = adstock_features_params, 
media_features=media_channels, 
tscv = tscv, 
is_multiobjective=False)

Wynik RMSE dla każdej części najlepszej próby:

experiment.best_trial.user_attrs["scores"]#[162390.01010327024, 114089.35799374945, 79415.8649240292]

Parametry Adstock odpowiadające najlepszej wersji próbnej:

experiment.best_trial.user_attrs["adstock_alphas"]#{'tv_S': 0.5343389820427953,
# 'ooh_S': 0.21179063584028718,
# 'print_S': 0.27877433150946473,
# 'facebook_S': 0.3447366707231967,
# 'search_S': 0.11609804659096469,
# 'newsletter': 0.2559060243894163}

Parametry modelu odpowiadające najlepszej próbie:

experiment.best_trial.user_attrs["params"]#{'n_estimators': 17,
# 'min_samples_leaf': 2,
# 'min_samples_split': 4,
# 'max_depth': 7,
# 'ccp_alpha': 0.19951653203058856,
# 'bootstrap': True,
# 'criterion': 'squared_error'}

Ostateczny model

Buduję finalny model wykorzystując zoptymalizowane parametry, podając do analizy okres początkowy i końcowy. Model budowany jest najpierw na wszystkich danych aż do końca okresu analizy. Prognozy i wartości SHAP są pobierane tylko dla okresu analizy.

best_params = experiment.best_trial.user_attrs["params"]
adstock_params = experiment.best_trial.user_attrs["adstock_alphas"]
result = model_refit(data = data, 
target = target,
features = features, 
media_channels = media_channels, 
organic_channels = organic_channels, 
model_params = best_params, 
adstock_params = adstock_params, 
start_index = START_ANALYSIS_INDEX, 
end_index = END_ANALYSIS_INDEX)

Pierwszym wykresem, który należy sprawdzić, jest dopasowanie modelu do danych z okresu analizy wynoszącego 92 tygodnie:

MAPE poprawiło się o 40%, a NRMSE o 17% w porównaniu z podejściem bayesowskim.

Następnie narysujmy stosunek wydatków do udziału w efekcie:

Udział efektu oblicza się na podstawie bezwzględnej sumy wartości SHAP dla każdego kanału medialnego w przedziale analizy i normalizuje się przez całkowitą sumę wartości SHAP wszystkich kanałów medialnych.

Udział efektu jest prawie zgodny z udziałem efektu z poprzedniego artykułu. Zauważam tylko jedną niespójność pomiędzy udziałem efektu szukaj kanał.

Zmniejszający się efekt zwrotu/nasycenia

Nie zastosowałem żadnych transformacji nieliniowych, aby jawnie modelować malejące zyski. Sprawdźmy więc, czy Random Forest może uchwycić jakiekolwiek nieliniowości.

Osiąga się to za pomocą wykresu punktowego, który pokazuje wpływ pojedynczego kanału medialnego na przewidywania dokonane przez model, gdzie oś x to wydatki na media, oś y to wartość SHAP dla tego kanału medialnego, która reprezentuje, ile znajomość konkretnego wydatku zmienia wynik modelu dla tej konkretnej prognozy. Linia pozioma odpowiada wartości SHAP wynoszącej 0. Linia pionowa odpowiada średnim wydatkom w kanale. Zielona linia to NAJNIŻSZA krzywa wygładzania.

Patrząc na wszystkie kanały medialne widzimy, że wyższe wydatki wiążą się ze wzrostem przychodów. Ale zależności nie zawsze są liniowe.

Nabierający drukuj_S obserwujemy niewielki spadek przychodów przy wydatkach do 25 tys. Następnie zaczyna rosnąć aż do ok. 90 tys., gdzie następuje wyhamowanie wzrostu przychodów.

Nabierający facebook_S praktycznie nie obserwujemy zmian w przychodach przy wydatkach do 90 tys. i po 250 tys. Wydatki od 90 tys. do 250 tys. to prawdopodobnie najbardziej optymalne wydatki.

Niektóre kanały medialne, np facebook_S, print_S i search_S charakteryzują się dużą rozbieżnością między wartościami SHAP przy tych samych wydatkach. Można to wyjaśnić interakcją z innymi kanałami medialnymi i należy to dokładniej zbadać.

To rozwiązanie może zarządzać optymalizacją wielocelową. Pomysł wyszedł od Robyn, która wprowadziła drugą metrykę optymalizacji RSSD (odległość pierwiastkowa rozkładu)

Odległość uwzględnia relację pomiędzy udziałem w wydatkach a udziałem w dekompozycji współczynnika kanału. Jeśli odległość jest zbyt duża, jej wynik może być zbyt nierealistyczny – np. największy efekt daje działalność medialna przy najmniejszych wydatkach

W przypadku optymalizacji wieloobiektowej Optuna wyznaczy tzw. Front Pareto, czyli zbiór wszystkich optymalnych rozwiązań. Procedura będzie taka sama, jak w przypadku pojedynczego przypadku optymalizacji: dla każdego modelu należącego do Frontu Pareto pobieramy jego parametry, budujemy finalny model i wizualizujemy wyniki.
Na poniższym wykresie wszystkie punkty o kolorze czerwonawym należą do zbioru optymalnego.

W tym artykule kontynuowałem badanie sposobów ulepszenia modeli Marketing Mix przy użyciu bardziej złożonych algorytmów, które są w stanie wychwycić nieliniowość i zmienne interakcje. W rezultacie cały potok zostaje uproszczony poprzez pominięcie kroku transformacji nieliniowej, który jest zawsze stosowany przy stosowaniu regresji liniowej. Zastosowanie wartości SHAP umożliwiło dalszą analizę krzywych udziału efektu i odpowiedzi. Moim drugim celem było osiągnięcie spójnych wyników w przypadku różnych podejść. Porównanie wyników z poprzedniego artykułu, w którym zastosowałem modelowanie bayesowskie, z wynikami tego artykułu wykazało wysoki stopień spójności w rozłożonych efektach na kanał medialny.

Cały kod można pobrać z mojego repozytorium Github

Dziękuje za przeczytanie!