Beta regression,Rejecting initial value,Error evaluating the log probability at the initial value

Modeling
1

Hello,I’m trying to fit the following beta regression in Stan :
Y_{i}\sim beta(p_{i},q_{i})
p_{i}=((1-\mu_{i})\mu_{i}^{2}-\mu_{i}\sigma^{2})/\sigma^{2}
q_{i}=(1-\mu_{i})(\mu_{i}-\mu_{i}^{2}-\sigma^{2})/\sigma^{2}
\mu_{i}=w_{1}x_{i,1}+...w_{k}x_{i,k}+w_{k+1}z_{i}

and make inference for weights w.

Where w\sim Dirichlet (1,...1)
, \sigma^{2}\sim Uniform(0,min((1-\mu)\mu))
and z\sim Beta(0.5,0.5)

So in order to do that I implement the following code:

data{
int<lower=0> N;
int<lower=0> K;
vector[K] C;
vector<lower=0,upper=1>[N] Y;
matrix<lower=0,upper=1>[N,K] X;
}

parameters{
vector<lower=0.0,upper=1.0>[K+1] W;
}
transformed parameters {
vector<lower=0, upper=1>[N] mu;
vector<lower=0>[N] p;
vector<lower=0>[N] q;
real<lower=0,upper=1> Z;
real<lower=0> sig;
for (i in 1:N) {
mu[i] = X[i,]*W[1:K]+W[K+1]*Z;
p[i] = ((1-mu[i])*mu[i]*mu[i]-mu[i]*sig*sig)/(sig*sig);
q[i] = (1-mu[i])*(mu[i]-mu[i]*mu[i]-sig*sig)/(sig*sig);
}
}
model{
target += beta_lpdf(Z|0.5,0.5);
target += uniform_lpdf(sig|0,min((1-mu).*mu));
target += dirichlet_lpdf(W|C);
target += beta_lpdf(Y|p,q);
}

However, I get the following message multiple times

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: model31f066875_0b43718f5170bcc3da95f3b8a5114180_namespace::write_array: mu[1] is nan, but must be greater than or equal to 0  (in 'model31f066875_0b43718f5170bcc3da95f3b8a5114180' at line 14)

I’ve read that in order to solve these problems Rejecting initial value,Error evaluating the log probability at the initial value. you have to change the constraints of the variables.I made some attempts but I still get the same message.

The way that the model is constructed it seems straightforward , thus I don’t know what I’m doing wrong.

2

Since you want to infer sigma and Z, those need to be parameters. Transformed parameters are not inferred.

Also coding that upper bound on sigma is going to be nasty. Are you sure you need it? Is there another way to do the prior maybe? I don’t think it’s differentiable either, so it would be problematic either way.

Also, welcome to the forums!

3

Maybe you have already done this, but it can’t hurt to ask.

Did you check that Y does not contain any zeros or ones?
If Y does contain zeros or ones, this could explain your problem, because beta_lpdf(Y|p,q) is not finite if p \neq 1 and Y is 1 or 0. In this case you would need to run a zero and/or one inflated beta regression.

4

In my notes I have that the prior that I used for \sigma guarantees that the parameters p and q will be non-negative.But in Stan we delcare the domain for each parameter, so I suppose that this kind of prior is not needed.However, there is a second prior for \sigma which I’m going to use ,\sigma =1/tau where tau\sim gamma(0.1,0.1) and check if the problem is fixed.

5

Thank you for mentioning this , because I haven’t checked it.Luckily all the values of Y belong in (0,1).

6

I think the key thing @bbbales2 is saying is that you need move the declaration of Z from the transformed parameters block to the parameters block.

7

@bbbales2 @Guido_Biele I did it and still doesn’t work.However,I think the problem is the way that I define the weights W.Because , I get that that mu_{i} is nan or that it is greater than 1.Also,we know that mu[i] = X[i,]*W[1:K]+W[K+1]*Z[i]; ,consequently the only parameter that can “harm” the validity of the calculation is W.With all these in mind I changed my code to the following(I also changed the prior for sigma):

data{
int<lower=0> N;
int<lower=0> K;
vector[K+1] C;
vector<lower=0,upper=1>[N] Y;
matrix<lower=0,upper=1>[N,K] X;
}

parameters{
simplex[K+1] W;
real<lower=0> tau;
vector<lower=0.0,upper=1.0>[N] Z;
 }

transformed parameters {
vector<lower=0.0,upper=1.0>[N] mu;
vector<lower=0.0>[N] p;
vector<lower=0.0>[N] q;
real<lower=0.0> sig;
sig = 1/tau;
for (i in 1:N) {
mu[i] = X[i,]*W[1:K]+W[K+1]*Z[i];
p[i] = ((1-mu[i])*mu[i]*mu[i]-mu[i]*sig*sig)/(sig*sig);
q[i] = (1-mu[i])*(mu[i]-mu[i]*mu[i]-sig*sig)/(sig*sig);
}
}

model{

W ~ dirichlet(C);
Z ~ beta(0.5,0.5);
tau ~ gamma(0.1,0.1);
target += beta_lpdf(Y|p,q);
}

'

which still rejects a lot of initial values, but this time not because mu behaves bad but because p is less than 0 some times.

8

It might be easier to just switch do a different parameterization.

Your outcome is a proportion of some sort, and you have covariates to describe that. Here’s a way of fitting a model like that with rstanarm (and you could write an equivalent Stan model too): https://cran.r-project.org/web/packages/rstanarm/vignettes/betareg.html