# Managing indices/data shape in a multi-level linear regression for incremental model updating

I was hoping to do something similar to this https://github.com/pymc-devs/pymc3/blob/master/docs/source/notebooks/updating_priors.ipynb but for the model defined below (toy data is provided and variables have been semi-anonymized). I’m realizing that indexing etc is going to be kind of a headache though. Any suggestions?

``````import numpy as np
import pymc3 as pm

c_count = 2
coh_count = 58
c_c_indices = np.asarray([0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1,
0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1,
0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1])
coh_index = sorted(np.tile(np.arange(0,58), [52,]))
exposure = np.random.randint(10, 200, len(coh_index))
observed = np.random.randint(30, 400, len(coh_index))
x = np.tile(np.log(np.arange(1,53)), [58,])

with pm.Model() as model:
sigma_1 = pm.Exponential('sigma_1', .1)
sigma_2 = pm.Exponential('sigma_2', .1)
sigma_3 = pm.Exponential('sigma_3', .1)
sigma_4 = pm.Exponential('sigma_4', .1)

global_m = pm.Normal('global_m',
mu=-1,
sigma=1)
global_b = pm.Normal('global_b',
mu=-1,
sigma=1)
c_m = pm.Normal('c_m',
mu=global_m,
sigma=sigma_1,
shape=c_count)
c_b = pm.Normal('c_b',
mu=global_b,
sigma=sigma_2,
shape=c_count)
c_c_m = pm.Normal('c_c_m',
mu=c_m[c_c_indices],
sigma=sigma_3,
shape=coh_count)
c_c_b = pm.Normal('c_c_b',
mu=c_b[c_c_indices],
sigma=sigma_4,
shape=coh_count)

theta = (c_c_b[coh_index] + x * c_c_m[coh_index] + pm.math.log(exposure))
txn_like = pm.Poisson('txn_like', mu=np.exp(theta), observed=observed)

with model:
trace = pm.sample(2000)
``````