ODE non-converging; Rhat beyond 1

Modeling
1

Hi,

I am running the following model, where nobs is defined as a vector of the dimension 698 participants, such as (3, 3, …, 2, …) indicating the number of observations per individual.
I do not know if there is something wrong with the data, or the initial values, or if the problem belongs to the code.
The Rhat goes beyond 1, especially for the sigma parameters. What could be causing this?

functions {
  real[] sldecay(real t, real[] y, real[] theta,  real[] x_r, int[] x_i) {
    real dydt[3];
    dydt[1] = -theta[1]*y[1];
    dydt[2] = -theta[2]*y[2];
    dydt[3] = theta[1]*y[1] + theta[2]*y[2] - theta[3]*y[3];
    
    return dydt;
  }
}
data {
  int<lower=1> N;
  int<lower=1> nsubjects;
  int<lower=1> subj_ids[N];
  real <lower=0> y[N];
  int <lower=1> ntimes;
  int <lower=1> nobs[nsubjects];
  real ts[ntimes];
  real t0;
}
transformed data {
    real x_r[0];
    int x_i[0]; 
}
parameters {
  real logps;
  real logpl;
  real logmuab;
  real sigmalogmuab;
  
  real logAB0;
  real sigmalogAB0;
  real logps0;
  real logpl0;
  real <lower = 0> sigma;
  
  vector[2] r[nsubjects];
}
model{
  real mu[ntimes, 3];
  real inits[3];
  real theta[3];
  real new_mu[ntimes, 3];
  vector[nsubjects] rAB0;
  vector[nsubjects] rmuab;
  // transformed parameters 
  real ps = exp(logps);
  real pl = exp(logpl);
  real muab = exp(logmuab);
  real sigma_muab = exp(sigmalogmuab);
  real AB0 = exp(logAB0);
  real sigma_AB0 = exp(sigmalogAB0);
  real ps0 = exp(logps0);
  real pl0 = exp(logpl0);
  
  // defining the correlation 
  vector[2] mean_vec;
  matrix[2,2] Sigma; 
  mean_vec[1] = -pow(sigma_AB0, 2.0)/2;
  mean_vec[2] = -pow(sigma_muab, 2.0)/2;
  
  Sigma[1,1] = pow(sigma_AB0, 2.0);
  Sigma[2,2] = pow(sigma_muab, 2.0);
  Sigma[1,2] = 0;
  Sigma[2,1] = 0;
  
  //prior distributions
  logps ~ normal(0, 1);
  logpl ~ normal(0, 1);
  logmuab ~ normal(0, 10);
  sigmalogmuab ~ normal(0, 1);
  logAB0 ~ uniform(5, 10);
  sigmalogAB0 ~ normal(0, 1);
  logps0 ~ normal(0, 1);
  logpl0 ~ normal(0, 1);
  sigma ~ normal(1, 10);
  
  for (j in 1:nsubjects) {
    r[j] ~ multi_normal(mean_vec, Sigma);
    rAB0[j] = r[j,1];
    rmuab[j] = r[j,2];
    //defining the initial conditions
    inits[3] = AB0*exp(rAB0[j]);
    inits[1] = inits[3]/ps0;
    inits[2] = inits[3]/pl0;
    //defining the parameters
    theta[1] = ps;
    theta[2] = pl;
    theta[3] = muab*exp(rmuab[j]);
    // ode solver
    mu[1:nobs[j]] = integrate_ode_rk45(sldecay, inits, t0, ts[1:nobs[j]], theta, x_r, x_i);}
    for (t in 1:nobs[nsubjects]){
      new_mu[t,1] = mu[t,1];
      new_mu[t,2] = mu[t,2];
      new_mu[t,3] = log(mu[t,3]) - pow(sigma, 2.0)/2;
      y[t] ~ lognormal(new_mu[t], sigma);  
    }
}
2

The final loop with y ~ lognormal(...) only visits the one participant whose number is nsubjects (ie. the last participant in the dataset) so all other participants are ignored. Hierarchical parameters can’t be inferred with only one data point and I think that’s why the sampler struggles.

You probably want to concatenate all the participants into one array.

model {
  array[N, 3] real mu;
  array[3] real inits;
  array[3] real theta;
  array[N, 3] real new_mu;
  vector[nsubjects] rAB0;
  vector[nsubjects] rmuab;
  ...
  logps0 ~ normal(0, 1);
  logpl0 ~ normal(0, 1);
  sigma ~ normal(1, 10);

  // n_index is the 
  int n_index = 1;
  for (j in 1 : nsubjects) {
    if (subj_id[n_index] != j) {
      reject("data must be sorted by subject!");
    }
    r[j] ~ multi_normal(mean_vec, Sigma);
    rAB0[j] = r[j, 1];
    rmuab[j] = r[j, 2];
    //defining the initial conditions
    inits[3] = AB0 * exp(rAB0[j]);
    inits[1] = inits[3] / ps0;
    inits[2] = inits[3] / pl0;
    //defining the parameters
    theta[1] = ps;
    theta[2] = pl;
    theta[3] = muab * exp(rmuab[j]);
    // ode solver
    mu[n_index : n_index+nobs[j]-1] = integrate_ode_rk45(sldecay, inits, t0, ts[1 : nobs[j]],
                                         theta, x_r, x_i);
    n_index = n_index + nobs[j];
  }
  for (t in 1 : N) {
    new_mu[t, 1] = mu[t, 1];
    new_mu[t, 2] = mu[t, 2];
    new_mu[t, 3] = log(mu[t, 3]) - pow(sigma, 2.0) / 2;
    y[t] ~ lognormal(new_mu[t], sigma);
  }
}

Note that the model assumes the times at which participants were sampled are always consecutive. This seems unlikely if the typical number of samples is 3 or 2 but the range of ts is 1 to 96 as you mentioned earlier. (in particular, there would have to be a subject with 96 samples to make full use of that range)

3

Thanks!

I am not sure if I understand this correctly, but ts ranges from 1 to 96 but has the dimension of N, which according to the number of observations per subject we would end up with the same number of observations.

4

Do you mean ntimes (the length of ts) is equal to N (the length of y)? If so the indexing in the integration should presumably be

mu[n_index : n_index+nobs[j]-1] = integrate_ode_rk45(sldecay, inits, t0,
                                         ts[n_index:n_index+nobs[j]-1],
                                         theta, x_r, x_i);
5

Do you mean in the code as subj_id the subj_ids from the data?
When I define the array the stan code does not recognize array and it says that it can not recognize that variable. Why is that and does it have an impact on the code?
I have run the cmdstan with the first code written here, and I have gotten the following issue, (whereas in Stan it could run but did not converge):

cmdstanerror
cmdstanerror988×367 18.9 KB

6

I still get the same error message after trying out your code:

cmdstanerror2
cmdstanerror21105×347 26.9 KB

7

You mean it does not recognize the keyword array? Sorry, I used the latest CmdStan syntax for array declaration but RStan is few versions behind and the array declarations should have been

model {
  real mu[N, 3];
  real inits[3];
  real theta[3];
  real new_mu[N, 3];

This makes no difference to the model behaviour, it’s just a style change.

It’s hitting the reject() statement in my code which means that I don’t understand how your data is organized. The code assumes the data is sorted so that subj_ids are in increasing order, and there are exactly nobs[k] consecutive ks in the subj_ids array.

8

Hi,

The data is indeed sorted by subj_ids, and the nobs vector is a vector specifying the number of observations per subject (has dimension 2081, and goes like this: {3, 3, 3, 2, 3…})

9

So if nobs={3, 3, 3, 2, 3…} then subj_ids={1,1,1,2,2,2,3,3,3,4,4,5,5,5,6…}
The reject statement is just checking that subj_ids[sum(nobs[1:k])] == k.

10

Just to be sure I have changed the dimension of the loop to subj_ids[N], which calls all the subjects in the dataset, with repeated measures, and these are sorted.
The ts on the other hand that goes inside the function integrate_ode, goes from 1 to 96 (dimension 2081 as subj_ids). I start to think that there is something wrong with the times defined. I really do not see where the error comes from anymore…

cmderror3
cmderror31108×281 21.3 KB 

functions {
  real[] sldecay(real t, real[] y, real[] theta,  real[] x_r, int[] x_i) {
    real dydt[3];
    dydt[1] = -theta[1]*y[1];
    dydt[2] = -theta[2]*y[2];
    dydt[3] = theta[1]*y[1] + theta[2]*y[2] - theta[3]*y[3];
    
    return dydt;
  }
}
data {
  int<lower=1> N;
  int<lower=1> nsubjects;
  int<lower=1> subj_ids[N];
  real <lower=0> y[N];
  int <lower=1> ntimes;
  int <lower=1> nobs[nsubjects];
  real ts[ntimes];
  real t0;
}
transformed data {
    real x_r[0];
    int x_i[0]; 
}
parameters {
  real logps;
  real logpl;
  real logmuab;
  real sigmalogmuab;
  
  real logAB0;
  real sigmalogAB0;
  real logps0;
  real logpl0;
  real <lower = 0> sigma;
  
  vector[2] r[nsubjects];
}
model{
  real mu[ntimes, 3];
  real inits[3];
  real theta[3];
  real new_mu[ntimes, 3];
  vector[nsubjects] rAB0;
  vector[nsubjects] rmuab;
  //real eval_time[1];
  // transformed parameters 
  real ps = exp(logps);
  real pl = exp(logpl);
  real muab = exp(logmuab);
  real sigma_muab = exp(sigmalogmuab);
  real AB0 = exp(logAB0);
  real sigma_AB0 = exp(sigmalogAB0);
  real ps0 = exp(logps0);
  real pl0 = exp(logpl0);
  // defining the correlation 
  vector[2] mean_vec;
  matrix[2,2] Sigma; 
  mean_vec[1] = -pow(sigma_AB0, 2.0)/2;
  mean_vec[2] = -pow(sigma_muab, 2.0)/2;
  
  Sigma[1,1] = pow(sigma_AB0, 2.0);
  Sigma[2,2] = pow(sigma_muab, 2.0);
  Sigma[1,2] = 0;
  Sigma[2,1] = 0;
  
  //prior distributions
  logps ~ normal(0, 1);
  logpl ~ normal(0, 1);
  logmuab ~ normal(0, 10);
  sigmalogmuab ~ normal(0, 1);
  logAB0 ~ uniform(5, 10);
  sigmalogAB0 ~ normal(0, 1);
  logps0 ~ normal(0, 1);
  logpl0 ~ normal(0, 1);
  sigma ~ normal(1, 10);
  
  for (j in 1:subj_ids[N]) {
    r[j] ~ multi_normal(mean_vec, Sigma);
    rAB0[j] = r[j,1];
    rmuab[j] = r[j,2];
    //defining the initial conditions
    inits[3] = AB0*exp(rAB0[j]);
    inits[1] = inits[3]/ps0;
    inits[2] = inits[3]/pl0;
    //defining the parameters
    theta[1] = ps;
    theta[2] = pl;
    theta[3] = muab*exp(rmuab[j]);
    // ode solver
     mu[1: nobs[j]] = integrate_ode_rk45(sldecay, inits, t0,
                                         ts[1: nobs[j]],
                                         theta, x_r, x_i);
      new_mu[j,1] = mu[j,1];
      new_mu[j,2] = mu[j,2];
      new_mu[j,3] = log(mu[j,3]) - pow(sigma, 2.0)/2;
      y[j] ~ lognormal(new_mu[j, 3], sigma);  
    }
}
11

The indexing in that model does not make sense—it has only one y[j] ~ lognormal(...) per subject.

12

What should then be the index for y[j]?
I have tried now the following code but I still get the same error message.

{
  real mu[ntimes, 3];
  real inits[3];
  real theta[3];
  real new_mu[ntimes, 3];
  vector[nsubjects] rAB0;
  vector[nsubjects] rmuab;
  //real eval_time[1];
  // transformed parameters 
  real ps = exp(logps);
  real pl = exp(logpl);
  real muab = exp(logmuab);
  real sigma_muab = exp(sigmalogmuab);
  real AB0 = exp(logAB0);
  real sigma_AB0 = exp(sigmalogAB0);
  real ps0 = exp(logps0);
  real pl0 = exp(logpl0);
  // defining the correlation 
  vector[2] mean_vec;
  matrix[2,2] Sigma; 
  mean_vec[1] = -pow(sigma_AB0, 2.0)/2;
  mean_vec[2] = -pow(sigma_muab, 2.0)/2;
  
  Sigma[1,1] = pow(sigma_AB0, 2.0);
  Sigma[2,2] = pow(sigma_muab, 2.0);
  Sigma[1,2] = 0;
  Sigma[2,1] = 0;
  
  //prior distributions
  logps ~ normal(0, 1);
  logpl ~ normal(0, 1);
  logmuab ~ normal(0, 10);
  sigmalogmuab ~ normal(0, 1);
  logAB0 ~ uniform(5, 10);
  sigmalogAB0 ~ normal(0, 1);
  logps0 ~ normal(0, 1);
  logpl0 ~ normal(0, 1);
  sigma ~ normal(1, 10);
  
  for (j in 1:nsubjects) {
    r[j] ~ multi_normal(mean_vec, Sigma);
    rAB0[j] = r[j,1];
    rmuab[j] = r[j,2];
    //defining the initial conditions
    inits[3] = AB0*exp(rAB0[j]);
    inits[1] = inits[3]/ps0;
    inits[2] = inits[3]/pl0;
    //defining the parameters
    theta[1] = ps;
    theta[2] = pl;
    theta[3] = muab*exp(rmuab[j]);
    // ode solver
     mu[1: nobs[j]] = integrate_ode_rk45(sldecay, inits, t0,
                                         ts[1: nobs[j]],
                                         theta, x_r, x_i);
      for (i in 1:N){
      new_mu[i,1] = mu[i,1];
      new_mu[i,2] = mu[i,2];
      new_mu[i,3] = log(mu[i,3]) - pow(sigma, 2.0)/2;
      y[i] ~ lognormal(new_mu[i, 3], sigma);  
    }
  }
}
13

The index at which the data for subject j is stored; i guess it’s something like

    sum(nobs[1:j]) : sum(nobs[1:j+1])

It would be easier to integrate each data point individually, without trying group them by subbject

model {
  ...
  sigma ~ normal(1, 10);
  
  for (j in 1:nsubjects) {
    r[j] ~ multi_normal(mean_vec, Sigma);
    rAB0[j] = r[j,1];
    rmuab[j] = r[j,2];
  }
  for (j in 1:N) {
    //defining the initial conditions
    inits[3] = AB0*exp(rAB0[subj_ids[j]]);
    inits[1] = inits[3]/ps0;
    inits[2] = inits[3]/pl0;
    //defining the parameters
    theta[1] = ps;
    theta[2] = pl;
    theta[3] = muab*exp(rmuab[j]);
    // ode solver
    mu[j:j] = integrate_ode_rk45(sldecay, inits, t0,
                                         ts[j:j],
                                         theta, x_r, x_i);
    new_mu[j,1] = mu[j,1];
    new_mu[j,2] = mu[j,2];
    new_mu[j,3] = log(mu[j,3]) - pow(sigma, 2.0)/2;
    y[j] ~ lognormal(new_mu[j, 3], sigma);  
  }
}

But I think last time it was too slow that way. Solving the differential equation in closed form would help with that. Simplifying the result from WolframAlpha, I get

y_{1}=C_{1}e^{-\theta_{1}t}
y_{2}=C_{2}e^{-\theta_{2}t}
y_{3}=\frac{C_{1}\theta_{1}}{\theta_{3}-\theta_{1}}e^{-\theta_{1}t}+\frac{C_{2}\theta_{2}}{\theta_{3}-\theta_{2}}e^{-\theta_{2}t}+C_{3}e^{-\theta_{3}t}

and you can check it does satisfy the equation.
So maybe this will work

functions {
  real[] analytic_int(real[] inits, real t0, real ts, real[] theta) {
    real y[3];
    real dt = ts - t0;
    y[1] = inits[1]*exp(-theta[1]*dt);
    y[2] = inits[2]*exp(-theta[2]*dt);
    real C = inits[3] + inits[1]*theta[1]/(theta[1]-theta[3]) + inits[2]*theta[2]/(theta[2]-theta[3]);
    y[3] = theta[1]*y[1]/(theta[3] - theta[1]) + theta[2]*y[2]/(theta[3] - theta[2]) + C*exp(theta[3]*dt);
    return inits;
  }
}
...
model {
  ...
    // ode solver
    //mu[j:j] = integrate_ode_rk45(sldecay, inits, t0,
    //                                     ts[j:j],
    //                                     theta, x_r, x_i);
    mu[j] = analytic_int(inits, t0, ts[j], theta);
    ...
}
1 Like
14

Unfortunately none of the solutions are working :(
With the code below I get the following error message:

cmdstanjj
cmdstanjj1104×223 16.4 KB 

functions {
  real[] sldecay(real t, real[] y, real[] theta,  real[] x_r, int[] x_i) {
    real dydt[3];
    dydt[1] = -theta[1]*y[1];
    dydt[2] = -theta[2]*y[2];
    dydt[3] = theta[1]*y[1] + theta[2]*y[2] - theta[3]*y[3];
    
    return dydt;
  }
}
data {
  int<lower=1> N;
  int<lower=1> nsubjects;
  int<lower=1> subj_ids[N];
  real <lower=0> y[N];
  int <lower=1> ntimes;
  int <lower=1> nobs[nsubjects];
  real ts[ntimes];
  real t0;
}
transformed data {
    real x_r[0];
    int x_i[0]; 
}
parameters {
  real logps;
  real logpl;
  real logmuab;
  real sigmalogmuab;
  
  real logAB0;
  real sigmalogAB0;
  real logps0;
  real logpl0;
  real <lower = 0> sigma;
  
  vector[2] r[nsubjects];
}
model{
  real mu[ntimes, 3];
  real inits[3];
  real theta[3];
  real new_mu[ntimes, 3];
  vector[nsubjects] rAB0;
  vector[nsubjects] rmuab;
  real eval_time[1];
  // transformed parameters 
  real ps = exp(logps);
  real pl = exp(logpl);
  real muab = exp(logmuab);
  real sigma_muab = exp(sigmalogmuab);
  real AB0 = exp(logAB0);
  real sigma_AB0 = exp(sigmalogAB0);
  real ps0 = exp(logps0);
  real pl0 = exp(logpl0);
  // defining the correlation 
  vector[2] mean_vec;
  matrix[2,2] Sigma; 
  mean_vec[1] = -pow(sigma_AB0, 2.0)/2;
  mean_vec[2] = -pow(sigma_muab, 2.0)/2;
  
  Sigma[1,1] = pow(sigma_AB0, 2.0);
  Sigma[2,2] = pow(sigma_muab, 2.0);
  Sigma[1,2] = 0;
  Sigma[2,1] = 0;
  
  //prior distributions
  logps ~ normal(0, 1);
  logpl ~ normal(0, 1);
  logmuab ~ normal(0, 10);
  sigmalogmuab ~ normal(0, 1);
  logAB0 ~ uniform(5, 10);
  sigmalogAB0 ~ normal(0, 1);
  logps0 ~ normal(0, 1);
  logpl0 ~ normal(0, 1);
  sigma ~ normal(1, 10);
  
  for (j in 1:nsubjects) {
    r[j] ~ multi_normal(mean_vec, Sigma);
    rAB0[j] = r[j,1];
    rmuab[j] = r[j,2];
  }
  for (i in 1:N){
    //defining the initial conditions
    inits[3] = AB0*exp(rAB0[subj_ids[i]]);
    inits[1] = inits[3]/ps0;
    inits[2] = inits[3]/pl0;
    //defining the parameters
    theta[1] = ps;
    theta[2] = pl;
    theta[3] = muab*exp(rmuab[subj_ids[i]]);
    eval_time[1] = ts[i];
    if (eval_time[1] == 0){
      new_mu[i,1] = inits[1];
      new_mu[i,2] = inits[2];
      new_mu[i,3] = log(inits[3]) - pow(sigma, 2.0)/2;}
    else{
    // ode solver
     mu = integrate_ode_rk45(sldecay, inits, t0,
                                         eval_time,
                                         theta, x_r, x_i);
      new_mu[i,1] = mu[1,1];
      new_mu[i,2] = mu[1,2];
      new_mu[i,3] = log(mu[1,3]) - pow(sigma, 2.0)/2;
    }
      y[i] ~ lognormal(new_mu[i, 3], sigma);  
    }
}

And solving it analytically, with the code below, I get the following error message:
log0

functions {
  real[] analytic_int(real[] inits, real t0, real ts, real[] theta) {
    real y[3];
    real dt = ts - t0;
    real C;
    y[1] = inits[1]*exp(-theta[1]*dt);
    y[2] = inits[2]*exp(-theta[2]*dt);
    C = inits[3] + inits[1]*theta[1]/(theta[1]-theta[3]) + inits[2]*theta[2]/(theta[2]-theta[3]);
    y[3] = theta[1]*y[1]/(theta[3] - theta[1]) + theta[2]*y[2]/(theta[3] - theta[2]) + C*exp(theta[3]*dt);
    return inits;
  }
}
data {
  int<lower=1> N;
  int<lower=1> nsubjects;
  int<lower=1> subj_ids[N];
  real <lower=0> y[N];
  int <lower=1> ntimes;
  int <lower=1> nobs[nsubjects];
  real ts[ntimes];
  real t0;
}
transformed data {
    real x_r[0];
    int x_i[0]; 
}
parameters {
  real logps;
  real logpl;
  real logmuab;
  real sigmalogmuab;
  
  real logAB0;
  real sigmalogAB0;
  real logps0;
  real logpl0;
  real <lower = 0> sigma;
  
  vector[2] r[nsubjects];
}
model{
  real mu[ntimes, 3];
  real inits[3];
  real theta[3];
  real new_mu[ntimes, 3];
  vector[nsubjects] rAB0;
  vector[nsubjects] rmuab;
  real eval_time[1];
  // transformed parameters 
  real ps = exp(logps);
  real pl = exp(logpl);
  real muab = exp(logmuab);
  real sigma_muab = exp(sigmalogmuab);
  real AB0 = exp(logAB0);
  real sigma_AB0 = exp(sigmalogAB0);
  real ps0 = exp(logps0);
  real pl0 = exp(logpl0);
  // defining the correlation 
  vector[2] mean_vec;
  matrix[2,2] Sigma; 
  mean_vec[1] = -pow(sigma_AB0, 2.0)/2;
  mean_vec[2] = -pow(sigma_muab, 2.0)/2;
  
  Sigma[1,1] = pow(sigma_AB0, 2.0);
  Sigma[2,2] = pow(sigma_muab, 2.0);
  Sigma[1,2] = 0;
  Sigma[2,1] = 0;
  
  //prior distributions
  logps ~ normal(0, 1);
  logpl ~ normal(0, 1);
  logmuab ~ normal(0, 10);
  sigmalogmuab ~ normal(0, 1);
  logAB0 ~ uniform(5, 10);
  sigmalogAB0 ~ normal(0, 1);
  logps0 ~ normal(0, 1);
  logpl0 ~ normal(0, 1);
  sigma ~ normal(1, 10);
  
  for (j in 1:nsubjects) {
    r[j] ~ multi_normal(mean_vec, Sigma);
    rAB0[j] = r[j,1];
    rmuab[j] = r[j,2];
  }
  for (i in 1:N){
    //defining the initial conditions
    inits[3] = AB0*exp(rAB0[subj_ids[i]]);
    inits[1] = inits[3]/ps0;
    inits[2] = inits[3]/pl0;
    //defining the parameters
    theta[1] = ps;
    theta[2] = pl;
    theta[3] = muab*exp(rmuab[subj_ids[i]]);
    // ode solver
      mu[i] = analytic_int(inits, t0, ts[i], theta);
      new_mu[i,1] = mu[1,1];
      new_mu[i,2] = mu[1,2];
      new_mu[i,3] = log(mu[1,3]) - pow(sigma, 2.0)/2;
      y[i] ~ lognormal(new_mu[i, 3], sigma);  
    }
}
15

These should be

      new_mu[i,1] = mu[i,1];
      new_mu[i,2] = mu[i,2];
      new_mu[i,3] = log(mu[i,3]) - pow(sigma, 2.0)/2;

That’s not what’s causing the error though, and nothing else stands out to me.

Are there any j such that y[j]==0. Lognormal is not compatible with exact zeros. If the measurement is rounded to zero for small values you have to change the likelihood, for example something like this may work:

      if (y[i] > 0.001) {
          y[i] ~ lognormal(new_mu[i, 3], sigma);
      } else {
          target += lognormal_lcdf(0.001, new_mu[i,3], sigma);
      }

where I have chosen 0.001 rather arbitrarily to represent the smallest measurable nonzero value.

1 Like
16

this needs a <lower = 5, upper = 10> bound to match its prior:

otherwise the initial values will never be feasible. This is also a very strange prior for a logarithmic quantity.

2 Likes
17

Hi,

Thanks for all the help. Unfortunately, I do not have zero values in my response. As a matter of fact, the lowest is 36 (therefore I have chosen those priors). Also, cause when I ran the model without any random effects the value for logAB0 was converging around 5. Do you have any better suggestions on this?

18

I just noticed the code I provided was missing a minus in front of theta[3]. That should have been:

   y[3] = theta[1]*y[1]/(theta[3] - theta[1]) + theta[2]*y[2]/(theta[3] - theta[2]) + C*exp(-theta[3]*dt);

It’s also possible this is too inaccurate when \theta_1 \approx \theta_3 or \theta_1 \approx \theta_3. Sampling with init=0 argument wont work because of that but you could still try init=0.001 to get more predictable initialization.

Is there a specific reason for uniform(5,10) prior for logAB0 other than that it looks consistent with the data? If the posterior concentrates near prior boundaries you’d better have a good reason for the boundary. Maybe you could try normal(5,10) instead.

The problem might be limited to a few subjects so you could try fitting only a subset of the data.

19

Currently with the code below, and adding initial values to the function sample;
rAB0 = rep(0, participants)
rmuab = rep(0, participants)
r = cbind(rAB0, rmuab)
initebl = function(){
list(logps = 1, logpl = 1, logmuab = 1,
sigmalogmuab = 0.5, logAB0 = 5, sigmalogAB0 = 0.6,
logps0 = 1, logpl0 = 1, sigma = 0.5, r = r)
}
I get the following error message, which I do not see in which line is the issue. I have also tried different initial values but that hasn’t worked either.

cmdinit
cmdinit1062×212 8.99 KB 

functions {
  real[] analytic_int(real[] inits, real t0, real ts, real[] theta) {
    real y[3];
    real dt = ts - t0;
    real C;
    y[1] = inits[1]*exp(-theta[1]*dt);
    y[2] = inits[2]*exp(-theta[2]*dt);
    C = inits[3] + inits[1]*theta[1]/(theta[1]-theta[3]) + inits[2]*theta[2]/(theta[2]-theta[3]);
    y[3] = theta[1]*y[1]/(theta[3] - theta[1]) + theta[2]*y[2]/(theta[3] - theta[2]) + C*exp(-theta[3]*dt);
    return inits;
  }
}
data {
  int<lower=1> N;
  int<lower=1> nsubjects;
  int<lower=1> subj_ids[N];
  real <lower=0> y[N];
  int <lower=1> ntimes;
  int <lower=1> nobs[nsubjects];
  real ts[ntimes];
  real t0;
}
transformed data {
    real x_r[0];
    int x_i[0]; 
}
parameters {
  real  logps;
  real  logpl;
  real  logmuab;
  real   sigmalogmuab;
  
  real <lower= 0> logAB0;
  real  sigmalogAB0;
  real  logps0;
  real  logpl0;
  real  sigma;
  
  vector[2] r[nsubjects];
}
model{
  real mu[ntimes, 3];
  real inits[3];
  real theta[3];
  real new_mu[ntimes, 3];
  vector[nsubjects] rAB0;
  vector[nsubjects] rmuab;
  real eval_time[1];
  // transformed parameters 
  real ps = exp(logps);
  real pl = exp(logpl);
  real muab = exp(logmuab);
  real sigma_muab = exp(sigmalogmuab);
  real AB0 = exp(logAB0);
  real sigma_AB0 = exp(sigmalogAB0);
  real ps0 = exp(logps0);
  real pl0 = exp(logpl0);
  real sigmap = exp(sigma);
  // defining the correlation 
  vector[2] mean_vec;
  matrix[2,2] Sigma; 
  mean_vec[1] = -pow(sigma_AB0, 2.0)/2;
  mean_vec[2] = -pow(sigma_muab, 2.0)/2;
  
  Sigma[1,1] = pow(sigma_AB0, 2.0);
  Sigma[2,2] = pow(sigma_muab, 2.0);
  Sigma[1,2] = 0;
  Sigma[2,1] = 0;
  
  //prior distributions
  logps ~ normal(0, 1);
  logpl ~ normal(0, 1);
  logmuab ~ normal(0, 10);
  sigmalogmuab ~ normal(0, 1);
  logAB0 ~ normal(5, .5);
  sigmalogAB0 ~ normal(0, 1);
  logps0 ~ normal(0, 1);
  logpl0 ~ normal(0, 1);
  sigma ~ normal(0, 1);
  
  for (j in 1:nsubjects) {
    r[j] ~ multi_normal(mean_vec, Sigma);
    rAB0[j] = r[j,1];
    rmuab[j] = r[j,2];
  }
  for (i in 1:N){
    //defining the initial conditions
    inits[3] = AB0*exp(rAB0[subj_ids[i]]);
    inits[1] = inits[3]/ps0;
    inits[2] = inits[3]/pl0;
    //defining the parameters
    theta[1] = ps;
    theta[2] = pl;
    theta[3] = muab*exp(rmuab[subj_ids[i]]);
    // ode solver
      mu[i] = analytic_int(inits, t0, ts[i], theta);
      new_mu[i,1] = mu[i,1];
      new_mu[i,2] = mu[i,2];
      new_mu[i,3] = log(mu[i,3]) - pow(sigma, 2.0)/2;
      y[i] ~ lognormal(new_mu[i, 3], sigma);
    }
}
20

Those initial values make all theta equal and division by (theta[3] - theta[1]) cannot work.
Have you tried initial values where logpl != logmuab != logps?
Alternatively I think you can special-case the “all theta are equal” scenario:

if (theta[1] == theta[3] && theta[2] == theta[3]) {
  y[3] = theta[1]*t*y[1] + theta[2]*dt*y[2] + inits[3]*exp(-theta[3]*dt)
} else {
   y[3] = theta[1]*y[1]/(theta[3] - theta[1]) + theta[2]*y[2]/(theta[3] - theta[2]) + C*exp(-theta[3]*dt);
}