# Beta Regression in pyMC3

**URL:** <https://discourse.pymc.io/t/beta-regression-in-pymc3/301>\
**Category:** Questions\
**Created:** [August 29, 2017, 11:55pm UTC](https://discourse.pymc.io/t/beta-regression-in-pymc3/301 "2017-08-29T23:55:49Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![bdyetton](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bdyetton/32/208_2.png) [@bdyetton](https://discourse.pymc.io/u/bdyetton)\
**Post date:** [August 29, 2017, 11:55pm UTC](https://discourse.pymc.io/t/beta-regression-in-pymc3/301/1 "2017-08-29T23:55:49Z")

</div>

I am trying to fit a regression model to estimate a percentage, `'sleepEf'` (scaled to 0:1) from a single predictor variable `'tst'` in pyMC3. Given that a percentage is bounded, I believe I should use a beta distributed outcome variable with a logit transform. Here is my model in pyMC3:

```
with pm.Model() as pooled_model:
    b_intr = pm.Normal('b_intr', mu=0.8, sd=100 ** 2)
    b_tst = pm.Normal('b_tst', mu=0.0, sd=100 ** 2)
    
    model_err = pm.HalfNormal('model_err', sd=50)

    # Expected value
    y_est = pm.math.invlogit(b_intr + \
            b_tst * data['tst'])

    # Data likelihood
    y_like = pm.Beta('y_like', mu=y_est, sd=model_err, observed=data['sleepEf'])

```

I have made sure that the `sleepEf` variable is never exactly 0 or 1. Here is a histogram of it:  
[![Histogram for sleepEf](https://yyz2.discourse-cdn.com/flex036/uploads/pymc3/original/1X/28b585bf4e9134f043b225c89246ebd95519f23d.png)](https://discourse.pymc.io/uploads/pymc3/original/1X/28b585bf4e9134f043b225c89246ebd95519f23d.png)

When I try to fit the MAP estimate and then sample I get errors:

```
with pooled_model:
    start = pm.find_MAP(fmin=sp.optimize.fmin_powell)
    hierarchical_trace = pm.sample(2000, step=pm.Metropolis(), start=start, tune=1000)

```

> ValueError: Optimization error: max, logp or dlogp at max have  
> non-finite values. Some values may be outside of distribution support.  
> max: {‘b\_intr’: array(3.3879289615458417), ‘b\_tst’:  
> array(2.587928961545842)} logp: array(-inf) dlogp: array([  
> -2.58792896e-08, -2.58792896e-08])Check that 1) you don’t have hierarchical parameters, these will lead to points with infinite  
> density. 2) your distribution logp’s are properly specified. Specific  
> issues:

Have I incorrectly specified the model?

Also note that if I do not find\_MAP before sampling, my parameter values never change from their prior values.

Thanks!

I have also asked this question [here](https://stats.stackexchange.com/questions/300426/beta-regression-model-to-estimate-a-percentage-in-pymc3):

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [August 30, 2017, 6:31am UTC](https://discourse.pymc.io/t/beta-regression-in-pymc3/301/2 "2017-08-30T06:31:33Z")

</div>

The problem is that the model\_err is too large - after internal transform the alpha and beta in the `pm.Beta` become negative, which is outside of the range of support.

You can check the value of `mu` and `sd` by doing:

```auto
y_est.tag.test_value
model_err.tag.test_value

```

And transform the `mu` and `sd` to `alpha` and `beta` by doing:

```auto
kappa = mu * (1 - mu) / sd**2 - 1
alpha = mu * kappa
beta = (1 - mu) * kappa

```

This is what pymc3 does internally.

---

<div class="post-metadata">

**Author:** ![aseyboldt](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/aseyboldt/32/5795_2.png) [@aseyboldt](https://discourse.pymc.io/u/aseyboldt)\
**Post date:** [August 30, 2017, 7:51am UTC](https://discourse.pymc.io/t/beta-regression-in-pymc3/301/3 "2017-08-30T07:51:56Z")

</div>

Parameterizing the beta distribution via mu and sd is always a bit of a mess. It might work much better if you logit transform your dataset instead, and then use a normal likelihood.

---

<div class="post-metadata">

**Author:** ![bdyetton](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bdyetton/32/208_2.png) [@bdyetton](https://discourse.pymc.io/u/bdyetton)\
**Post date:** [August 30, 2017, 5:45pm UTC](https://discourse.pymc.io/t/beta-regression-in-pymc3/301/4 "2017-08-30T17:45:02Z")

</div>

Thanks guys! i got it to work with a smaller sd.

@junpenglao, maybe there should be something in the documentation on the limits of sd, (currently it only states `sd > 0`). Is this the kind of thing i would open a github issue for?

@aseyboldt, i think im going to stick with the beta regression because the units of my data are readily interpretable and its working now. However, im curious: i believe log transforming the data would solve the bounds issue on the mean, but would the model error be accurate when _Normal_(mu\>0.8 or mu\<0.2, sd=model\_err), i.e. when mu gets close to the bounds?

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [August 30, 2017, 8:33pm UTC](https://discourse.pymc.io/t/beta-regression-in-pymc3/301/5 "2017-08-30T20:33:59Z")

</div>

Yeah: a beta distribution with mean mu, the var should be smaller than mu(1 - mu). A PR to edit the docstring and also a better check of this would be great.

---

<div class="post-metadata">

**Author:** ![bdyetton](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bdyetton/32/208_2.png) [@bdyetton](https://discourse.pymc.io/u/bdyetton)\
**Post date:** [August 30, 2017, 9:59pm UTC](https://discourse.pymc.io/t/beta-regression-in-pymc3/301/6 "2017-08-30T21:59:22Z")

</div>

PR is here: [https://github.com/pymc-devs/pymc3/pull/2534](https://github.com/pymc-devs/pymc3/pull/2534)  
I just edited the docstring.

---

<div class="post-metadata">

**Author:** ![aseyboldt](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/aseyboldt/32/5795_2.png) [@aseyboldt](https://discourse.pymc.io/u/aseyboldt)\
**Post date:** [August 31, 2017, 3:46pm UTC](https://discourse.pymc.io/t/beta-regression-in-pymc3/301/7 "2017-08-31T15:46:08Z")

</div>

@bdyetton Thanks for the PR!

Just for clarification: I didn’t mean a log transform, but a logit transform.  
I’m not sure what “accurate” means in this context, a logit-normal model is a bit different from the beta distribution, but unless you have a good reason why the beta would somehow be the right one, I don’t see a clear advantage to using it. They are in fact quite similar, a normal distribution and the logit of a beta aren’t _that_ different, the beta has a bit longer tails. But the interpretation of mean and std in the logit-normal case are I think a bit more natural.

---

<div class="post-metadata">

**Author:** ![ericmjl](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ericmjl/32/196_2.png) [@ericmjl](https://discourse.pymc.io/u/ericmjl)\
**Post date:** [September 18, 2017, 3:21pm UTC](https://discourse.pymc.io/t/beta-regression-in-pymc3/301/8 "2017-09-18T15:21:48Z")

</div>

Just chipping in with a “hooray!!!” for a PR to PyMC3, @bdyetton!