An intercept is an additive regression parameter that shifts predicted values before the effect of input variables is applied. In a hierarchical model, each county receives its own intercept, while those intercepts are sampled from a shared Gaussian distribution. section 6.3 section 1.4
The treatment connects intercepts across two modeling settings:
- Bayesian regression estimates an intercept for each county alongside a county-specific slope, allowing groups to share statistical structure without forcing identical parameters. section 6.3
- Deep neural networks include bias parameters with weight matrices, making the bias the layer-level counterpart of an additive intercept. section 1.4
County-level radon regression
A radon model predicts house radon levels from floor information, including whether a house has a basement. Each county receives an intercept and slope; hierarchical priors relate these parameters through shared distributions. section 6.3
Listing 6.2 Hierarchical Bayesian regression
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import pymc3 as pm
def main():
data = pd.read_csv('./data/radon.txt') #1
county_names = data.county.unique()
county_idx = data['county_code'].values
with pm.Model() as hierarchical_model:
mu_a = pm.Normal('mu_alpha', mu=0., sd=100**2) #2
sigma_a = pm.Uniform('sigma_alpha', lower=0, #2
upper=100) #2
mu_b = pm.Normal('mu_beta', mu=0., sd=100**2) #2
sigma_b = pm.Uniform('sigma_beta', lower=0, #2
upper=100) #2
a = pm.Normal('alpha', mu=mu_a, sd=sigma_a,
shape=len(data.county.unique())) #3
b = pm.Normal('beta', mu=mu_b, sd=sigma_b,
shape=len(data.county.unique())) #4
eps = pm.Uniform('eps', lower=0, upper=100) #5
radon_est = a[county_idx] + b[county_idx] *
data.floor.values #6
y_like = pm.Normal('y_like', mu=radon_est, sd=eps,
observed=data.log_radon) #7
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)
pm.traceplot(hierarchical_trace[500:])
plt.show()
if __name__ == "__main__":
main()
The implementation defines hyperpriors for the intercept distribution and creates one Gaussian alpha parameter for each county. section 6.3
