## Hierarchical Regression

Hierarchical models enable sharing of features among groups. The parameters of the model are assumed to be sampled from a common distribution that models similarity between groups. The figure below shows three different scenarios that illustrate the benefit of hierarchical modeling. In the figure on the left, we have a single set of parameters $\theta$ that model the entire sequence of observations referred to as a pooled model. Here any variation in data is not modelled explicitly since we are assuming a common set of parameters that give rise to the data. On the other hand, we have an unpooled scenario where we model a different set of parameters for each observation. In the unpooled case, we are assuming that there is no sharing of parameters between groups of observations and that each parameter is independent. The hierarchical model combines the best of both worlds: it assumes that there's a common distribution from which individual parameters are sampled and therefore captures similarities between groups.



In Bayesian Hierarchical Regression, we can assign priors on model parameters and use MCMC sampling to infer posterior distribution. We use the radon dataset to regress radon gas levels in houses of different counties based on the floor number (in particular if there's a basement or not). Thus our regression model will look like the following:
\begin{equation}
    \alpha_c \sim N(\mu_a, \sigma_{a}^{2})\\
    \beta_c \sim N(\mu_{\beta}, \sigma_{\beta}^{2})\\
    \mathrm{radon}_c = \alpha_c + \beta_c \times \mathrm{floor}_{i,c} + \epsilon_c
\end{equation}

Notice that subscript $c$ indicates a county, thus we are learning an intercept and a slope for each county sampled from a shared Gaussian distribution. Thus, we are assuming a hierarchical model in which our parameters ($\alpha_c$ and $\beta_c$) are sampled from a common distribution. Having specified the model, we can implement it using PyMC3.


```python
%matplotlib inline

import numpy as np
import pandas as pd

import seaborn as sns
import matplotlib.pyplot as plt

import pymc3 as pm 
```


```python
#load data
data = pd.read_csv('./radon.txt')
data.head()
```




<div>
<style>
    .dataframe thead tr:only-child th {
        text-align: right;
    }

    .dataframe thead th {
        text-align: left;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Unnamed: 0</th>
      <th>idnum</th>
      <th>state</th>
      <th>state2</th>
      <th>stfips</th>
      <th>zip</th>
      <th>region</th>
      <th>typebldg</th>
      <th>floor</th>
      <th>room</th>
      <th>...</th>
      <th>pcterr</th>
      <th>adjwt</th>
      <th>dupflag</th>
      <th>zipflag</th>
      <th>cntyfips</th>
      <th>county</th>
      <th>fips</th>
      <th>Uppm</th>
      <th>county_code</th>
      <th>log_radon</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0</td>
      <td>5081.0</td>
      <td>MN</td>
      <td>MN</td>
      <td>27.0</td>
      <td>55735</td>
      <td>5.0</td>
      <td>1.0</td>
      <td>1.0</td>
      <td>3.0</td>
      <td>...</td>
      <td>9.7</td>
      <td>1146.499190</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>1.0</td>
      <td>AITKIN</td>
      <td>27001.0</td>
      <td>0.502054</td>
      <td>0</td>
      <td>0.832909</td>
    </tr>
    <tr>
      <th>1</th>
      <td>1</td>
      <td>5082.0</td>
      <td>MN</td>
      <td>MN</td>
      <td>27.0</td>
      <td>55748</td>
      <td>5.0</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>4.0</td>
      <td>...</td>
      <td>14.5</td>
      <td>471.366223</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>1.0</td>
      <td>AITKIN</td>
      <td>27001.0</td>
      <td>0.502054</td>
      <td>0</td>
      <td>0.832909</td>
    </tr>
    <tr>
      <th>2</th>
      <td>2</td>
      <td>5083.0</td>
      <td>MN</td>
      <td>MN</td>
      <td>27.0</td>
      <td>55748</td>
      <td>5.0</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>4.0</td>
      <td>...</td>
      <td>9.6</td>
      <td>433.316718</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>1.0</td>
      <td>AITKIN</td>
      <td>27001.0</td>
      <td>0.502054</td>
      <td>0</td>
      <td>1.098612</td>
    </tr>
    <tr>
      <th>3</th>
      <td>3</td>
      <td>5084.0</td>
      <td>MN</td>
      <td>MN</td>
      <td>27.0</td>
      <td>56469</td>
      <td>5.0</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>4.0</td>
      <td>...</td>
      <td>24.3</td>
      <td>461.623670</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>1.0</td>
      <td>AITKIN</td>
      <td>27001.0</td>
      <td>0.502054</td>
      <td>0</td>
      <td>0.095310</td>
    </tr>
    <tr>
      <th>4</th>
      <td>4</td>
      <td>5085.0</td>
      <td>MN</td>
      <td>MN</td>
      <td>27.0</td>
      <td>55011</td>
      <td>3.0</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>4.0</td>
      <td>...</td>
      <td>13.8</td>
      <td>433.316718</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>3.0</td>
      <td>ANOKA</td>
      <td>27003.0</td>
      <td>0.428565</td>
      <td>1</td>
      <td>1.163151</td>
    </tr>
  </tbody>
</table>
<p>5 rows × 30 columns</p>
</div>




```python
county_names = data.county.unique()
county_idx = data['county_code'].values
```


```python
with pm.Model() as hierarchical_model:
        
    # Hyperpriors
    mu_a = pm.Normal('mu_alpha', mu=0., sd=100**2)
    sigma_a = pm.Uniform('sigma_alpha', lower=0, upper=100)
    mu_b = pm.Normal('mu_beta', mu=0., sd=100**2)
    sigma_b = pm.Uniform('sigma_beta', lower=0, upper=100)
    
    # Intercept for each county, distributed around group mean mu_a
    a = pm.Normal('alpha', mu=mu_a, sd=sigma_a, shape=len(data.county.unique()))
    # Intercept for each county, distributed around group mean mu_a
    b = pm.Normal('beta', mu=mu_b, sd=sigma_b, shape=len(data.county.unique()))
    
    # Model error
    eps = pm.Uniform('eps', lower=0, upper=100)
    
    # Expected value
    radon_est = a[county_idx] + b[county_idx] * data.floor.values
    
    # Data likelihood
    y_like = pm.Normal('y_like', mu=radon_est, sd=eps, observed=data.log_radon)
```

Having specified the graphical model, we can run inference using the state of the art No-U-Turn Sampler (NUTS) initialized with ADVI.


```python
with hierarchical_model:
    # Use ADVI for initialization
    mu, sds, elbo = pm.variational.advi(n=100000)
    step = pm.NUTS(scaling=hierarchical_model.dict_to_array(sds)**2, is_cov=True)
    hierarchical_trace = pm.sample(5000, step, start=mu)
```

    /usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:3: DeprecationWarning: Old ADVI interface and sample_vp is deprecated and will be removed in future, use pm.fit and pm.sample_approx instead
      This is separate from the ipykernel package so we can avoid doing imports until
    Average ELBO = -1,118.2:  49%|████▉     | 48907/100000 [00:15<00:16, 3116.50it/s]Median ELBO converged.
    Finished [100%]: Average ELBO = -1,103.6
    
    WARNING (theano.tensor.blas): We did not found a dynamic library into the library_dir of the library we use for blas. If you use ATLAS, make sure to compile it with dynamics library.
    100%|█████████▉| 5494/5500 [00:37<00:00, 143.16it/s]/usr/local/lib/python2.7/dist-packages/pymc3/step_methods/hmc/nuts.py:440: UserWarning: The acceptance probability in chain 0 does not match the target. It is 0.667219685913, but should be close to 0.8. Try to increase the number of tuning steps.
      % (self._chain_id, mean_accept, target_accept))
    100%|██████████| 5500/5500 [00:37<00:00, 144.83it/s]


Let's examine the trace plots of all our latent variables in the hierarchical regression model:



```python
pm.traceplot(hierarchical_trace[500:])
plt.show()
```

From the traceplots, we can see convergence in our posterior distributions for $\alpha_c$ and $\beta_c$ indicating different intercepts and slopes for different counties. In addition, we also recover the posterior distribution of the shared parameters. $\mu_{\alpha}$ tells us that the group mean of log radon levels is close to $1.5$, while $\mu_{\beta}$ tells us that the slope is negative with a mean of $-0.65$ and thus having no basement decreases radon levels.


```python

```

23e1142b6897a9cd042e46d053536aa6c8e5902e
ipynb
Jupyter Notebook
chp01/hierarchical_regression.ipynb
gerket/experiments_with_python
5dd6dbd69deaaa318bfa7d2c3c9f7fae6220c460
MIT
2017-08-22T13:14:54.000Z
2022-03-28T17:56:59.000Z
chp01/hierarchical_regression.ipynb
gerket/experiments_with_python
5dd6dbd69deaaa318bfa7d2c3c9f7fae6220c460
MIT
2017-07-31T00:52:36.000Z
2018-10-01T14:29:51.000Z
chp01/hierarchical_regression.ipynb
gerket/experiments_with_python
5dd6dbd69deaaa318bfa7d2c3c9f7fae6220c460
MIT
2017-08-23T08:08:32.000Z
2022-03-09T07:04:01.000Z
Qwen/Qwen-72B
1. YES
2. YES