Creating a three-level logistic regression model in pymc3

Few changes (note I changed yhat to theta):

theta = pm.Deterministic( 'theta', pm.invlogit(sum(mid_level[i] for i in mid_to_bot_idx)+intercept) )
yact = pm.Bernoulli( 'yact', p=theta, observed=y )

Turns out the solution was simple: it appears that any distribution (like pm.Normal) can accept a vector of means as an argument, so replacing this line

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)]

with this

mid_level = pm.Normal('mid_level',
                       mu=top_level[mid_to_top_idx],
                       tau=mid_level_tau,
                       shape=k_mid)

works. The same method can also be used to specify individual standard deviations for each of the mid-level groups.

EDIT: Fixed typo