# Poisson regression for market media modeling?

I put together a related question linked for added context; it generally overviews a model proposed by Google in 2017 for (Marketing) Media Mix Modelling (MMM.) The model updates the pretense of MMM in that it accounts for delay (advertising spend today might have peak influential effect on customers 2-3 days later) and saturation (diminishing returns in target variable after a certain amount of spend.)

As detailed in the question, the model is extremely sensitive to variable scale; finding appropriate priors is a non-trivial task. In best case scenario, the effect of inputs is minimized and the intercept is effectively predicted for all observations. At worst, the sampler fails to converge or fails to compute gradients and quits.

In response, I've created an architecture of my own that has promising results: Namely, I still use the decay function Google proposed; however, I've multiplied this output per channel by 0.5^{Power} where power is inferred per channel parameter. This mitigates the effect of scale disparity. Using pyMC3 code:

import arviz as az
import pymc3 as pm

with pm.Model() as train_model:
#var,      dist, pm.name,          params,  shape
p     = pm.Gamma('power',          2 , 2,   shape=X.shape[1]) # raises shrinkage to power

alpha = pm.Beta('alpha',           3 , 3,   shape=X.shape[1]) # retain rate in adstock
theta = pm.Uniform('theta',        0 , 12,  shape=X.shape[1]) # delay in adstock

tau    = pm.Normal('intercept',    0,  5                    ) # model intercept
noise = pm.Gamma('noise',          3,  1                    ) # variance about y

computations = []
for idx,col in enumerate(train.columns):
delay = geometric_adstock(x=train[col].values, alpha=alpha[idx], theta=theta[idx], L=12)
comp = 0.5**p[idx] * delay
computations.append(comp)

y_hat = pm.Normal('y_hat', mu= tau + sum(computations),
sigma=noise,
observed=train_y)

trace_train = pm.sample(chains=4)


And the trained model's performance on unseen test data:

So far, I've just built up background context for my actual question.

This approach works well in this situation; however, I would like to find a way to capture saturation in this model by other means. Of note, the scale of my data is such that the inputs (marketing spend) are much higher than the output (new customers gained that interval.) Because of the events per interval nature of my output, would Poisson regression be appropriate?

It makes little sense for my output to take on a negative number (very unlikely that a marketing campaign would cause customers to leave, though I suppose possible). Likewise, I think that Poisson sort of has a saturation effect built into it, whereas linear regression inherently doesn't taper off expectations at higher level inputs. (Could someone verify my thinking?)

Lastly, due to the 0.5^{power} transformation, all inputs are drawn towards the same scale as the target variable. I might not be able to infer a saturation rate per channel (as Google's approach tried to do) but once all channels are similarly scaled, a single saturation effect, captured by the Poisson likelihood should function well, regardless. (Is this line of thought reasonable?)

lastly, I have not implemented this yet, as I wanted to explore its theoretical justifications/shortcomings first. However, I'm wondering if the exponential nature of the Poisson PDF and the power parameters "fight" each other; it might produce bizarre curvature and tricky geometry for a gradient based sampler. Any thoughts?

Edit:

See the below distribution of sales data. Note there are only 200 examples in my dataset. So my data might be normally distributed with skewness by chance or Poisson distributed and properly represented.

• It would help if you provided information about the target variable, e.g., how are the counts distributed, i.e., what's their pdf? And, more importantly, are they aggregated from multiple sources, locations, channels, etc...? If the counts per interval are large and come from more than one source then the lognormal transformation would be a more appropriate assumption than the poisson. Next, there are many ways to capture saturation. One of the best books about MMMs is Lee Cooper's Market Share Analysis. In it, he describes several different functional forms for saturation. Dec 14, 2020 at 20:34
• @MikeHunter, thanks for the question, comment, and book recommendation! I've added a graphic for target metric. (I wrote sales data but meant customer acquisitions- absent minded mistake.) Yes, the target variable is an aggregation of latent acquisitions by channel. Dec 14, 2020 at 21:31
• Thanks, it's much clearer now. Poisson or negative binomial regression modeling would be justified. One more question: why not break out the counts by channel? In my previous work with data similar to this, channel proved to be a significant predictor of, in one case, the likelihood of attrition. Dec 15, 2020 at 17:06
• @MikeHunter, if you follow the link to my related question, there's more info on the kaggle data. I have 3 marketing channel spend data and only one target variable. My model does infer latent per channel parameters, but I think that's the best I could do, unfortunately. Dec 15, 2020 at 17:45
• I am not sure the 0.5**p makes sense, do you just want a positive number? You can use the HalfNormal distribution. I generally scale all my inputs by dividing by the mean. Also check github.com/google/lightweight_mmm May 27 at 10:57

1. Is Poisson regression appropriate given data's shape and information supplied?

The long answer is maybe but not without some transformations of the data. I took a peak at the data set you referenced and your sales variable is not an integer. The reason the longer answer is a "maybe" is that it is at least conceivable that sales data would count something like the number of sales, which would be an integer and always $$\ge 0$$.

1. Would a continuous PDF be more appropriate and/or sample better due to gradient computation reasons (ex: Gamma, exponential, lognormal)?

1. Any other recommendations?

I am a little curious as to the need to re-scale. I am not super familiar with PyMC's backend. You mentioned Stan. I have a simple lagged time-series regression to demonstrate that the lognormal using a default set of mildly informative priors behaves reasonably well.

Apologies for taking you the R direction here, but we work with the tools we know :)

library(brms)
library(glue)
library(tidybayes)
library(tidyverse)

ROOT_DIR <- '~/github/kaggle_fun'
DATA_DIR <- '{ROOT_DIR}/data' %>% glue()

dat %>%
select(-X) %>%
psych::pairs.panels()

# Create lagged variables
dat <- dat %>%
mutate(
TV_lag1 = lag(TV),
newspaper_lag1 = lag(newspaper),
sales_lag1 = lag(sales),
TV_lag2 = lag(TV, 2),
newspaper_lag2 = lag(newspaper, 2),
sales_lag2 = lag(sales, 2),
TV_lag3 = lag(TV, 3),
newspaper_lag3 = lag(newspaper, 3),
sales_lag3 = lag(sales, 3)
)

yvar <- 'sales'
n_lags <- 3

model_formula <- '
sales ~
{paste(xvars, collapse = " + ")} +
{paste(
paste0(c(yvar, xvars), "_lag", rep(1:n_lags, each = length(c(yvar, xvars)))),
collapse = " + "
)
}
' %>% glue() %>%
brmsformula()

# Going to just simplify by dropping cases missing lagged values (at the start of the collection window)
dat <- na.omit(dat)
set.seed(seed = 16100107)

# since this is a time series we'll split the data but maintain continuity

# Let's start with a Poisson distribution here using default priors (we can see them later)
fit_lgnrml <- brm(
model_formula + lognormal(),
data = dat_train,
chains = 3,
iter = 2000,
warmup = 1000,
future = TRUE
)

summary(fit_lgnrml)
pp_check(fit_lgnrml, nsamples = 100)

dat_test %>%
mutate(
row_idx = X,
y_hat = .value
) %>%
ungroup() %>%
select(row_idx, y_hat) %>%
ggplot() +
stat_lineribbon(aes(x = row_idx, y = y_hat), color = "#08519C",
fill = '#426EBD', lwd = .5, .width = .95, alpha = .7) +
geom_point(data = dat_test, aes(y = sales, x = X), alpha = .5) +
geom_line(data = dat_test, aes(y = sales, x = X), lty = 'dashed', alpha = .5) +
theme_bw() +
labs(x = "", y = "Sales", caption = "Dashed line and connected points represent observed timeseries")


which will get you this plot:

And because I added a little abstraction to the creation of the model code to make it easier for me to play around with the model a bit, I have included a copy of the regression model below (with some manual editing for clarity).

sales ~ TV + radio + newspaper +
sales_lag1 + TV_lag1 + radio_lag1 + newspaper_lag1 +
sales_lag2 + TV_lag2 + radio_lag2 + newspaper_lag2 +
sales_lag3 + TV_lag3 + radio_lag3 + newspaper_lag3


UPDATE

And here is the Stan code that the model and training data generate:

// generated with brms 2.14.4
functions {
}
data {
int<lower=1> N;  // total number of observations
vector[N] Y;  // response variable
int<lower=1> K;  // number of population-level effects
matrix[N, K] X;  // population-level design matrix
int prior_only;  // should the likelihood be ignored?
}
transformed data {
int Kc = K - 1;
matrix[N, Kc] Xc;  // centered version of X without an intercept
vector[Kc] means_X;  // column means of X before centering
for (i in 2:K) {
means_X[i - 1] = mean(X[, i]);
Xc[, i - 1] = X[, i] - means_X[i - 1];
}
}
parameters {
vector[Kc] b;  // population-level effects
real Intercept;  // temporary intercept for centered predictors
real<lower=0> sigma;  // residual SD
}
transformed parameters {
}
model {
// likelihood including all constants
if (!prior_only) {
// initialize linear predictor term
vector[N] mu = Intercept + Xc * b;
target += lognormal_lpdf(Y | mu, sigma);
}
// priors including all constants
target += student_t_lpdf(Intercept | 3, 2.6, 2.5);
target += student_t_lpdf(sigma | 3, 0, 2.5)
- 1 * student_t_lccdf(0 | 3, 0, 2.5);
}
generated quantities {
// actual population-level intercept
real b_Intercept = Intercept - dot_product(means_X, b);
}


We're basically saying that $$y$$ is the result of a process that is normally distributed in log-space. The linear model predicts in the log-space and then you just need to exponentiate back across the link function to get to the original scale. The add_fitted_draws function I used to create the second plot above does just that.

The crux of the model is here in a matrix-notated form:

vector[N] mu = Intercept + Xc * b;
target += lognormal_lpdf(Y | mu, sigma);


And there are some transformations (mainly centering of predictors) that make the estimation steps a bit easier/more efficient (this happens by default when using brms). The intercept is then transformed back at the end.

generated quantities {
// actual population-level intercept
real b_Intercept = Intercept - dot_product(means_X, b);
}

• Thanks for the detailed response, I'm not very fluent in R so I have some follow-ups: (1) You went with the lognormal likelihood; did you need to transform your variables to log space for this work? (I couldn't tell either way.) Dec 23, 2020 at 16:23
• And (2) Yes, to make use of Google's model, I had to create a power parameter that scales inputs (but discarded their saturation function.) One thing I liked about their model was the geometric function, it infers a decay and delay parameters on advertising spend. This doesn't improve model accuracy (as yours is just as good) but it does offer some real-world insights, which assist to better understanding optimal spend strategy. Dec 23, 2020 at 16:30
• No the variables all enter the the brms call (my preferred Stan wrapper in R) as is and they are passed along to Stan without transformations. I am going to add the Stan code above and some additional commentary so that this could be ported over time python. Dec 23, 2020 at 16:41
• That is a fair point and it is always nice when you have a model parameter that has a clear mapping to relevant phenomena. Dec 23, 2020 at 17:00