Seeds: Random effect logistic regression

Questions
14

Folks I am facing similar issue, my code (@cluhmann and @ckrapu

import pymc3 as pm

import arviz as az

# ri~ Binomial(pi, ni)

#

# logit(pi) = α0 + α1x1i + α2x2i + α12x1ix2i + bi

#

# bi~ Normal(0, τ)

import numpy as np, pymc3 as pm, theano.tensor as T

# germinated seeds

r =  np.array([10, 23, 23, 26, 17, 5, 53, 55, 32, 46, 10, 8, 10, 8, 23, 0, 3, 22, 15, 32, 3])

# total seeds

n =  np.array([39, 62, 81, 51, 39, 6, 74, 72, 51, 79, 13, 16, 30, 28, 45, 4, 12, 41, 30, 51, 7])

# seed type

x1 = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])

# root type

x2 = np.array([0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1])

# number of plates

N =  x1.shape[0]

import math

with pm.Model() as m:

    ### hyperpriors

    sigma = pm.Uniform('sigma',lower=0.,upper=10.)

    tau =1/pow(sigma,2)

 

    ### parameters

    # fixed effects

    alpha_0 = pm.Normal('alpha_0', mu=0.,tau=1e-6)

    alpha_1 = pm.Normal('alpha_1', mu=0.,tau=1e-6)

    alpha_2 = pm.Normal('alpha_2', mu=0.,tau=1e-6)

    alpha_12 = pm.Normal('alpha_12', mu=0.,tau=1e-6)

    # random effect

    b = pm.Normal('b', 0., tau, shape=(N,))

    # expected parameter

    logit_p = (alpha_0 + alpha_1 * x1 + alpha_2 * x2 + alpha_12 * x1 * x2 + b)

    p = T.exp(logit_p) / (1 + T.exp(logit_p))

    ### likelihood

    obs = pm.Binomial('obs', n, p, observed=r)

n = 10000

with m:

    start = pm.find_MAP({'sigma':8., 'alpha_0': 0., 'alpha_1': 0., 'alpha_2': 0., 'alpha_12': 0.})

    step = pm.HamiltonianMC(scaling=start)

    ptrace = pm.sample(n, step, start, progressbar=False)

    start = pm.find_MAP({'sigma':1., 'alpha_0': 0., 'alpha_1': 0., 'alpha_2': 0., 'alpha_12': 0.,'b': np.zeros(N)})

    step = pm.HamiltonianMC(scaling=start)

    ptrace = pm.sample(n, step, start, progressbar=False)

    

burn = 1000

pm.summary(ptrace[burn:])

It is nowhere close to actual params. Any help appreciated