# How to specify this model directly in Pymc3?

**URL:** <https://discourse.pymc.io/t/how-to-specify-this-model-directly-in-pymc3/4517>\
**Category:** Questions\
**Created:** [February 17, 2020, 12:47pm UTC](https://discourse.pymc.io/t/how-to-specify-this-model-directly-in-pymc3/4517 "2020-02-17T12:47:34Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![tankenot](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/tankenot/32/2491_2.png) [@tankenot](https://discourse.pymc.io/u/tankenot)\
**Post date:** [February 17, 2020, 12:47pm UTC](https://discourse.pymc.io/t/how-to-specify-this-model-directly-in-pymc3/4517/1 "2020-02-17T12:47:34Z")

</div>

Hello,

I would like to model some repeated-measures data with Pymc3. The data comes from an experiment where n=17 subjects experienced three conditions [baseline, treatment1, treatment2] and a reaction time was measured for each condition. The final data set has 51 reaction times (3 conditions \* 17 subjects).

1. I have managed to fit the following varying-intercepts model with the bambi-package:

> b\_model = bambi.Model(data)  
> results = b\_model.fit(‘reaction\_time ~ condition’, random = [‘1|subject\_id’], tune = 2000,  
> samples = 3000, target\_accept = 0.9)

1. I would like to create the same model (or a replica of the bambi-model) directly with Pymc3, but cannot get it to work. As I’ve understood it a hierarchical model with varying intercept can be specified in Pymc3 using hyperpriors and priors (and in my case only for the intercept). The model I’ve specified so far is:

> with pm.Model() as bambi\_replica:  
> # Hyperpriors  
> mu\_intercept = pm.Normal(‘mu\_intercept’, mu=0, sd=1)  
> sd\_intercept = pm.HalfNormal(‘sd\_intercept’, sd=1)  
> # Priors  
> Intercept = pm.Normal(‘Intercept’, mu = mu\_intercept, sd=sd\_intercept, shape=n\_subjects)  
> reaction\_time\_sd = pm.HalfNormal(‘reaction\_time\_sd’, sd = data.reaction\_time.std()\*2)  
> condition = pm.Normal(‘condition’, mu = 0, sd = 1, shape=2)  
> # mu  
> mu = Intercept[tp\_idx] + condition[0]\*x\_1 + condition[1]\*x\_2  
> # Likelihood (sampling distribution) of observations  
> reaction\_time = pm.Normal(‘Y\_obs’, mu=mu, sd = reaction\_time\_sd, observed=reaction\_time\_data)

Here:

- n\_subjects = 17 (since I have 17 subjects)
- x\_1 = array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1  
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,  
0, 0, 0, 0, 0, 0, 0], dtype=int64) (to identify indexes for treatment 1)
- x\_2 = array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,  
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  
1, 1, 1, 1, 1, 1, 1], dtype=int64) (to identify indexes for treatment 2)
- reaction\_time\_data = array([2.918456, 3.124924, 3.013115, 3.117612, 2.91747 , 3.157323,  
3.101546, 3.595617, 2.961656, 2.960227, 2.537308, 2.908257,  
2.971608, 3.074054, 2.870034, 2.867201, 3.149056, 3.186297,  
3.065235, 4.028026, 3.068959, 2.902312, 3.276905, 3.398952,  
3.286361, 3.105795, 3.748059, 2.893177, 4.078607, 3.412715,  
3.849621, 3.440311, 3.416428, 3.169723, 3.256486, 3.014277,  
2.745785, 4.160992, 3.226214, 4.321797, 3.05727 , 3.231891,  
3.611861, 3.229321, 2.851948, 2.622444, 2.844787, 3.229938,  
3.055713, 2.722098, 3.158651])
- tp\_idx = array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,  
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,  
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16],  
dtype=int64)

**Questions:**

1. Would this approach correspond to the varying-intercept model specified with bambi (i.e., ‘reaction\_time ~ condition’, random = [‘1|subject\_id’])? I understand that the priors may not be specified exactly the same, but I’m mainly interested in if the overall structure of the model corresponds to a varying-intercepts model that could be used to model the data I have.
2. With bambi it seems to work fine to fit the model (i.e., the model seems to converge, traceplots look good). The model I specified directly with Pymc3, however, does not seem to converge and I get error-messages looking like below. I can’t understand why this is. Any clues?

> There were 104 divergences after tuning. Increase `target_accept` or reparameterize.  
> The acceptance probability does not match the target. It is 0.7975575845495201, but should be close to 0.95. Try to increase the number of tuning steps.

The estimated number of effective samples is smaller than 200 for some parameters.

All suggestions are appreciated!

---

<div class="post-metadata">

**Author:** ![AlexAndorra](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/alexandorra/32/9142_2.png) [@AlexAndorra](https://discourse.pymc.io/u/AlexAndorra)\
**Post date:** [February 18, 2020, 9:38am UTC](https://discourse.pymc.io/t/how-to-specify-this-model-directly-in-pymc3/4517/2 "2020-02-18T09:38:10Z")

</div>

Hi Linda,  
At first glance I’d say your model seems right for a varying-intercepts model, and I’d expect it to run 🤔

- Have you tried standardizing your `reaction_time_data`? I’m not sure it will have an impact, as your data don’t seem to be on a huge scale, but you never know…

- Which parameters have a small effective sample size? And which don’t look good in general? This could come from `sd_intercept` being very close to zero (i.e subjects being very similar to each other), in which case a non-centered parametrization could do the trick.

In any case, if you haven’t seen it already, [here is the port to PyMC](https://github.com/pymc-devs/resources/blob/master/Rethinking/end-of-chapter-practice-problems/ch-13.ipynb) of end-of-chapter exercices for the chapter on varying-intercepts of McElreath’s book – it’s really useful!  
Hope this helps 🖖

---

<div class="post-metadata">

**Author:** ![tankenot](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/tankenot/32/2491_2.png) [@tankenot](https://discourse.pymc.io/u/tankenot)\
**Post date:** [February 24, 2020, 3:20pm UTC](https://discourse.pymc.io/t/how-to-specify-this-model-directly-in-pymc3/4517/3 "2020-02-24T15:20:43Z")

</div>

Thanks @AlexAndorra 🙂

- I standardized the reaction time

> standardized\_rt = (data.reaction\_time - reaction\_time\_data.mean())/reaction\_time\_data.std()).

That only did not seem to make a difference.

- I think I tried a non-centered parametrization:

> mu\_intercept = pm.Normal(‘mu\_intercept’, mu=reaction\_time\_data.mean(), sd = reaction\_time\_data.std()\*2)).

That did also not make a difference.

- I’ve noticed that the bambi-model uses a uniform prior on the reaction\_time\_sd. A combination of: 1) standardizing the reaction time and changing the prior from half-normal to a uniform for the reaction\_time\_sd with limits: lower = 0, upper = 0.4 gives me a model that converges. However, the traceplot for the standard-deviation looks cut (result of the cut in the uniform prior I assume):

 ![trace_sd](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/f/f2416f23d7826daa84e3ac8a91bd49a046ceec61.png)

It seems that the model-convergence is very dependent on this prior choice. Any clues why?

- Thanks for the suggestion on further readings about varying-intercepts. I had not seen that before. Will check it out! Your suggestions have been very helpful indeed!

---

<div class="post-metadata">

**Author:** ![nkaimcaudle](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/nkaimcaudle/32/2142_2.png) [@nkaimcaudle](https://discourse.pymc.io/u/nkaimcaudle)\
**Post date:** [February 25, 2020, 1:36am UTC](https://discourse.pymc.io/t/how-to-specify-this-model-directly-in-pymc3/4517/4 "2020-02-25T01:36:24Z")

</div>

You only have three data points per a subject so the choice of priors is important. I got a model to run using a non-centered parameterisation and increasing the target\_accept to 0.95. By non-centered Alex means that you sample the subjects as N(0,1) and then add the population mean and scale by standard deviation afterwards.

```
with pm.Model() as bambi_replica:
    # Hyperpriors
    mu_intercept = pm.Normal('mu_intercept', mu=0, sd=1)
    sd_intercept = pm.Exponential('sd_intercept', 1.)
    
    # Priors
    #Intercept = pm.Normal('Intercept', mu=mu_intercept, sd=sd_intercept, shape=n_subjects)
    Intercept = pm.Normal('Intercept', mu=0, sd=1, shape=n_subjects)
    reaction_time_sd = pm.Exponential('reaction_time_sd', reaction_time_data.std()*2)
    condition = pm.Normal('condition', mu=0, sd=1, shape=2)
    
    # mu
    #mu = Intercept[tp_idx] + condition[0]*x_1 + condition[1]*x_2
    mu = (mu_intercept + sd_intercept*Intercept[tp_idx]) + condition[0]*x_1 + condition[1]*x_2
    
    # Likelihood (sampling distribution) of observations
    reaction_time = pm.Normal('Y_obs', mu=mu, sd=reaction_time_sd, observed=reaction_time_data)
    
    trace = pm.sample(nuts_kwargs={'target_accept':0.95})
pm.summary(trace)

```

 ![2020-02-25_09-34](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/f/f1dc4e25558be3ac9585086312104671c4a29081.png)

In general flat or uniform priors on scale (standard deviation) parameters is not a good idea, you want to have more regularisation towards zero. Your trace of reaction\_time\_sd is clearly pushing against the upper limit you’ve imposed of 0.4

---

<div class="post-metadata">

**Author:** ![AlexAndorra](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/alexandorra/32/9142_2.png) [@AlexAndorra](https://discourse.pymc.io/u/AlexAndorra)\
**Post date:** [February 25, 2020, 12:13pm UTC](https://discourse.pymc.io/t/how-to-specify-this-model-directly-in-pymc3/4517/5 "2020-02-25T12:13:50Z")

</div>

Thanks Nicholas – couldn’t have said it better!  
PS: note that in the newest PyMC version, you can just add this arg to `pm.sample`: `target_accept=0.95` – more readable than the `nuts_kwargs` dict 😉

@tankenot, as you’ll see in Nicholas’s trace, a lot of the mass of `sd_intercept` is very near zero, as I suspected above. This is typically the case when a non-centered parametrization will work a lot better.  
Hope this helps 🖖

---

<div class="post-metadata">

**Author:** ![tankenot](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/tankenot/32/2491_2.png) [@tankenot](https://discourse.pymc.io/u/tankenot)\
**Post date:** [March 3, 2020, 8:12am UTC](https://discourse.pymc.io/t/how-to-specify-this-model-directly-in-pymc3/4517/6 "2020-03-03T08:12:15Z")

</div>

Thanks a lot to both of you! @AlexAndorra, @nkaimcaudle  
It’s been very helpful! 🤩