A problem when I try to combine mutiple logistic growth models together

Modeling
1

My final goal is to build a hierarchical model of logistic population growth models for multiple fish species. The problem I have now is that the model work “ok” independently, but can not work if I combine them together.
The error I get is as follows:
"Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: model474863befb1__namespace::log_prob: P[sym1__, sym2__] is 0.0001, but must be greater than or equal to 0.001 (in ‘string’, line 24, column 1 to column 47)
Chain 1:
Chain 1: Initialization between (-2, 2) failed after 100 attempts.
Chain 1: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”
[1] “error occurred during calling the sampler; sampling not done” "

Notice: The time-series data have missing entries, which are labeled as “0” and filtered out in the model section in stan. The demo data are provided.

Independent stan model is as follows

data{
  int<lower=1> N;
  vector[N] Catch;
  vector[N] CPUE;
  vector[N] Stocking;
  vector[N] Environment;
  real area;
}

parameters{
  real logK;
  real logR;
  real logQ;
  real logSigma;
  real impact_E;
  real<lower=0.001,upper=1.5> P0;
}

transformed parameters{
 real<lower=0> k;
 real<lower=0> r;
 real<lower=0> q;
 real<lower=0> sigma;
 vector<lower=0>[N] P;
 vector[N] CPUE_hat;

 k = exp(logK);
 r = exp(logR);
 q = exp(logQ);
 sigma = exp(logSigma);

 P[1] = P0;
 CPUE_hat[1] = k*q*P[1];

 for (t in 2:N){
      P[t]=P[t-1]+(r*exp(impact_E*Environment[t-1]))*P[t-1]*(1-P[t-1])-Catch[t-1]/(k*area)+Stocking[t-1]/(k*area);
      if(P[t]<0.0001){
        P[t]=0.0001; // prevent negative P
      } else {
        P[t]=P[t];
 }
  CPUE_hat[t] = k*area*q*P[t];
}
}

model{
  logR ~ uniform(log(0.1),log(0.5));
  logK ~ uniform(log(1),log(1000));
  logQ ~ uniform(log(0.00001),log(1));
for (t in 1:N){
  if (CPUE[t]>0){
    log(CPUE[t]) ~ normal(log(CPUE_hat[t]),sigma); // filter out the missing entries
  }
 }
}

Combined model is as follows.

data{
  int Nstocks;
  int Nyears;
  matrix[Nyears,Nstocks] CPUE;
  matrix[Nyears,Nstocks] Catch;
  matrix[Nyears,Nstocks] Stocking;
  matrix[Nyears,Nstocks] Environment;
  real Area [Nstocks];
  int Nsurvey_year[Nstocks];
}

parameters{
  vector[Nstocks] logK;
  vector[Nstocks] logR;
  vector[Nstocks] logQ;
  vector[Nstocks] logSigma;
  vector[Nstocks] impact_E;
  vector<lower=0.001,upper=2>[Nstocks] P0;
}

transformed parameters{
 vector[Nstocks] k;
 vector[Nstocks] r;
 vector[Nstocks] q;
 vector[Nstocks] sigma;
 matrix<lower=0.001,upper=2>[Nyears,Nstocks] P;
 matrix[Nyears,Nstocks] CPUE_hat;

 for (i in 1:Nstocks){
   k[i] = exp(logK[i]);
   r[i] = exp(logR[i]);
   q[i] = exp(logQ[i]);
   sigma[i] = exp(logSigma[i]);
   P[1,i] = P0[i];
   for (o in 2:Nsurvey_year[i]){
     P[o,i] = P[o-1,i]+(r[i]*exp(impact_E[i]*Environment[o-1,i]))*P[o-1,i]*(1-P[o-1,i])-Catch[o-1,i]/(k[i])+Stocking[o-1,i]/(k[i]);
     if (P[o,i]<0.0001){
       P[o,i] = 0.0001;
     } else {
       P[o,i] = P[o,i];
     }
     CPUE_hat[o,i] = k[i]*P[o,i]*q[i];
   }
  }
 }

model{
  for (i in 1:Nstocks){
    logR[i] ~ uniform(log(0.1),log(0.5));
    logK[i] ~ uniform(log(1),log(1000));
    logQ[i] ~ uniform(log(0.00001),log(1));
    for (o in 1:Nsurvey_year[i]){
      if (CPUE[o,i]>0){
        log(CPUE[o,i]) ~ normal(log(CPUE_hat[o,i]),sigma[i]);
      }
    }
  }
}

combine_data.RData (1.3 KB)

2

The error message is:

P[sym1__, sym2__] is 0.0001, but must be greater than or equal to 0.001 (in ‘string’, line 24, column 1 to column 47)

Which is occurring because you’ve declared the lower-bound for P as 0.001:

 matrix<lower=0.001,upper=2>[Nyears,Nstocks] P;

But you’re setting values to 0.0001:

     if (P[o,i]<0.0001){
       P[o,i] = 0.0001;
     }

You need to update either the lower-bound or your default value here