Bayesian analysis of Multi-cell experiment

Questions v3
1

I have an experiment running with four cells divided by two axes; the axes are (1) experiment group (test, control) and demand period (high, low). There are four total combinations therein, test-high, test-low, control-high & control-low.

The unit is sales per hour. The test group has been exposed to a marketing strategy that nudges consumers towards the high-demand period, whereas the control has a time-agnostic marketing stimulus.

To complicate matters, the high demand periods are a small fraction of time in comparison to the low-demand periods. (ex low: Monday-Friday, high: Saturday-Sunday.)

Iā€™m interested in comparing the lift test (high_demand_test / low_demand_test) to lift control (high_demand_control / low_demand_control).

Below is some simulated data:

def gen_data(test=True, high_demand=True, size=10):
  def gen_single(test=test, high_demand=high_demand):
    val = np.random.normal(25, 3)
    if test:
      val += np.random.normal(5, 1)
    if high_demand: 
      val += np.random.normal(3, 1)
    return val 
  return [gen_single() for s in range(size)]

config = [# test, high_demand, size
          (True,  True, 10),
          (True, False, 30),
          (False, True, 10),
          (False, False, 30)
          ]

data = []
for test, high_demand, size in config:
   data.extend(gen_data(test=test, high_demand=high_demand, size=size))

mu_pooled = np.mean(data)
std_pooled = np.std(data)

And the analysis:

with pm.Model() as model:
  mu_arr = pm.Normal('mu_i', mu=mu_pooled, sigma=std_pooled, shape=4)
  error = pm.InverseGamma('error',3,2)
  
  base = 0
  obs_arr = []
  for idx, size in enumerate(config):
    elements = data[base : base+size[2] -1]
    base += size[2]
    obs_arr.append(np.mean(elements))
  
  test_ratio = pm.Deterministic('test_ratio', mu_arr[0] / mu_arr[1])
  control_ratio = pm.Deterministic('control_ratio', mu_arr[2] / mu_arr[3])
  delta = pm.Deterministic('delta', test_ratio - control_ratio)

  lik = Normal('likelihood', mu=mu_arr, sigma=error, observed=obs_arr)
  trace = pm.sample(chains=4 )

The divergences are in the range of 1-5 on two chains using NUTS. (See link for simplicity.)

Can you recommend some advice on parametrizing this model?

2

The inverse-Gamma is usually used as a variance prior instead of a sigma prior; when I changed sigma=error to sigma=error**0.5, the divergences went away.

4 Likes