Random intercept binomial logistic model in Python using Stan

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.11 Random intercept binomial logistic model in Python using Stan

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

import numpy as np
import statsmodels.api as sm
import pystan

from scipy.stats import norm, uniform, bernoulli

y = [6,11,9,13,17,21,8,10,15,19,7,12,8,5,13,17,5,12,9,10]
m = [45,54,39,47,29,44,36,57,62,55,66,48,49,39,28,35,39,43,50,36]
x1 = [1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0]
x2 = [1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0]

Groups = range(len(y))

X = sm.add_constant(np.column_stack((x1,x2)))
K = X.shape[1]

model_data = {}
model_data['Y'] = y # response
model_data['X'] = X # covariates
model_data['K'] = K # num. betas
model_data['m'] = m # binomial denominator
model_data['N'] = len(y) # sample size
model_data['re'] = Groups # random effects
model_data['b0'] = np.repeat(0, K)
model_data['B0'] = np.diag(np.repeat(100, K))
model_data['a0'] = np.repeat(0, len(y))
model_data['A0'] = np.diag(np.repeat(1, len(y)))

# Fit
stan_code = """
data{
int<lower=0> N;
int<lower=0> K;
matrix[N, K] X;
int<lower=0> Y[N];
int re[N];
int m[N];
vector[K] b0;
matrix[K, K] B0;
vector[N] a0;
matrix[N, N] A0;
}
parameters{
vector[K] beta;
vector[N] a;
real<lower=0> sigma;
}
transformed parameters{
vector[N] eta;
vector[N] p;

eta = X * beta;
for (i in 1:N){
p[i] = eta[i] + a[re[i]+1];
}
}
model{
sigma ~ cauchy(0, 25);

beta ~ multi_normal(b0, B0);
a ~ multi_normal(a0, sigma * A0);

Y ~ binomial_logit(m, p);
}
"""

fit = pystan.stan(model_code=stan_code, data=model_data, iter=5000, chains=3, thin=10,
warmup=4000, n_jobs=3)

# Output
nlines = 29 # number of lines in screen output

output = str(fit).split('\n')

for item in output[:nlines]:
print(item)

GET SOURCE

Output on screen:

Inference for Stan model: anon_model_617c4756be1da3f12c74708416b43501.
3 chains, each with iter=5000; warmup=4000; thin=10;
post-warmup draws per chain=100, total post-warmup draws=300.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
beta[0] -1.11    0.02 0.29 -1.72 -1.28 -1.09 -0.92 -0.54 260 1.0
beta[1] 0.22 0.02 0.34 -0.46 7.0e-3 0.24 0.44 0.82    300 1.0
beta[2] -0.26    0.02   0.35 -1.07 -0.45 -0.25 -0.05 0.4    272 1.0
a[0] -0.57 0.02 0.43 -1.54 -0.84 -0.55 -0.27 0.21    300 0.99
a[1] -0.41 0.02 0.39 -1.19 -0.66 -0.43 -0.14 0.31 254 1.01
a[2] -0.07 0.02 0.39 -0.86 -0.31 -0.08 0.21 0.67     300 1.0
a[3]    -0.1 0.02 0.37 -0.79 -0.33 -0.1 0.13 0.65 300 1.01
a[4] 1.1    0.02 0.42   0.35 0.8 1.06 1.35   2.06 294   0.99
a[5] 0.64 0.02 0.36 -0.09 0.42 0.64 0.87 1.41 300 1.01
a[6] -0.13 0.02 0.43 -0.97 -0.41 -0.13 0.17 0.75 300 1.0
a[7] -0.53 0.02 0.41 -1.45 -0.79 -0.51 -0.25 0.25 300 1.0
a[8] 8.9e-3 0.02 0.39 -0.81   -0.25 0.02 0.27 0.77 272 1.0
a[9] 0.22 0.02 0.35 -0.39 -0.02 0.21 0.44 0.96 300   1.01
a[10] -0.59 0.02 0.39 -1.39 -0.86 -0.58   -0.31 0.13 269 0.99
a[11] -6.2e-3 0.02 0.4 -0.78 -0.28 5.1e-3 0.26 0.84 300 1.0
a[12] -0.23 0.03 0.41 -1.01 -0.48 -0.23 0.05 0.67    259 1.01
a[13] -0.56    0.02 0.41 -1.43 -0.82 -0.51 -0.29 0.19 288 1.0
a[14] 0.92 0.02 0.41 0.11 0.62 0.92 1.21 1.78 300 1.0
a[15] 0.84 0.02 0.37 0.18 0.56 0.88 1.09 1.55 300   0.99
a[16] -0.41     0.03    0.43   -1.33 -0.69 -0.41 -0.09 0.41 226 0.99
a[17] 0.09 0.02 0.36 -0.62 -0.15 0.07 0.31 0.87 278 1.0
a[18]    -0.17 0.03 0.43 -1.03 -0.44 -0.17 0.11 0.74 280 1.0
a[19] 0.11    0.02 0.35 -0.61 -0.14 0.11 0.36 0.78 300 1.0
sigma 0.5 0.01 0.23 0.19 0.33   0.45 0.62 1.17 300 0.99

HSI

HSI