title: Creating a three-level logistic regression model in pymc3
body: <p>I'm attempting to create a three-level logistic regression model in pymc3. There is a top level, mid level, and an individual level, where the mid-level coefficients are estimated from top-level coefficients. I'm having difficulty specifying the proper data structure for the mid level, however.</p>
  
  <p>Here's my code:</p>
  
  <pre><code>with pm.Model() as model:
      # Hyperpriors
      top_level_tau = pm.HalfNormal('top_level_tau', sd=100.)
      mid_level_tau = pm.HalfNormal('mid_level_tau', sd=100.)    
  
      # Priors
      top_level = pm.Normal('top_level', mu=0., tau=top_level_tau, shape=k_top)
      mid_level = [pm.Normal('mid_level_{}'.format(j),
                             mu=top_level[mid_to_top_idx[j]],
                             tau=mid_level_tau)
                   for j in range(k_mid)]
  
      intercept = pm.Normal('intercept', mu=0., sd=100.)
  
      # Model prediction
      yhat = pm.invlogit(mid_level[mid_to_bot_idx] + intercept)
  
      # Likelihood
      yact = pm.Bernoulli('yact', p=yhat, observed=y)
  </code></pre>
  
  <p>I'm getting the error <code>"only integer arrays with one element can be converted to an index"</code> (on line 16), which I think is related to the fact that the <code>mid_level</code> variable is a list, not a proper pymc Container. (I don't see the Container class in the pymc3 source code, either.)</p>
  
  <p>Any help would be appreciated.</p>
  
  <p>Edit: Adding some mock data</p>
  
  <pre><code>y = np.array([0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0])
  mid_to_bot_idx = np.array([0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 3, 2, 3, 3, 3, 3, 2, 2, 2, 3, 3, 3, 3, 2, 3, 2, 3, 3, 3, 3, 2, 3, 2, 3, 3, 3, 3, 2, 3, 2, 3, 3, 2, 2, 3, 2, 2, 3, 3, 3, 3, 2, 2, 2, 3, 2, 3, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 2, 3, 2, 2, 3, 3, 2, 2, 3, 2])
  mid_to_top_idx = np.array([0, 0, 1, 1])
  k_top = 2
  k_mid = 4
  </code></pre>
  
  <p>Edit #2:</p>
  
  <p>There seem to be a few different ways to solve this issue, although none are completely satisfactory:</p>
  
  <p>1) One can reframe the model as:</p>
  
  <pre><code>with pm.Model() as model:
      # Hyperpriors
      top_level_tau = pm.HalfNormal('top_level_tau', sd=100.)
      mid_level_tau = pm.HalfNormal('mid_level_tau', sd=100.)    
  
      # Priors
      top_level = pm.Normal('top_level', mu=0., tau=top_level_tau, shape=k_top)
      mid_level = pm.Normal('mid_level', mu=0., tau=mid_level_tau, shape=k_top)
      intercept = pm.Normal('intercept', mu=0., sd=100.)
  
      # Model prediction
      yhat = pm.invlogit(top_level[top_to_bot_idx] + mid_level[mid_to_bot_idx] + intercept)
  
      # Likelihood
      yact = pm.Bernoulli('yact', p=yhat, observed=y)
  </code></pre>
  
  <p>This seems to work, although I can't figure out how to extend it to the case where the mid-level variance is not constant for all of the mid-level groups.</p>
  
  <p>2) One can wrap the mid-level coefficients into a Theano tensor using theano.tensor.stack: i.e.,</p>
  
  <pre><code>import theano.tensor as tt
  mid_level = tt.stack([pm.Normal('mid_level_{}'.format(j),
                             mu=top_level[mid_to_top_idx[j]],
                             tau=mid_level_tau)
                   for j in range(k_mid)])
  </code></pre>
  
  <p>But this seems to run very slowly on my actual data set (30k observations), and it makes plotting inconvenient (each of the mid_level coefficients gets its own trace using <code>pm.traceplot</code>).</p>
  
  <p>Anyway, some advice/input from the developers would be appreciated.</p>