Rejecting intial value: location parameter is - Inf but must be finite!

Modeling
1

Hi all,

I have been facing some issues when coding this exponential decay with random effects. I have a longitudinal data of 7 measurements per patient.
When I run the following code with a sub-sample of 3 measurements per individual I don’t get any errors. However, when I run the same with the complete dataset I get the following error message:

rstanerror
rstanerror1133×208 9.73 KB

I have also defined initial values as:
initebl = function(){
nsubjects ← 700
pop_var = c(0.5, 0.5)
mua_pop = rlnorm(1, 0, pop_var[1])
AB0_pop = exp(rnorm(1, log(37), pop_var[2]))
omega ← abs(rnorm(2, 0, pop_var))
theta_pop ← c(mua_pop, AB0_pop)
theta ← matrix(NA, nsubjects, length(theta_pop))
for (i in 1:nsubjects){
theta[i, ] = rlnorm(length(theta_pop), log(theta_pop), omega)
}
list(mua_pop = mua_pop, AB0_pop = AB0_pop, realsigmaA = abs(rnorm(1, 0, 1)),
sigma = abs(rnorm(1, 0, 1)), omega = omega, theta = theta)
}

Why am I getting this initialization error message? How could I avoid it?

functions{
  real[] numerical_int(real t0, real ts, real[] theta){
    real dt = ts - t0;
    real y[1];
    y[1] = theta[2]*exp(-theta[1]*dt);
    return y;
  }
}
data{
  int <lower=1> N;
  int <lower=1> nsubjects;
  int <lower=1> subj_ids[N];
  real <lower=0> y[N];
  real ts[N];
  int <lower=1> iobs[N];
  real t0;
}
transformed data{
  int nTheta = 2;
  int niv = 2;
  //prior sd
  real prior_sd[niv] = {0.5, 0.5};
}
parameters{
  //population parameters
  real <lower=0> mua_pop;
  real <lower=0> AB0_pop;
  //residual error
  real <lower=0> realsigmaA;
  real <lower=0> sigma;
  //inter-individual variability
  vector <lower=0>[niv] omega;
  real <lower=0> theta[nsubjects, nTheta];
}
transformed parameters{
  vector <lower=0>[nTheta] theta_pop = to_vector({mua_pop, AB0_pop});
  real mu[N, 1];
  //vector [new_mu];
  //individual parameters
  vector <lower=0> [nsubjects] mua = to_vector(theta[, 1]);
  vector <lower=0> [nsubjects] AB0 = to_vector(theta[, 2]);
  
  for (j in 1:N){
    mu[j, ] = numerical_int(t0, ts[j], {mua[subj_ids[j]], AB0[subj_ids[j]]});
  }
}
model{
  //priors 
  mua_pop ~ lognormal(0, prior_sd[1]);
  AB0_pop ~ lognormal(log(37), prior_sd[2]);
  realsigmaA ~ normal(0, 1);
  sigma ~ normal(0, realsigmaA);
  omega ~ normal(0, 1);
  
   // interindividual variability
  for (j in 1:nsubjects){
    theta[j, ] ~ lognormal(log(theta_pop), omega);
  }

  // likelihood
  y ~ lognormal(log(mu[iobs, 1]), sigma);
} 
generated quantities {
  real pred_antibody[N];
  for (i in 1:N) {
  pred_antibody[i] = exp(normal_rng(log(y[i]), sigma));
  }
}
2

The error is due to getting a -infinity log density (0 density), which means either you’ve given it parameters out of support (the usual problem) or that there are floating point issues somewhere. The usual solution for parameters out of support is to add appropriate constraints to the parameters so that any value of parameters satisfying the constraints has a finite log density.

You can insert print("location 1234 target() = ", target()) statements to see when things go off the rails—it prints the accumulated log density up to that point.

Why is mu[iobs, 1] the only value of mu used in the program? You can just compute that value rather than all N in the numerical_int loop.

Finally, the generated quantities can be done as a one-liner (and brought up to our current array syntax—the old style you used is being eliminated next release):

array[N] real pred_antibody = exp(normal_rng(log(y), sigma);

Rather than to_vector on an array, you can use the row vector constructor and transpose:

vector<lower=0>[nTheta] theta_pop = [mua_pop, AB0_pop]';

If you are assuming nTheta = 2 and hard coding that as you seem to do for mua and AB0, then it’s a lot easier just to declare those two vectors and write the sampling for theta on two lines.