I am trying to implement the suggestion by Gelman et al. (2008) to use weakly informative priors for bayesian logistic regression.


I am trying to achieve this with pymc3, and I am unsure on how to handle the standardization Gelman proposes for binary input variables:


Standardizing the data is in my case simply data['score'] - data['score'].mean(), but how to deal with that in the model specification?

Here is some working example code, but without the standardization.

import scipy as sp
from scipy import stats
import pandas as pd
from patsy import dmatrix
import pymc3 as pm
import theano.tensor as Tht
import matplotlib.pyplot as plt
import seaborn.apionly as sns

# data generation
labels = list('ABCDE')
true_freqs = [0.9,0.8,0.5,0.8,0.4]
Neach = 20

data = []
for label,freq in zip(labels,true_freqs):
    d = sp.rand(Neach) < freq

data = pd.concat(data)
Nobs = data.shape[0]
data.index = range(Nobs)

# bayesian logistic regression with pymc3
nSamples = 2000
burn = 1000
nsim = nSamples - burn

model_matrix = dmatrix(' ~ label', data=data)
x0,x1,x2,x3,x4 = [sp.array(model_matrix[:,i]).astype('int64') for i in range(5)]

with pm.Model() as model:
    # betas
    beta0 = pm.Cauchy('beta0', 0., 10.0)
    beta1 = pm.Cauchy('beta1', 0., 2.5)
    beta2 = pm.Cauchy('beta2', 0., 2.5)
    beta3 = pm.Cauchy('beta3', 0., 2.5)
    beta4 = pm.Cauchy('beta4', 0., 2.5)

    # logit
    logit_p =  (beta0*x0 + beta1*x1 + beta2*x2 + beta3*x3 + beta4*x4)
    p = Tht.exp(logit_p) / (1 + Tht.exp(logit_p))

    # likelihood
    likelihood = pm.Binomial('likelihood',n=1,p=p,observed=data['score'])

    # inference
    start = pm.find_MAP()
    step = pm.NUTS(scaling=start)
    trace = pm.sample(nSamples, step, progressbar=True)

### inspect

# credible intervals from posterior
# get simulated betas
varnames = ['beta0','beta1','beta2','beta3','beta4']
betas_mc = sp.zeros((len(varnames),nsim))
for i,varname in enumerate(varnames):
    betas_mc[i,:] = trace.get_values(varname)[nSamples-burn:]
betas_mc = sp.matrix(betas_mc)

# calculate posteriors
ps_sim = sp.zeros((Nobs,nsim))
X = sp.matrix(model_matrix)

for i in range(nsim):
    ps = X * betas_mc[:,i]
    ps_sim[:,i] = stats.logistic.cdf(ps).flatten()

# plotting
fig, axes = plt.subplots()

# calculate and add CrIs
for i, (label, group) in enumerate(data.groupby('label')):
    posterior_mc = ps_sim[group.index][0,:]
    CrIs = sp.percentile(posterior_mc,(2.5,97.5))

Here is the output: enter image description here enter image description here

If I simply plug in the centered values as such

likelihood = pm.Binomial('likelihood',n=1,p=p,observed=data['score'] - data['score'].mean())

then I get oddly looking chains and generally weird results, I guess because this does something completely different than I expected and leads to some numerical problems.

Long story short: How to adjust a pmc3 model for binary input, but not ranging from 0 to 1?


1 Answer 1


Gelman is talking about centering the binary predictors ("input"), not the binary response ("output"). You centered the binary response.

  • $\begingroup$ can you please elaborate? I am not entirely sure how the terms input and ouput are used here. The only other thing I can see that is binary and could be centered would be the values in the model matrix, but that would either make mathematically no sense at all or I don't get the point. $\endgroup$
    – grg rsr
    Aug 11, 2017 at 10:36
  • $\begingroup$ @grgrsr Centering values in the model matrix is indeed exactly what Gelman is talking about. As for what the point is, there is a huge literature on this in statistics: google.com/search?q=regression+center+predictors $\endgroup$ Aug 11, 2017 at 14:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.