# Hierarchical Logistic Regression. How to reparameterize?

**URL:** <https://discourse.pymc.io/t/hierarchical-logistic-regression-how-to-reparameterize/8215>\
**Category:** Questions\
**Created:** [October 25, 2021, 3:39pm UTC](https://discourse.pymc.io/t/hierarchical-logistic-regression-how-to-reparameterize/8215 "2021-10-25T15:39:07Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![vanko](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/vanko/32/4509_2.png) [@vanko](https://discourse.pymc.io/u/vanko)\
**Post date:** [October 25, 2021, 3:39pm UTC](https://discourse.pymc.io/t/hierarchical-logistic-regression-how-to-reparameterize/8215/1 "2021-10-25T15:39:08Z")

</div>

I have a simple case with 30 groups(NETWORKS) and only 1 RV(EXP\_FLAG) where dependent variable binary class (0 or 1). 65100 rows × 3 columns  
 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/4/4eec42772caa780693ed4bebb66495a211147ed0.png)

**#split variables and define network idx**  
x, y = data\_tv\_s.EXP\_FLAG, data\_tv\_s.CONV\_FLAG  
network\_idx = data\_tv\_s[‘NETWORK\_IDX’].values

**Building Individual (unpooled) model:**  
with pm.Model() as unpooled\_model:

```
# priors
alpha = pm.Normal('alpha', mu=0, sd=10, shape =len(data_tv_s.NETWORK_IDX.unique()))
beta = pm.Normal('beta', mu=0, sd=10, shape =len(data_tv_s.NETWORK_IDX.unique()))

m = alpha[network_idx] + pm.math.dot(x, beta[network_idx])

# inverse link function with alias for the Theano function
theta = pm.Deterministic('theta', pm.math.sigmoid(m)) 

# boundary decision, which is the value used to separate class of the target
bd = pm.Deterministic('bd', -alpha[network_idx]/beta[network_idx])

conv = pm.Bernoulli('conv', p=theta, observed=y)

priors_unpooled = pm.sample_prior_predictive()

# posterior/create the trace
trace_unpooled = pm.sample(chains= 2, target_accept = 0.95)

```

Takes 10hr to run and results are following:

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/2/28bfa602f1e66ca838a06c817381f8fd048ce1e3.png)

Then I’m building **HIERARCHICAL** with 4chains:  
with pm.Model() as hr\_model:

```
# hyper-priors
alpha_m = pm.Normal('alpha_m', mu=0, sd=10)
alpha_sigma = pm.HalfNormal('alpha_sigma', sd= 10)
beta_m = pm.Normal('beta_m', mu=0, sd=10)
beta_sigma = pm.HalfNormal('beta_sigma', sd=10)

# priors
alpha = pm.Normal('alpha', 
                  mu=alpha_m, 
                  sd=alpha_sigma, 
                  shape =len(data_tv_s.NETWORK_IDX.unique()))
beta = pm.Normal('beta', 
                 mu=beta_m, 
                 sd=beta_sigma, 
                 shape =len(data_tv_s.NETWORK_IDX.unique()))

m = alpha[network_idx] + pm.math.dot(x, beta[network_idx])

# inverse link function with alias for the Theano function
theta = pm.Deterministic('theta', pm.math.sigmoid(m)) 

# boundary decision, which is the value used to separate class of the target
bd = pm.Deterministic('bd', -alpha[network_idx]/beta[network_idx])

conv = pm.Bernoulli('conv', p=theta, observed=y)

priors_hr = pm.sample_prior_predictive()

# posterior/create the trace
trace_hr = pm.sample(chains= 4, target_accept = 0.95)

```

Takes for me 20 Hours to run with results:

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/5/5f3bed2d0716483b269f3313a4442d09d182c1ea.png)

**r\_hat is so high and each chain reached the maximum tree depth. What do I do wrong in here ?**  
My forest for **betas** in both unpooled and Hierarchical models for each group has wide HDI interval. Following images are from Hierarchical models:

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/b/b130c12bd8684d520d11b279a718d1648fa2160c.png)

Intercept and Variance in between groups:  
 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/1/1838b6d3ba75afdf13f38eff44a3e06d6ce68070.png)

---

<div class="post-metadata">

**Author:** ![tcapretto](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/tcapretto/32/4902_2.png) [@tcapretto](https://discourse.pymc.io/u/tcapretto)\
**Post date:** [October 26, 2021, 9:10pm UTC](https://discourse.pymc.io/t/hierarchical-logistic-regression-how-to-reparameterize/8215/2 "2021-10-26T21:10:26Z")

</div>

You could try to fit this model using [Bambi](https://bambinos.github.io/bambi/main/index.html).

The non-hierarchical version

```python
model = bmb.Model(
  "CONV_FLAG ~ EXP_FLAG + NETWORK_IDX", 
  data_tv_s,
  family="bernoulli"
)
idata = model.fit()

```

The hierarchical version (i.e. varying intercept and varying slopes)

```python
model = bmb.Model(
  "CONV_FLAG ~ EXP_FLAG + NETWORK_IDX + (EXP_FLAG|NETWORK_IDX)", 
  data_tv_s,
  family="bernoulli"
)
idata = model.fit()

```

It can also be good to select a sub-sample of the data in the first attempts to fit the model to detect possible improvements faster.