first of all, thank you very much for this awesome library. It is really great and I enjoy using it. I am still pretty new to Bayesian Modelling so I hope my questions are not too stupid.
I have some panel data about monthly sales (dependent variable) in 20 cities for 5 years as well as marketing and other activities done in those (independent variables/control variables). For each city I also have the corresponding county and state information as an additional variable. The sales variable represents the number of packages of a specific product sold.
I would like to use this data for marketing mix modelling and based on my research hierarchical bayesian regression is the way to go. So I started implementing the model and currently my model code looks as follows:
# Hyperprios mu_0 = pm.Normal('mu_intercept', mu=5, sigma=1) mu_1 = pm.HalfCauchy('mu_act1', beta=0.5) mu_2 = pm.HalfCauchy('mu_act2', beta=0.5) mu_3 = pm.HalfCauchy('mu_act3', beta=0.5) mu_4 = pm.HalfCauchy('mu_act4', beta=0.5) mu_23 = pm.HalfCauchy('mu_sales_lag_1', beta=1) mu_26 = pm.HalfCauchy('mu_ext', beta=1) mu_27 = pm.Normal('mu_market', mu=5, sigma=1) # Priors, the coefficients are expected to be rather small, that's why the HalfNormal distribution was chosen beta_0 = pm.HalfNormal('intercept', sigma=mu_0, shape=n_cities) # marketing activities beta_1 = pm.HalfNormal('act1', sigma= mu_1, shape=n_cities) beta_2 = pm.HalfNormal('act2', sigma=mu_2, shape=n_cities) beta_3 = pm.HalfNormal('act3', sigma=mu_3, shape=n_cities) beta_4 = pm.HalfNormal('act4', sigma=mu_4, shape=n_cities) # control variables beta_23 = pm.HalfNormal('sales_lag_1', sigma=mu_23, shape=n_cities) beta_26 = pm.HalfNormal('ext', sigma=mu_26, shape=n_cities) beta_27 = pm.HalfNormal('market', sigma=mu_27, shape=n_cities) # data & hierarchical index of city level mu = beta_0[city_idx] + beta_1[city_idx] * X['act1']+ beta_2[city_idx] * X['act2'] + beta_3[city_idx] * X['act3']+ beta_4[city_idx] * X['act4'] + beta_23[city_idx] * X['sales_lag_1'] + beta_26[city_idx] * X['ext'] + beta_27[city_idx] * X['market'] Y_obs = pm.Poisson('Y_obs', mu=mu, observed=y)
The model works and I get the coefficients for all cities in the data, but I would like to report the coefficients on different hierarchical levels as well, e.g. per county and state, and also for the whole country. Could you point me into the right direction on how to best do that? Also, how could I get the coefficients for the national level?
Also, of course, any feedback on the code above is highly appreciated
Thank you very much in advance!