data {
 int<lower=1> N; //number of reserve data points 
 int<lower=1> RE;//number of remote data points
 real b[N]; //response variable (y axes) in reserves(log-transformed)
 real br[RE]; //response variable (y axes) in remote(log-transformed)
 //explanatory variables reserves
 int<lower=0> ag[N]; //age reserve
 int rh_c[N]; // sampling covariate 1 
 int rh_l[N]; //sampling covariate 2
 int sm_p[N]; //sampling covariate 3 
 //explanantory variables remote
 int rh_c2[RE]; // sampling covariate 1 
 int rh_l2[RE];  // sampling covariate 2 
 int sm_d[RE];  // sampling covariate 4 
 //random effects
 int<lower=1> R; //number of data regions for reserves(groups)
 int<lower=1, upper=R> pr[N]; //region id for reserves
 int<lower=1> R2; //number of data regions for remote(groups)
 int<lower=1, upper=R2> pr2[RE]; //region id for remote
}
parameters {
 vector[4] beta; //slopes of explanatory variables 
 real<lower=0> sigma_e; //error sd reserves
 real<lower=0> sigma_r; //error sd remote
 vector[R] u_tilde; //for non centered
 real<lower=0> sigma_u; //group sd reserves
 vector[R2] u2_tilde; //for non centered
 real<lower=0> sigma_u2; //group sd remote
 real<lower=0> r;//intrinsic growth rate
 real<lower=0> bmin; //response at reserve age 0
 real<lower=0> B0; //carrying capacity
 real<lower=0,upper=1> p; //portion of pop growth exported
}
transformed parameters {
 vector[N] mu;//mean response variable reserves (logistic growth)
 vector[N] k; //site-specific reserve carrying capacity given covariates and exports
 vector[RE] mu2;//mean response variable remote (linear)
 vector[RE] k2;//site-specific remote carrying capapcity given enevironmental factors
 vector[R] u; //region intercepts reserves
 vector[R2] u2; //region intercepts remote
 real reserve_eq;// reserve carrying capacity at average conditions
 real r_apparent; //reserve recovery rates
 //reserve submodel
 r_apparent=r*(1-p); //recovery in reserves a function of "true" intrinsic growth rate and portion exports
 reserve_eq=r*(1-p)/(r/B0); //or B0*(1-p), reserve carrying capacity (if p=0, reserve_eq=B0)
 for(i in 1:R){
    u[i] = sigma_u .* u_tilde[i];//non-centered
  }
 for (i in 1:N) {
  k[i] = log(reserve_eq) +beta[1] * rh_c[i] +beta[2]*rh_l[i]; 
  mu[i] = log(exp(k[i])/(1+((exp(k[i])-bmin)/bmin)*exp(-(r_apparent)*ag[i])))+ beta[3]*sm_p[i]+u[pr[i]];
 }
 //remote submodel
  for(i in 1:R2){
    u2[i] = sigma_u2 .* u2_tilde[i];
  }
 for (i in 1:RE) {
  k2[i]=log(B0)+beta[1]*rh_c2[i]+ beta[2]*rh_l2[i]+beta[4]*sm_d[i];
  mu2[i] = log(exp(k2[i]))+u2[pr2[i]];
 }
}
model {
 //priors
 beta[1] ~ normal (0,10); //prior slope
 beta[2] ~ normal (0,10); //prior slope
 beta[3] ~ normal (0,10); //prior slope
 beta[4] ~ normal (0,10); //prior slope
 sigma_u2 ~ cauchy (0,1); //prior sd for varying intercept(strong based on Gelman)
 u2_tilde ~ normal (0, 1); //prior noncenteredversion
 sigma_u ~ cauchy (0,1); //prior sd for varying intercept (strong based on Gelman)
 u_tilde ~ normal(0, 1); //prior noncenteredversion
 r ~ normal (0.2,1); //weekly informative prior log biomass intrinsic growth rate
 sigma_e ~ cauchy(0,2.5); //uninformative prior sd
 sigma_r ~ cauchy(0,2.5); //uninformative prior sd
 bmin ~ normal (0,200); //weakly informative prior reserve biomass at age 0
 B0 ~ normal(120,200);//weakly informative prior for unfished biomass
 p~ uniform(0,1);// prior for portion of growth exported
//likelihood  
 for(n in 1:N){
            target += normal_lpdf(b[n] | mu[n], sigma_e);
    }
  for(n in 1:RE){
            target += normal_lpdf(br[n] | mu2[n], sigma_r);
    }
}
generated quantities {
 vector[N+RE] log_lik;
 for (n in 1:N) {
 log_lik[n] = normal_lpdf(b[n]| mu, sigma_e);
 }
 for (n in (N+1):(N+RE)) {
 log_lik[n] = normal_lpdf(br[(n-N)]| mu2, sigma_r);
 }
}