Beta Regression - Needs more restrictions but why, where, and how?

[Ilan Strauss]

Hi,

I am modelling a proportion (a real number between 0 and 1) - my y[N] - using, at the moment, a single regressor (my X) - which is also a real number between 0 and 200.  Im trying to model this using a Bayesian beta regression in Stan. Eventually we will want to make it a more complex model (hierarcichal and quadratic?) But this simple version keeps giving me grief. It doesn't run unless I put constraints on everything (which requires moving things to be parameters even when in theory they shouldn't really). I can't see why this model below would need more restrictions. And if so where and how? Many thanks in advance for any assistance.

Current error message:

R CODE AND ERROR OUTPUT
---

N <- nrow(ds)               \\we have 100 data points (N=100)

y <- as.numeric(ds$usage_p)

Xvar <- as.matrix(ds$distance)

x1 = as.numeric(Xvar)

K <- ncol(x1)

data2 <- list(N=N,x1=x1,y=y)

model_2 <- stan(file = "beta6.stan", data = data2, iter = 2000, chains = 5)

Exception thrown at line 28: beta_log: Second shape parameter is 0, but must be > 0!    38

  Error evaluating the log probability at the initial value.                            32

Rejecting initial value:                                                                32

---

STAN CODE:
---

data {

  int<lower=1> N;                   // sample size

  real<lower=0,upper=1> y[N];     // response 

  vector[N] x1;                    // predictor matrix

}

parameters {

  real<lower=0.1> phi; 

  real<lower=0> a;

  real<lower=0> b;

    }

model {

   real beta_a;

   real beta_b;

   real mu; 

   

    for (i in 1:N) {

        mu = inv_logit(a + b*x1[i]);

        beta_a = mu * phi;

        beta_b = (1-mu) * phi;

       y[i] ~ beta(beta_a, beta_b);

    }

   phi ~ cauchy(0, 2.5); // half-cauchy priors

    a ~ cauchy(0, 2.5);

    b ~ cauchy(0, 2.5);

}

---

I know it would be more efficient to increment the log posterior myself. I've tried to look for posts on similar(ish) issues. Found a few but nothing which helped me sufficiently on this specific issue.

Many thanks!

Ilan Strauss

[Ben Goodrich]

On Monday, October 10, 2016 at 10:57:33 PM UTC-4, Ilan Strauss wrote:

I can't see why this model below would need more restrictions. And if so where and how? Many thanks in advance for any assistance.

Possibly just need to squeeze the initial values by specifying the init_r argument to stan() to be something less than 2. Or maybe your predictor has crazy units?

Ben

[Bob Carpenter]

Why are you constraining a and b to be positive here? They're
inside the inverse logit transform, so it'd be more typical to let
them be unconstrained. The default random inits for positive constrained
values will be exp(uniform(-2, 2)) or roughly in (1/9, 9).
If these all pile up to produce large mu values, then it's possible
beta_b is just being driven to 0. Reducing the range may help, but
if you choose 1, the inits will be in (1/3, 3) or thereabouts, which may
still be a problem.

You could always try custom inits, which shouldn't be too hard
for this problem.

And why lower=0.1 on phi? We only recommend constraints where
physically/logically necessary. Otherwise, they can cause probability
mass to pile up at the boundary, which will cause computational
issues because of the logit transform to unconstrained space taking
these to large absolute negative values. That may not be a problem
here, though.

I also embedded some comments on efficiency of the program.

Not sure what you mean by that. The sampling statements are at least
as efficient as target += ... It'd be a bit more efficient to collect
the beta_a and beta_b up as vectors beta_as and beta_bs
and write y ~ beta(beta_as, beta_bs); But given there are no shared
computations, it might not help much.

> I've tried to look for posts on similar(ish) issues. Found a few but nothing which helped me sufficiently on this specific issue.
>
> Many thanks!
>
> Ilan Strauss
>

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Ilan Strauss]

Hi Ben and Bob,

Many thanks for the helpful replies. It is appreciated.

I removed the restrictions on my regression coefficients. Does help much. Phi has a restriction since for a beta regression it has to be > 0.     I tried removing the restriction anyways but doesn't help. I run Stan() with different initial values of less than 2  (0.5, 1, 1.5, 2) doesn't help that much.

When I print my output I get the following - most notably, the regression parameters are being pushed to zero.

OUTPUT

The following numerical problems occured the indicated number of times after warmup on chain 4

                                                                                         count

Exception thrown at line 20: beta_log: Second shape parameter is 0, but must be > 0!        17

Exception thrown at line 20: beta_log: First shape parameter is inf, but must be finite!     1

      mean se_mean  sd  2.5%   25%   50%   75% 97.5% n_eff Rhat

phi   21.8     0.1 2.9  16.3  19.7  21.8  23.7  27.6  2167    1

a      0.0     0.0 0.1  -0.2  -0.1   0.0   0.0   0.1  1791    1

b      0.0     0.0 0.0   0.0   0.0   0.0   0.0   0.0  2297    1

lp__ 135.7     0.0 1.1 132.7 135.1 135.9 136.5 136.9  2177    1

I wasn't quite sure about your one comment Bob and if it was different to Ben's suggestion: "You could always try custom inits, which shouldn't be too hard for this problem. "

Way forward?

Thanks again,

Ilan

On Monday, October 10, 2016 at 10:57:33 PM UTC-4, Ilan Strauss wrote:

[Bob Carpenter]

I know it has to be > 0, but didn't understand why it was
originally constrained to > 0.1.

Especially given that output I'm guessing your problem
is constraining the regression coefficients > 0. I'm guessing
if you remove that constraint, it'll work. When you add
a constraint and the value wants to be outside the constraint,
the draws pile up at the boundary.

Our recommendation is *always* to fit simulated data first to
see if the model works as expected. If the model fits simulated
data but nor real data, it typically indicates model misspecification.

I'm curious as to why you'd want to use beta regression rather
than a linear regression on logit(y) vs. x.

I was curious about beta regression (Ben and crew are adding it
to RStanArm, I believe), so I went and did the simulation myself
and realized that this model is very finicky.

data {
int<lower=1> N;
real<lower=0,upper=1> y[N];
vector[N] x;
}
parameters {
real<lower=0> phi;
real alpha;
real beta;
}
model {
for (n in 1:N) {
real mu;
mu = inv_logit(alpha + beta * x[n]);
y[n] ~ beta(mu * phi, (1 - mu) * phi);
}
}

I just used flat uniform priors here; if you do use a Cauchy and
expect high phi values, you should adjust the scale accordingly.

Here's the R simulation and fit code; I added some control to avoid
divergences.

inv_logit <- function (u) 1 / (1 + exp(-u));
logit <- function(v) log(v / (1 - v));

N <- 500;
alpha <- -0.5;
beta <- 1.9;
phi <- 57.4;
x <- rnorm(N);
y <- rep(NA, N);
for (n in 1:N) {
mu <- inv_logit(alpha + beta * x[n]);
y[n] <- rbeta(1, mu * phi, (1 - mu) * phi);
}

library(rstan);
fit <- stan("beta-reg.stan", data=c("N", "y", "x"),
control=list(stepsize=0.01, adapt_delta=0.95));


You can see that it fits just fine:

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

phi 57.12 0.06 3.45 50.63 54.83 57.01 59.30 64.02 2864 1
alpha -0.50 0.00 0.01 -0.53 -0.51 -0.50 -0.49 -0.47 2733 1
beta 1.89 0.00 0.02 1.85 1.88 1.89 1.90 1.93 2692 1
lp__ 902.42 0.03 1.22 899.11 901.92 902.74 903.29 903.73 2093 1

Recall that the simulated params were:

alpha <- -0.5;
beta <- 1.9;
phi <- 57.4;

Hope that helps.

- Bob

[Ilan Strauss]

After ensuring that we have a mu[i] - instead of mu - and adjusting the code accordingly seems to be working now. Also changing things to real away from vector etc.

On Monday, October 10, 2016 at 10:57:33 PM UTC-4, Ilan Strauss wrote:

[Ilan Strauss]

Thanks so much! Just saw your response now. Very helpful.

Best,

Ilan

On Monday, October 10, 2016 at 10:57:33 PM UTC-4, Ilan Strauss wrote: