Bayesian random intercept NB mixed model in Python using pymc3

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.21 Bayesian random intercept negative binomial mixed model in Python using pymc3

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

import numpy as np
import pymc3 as pm
import statsmodels.api as sm

from scipy.stats import norm, uniform, nbinom

# Data
np.random.seed(1656) # set seed to replicate example
N = 2000 # number of obs in model
NGroups = 10

x1 = uniform.rvs(size=N)
x2 = uniform.rvs(size=N)

Groups = np.array([200 * [i] for i in range(NGroups)]).flatten()
a = norm.rvs(loc=0, scale=0.5, size=NGroups)
eta = 1 + 0.2 * x1 - 0.75 * x2 + a[list(Groups)]
mu = np.exp(eta)

y = nbinom.rvs(mu, 0.5)

with pm.Model() as model:
# Define priors
alpha = 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.NegativeBinomial('y', mu=np.exp(eta), alpha=alpha, 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)

# Print summary to screen
pm.summary(trace)

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

GET SOURCE

Output on screen:

beta1:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
1.047 0.171 0.005 [0.720, 1.404]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
0.707 0.937 1.048 1.153 1.397

beta2:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
0.216 0.072 0.001 [0.082, 0.364]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
0.076 0.168 0.216 0.264 0.359

beta3:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
-0.822 0.073 0.001 [-0.959, -0.670]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
-0.968 -0.870 -0.821 -0.773 -0.678

a_param:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
0.343 0.172 0.005 [-0.009, 0.675]
0.238 0.173 0.005 [-0.096, 0.593]
0.103 0.175 0.005 [-0.243, 0.458]
-0.984 0.182 0.005 [-1.350, -0.631]
0.195 0.173 0.005 [-0.144, 0.544]
-0.328 0.176 0.005 [-0.682, 0.030]
0.224 0.172 0.005 [-0.126, 0.556]
0.452 0.172 0.005 [0.090, 0.784]
0.081 0.173 0.005 [-0.279, 0.412]
-0.302 0.174 0.005 [-0.671, 0.024]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
-0.004 0.235 0.342 0.453 0.682
-0.110 0.128 0.239 0.347 0.582
-0.247 -0.007 0.103 0.212 0.455
-1.353 -1.099 -0.981 -0.865 -0.632
-0.154 0.085 0.195 0.304 0.537
-0.688 -0.440 -0.328 -0.215 0.026
-0.120 0.115 0.223 0.334 0.568
0.101 0.342 0.452 0.560 0.797
-0.269 -0.029 0.079 0.192 0.424
-0.652 -0.414 -0.300 -0.192 0.048

sigma:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
2.456 0.183 0.006 [2.118, 2.824]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
2.134 2.329 2.446 2.570 2.854

sigma_a:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
0.507 0.151 0.003 [0.273, 0.807]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
0.306 0.404 0.476 0.575 0.893

HSI

HSI