# Seeds: Random effect logistic regression

**URL:** <https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094>\
**Category:** Questions\
**Created:** [March 24, 2022, 9:33pm UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094 "2022-03-24T21:33:16Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Mohan\_Nayak](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mohan_nayak/32/4987_2.png) [@Mohan\_Nayak](https://discourse.pymc.io/u/Mohan_Nayak)\
**Post date:** [March 24, 2022, 9:33pm UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/1 "2022-03-24T21:33:17Z")

</div>

[https://www.multibugs.org/examples/latest/Seeds.html](https://www.multibugs.org/examples/latest/Seeds.html)

I am trying to solve this example using PyMC3 but getting mean, SD 0. For Sampling I have tried HamiltonMC and Metropolis.

Any hint will be appreciated ?

---

<div class="post-metadata">

**Author:** ![cluhmann](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/cluhmann/32/3083_2.png) [@cluhmann](https://discourse.pymc.io/u/cluhmann)\
**Post date:** [March 25, 2022, 12:23am UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/2 "2022-03-25T00:23:42Z")

</div>

Welcome!

I would suggest checking out this [notebook](https://docs.pymc.io/projects/examples/en/latest/generalized_linear_models/GLM-logistic.html) for an example of logistic regression and [this notebook](https://docs.pymc.io/projects/examples/en/latest/generalized_linear_models/GLM-hierarchical.html) or possibly [this one](https://docs.pymc.io/projects/examples/en/latest/case_studies/multilevel_modeling.html?highlight=radon#varying-intercept-model) for an example of how to build in “random effects”. If you have specific questions, feel free to post them.

---

<div class="post-metadata">

**Author:** ![Mohan\_Nayak](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mohan_nayak/32/4987_2.png) [@Mohan\_Nayak](https://discourse.pymc.io/u/Mohan_Nayak)\
**Post date:** [March 25, 2022, 4:02am UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/3 "2022-03-25T04:02:22Z")

</div>

Thanks @cluhmann for your response.

Basically I have tried HamiltonMC and Variational Inference. Below is the code, but I am not getting satisfactory result. Probably I am not using correct hyperparameter ?

I am trying to solve seeds problem [Seeds](https://www.multibugs.org/examples/latest/Seeds.html). In my case, mean and sd are 0 and no where close to the numbers it should have.

[Seeds](https://github.com/pyare007/Data-Science/blob/main/Seeds.ipynb)

---

<div class="post-metadata">

**Author:** ![cluhmann](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/cluhmann/32/3083_2.png) [@cluhmann](https://discourse.pymc.io/u/cluhmann)\
**Post date:** [March 25, 2022, 4:41am UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/4 "2022-03-25T04:41:20Z")

</div>

I think the biggest issue is that the SDs of your priors are very small. Things seem to run ok once setting them to 1.0.

As a side note, you can also use use the built in `invlogit()` if you want to skip doing some of the intermediate calculations by hand:

```python
obs = pm.Binomial('obs',
                  n,
                  pm.math.invlogit(logit_p),
                  observed=r)

```

---

<div class="post-metadata">

**Author:** ![Mohan\_Nayak](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mohan_nayak/32/4987_2.png) [@Mohan\_Nayak](https://discourse.pymc.io/u/Mohan_Nayak)\
**Post date:** [March 25, 2022, 5:01am UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/5 "2022-03-25T05:01:11Z")

</div>

Yes, I tried changing sd before and was not getting 0 but we have consider this while using model.

```
     alpha0 ~ dnorm(0.0,1.0E-6)
     alpha1 ~ dnorm(0.0,1.0E-6)
     alpha2 ~ dnorm(0.0,1.0E-6)
     alpha12 ~ dnorm(0.0,1.0E-6)
     sigma ~ dunif(0,10)
     tau <- 1 / pow(sigma, 2)

```

The example considered here [Seeds](https://www.multibugs.org/examples/latest/Seeds.html) is done in R and tried with 2 initialization values and never got 0. So I am thinking I am doing some mistake or missing some hyperparameter in sample function.

I am not sure if R does it better.

---

<div class="post-metadata">

**Author:** ![cluhmann](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/cluhmann/32/3083_2.png) [@cluhmann](https://discourse.pymc.io/u/cluhmann)\
**Post date:** [March 25, 2022, 5:22am UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/6 "2022-03-25T05:22:44Z")

</div>

> [@Mohan\_Nayak](#):
>
> …was not getting 0…

I’m not sure what you mean by “getting zero”.

I saw that the small (1.0E-6) SDs in the model you referenced. Maybe try increasing the number of tuning steps would help (maybe 5,000 or 10,000?).

---

<div class="post-metadata">

**Author:** ![Mohan\_Nayak](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mohan_nayak/32/4987_2.png) [@Mohan\_Nayak](https://discourse.pymc.io/u/Mohan_Nayak)\
**Post date:** [March 25, 2022, 6:11am UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/7 "2022-03-25T06:11:35Z")

</div>

oh I am sorry for not being specific. I meant, mean and variance in my case was always zero. When you said, try changing variance. That’s what my answer was that by changing sd to 1, I stopped getting zero as mean and variance.

As per you advice, I tried with 5000 & 10000 but still no luck.

Luckily I found a sample python implementation [https://github.com/aflaxman/pymc-examples/blob/master/seeds\_re\_logistic\_regression\_pymc3.ipynb](https://github.com/aflaxman/pymc-examples/blob/master/seeds_re_logistic_regression_pymc3.ipynb) using pymc3. In this, the author used almost same mean and variance as defined in multibugs seeds example and got non-zero mean and variance. But when I ran same code, I am getting mean and variance as 0.

I am doing something wrong somewhere for sure.

---

<div class="post-metadata">

**Author:** ![Mohan\_Nayak](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mohan_nayak/32/4987_2.png) [@Mohan\_Nayak](https://discourse.pymc.io/u/Mohan_Nayak)\
**Post date:** [March 25, 2022, 9:25pm UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/8 "2022-03-25T21:25:12Z")

</div>

@cluhmann - As per your advice, I have changed prior sd to 1 and can see that the mean & variance are approximately matching with Multibugs SEEDS example.  
What reason I should write while modifying standard deviation ? As I am completely using a different model from original setup defined now.  
What I meant is how to find if we should use this distribution and initial values for alphas, beta, Sigma etc. Pls advice.

---

<div class="post-metadata">

**Author:** ![cluhmann](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/cluhmann/32/3083_2.png) [@cluhmann](https://discourse.pymc.io/u/cluhmann)\
**Post date:** [March 25, 2022, 9:31pm UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/9 "2022-03-25T21:31:24Z")

</div>

The results on that page suggest that some of the posterior estimates (e.g., `alpha2`) are very far from zero. So constructing a prior with a mean of 0 and an SD of 1e-6 suggest that that those results are not very credible a priori. Another way to think about it is that those priors suggest that the observed data is (nearly) impossible. Once you set the SD to 1, the data is much more plausible as far as the model is concerned.

---

<div class="post-metadata">

**Author:** ![Mohan\_Nayak](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mohan_nayak/32/4987_2.png) [@Mohan\_Nayak](https://discourse.pymc.io/u/Mohan_Nayak)\
**Post date:** [March 25, 2022, 9:48pm UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/10 "2022-03-25T21:48:04Z")

</div>

That’s absolutely make sense. Thanks a zillion.

I am still in dilema how in R, they are getting these values with small sd.

Sorry, one last question, in multibugs seed example, I see there is one point where values are getting initilized as below, I m planning to try this now. How we do that in pymc3 ?

##### Inits for chain 1

`list(alpha0 = 0, alpha1 = 0, alpha2 = 0, alpha12 = 0, sigma = 8)`

##### Inits for chain 2

`list(alpha0 = 0, alpha1 = 0, alpha2 = 0, alpha12 = 0, sigma = 1,beta = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))`

---

<div class="post-metadata">

**Author:** ![cluhmann](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/cluhmann/32/3083_2.png) [@cluhmann](https://discourse.pymc.io/u/cluhmann)\
**Post date:** [March 25, 2022, 10:09pm UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/11 "2022-03-25T22:09:32Z")

</div>

You should be able to do something like this:

```python
pm.sample(start={'alpha_0':0, 
                 'alpha_1':0, 
                 'alpha_2':0, 
                 'alpha_12':0, 
                 # etc
                }
         )

```

---

<div class="post-metadata">

**Author:** ![ckrapu](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ckrapu/32/1366_2.png) [@ckrapu](https://discourse.pymc.io/u/ckrapu)\
**Post date:** [March 25, 2022, 11:57pm UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/12 "2022-03-25T23:57:06Z")

</div>

My understanding was that the parameters to dnorm are in terms of mean and precision instead of variance/as.

---

<div class="post-metadata">

**Author:** ![Mohan\_Nayak](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mohan_nayak/32/4987_2.png) [@Mohan\_Nayak](https://discourse.pymc.io/u/Mohan_Nayak)\
**Post date:** [March 26, 2022, 12:35am UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/13 "2022-03-26T00:35:20Z")

</div>

@ckrapu - wow, thanks for the pointer. I just went through the MULTIBUG document and found that it is not sd but tau. you are absolutely correct and thanks for jumping and being a savior.

In the BUGS language it is used as  
x ~ dnorm(mu, tau)

@cluhmann - I am sorry, I got confused and thought this is R code where as it is Multibug specific code. I came to know after running the example in MULTIBUG IDE([Tutorial](https://www.multibugs.org/documentation/latest/Tutorial.html))

As soon as I changed my code to below in pymc3, I can see value approximately matching to the result in MULTIBUGS.  
alpha\_0 = pm.Normal(‘alpha\_0’, mu=0., tau=1e-6)

Let me try few things and might seek some help.

Thanks a zillion again @cluhmann and @ckrapu

---

<div class="post-metadata">

**Author:** ![Parikshit\_Mukherjee](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/parikshit_mukherjee/32/4937_2.png) [@Parikshit\_Mukherjee](https://discourse.pymc.io/u/Parikshit_Mukherjee)\
**Post date:** [March 26, 2022, 11:35am UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/14 "2022-03-26T11:35:58Z")

</div>

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

```python
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

---

<div class="post-metadata">

**Author:** ![ckrapu](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ckrapu/32/1366_2.png) [@ckrapu](https://discourse.pymc.io/u/ckrapu)\
**Post date:** [March 26, 2022, 2:38pm UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/15 "2022-03-26T14:38:07Z")

</div>

You shouldn’t be using the `scaling` keyword argument with the values you have. I’m surprised it lets you sample at all because that argument sets the diagonal terms of a preconditioning matrix. If you set some of those values to zero, you are forcing the sampler to make moves of size zero, going nowhere. I recommend getting rid that and using NUTS with no additional arguments.

---

<div class="post-metadata">

**Author:** ![cluhmann](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/cluhmann/32/3083_2.png) [@cluhmann](https://discourse.pymc.io/u/cluhmann)\
**Post date:** [March 26, 2022, 2:38pm UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/16 "2022-03-26T14:38:29Z")

</div>

> [@Parikshit\_Mukherjee](#):
>
> `b = pm.Normal('b', 0., tau, shape=(N,))`

This should probably be `tau=tau`. Right now I think this is using `tau` as the SD.

---

<div class="post-metadata">

**Author:** ![Parikshit\_Mukherjee](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/parikshit_mukherjee/32/4937_2.png) [@Parikshit\_Mukherjee](https://discourse.pymc.io/u/Parikshit_Mukherjee)\
**Post date:** [March 26, 2022, 6:31pm UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/17 "2022-03-26T18:31:40Z")

</div>

Hi @ckrapu@cluhmann  
Thanks a lot for the help. It works.  
One last query I wanted to use this model to predict ,I am following this

[https://docs.pymc.io/en/v3/pymc-examples/examples/diagnostics\_and\_criticism/posterior\_predictive.html](https://docs.pymc.io/en/v3/pymc-examples/examples/diagnostics_and_criticism/posterior_predictive.html)

but not able to identify clearly how to implement in my case

we are looking for a simple example

---

<div class="post-metadata">

**Author:** ![cluhmann](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/cluhmann/32/3083_2.png) [@cluhmann](https://discourse.pymc.io/u/cluhmann)\
**Post date:** [March 26, 2022, 6:33pm UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/18 "2022-03-26T18:33:48Z")

</div>

> [@Parikshit\_Mukherjee](#):
>
> I wanted to use this model to predict

Can you be more specific?

---

<div class="post-metadata">

**Author:** ![Parikshit\_Mukherjee](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/parikshit_mukherjee/32/4937_2.png) [@Parikshit\_Mukherjee](https://discourse.pymc.io/u/Parikshit_Mukherjee)\
**Post date:** [March 26, 2022, 6:38pm UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/19 "2022-03-26T18:38:10Z")

</div>

Hi @cluhmann

Basically I am referring to the **[Seeds: Random effect logistic regression]** model that I implemented using pymc3 I want to use it for prediction now. Basically I want to use the learnt posterior distribution to make prediction on either held out test data points or training data points and find the prediction error (if any)

---

<div class="post-metadata">

**Author:** ![cluhmann](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/cluhmann/32/3083_2.png) [@cluhmann](https://discourse.pymc.io/u/cluhmann)\
**Post date:** [March 26, 2022, 6:46pm UTC](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094/20 "2022-03-26T18:46:18Z")

</div>

For “predicting” the training data, you can use `sample_prior_predictive()` as illustrated in [this notebook](https://docs.pymc.io/en/latest/learn/examples/posterior_predictive.html). For predicting new data, you can check out the answer [in this thread](https://discourse.pymc.io/t/how-do-i-predict-on-new-unseen-real-data-using-pm-sample-posterior-predictive/6467/2) for how to use `sample_prior_predictive()` and `set_data()` to do so.

[Next page](https://discourse.pymc.io/t/seeds-random-effect-logistic-regression/9094.md?page=2)