Poisson regression for market media modeling?

Question

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.

Answer 1

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$.

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

I'd probably start with a lognormal prior for the likelihood.

3. 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 <- read.csv('{DATA_DIR}/Advertising.csv' %>% glue())

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

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

yvar <- 'sales'
xvars <- c('TV', 'radio', 'newspaper')
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 
training_mask <- 1:floor(nrow(dat) * .8)
dat_train <- dat[training_mask,]
dat_test <- dat[-training_mask,]

# 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, 
    control = list(adapt_delta = .95), 
    future = TRUE
)

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

dat_test %>% 
    add_fitted_draws(fit_lgnrml) %>% 
    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);
}