Question

This uses the Gelman radon data set and is based off this IPython notebook . Radon levels were measured in houses from all counties in several states. Here we want to know if the preence of a basement affects the level of radon, and if this is affected by which county the house is located in.

The data set provided is just for the state of Minnesota, which has 85 counties with 2 to 116 measurements per county. We only need 3 columns for this example county , log_radon , floor , where floor=0 indicates that there is a basement.

We will perfrom simple linear regression on log_radon as a function of county and floor.

radon = pd.read_csv('radon.csv')[['county', 'floor', 'log_radon']]
radon.head()

	county	floor	log_radon
0	AITKIN	1	0.832909
1	AITKIN	0	0.832909
2	AITKIN	0	1.098612
3	AITKIN	0	0.095310
4	ANOKA	0	1.163151

We will be creating lots of similar models, so it is worth wrapping definitions into a function to avoid repetition.

def make_model(x, y):
    # define priors
    a = pymc.Normal('slope', mu=0, tau=1.0/10**2)
    b = pymc.Normal('intercept', mu=0, tau=1.0/10**2)
    tau = pymc.Gamma("tau", alpha=0.1, beta=0.1)

    # define likelihood
    @pymc.deterministic
    def mu(a=a, b=b, x=x):
        return a*x + b

    y = pymc.Normal('y', mu=mu, tau=tau, value=y, observed=True)

    return locals()

Pooled model

If we pool the data across counties, this is the same as the simple linear regression model.

plt.scatter(radon.floor, radon.log_radon)
plt.xticks([0, 1], ['Basement', 'No basement'], fontsize=20);

m = pymc.Model(make_model(radon.floor, radon.log_radon))
mc = pymc.MCMC(m)
mc.sample(iter=1100, burn=1000)

[-----------------100%-----------------] 1100 of 1100 complete in 5.2 sec

abar = mc.stats()['slope']['mean']
bbar = mc.stats()['intercept']['mean']
radon.plot(x='floor', y='log_radon', kind='scatter', s=50);
xp = np.array([0, 1])
plt.plot(mc.trace('slope')()*xp[:, None] + mc.trace('intercept')(), c='red', alpha=0.1)
plt.plot(xp, abar*xp + bbar, linewidth=2, c='red');

Individual couty model

Inidividual couty models are done in the same way, except that we create a model for each county.

n = 0
i_as = []
i_bs = []
for i, group in radon.groupby('county'):

    m = pymc.Model(make_model(group.floor, group.log_radon))
    mc = pymc.MCMC(m)
    mc.sample(iter=1100, burn=1000)

    abar = mc.stats()['slope']['mean']
    bbar = mc.stats()['intercept']['mean']
    group.plot(x='floor', y='log_radon', kind='scatter', s=50);
    xp = np.array([0, 1])
    plt.plot(mc.trace('slope')()*xp[:, None] + mc.trace('intercept')(), c='red', alpha=0.1)
    plt.plot(xp, abar*xp + bbar, linewidth=2, c='red');
    plt.title(i)

    n += 1
    if n > 3:
        break

[-----------------100%-----------------] 1100 of 1100 complete in 3.0 sec

Hiearchical model

With a hierarchical model, there is an \(a_c\) and a \(b_c\) for each county \(c\) just as in the individual couty model, but they are no longer indepnedent but assumed to come from a common group distribution

\[\begin{split}a_c \sim \mathcal{N}(\mu_a, \sigma_a^2) \\ b_c \sim \mathcal{N}(\mu_b, \sigma_b^2)\end{split}\]

we furhter assume that the hyperparameters come from the following distributions

\[\begin{split}\mu_a \sim \mathcal{N}(0, 100^2) \\ \sigma_a \sim \mathcal{U}(0, 100) \\ \mu_b \sim \mathcal{N}(0, 100^2) \\ \sigma_b \sim \mathcal{U}(0, 100)\end{split}\]

county = pd.Categorical(radon['county']).codes

# County hyperpriors
mu_a = pymc.Normal('mu_a', mu=0, tau=1.0/100**2)
sigma_a = pymc.Uniform('sigma_a', lower=0, upper=100)
mu_b = pymc.Normal('mu_b', mu=0, tau=1.0/100**2)
sigma_b = pymc.Uniform('sigma_b', lower=0, upper=100)

# County slopes and intercepts
a = pymc.Normal('slope', mu=mu_a, tau=1.0/sigma_a**2, size=len(set(county)))
b = pymc.Normal('intercept', mu=mu_b, tau=1.0/sigma_b**2, size=len(set(county)))

# Houseehold priors
tau = pymc.Gamma("tau", alpha=0.1, beta=0.1)

@pymc.deterministic
def mu(a=a, b=b, x=radon.floor):
    return a[county]*x + b[county]

y = pymc.Normal('y', mu=mu, tau=tau, value=radon.log_radon, observed=True)

m = pymc.Model([y, mu, tau, a, b])
mc = pymc.MCMC(m)
mc.sample(iter=110000, burn=100000)

[-----------------100%-----------------] 110000 of 110000 complete in 235.1 sec

abar = a.stats()['mean']
bbar = b.stats()['mean']

xp = np.array([0, 1])
for i, (a, b) in enumerate(zip(abar, bbar)):
    plt.figure()
    idx = county == i
    plt.scatter(radon.floor[idx], radon.log_radon[idx])
    plt.plot(xp, a*xp + b, c='red');
    plt.title(radon.county[idx].unique()[0])
    if i >= 3:
        break