Байесовское моделирование маркетингового комплекса на Python с помощью PyMC

Автор: Дмитрий Иванов [Команда P9X]

~8 минут чтения

Фото Грега Ракози на Unsplash

В этой статье я хочу объединить две концепции, обсуждавшиеся в предыдущих постах: байесовское моделирование и моделирование маркетингового микса. Поскольку велика вероятность, что вам незнакомы обе темы, позвольте мне кратко познакомить вас с ними и предложить дополнительную литературу. Я расскажу:

  1. что такое моделирование маркетингового микса;
  2. что такое байесовское моделирование;
  3. почему имеет смысл объединить оба подхода.

Затем я покажу вам, как это сделать на практике, используя PyMC.

Предварительные сведения

Если вы жадно читаете мои статьи (спасибо!), вы можете пропустить несколько разделов, чтобы сразу перейти к коду. В противном случае продолжайте чтение.

Моделирование маркетингового микса

Фундаментальная проблема для каждого бизнеса — решить, на какие каналы тратить маркетинговый бюджет. Вы можете потратить 1 000 € на ТВ-рекламу, 2 000 € на радиорекламу и 3 000 € на веб-баннеры каждый день, полагаясь только на интуицию. Но так ли это хорошо?

Возможно, канал веб-баннеров уже насыщен, и потратить 1 500 € там так же хорошо, как 3 000 €. Тогда вы можете сэкономить 1 500 € или вложить их в другие, более эффективные каналы для увеличения продаж.

Или, возможно, какой-то канал имеет отрицательную рентабельность инвестиций — за каждый потраченный на рекламу Евро вы получаете обратно меньше одного Евро. Мы определённо не должны тратить слишком много денег на такой канал, по крайней мере, если он не является стратегически важным с точки зрения бизнеса.

Чтобы ответить на подобные вопросы, вы должны понимать, как различные расходы на медиа (ТВ, радио и т. д.) влияют на ваши продажи или другие ключевые показатели эффективности (KPI).

В моделировании маркетингового микса вы начинаете с набора данных о расходах на медиа. Обычно его дополняют некоторыми контрольными переменными, то есть дополнительной информацией обо всём, что может повлиять на целевой KPI, например, о праздниках, погоде, футбольных чемпионатах, локдаунах, цене продукта и многом другом. Чтобы сделать текст более кратким, мы опустим контрольные переменные здесь. Затем, конечно, вам нужен KPI, который вы хотите предсказать. Часто это продажи, количество новых клиентов и т. д. Таким образом, типичный набор данных может выглядеть следующим образом:

Байесовское моделирование

Многие оценки и модели основаны на подходе максимального правдоподобия. Например, представьте, что вы хотите оценить вероятность p выпадения орла при подбрасывании монеты. Вы подбрасываете её 10 раз и видите 8 орлов; что вы можете заключить? Естественной оценкой для вероятности будет p = 8 / 10 = 80%, что также является оценкой максимального правдоподобия. Вы также можете рассчитать доверительный интервал, чтобы увидеть, насколько эта оценка надёжна, но мы хотим пойти другим путём.

Представьте, что мы хотим учесть некоторые предварительные знания о вероятности p. Если вы, например, вытащили монету случайным образом из своего кошелька, нет причин думать, что монета смещена, поэтому p не должна быть слишком далеко от 50%, если вы, конечно, не фокусник.

С помощью байесовского моделирования вы можете учесть эти предварительные знания, чтобы получить оценку плотности p, то есть не одно значение, а целое распределение. Это распределение, скорее всего, будет иметь пик между оценкой максимального правдоподобия и предварительной оценкой, например, 65%.

Почему моделирование маркетингового микса с помощью Байеса?

Вы можете определить модели маркетингового микса с большим количеством гиперпараметров:

  • насыщенность;
  • сила переноса;
  • длина переноса и т. д.

Затем вы можете использовать метод оптимизации гиперпараметров, чтобы найти оптимальную комбинацию. Я сделал это в своей другой статье о моделировании маркетингового микса, Усовершенствованное моделирование маркетингового микса в Python.

Этот подход работает нормально, но есть одна вещь, которая мне не нравится:

Оценки гиперпараметров часто нестабильны.

Это означает, что совершенно разные наборы гиперпараметров могут давать одинаково хорошие модели.

Однако экстраполяция — это то, что нужно делать при оптимизации маркетингового бюджета. Если вы тратили от 0 € до 1 000 € в день на ТВ-рекламу в прошлом, для оптимизации вам нужно знать, что произойдёт, когда вы потратите 5 000 € или даже 10 000 €.

Вам нужна модель с отличными возможностями экстраполяции.

И обычно вам приходится выбирать между более чем двумя моделями. В этом случае вы можете действовать как минимум двумя способами:

  • Вы можете выбрать первую модель, которую создадите, потому что вы даже не осознаёте проблему. Этот подход прост, но опасен.
  • Вы можете выбрать модель, которая кажется вам правильной, некоторым экспертам в предметной области или заинтересованным сторонам. Для некоторых людей это нормально, но я предпочитаю не закладывать в модель ожидания по выходу, потому что возникает следующий вопрос:

Если кто-то уже знает ответ, зачем мне вообще строить модель, которая воспроизводит этот ответ?

Начало моделирования

Сначала давайте возьмём наш набор данных.

import pandas as pd

data = pd.read_csv(
  'https://raw.githubusercontent.com/Garve/datasets/4576d323bf2b66c906d5130d686245ad205505cf/mmm.csv',
  parse_dates=['Date'],
  index_col='Date'
)

X = data.drop(columns=['Sales'])
y = data['Sales']

Затем мы должны определить функции насыщения и переноса, аналогично предыдущей статье. На языке PyMC это может выглядеть так:

import pytensor.tensor as pt

def saturate(x, a):
    return 1 - pt.exp(-a*x)

def carryover(x, strength, length=21):
    w = pt.as_tensor_variable(
        [pt.power(strength, i) for i in range(length)]
    )

    x_lags = pt.stack(
        [pt.concatenate([
            pt.zeros(i),
            x[:x.shape[0]-i]
        ]) for i in range(length)]
    )

    return pt.dot(w, x_lags)

Функция насыщения должна быть простой для понимания. Перенос, однако, требует некоторой работы. В основном вы можете выразить преобразование переноса как умножение матрицы на вектор. Вам просто нужно сначала собрать матрицу x_lags и вектор w.

Анализ выходных данных модели

После этого мы можем посмотреть на обычные графики. Давайте начнём с posterior распределений. Выполнив

import arviz as az

az.plot_posterior(
    trace,
    var_names=['~contribution'],
    filter_vars='like'
)

вы получите

Здесь вы можете увидеть posterior для всех параметров. У всех них есть унимодальная (с одним пиком) форма. Вы также можете изучить, как пары переменных ведут себя вместе, с помощью

az.plot_pair(
    trace,
    var_names=['coef_TV', 'sat_TV'],
    marginals=True
)

Здесь вы можете увидеть, что параметры насыщенности и коэффициент регрессии не являются независимыми, а отрицательно коррелированы: чем выше коэффициент, тем ниже параметр насыщенности имеет тенденцию быть. Это имеет смысл, потому что более высокий коэффициент может компенсировать более медленное увеличение кривой насыщения (= более низкий sat_TV) и наоборот.

Давайте посмотрим на другой график:

az.plot_pair(
    trace,
    var_names=['car_TV', 'sat_TV'],
    marginals=True
)

Здесь мы видим, почему у оптимизации гиперпараметров могут быть проблемы. Каждая точка на этом графике — это потенциальная модель, которую вы могли бы найти с помощью оптимизации гиперпараметров.

Для действительно уникальной лучшей модели мы бы скорее увидели облако точек, плотно сконцентрированное вокруг одной точки (car_TV_true, sat_TV_true). Однако здесь мы видим, что сила переноса ТВ может иметь разумные значения от 0,4 до 0,5 в зависимости от параметра насыщенности.

Мы также должны проверить, хороша ли модель, прежде чем делать выводы. Выполнив

import matplotlib.pyplot as plt

with mmm:
    posterior = pm.sample_posterior_predictive(trace)

means = posterior.posterior_predictive['sales'].mean(dim=['chain', 'draw'])
stds = posterior.posterior_predictive['sales'].std(dim=['chain', 'draw'])

plt.figure(figsize=(20, 8))
plt.plot(y.values, linewidth=2, c='r', label='Observations')
plt.plot(means, linewidth=1, c='b', label='Mean prediction')
plt.fill_between(np.arange(len(y)), means - 2*stds, means + 2*stds, alpha=0.33)
plt.legend()

вы получите

так что похоже, что модель уловила что-то полезное.

Вклад каналов

Мы занимались распределениями, но для нашей любимой картинки с вкладом каналов давайте снова возьмём средние значения, чтобы получить одно значение. Поскольку мы ввели некоторые переменные вклада канала в код PyMC, мы можем легко извлечь их сейчас, используя короткую функцию compute_mean.

channels = ['Banners', 'Radio', 'TV']
unadj_contributions = pd.DataFrame(
    {'Base': trace.posterior['base'].mean().values},
    index=X.index
)

for channel in channels:
    unadj_contributions[channel] = trace.posterior[f'contribution_{channel}'].mean(dim=['chain', 'draw']).values

adj_contributions = (unadj_contributions
                     .div(unadj_contributions.sum(axis=1), axis=0)
                     .mul(y, axis=0)
                    )

ax = (adj_contributions
      .plot.area(
          figsize=(16, 10),
          linewidth=1,
          title='Predicted Sales and Breakdown',
          ylabel='Sales',
          xlabel='Date'
      )
     )

handles, labels = ax.get_legend_handles_labels()
ax.legend(
    handles[::-1], labels[::-1],
    title='Channels', loc="center left",
    bbox_to_anchor=(1.01, 0.5)
)

Резюме и перспективы

В этой статье мы обсудили, что моделирование маркетингового микса с использованием подхода максимального правдоподобия может быть затруднено из-за оценки гиперпараметров. Может быть много моделей, которые хорошо работают на тестовом наборе, но имеют совершенно разное поведение при экстраполяции — параметры довольно нестабильны.

Однако правильная экстраполяция — это ключ к оптимизации. Поэтому мы разработали базовую модель маркетингового микса с эффектами насыщения и переноса в байесовской постановке. Это оказалось полезным, потому что позволяет одновременно оценить все параметры, давая нам более стабильные оценки.

Мы реализовали это с помощью PyMC и создали красивую диаграмму вкладов, чтобы увидеть, какой вклад в продажи внесли расходы на каналы.

Расширения модели

Мы могли бы теперь взять байесовскую модель, которую мы создали, и расширить её ещё больше. Например, мы можем учесть изменение параметров во времени. Это может быть полезно, если сила переноса ТВ составляла 0,8 два года назад, но постепенно снизилась до 0,5 с течением времени, что является примером дрейфа концепции. Мы можем легко смоделировать это с помощью гауссовых случайных блужданий, как показано в моей другой статье о скользящей регрессии:

Скользящая регрессия в Python с помощью PyMC3

Таким же образом мы можем смоделировать изменение базовых показателей. До сих пор мы рассматривали базовый показатель как одно число, фиксированное на протяжении всего периода обучения. Но он мог также увеличиться со временем из-за наших стратегических усилий в области рекламы. Это трудно выяснить, используя наши старые модели максимального правдоподобия с оптимизацией гиперпараметров.