Bayesian Time Varying Coefficient Model with Applications to Marketing Mix Modeling

Expressing sales as a function of spending variables with diminishing marginal returns (Farris et al., 2015) is one of the fundamental properties in an attribution or marketing response model. In practice, Han and Gabor (2020), and Lewis and Wong (2018) adopted a similar strategy for their budget allocation, bidding, and attribution. In view of that, our model can be expressed in a multiplicative format as below

where $x_{t,p}$ are the regressors (i.e., the ads spending variables in our case), $\hat{y_{t}}$ is the marketing response, $g$ is a time-series process, $f$ is the cost curve function, $P$ is the number of regressors, and $T$ is the number of time points. Choice of $f$ is desired to have the following properties such that

• •

  $\hat{y}_{t}$ has an explainable structure to decompose into different driving factors,

• •

  temporal effects such as trend and seasonality of $\hat{y}_{t}$ are captured,

• •

  $f_{t,p}$ is differentiable and monotonically increasing,

• •

  $\hat{y}_{t}$ has diminishing marginal returns with respect to $x_{t,p}$.

Equation 1 has an intuitive form as

where $e^{l_{t}}$ is the trend component, $e^{s_{t}}$ is the seasonality, and $\beta_{t,p}$ are channel-specific time-varying coefficients.

With a log-log transformation, equation (2) can be re-written as

A natural idea is to use state-space models such as Dynamic Linear Model (DLM) (West and Harrison, 2006) or Kalman filter (Durbin and Koopman, 2012) to solve equation (3). However, there are some caveats associated with these approaches, especially given the goal we want to achieve with MMM:

• •

  DLM with Markov Chain Monte Carlo(MCMC) sampling is not efficient and could be costly especially for high-dimensional problems, which requires sampling for a large number of regressors and time steps.

• •

  Although the Kalman filter provides analytical solutions, it has limited room for further customization such as applying restrictions on coefficient signs (e.g., positive coefficient sign for marketing spend) or a t-distributed noise process that is more outlier robust.

Meanwhile, there have been parametric and non-parametric statistical methods proposed (Fan and Zhang, 2008). Wu and Chiang (2000) considered a nonparametric varying coefficient regression model with longitudinal dependent variable and cross-sectional covariates. Two kernel estimators based on componentwise local least squares criteria were proposed to estimate the time varying coefficients. Li et al. (2002) proposed a semiparametric smooth coefficient model as a useful yet flexible specification for studying a general regression relationship with time varying coefficients. It used a local least squares method with a kernel weight function to estimate the smooth coefficient function. Nonetheless, the frequentist approaches can be expensive when doing local estimates with respect to time dimension. There are also no straightforward ways to incorporate information from experimentation results. As such, we are motivated to develop a new approach to derive time varying coefficients under a Bayesian framework for our MMM applications.

In view of the increased complexity of regression problem in practical MMM, we propose a Bayesian Time Varying Coefficient (BTVC) model, as inspired by the Generalized Additive Models (GAM) (Hastie and Tibshirani, 1990) and kernel regression smoothing. The key idea behind BTVC is to express regression coefficients as a weighted sum of local latent variables.

First, we define a latent variable $b_{j,p}$ for the $p$-th regressor at time $t_{j}$, $p=1,\cdots,P$, $j=1,\cdots,J$, $t_{j}\in\{1,\cdots,T\}$. There are $J$ latent variables in total for each regressor. From the perspective of spline regression, $b_{j,p}$ can be viewed as a knot distributed at time $t_{j}$ for a regressor. $w$ is a time-based weighting function such that

It is intuitive to use a weighting function taking into account the time distance between $t_{j}$ and $t$,

where $k(\cdot,\cdot)$ is the kernel function, and the denominator is to normalize the weights across knots. In practice, we have different choices for the kernel functions, such as Gaussian kernel, quadratic kernel or any other custom kernels. In subsection 3.3, we will discuss this in more detail.

We can also rewrite Equation 4 into a matrix form

where $B$ is the $T\times P$ coefficient matrix with entries $\beta_{t,j}$, $K$ is the $T\times J$ kernel matrix with normalized weight entries $w_{j}(t)$, and $b$ is the $J\times P$ knot matrix with entries

where $B_{t}=(\beta_{t,1},\cdots,\beta_{t,p})$ and $X_{t}$ is the $t$-th row of regressor covariate matrix.

Besides the regression component, we can also apply Equation 4 to other components such as trend and seasonality in Equation 3. Specifically, for the trend component,

The trend component can be viewed as a dynamic intercept. For the seasonality component,

$X_{t,\text{seas}}$ is the $t$-th row of seasonality covariate matrix derived from Fourier series. In subsection 3.4, we will discuss the seasonality in more detail.

Instead of estimating the local knots directly (i.e., $b$, $b_{\text{lev}}$, and $b_{\text{seas}}$ in the above equations) by optimizing an objective function, we introduce the Bayesian framework along with customizable priors to conduct the posterior sampling.

To capture the sequential dynamics and cyclical patterns, we use the Laplace prior to model adjacent knots

The initial values ($b_{0,\text{lev}}$ and $b_{0,\text{seas}}$) can be sampled from a Laplace distribution with mean 0. A similar approach can be found in models implemented in Facebook’s Prophet package (Taylor and Letham, 2018), which uses Laplace prior to model adjacent change points of the trend component.

For the regression component, we introduce a two-layer hierarchy for more robust sampling due to the sparsity in the channel spending data,

where the superscript $+$ means a folded normal distribution (positive restriction on the coefficient signs).

In the hierarchy, the latent variable $\mu_{\text{reg}}$ depicts the overall mean of a set of knots of a single marketing lever. We can treat this as the overall estimate of a channel coefficient across time. This provides two favorable behaviors for the model: during a period with absence of spending for a channel, coefficient knot estimation of such a channel exhibits a shrinkage effect towards the overall estimate; it helps fight against over-fitting due to the volatile local structure.

The two-layer hierarchy is widely adopted in hierarchical Bayesian models and the shrinkage property is sometimes called the pooling effect on regression coefficients (Gelman et al., 2013). Figure 1 depicts the model flowchart of BTVC. Stochastic Variational Inference (Hoffman et al., 2013) is used to estimate the knot coefficient posteriors from which time varying coefficient estimates can be derived using the formulas above in Equation 7, Equation 8, and Equation 9.

For the kernel function used for trend and seasonality, we propose a customized kernel, i.e.,

For the kernel function used for regression, we adopt the Gaussian kernel, i.e.,

where $\rho$ is the scale parameter. Other kernels such as Epanechnikov kernel and quadratic kernel etc., could also be leveraged for the regression component.

Seasonality is a pattern that repeats over a regular period in a time series. To estimate seasonality, a standard approach is to decompose time-series into trend, seasonality and irregular components using Fourier analysis (De Livera et al., 2011). This method represents the time series by a set of elementary functions called basis such that all functions under study can be written as linear combinations of the elementary functions in the basis. These elementary functions involve the sine and cosine functions or complex exponential. The Fourier series approach describes the fluctuation of time series in terms of sinusoidal behavior at various frequencies.

Specifically, for a given period $S$ and a given order $k$, two series $cos(2k\pi t/S)$ and $sin(2k\pi t/S)$ will be generated to capture the seasonality. For example, with daily data, $S=7$ represents the weekly seasonality, while $S=365.25$ represents the yearly seasonality.

We implemented the BTVC model as a feature branch in our open sourced package Orbit (Ng et al., 2021) by Uber. Orbit is a software package aiming to simplify time series inferences and forecasting with structural Bayesian time series models for real-world cases and research. It provides a familiar and intuitive initialize-fit-predict interface for time series tasks, while utilizing the probabilistic programming languages such as Stan (Carpenter et al., 2017) and Pyro (Bingham et al., 2019) under the hood.

The BTVC deployment system is customized by leveraging Michelangelo (Hermann and Balso, 2017), a machine learning (ML) platform developed by Uber. Michelangelo provides centralized workflow management for end-to-end modeling process.With the help of Michelangelo, BTVC deployment system is able to automate data preprocessing, model training, validations, predictions, and measurement monitoring at scale.

The deployment workflow summarized in Figure 6 consists of three main components:

• •

  Hyperparameter management layer: this is to store and manage the various supplement data needed for BTVC model training such as normalization scalar, adstock (Jin et al., 2017), lift test based priors, as well as model specific hyperparameters.

• •

  Orchestration layer: this is to upload and trigger the model training job.

• •

  Model container: a docker container including all the essential modeling code to be integrated with Michelangelo’s ecosystem.