Bayesian random intercept binary logistic model in Python using pymc3

GET SOURCE

From: Bayesian Models for Astrophysical Data, Cambridge Univ. Press

you are kindly asked to include the complete citation if you used this material in a publication

Code 8.7 Bayesian random intercept binary logistic model in Python using pymc3

===============================================================

import numpy as np
import pymc3 as pm

from scipy.stats import norm, uniform, bernoulli

# Data
np.random.seed(13531) # set seed to replicate example
N = 4000 # number of obs in model
NGroups = 20

x1 = uniform.rvs(size=N)
x2 = uniform.rvs(size=N)
Groups = np.array([200 * [i] for i in range(20)]).flatten()

a = norm.rvs(loc=0, scale=0.5, size=NGroups)
eta = 1 + 0.2 * x1 - 0.75 * x2 + a[list(Groups)]
mu = 1.0/(1.0 + np.exp(-eta))
y = bernoulli.rvs(mu, size=N)

with pm.Model() as model:
# Define priors
sigma = pm.Uniform('sigma', 0, 100)
sigma_a = pm.Uniform('sigma_a', 0, 10)
beta1 = pm.Normal('beta1', 0, sd=100)
beta2 = pm.Normal('beta2', 0, sd=100)
beta3 = pm.Normal('beta3', 0, sd=100)

# priors for random intercept (RI) parameters
a_param = pm.Normal('a_param',
np.repeat(0, NGroups), # mean
sd=np.repeat(sigma_a, NGroups), # standard deviation
shape=NGroups) # number of RI parameters

eta = beta1 + beta2 * x1 + beta3 * x2 + a_param[Groups]

# Define likelihood
y = pm.Normal('y', mu=1.0/(1.0 + np.exp(-eta)), sd=sigma, observed=y)

# Fit
start = pm.find_MAP()    # Find starting value by optimization
step = pm.NUTS(scaling=start) # Initiate sampling
trace = pm.sample(7000, step, start=start, progressbar=False)

# Print summary to screen
pm.summary(trace)

====================================================================

Output on screen:

sigma_interval:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------

-5.397 0.011 0.001 [-5.419, -5.375]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
-5.418 -5.405 -5.398 -5.389 -5.374

sigma_a_interval:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
-2.842 0.205 0.012 [-3.293, -2.461]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
-3.243 -2.976 -2.845 -2.718 -2.399

beta1:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
1.140 0.160 0.004 [0.824, 1.450]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
0.831 1.032 1.137 1.245 1.461

beta2:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
0.223 0.121 0.001 [-0.015, 0.459]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
-0.014 0.143 0.224 0.305 0.461

beta3:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
-0.861 0.126 0.001 [-1.125, -0.629]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
-1.112 -0.945 -0.860 -0.776 -0.612

a_param:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
-0.343 0.183 0.004 [-0.690, 0.024]
-0.295 0.183 0.004 [-0.671, 0.042]
-0.192 0.187 0.004 [-0.552, 0.181]
-0.341 0.184 0.004 [-0.710, 0.014]
0.103 0.193 0.004 [-0.262, 0.489]
-0.410 0.181 0.004 [-0.766, -0.058]
0.004 0.193 0.004 [-0.382, 0.376]
-0.420 0.181 0.004 [-0.784, -0.066]
-1.063 0.180 0.004 [-1.432, -0.724]
-0.724 0.182 0.004 [-1.074, -0.359]
-0.047 0.187 0.004 [-0.379, 0.364]
0.345 0.205 0.004 [-0.053, 0.746]
0.395 0.210 0.004 [-0.006, 0.817]
0.890 0.250 0.005 [0.394, 1.385]
0.352 0.207 0.004 [-0.044, 0.754]
0.020 0.195 0.004 [-0.357, 0.409]
0.109 0.198 0.004 [-0.266, 0.510]
0.632 0.224 0.004 [0.201, 1.088]
0.859 0.249 0.005 [0.398, 1.377]
0.197 0.199 0.004 [-0.172, 0.601]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
-0.698 -0.464 -0.343 -0.220 0.018
-0.656 -0.418 -0.292 -0.170 0.062
-0.560 -0.320 -0.192 -0.065 0.175
-0.708 -0.459 -0.340 -0.220 0.018
-0.269 -0.026 0.100 0.229 0.486
-0.764 -0.531 -0.410 -0.289 -0.055
-0.367 -0.126 0.003 0.131 0.393
-0.787 -0.538 -0.417 -0.299 -0.067
-1.426 -1.181 -1.062 -0.941 -0.716
-1.083 -0.847 -0.721 -0.602 -0.362
-0.414 -0.172 -0.052 0.078 0.331
-0.042 0.205 0.338 0.477 0.765
-0.008 0.252 0.390 0.530 0.815
0.419 0.720 0.881 1.052 1.418
-0.034 0.211 0.348 0.487 0.770
-0.356 -0.109 0.018 0.147 0.410
-0.269 -0.022 0.104 0.239 0.509
0.206 0.480 0.625 0.778 1.100
0.401 0.688 0.845 1.017 1.385
-0.191 0.062 0.198 0.331 0.588

sigma:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
0.451 0.005 0.000 [0.441, 0.461]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|

0.441 0.447 0.451 0.455 0.461

sigma_a:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
0.561 0.112 0.007 [0.353, 0.777]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
0.376 0.485 0.549 0.619 0.833

HSI

HSI