REBA Model Technical Documentation

The REBA model is a Bayesian synthetic control model.

Given a target time series and a set of donor time series, REBA constructs a twin (a “synthetic control”) of the target time series from the donors using historical data, then measures the difference between the target and the synthetic control in another time period (the “test period”) during which there is some intervention in the target. Since the donors did not receive the intervention, the difference between the synthetic control and the target during the test period quantifies the effect of the intervention.

In this document we’ll discuss exactly how REBA constructs the synthetic control and how it calculates the intervention effect estimate, but first we need to give an overview of the data we use.

The data

In a geographic lift test, some regions receive an intervention and others do not. Examples of common interventions would be things like doubling spend in a marketing channel or turning off marketing across a swath of channels. Some KPI of interest, like revenue, is tracked in all regions.

Suppose there are three test regions that received the intervention (it could be any number of regions). If we call their revenue time series R1(t), R2(t), and R3(t), then we can define R(t) to be the average of their time series:

R(t) = ( R1(t) + R2(t) + R3(t) ) / 3

We use this R(t) as the target for which we construct the synthetic control.

We will also have a number of donor regions that did not receive the intervention. Call their revenue time series X1(t), X2(t), all the way up to XN(t). Sometimes we will refer to these as Xi(t), where i could be any number from 1 to N.

Finally, there is some time step t = T at which the intervention starts. The time after that point, t >= T, contains the intervention and any cooldown period. The time before that point, t < T, is the period we will use to construct a synthetic control from the donors Xi(t) to match R(t).

Generally we like to use around one year of historical data leading up to the test to construct the synthetic control.

The model and the effect estimate

REBA finds a combination of fixed weights W1, W2, …, WN (as many weights as there are donor regions) such that

W1*X1(t) + W2*X2(t) + ... + WN*XN(t) ≈ R(t)  for t < T

Note we used the approximately equals sign “≈”. The closest match will generally still not be a perfect fit.

We place some restrictions on the weights Wi:

1. The weights take the form Wi = M * Si

2. The multiplier M is positive: M > 0

3. The Si factors form a simplex:

   1. Every one must be positive or zero: Si >= 0

   2. They must sum to 1: S1 + S2 + ... + SN = 1

These restrictions basically mean that we are combining scaled-down versions of each donor in the synthetic control, then magnifying the result to reach the scale of R(t). The final synthetic control lies in the convex hull of the scaled donors M * Xi(t).

Supposing we’ve found a good combination of weights, define the synthetic control S(t) to be that weighted sum of donor regions:

S(t) = W1*X1(t) + W2*X2(t) + ... + WN*XN(t)

The quantified effect of the intervention is then the total difference between the target R(t) and the synthetic control S(t) from the test start at t = T until the end of the test and any included cooldown period. We also need to multiply that by the number of test regions to account for the fact that R(t) was an average of the actual test regions. So, in the end,

effect = (number of test regions) * sum_{t >= T} ( R(t) - S(t) )

When calculating this, we account for the imperfect match between the test and synthetic control by adding simulated errors into this S(t).

Implementation as a Bayesian model

We use a Gaussian likelihood to match R(t) and S(t).

The Stan code for the first part of the Bayesian model roughly looks like this:

data {
  int n_timesteps;  // pre-intervention timesteps
  int n_donor_regions;
  vector[n_timesteps] R;
  matrix[n_timesteps, n_donor_regions] X;

  real alpha;
}

parameters {
  simplex[n_donor_regions] S;
  real<lower=0> M;  // <lower=0> means this is bounded below by zero
  real<lower=0> sigma;
}

transformed parameters {
  // the weights
  vector[n_donor_regions] W = M * S;
  
  // the synthetic control
  vector[n_timesteps] S = X * W;
}

model {
  // priors
  S ~ dirichlet(alpha);

  // plus other priors on sigma and M

  // likelihood
  R ~ normal(S, sigma);
}

Note that we have a parameter sigma that measures the error between the intervened time series R and the synthetic control S over these pre-test time steps. This parameter gets a weak exponential prior that depends on the scale of the data.

There is also a hyperparameter alpha that controls the prior sparsity of the simplex S. This hyperparameter depends on the number of donor regions and is tuned so that every donor region has some chance of making up a large chunk (10% to 50%) of the synthetic control.

The prior on the multiplier parameter M is multiplicatively centered at 1 and has most of its density between 1/1.5 and 1.5.

Then, to calculate the effect estimate, we pass in the during- and post-intervention data, and simulate using the parameters we estimated in the previous step:

data {
  int n_test_timesteps;   // test timesteps
  int n_exposed_regions;  // number of regions with the intervention
  int n_donor_regions;
  vector[n_test_timesteps] R_test;
  matrix[n_test_timesteps, n_donor_regions] X_test;
}

generated quantities {
  // build synthetic control in the test timesteps
  vector[n_test_timesteps] S_test = X_test * W;

  // add simulated errors to the synthetic control
  vector[n_test_timesteps] S_test_sim = normal_rng(S_test, sigma);

  // calc the measured effect
  real effect = n_exposed_regions * sum(R_test - S_test_sim);
}

The quantity effect is the model’s estimate for the effect of the intervention.