Rejecting initial value: Chain 1: Log probability evaluates to log(0), i.e. negative infinity

Hi all,

I am running over and over into the same issue:

Even if I add restrictions to my parameters, change some prior distributions and set a list of initial values when I run the model in R as follows, I get the same issue:
nobs = nrow(data14NOMISSING)
t0 = 0.0
zAB0 ← data14NOMISSING$AB0
ts ← data14NOMISSING$Time
indivs = length(unique(data14NOMISSING$ID))
indivs
zASC0 ← rnorm(nobs, 654, 0.5)
antibodies = data14NOMISSING$AB
subj = data14NOMISSING$ID
zinitials = cbind(zAB0, zASC0)
datalist2 ← list(nobs = nobs, t0 = t0, y0 = zinitials,
ts = ts, indivs = indivs, antib = antibodies, subj = subj)
init = function(){
list(pAB0= -2, uAB = 0.053, uASC0 = -2, AB0 = 1766, ASC0 = 654, sigma = 1,
sigmapAB = 0.05, sigmauASC = 0.05, sigmaAB0 = 0.05,
sigmaASC0 = 0.05, rpAB = runif(nobs, 1, 1.5), ruASC = runif(nobs, 1, 1.5),
rAB0 = runif(nobs, 1, 1.5), rASC0 = runif(nobs, 1, 1.5))
}
Modelnosorted ← stan_model(“C:/Users/IGarciaFogeda/Documents/R/Rstan/HIERARCHICAL/Longitudinal/longitudinalmodel.stan”)
fit_modelnosorted<- sampling(Modelnosorted,
data = datalist2, init = init,
iter = 2000,
chains = 1,
seed = 0, control = list (stepsize = 0.1))

Where does the 0 of the log probability come from?

functions {
  vector model3(real t, vector y, real pAB, real uAB, real uASC) {
    vector[2] dydt;  
    dydt[1] = pAB*y[2]-uAB*y[1];
    dydt[2] = -uASC*y[2];
    
   return dydt;
  }
}
data {
  int <lower=1> nobs;
  real t0;
  vector[2] y0[nobs];
  real ts[nobs];
  int <lower=1> indivs;
  real <lower=0> antib[nobs];
  int <lower=1> subj[nobs];
}
parameters {
  real pAB0;
  real <lower=0> uAB;
  real uASC0;
  real <lower=5> AB0;
  real <lower=2> ASC0;
  real <lower=0> sigma;
  real <lower=0> sigmapAB;
  real <lower=0> sigmauASC;
  real <lower=0> sigmaAB0;
  real <lower=0> sigmaASC0;
  real <lower=0> rpAB[nobs];
  real <lower=0> ruASC[nobs];
  real <lower=0> rAB0[nobs];
  real <lower=0> rASC0[nobs];
}
model {
  vector[2] yhat[nobs];
  vector[2] yhatmu[nobs];
  real eval_time[1];
  //prior distributions
  pAB0 ~ normal(0, 2.5);
  uASC0 ~ normal(0, 2.5);
  uAB ~ uniform(0, 1);
  AB0 ~ uniform(5,10);
  ASC0 ~ uniform(2,4);
  sigmapAB ~ inv_gamma(0.01, 0.01);
  sigmauASC ~ inv_gamma(0.01, 0.01);
  sigmaAB0 ~ inv_gamma(0.01, 0.01);
  sigmaASC0 ~ inv_gamma(0.01, 0.01);
  sigma ~ inv_gamma(0.01, 0.01);
  for (j in 1:nobs) {
    real pAB = exp(pAB0)*exp(rpAB[subj[j]]);
    real uABt = uAB;
    real uASC = exp(uASC0)*exp(ruASC[subj[j]]);
    vector[2] zinitials;
    zinitials[1] = AB0 + rAB0[subj[j]];
    //are we sure this is not yielding to negative values? is it hierarchical?
    zinitials[2] = ASC0 + rASC0[subj[j]];
    // where mu of the r.e. is 0 or -sigmapAB/2
    rpAB[subj[j]] ~ normal(-sigmapAB/2, sigmapAB);
    ruASC[subj[j]] ~ normal(-sigmauASC/2, sigmauASC);
    rAB0[subj[j]] ~ normal(-sigmaAB0/2, sigmaAB0);
    rASC0[subj[j]] ~ normal(-sigmaASC0/2, sigmaASC0);
    eval_time[1] = ts[j];
    if (eval_time[1] == 0){
      yhatmu[j,1] = zinitials[1] - pow(sigma, 2.0)/2;
      yhatmu[j,2] = zinitials[2] - pow(sigma, 2.0)/2;}
    else{
      yhat = ode_rk45(model3, zinitials, t0, eval_time, pAB, uABt, uASC);
      yhatmu[j,2] = yhat[1,1] - pow(sigma, 2.0)/2;
    }  
  //likelihood
    antib[j] ~ lognormal(yhatmu[j,2], sigma); 
  }
}
generated quantities {
  real z_pred[nobs];
  for (t in 1:nobs){
    z_pred[t] = lognormal_rng(antib[t], sigma);
  }
}

I am not an expert, so treat my advice with caution.
For a few of your parameters, the prior does not provide support over the whole range of possible parameter values. For example, uAB has a lower limit of zero and no upper limit but its prior is uniform(0, 1). If the initial value is picked outside of [0, 1], the probability will be 0 and the log probability will be -Infinity, as in the error message. uAB, with a lower bound of zero, is transformed to a new variable equal to log(uAB) and the initial values are chosen on that variable. Any positive value of log(uAB) corresponds to uAB > 1, which will lead to the error. ASC0 and AB0 suffer from the same kind of problem. You need to adjust the prior to a distribution that covers the entire range of the parameter or change the constraint on the parameter to match the prior.
Also, for parameters with uniform priors, are there truly hard boundaries for the possible values? If not, I believe it is preferred to use a prior that has some support outside of the expected range. As I understand it, unnecessary hard boundaries can cause fitting problems.

1 Like