functions { vector growth(real t, vector y, vector theta) { vector[1] dydt; dydt[1] = (1-inv_logit(log(19)*(y[1] - theta[3]/theta[4]) / theta[5]))*(theta[1]-theta[2]*y[1]); return dydt; } /// PARTIAL SUM FOR RK45s real tags_sum_lpdf(real[] slice_obs_incr, int start, int end, real growth_obs_error_sd, real mat_scaler, real mat_p95, real expo, vector log_growth_as, vector log_growth_ks, vector Pred_Tag_Allocs, vector LENGTH_AT_RELEASE, real[] DAYS) { int N=end-start+1; vector[N] this_Pred_Tag_Inc; for (t in start:end) { vector[5] theta; // [1]anabolism [2] catabolsim [3] lmat [4] investment theta[1] = exp(log_growth_as[t]); theta[2] = exp(log_growth_ks[t]); theta[3] = Pred_Tag_Allocs[t]; theta[4] = mat_scaler; theta[5] = mat_p95/mat_scaler^expo; vector[1] er[1] = ode_rk45_tol(growth, LENGTH_AT_RELEASE[t:t], 0.0, DAYS[t:t], 1e-5, 1e-5, 100000, theta); real err = er[1,1]; real lar = LENGTH_AT_RELEASE[t]; this_Pred_Tag_Inc[t-start+1] = (err - lar) > 0 ? (err - lar) : 0; } return cauchy_lupdf(slice_obs_incr | this_Pred_Tag_Inc, growth_obs_error_sd); } /// partial sum for MAT obs real mat_sum_lpmf(int[] slice_n_mat, int start, int end, int[] N_MAT_OBS, real mat_site_sigma, real mat_p50, real mat_p95, vector unscaled_site_fxs, int[] MAT_SITE, vector MAT_SIZE) { return binomial_logit_lupmf(slice_n_mat | N_MAT_OBS[start:end], log(19)*(MAT_SIZE[start:end] - (mat_p50 + (mat_site_sigma*unscaled_site_fxs[MAT_SITE[start:end]])))/ mat_p95); } } data{ // this is over all maturation and growth sites int N_SITES; int N_INT; int N_QMAS; int L; vector[L] LENGTHS; int grainsize; // maturity int N_MAT; int MAT_SITE[N_MAT]; int MAT_QMA[N_MAT]; int OBS_MAT[N_MAT]; int N_MAT_OBS[N_MAT]; vector[N_MAT] MAT_SIZE; // maturity priors real mat_p50_PRIOR_MEAN; real mat_p50_PRIOR_SD; real mat_p95_PRIOR_MEAN; real mat_p95_PRIOR_SD; real mat_scaler_PRIOR_A; real mat_scaler_PRIOR_B; real mat_site_sigma_PRIOR_SD; real expo_PRIOR; // Growth / tag-recapture data int N_TAGS; //vector[N_TAGS] METHOD; int GROWTH_SITE[N_TAGS]; int GROWTH_QMA[N_TAGS]; vector[N_TAGS] LENGTH_AT_RELEASE; real DAYS[N_TAGS]; vector[N_QMAS] QTEMP; vector[N_SITES] TEMP; vector[N_SITES] FAST; vector[N_SITES] AVG; int SITE_QMA[N_SITES]; real OBS_INCR[N_TAGS]; // Growth / tag-recapture priors real growth_k_mean_PRIOR_LOG_MEAN; real growth_k_mean_PRIOR_LOG_SD; real growth_k_site_sd_PRIOR_SD; real growth_k_qma_sd_PRIOR_SD; real growth_alpha_PRIOR_MEAN; real growth_alpha_PRIOR_SD; real growth_a_site_sd_PRIOR_SD; real growth_a_qma_sd_PRIOR_SD; real growth_a_ind_sd_PRIOR_SD; real growth_obs_error_sd_PRIOR; } transformed data { real rel_tol = 1e-5; real abs_tol = 1e-5; int max_num_steps = 10000000; } parameters{ // maturation real mat_scaler; real mat_p50; real mat_p95; real mat_site_sigma; //locoation variability sd // real mat_qma_sigma; //qma variability sd // vector[N_QMAS] mat_qma_fx; //qma variability // growth real growth_k_mean; real growth_alpha; cholesky_factor_corr[2] site_Omega; matrix[2,N_SITES] unscaled_site_fx; vector[2] site_scales; cholesky_factor_corr[2] qma_Omega; matrix[2,N_QMAS] unscaled_qma_fx; vector[2] qma_scales; //real growth_k_ind_sd; //vector[N_TAGS] growth_k_ind_fx; real growth_a_ind_sd; vector[N_TAGS] growth_a_ind_fx; real Ea; real growth_obs_error_sd; real expo; real stratum_avg; real stratum_fast; } transformed parameters{ // Maturation vector[N_MAT] Pred_Mat_Prop_Mat; // predicted maturity proportion for mat data vector[N_TAGS] Pred_Tag_Alloc; // predicted maturity investment for tags vector[N_TAGS] log_growth_a; //predicted increments for tags vector[N_TAGS] log_growth_k; //predicted increments for tags matrix[N_SITES, 2] site_fx; matrix[N_QMAS, 2] qma_fx; site_fx = (diag_pre_multiply(site_scales,site_Omega) * unscaled_site_fx)'; qma_fx = (diag_pre_multiply(qma_scales,qma_Omega) * unscaled_qma_fx)'; log_growth_k = growth_k_mean + site_fx[GROWTH_SITE,1] + qma_fx[GROWTH_QMA,1] + Ea*TEMP[GROWTH_SITE]; log_growth_a = growth_alpha + site_fx[GROWTH_SITE,2] + growth_a_ind_fx + qma_fx[GROWTH_QMA,2] + stratum_avg * AVG[GROWTH_SITE] + stratum_fast * FAST[GROWTH_SITE]; Pred_Tag_Alloc = mat_p50 - (mat_site_sigma*unscaled_site_fx[1,GROWTH_SITE])'; Pred_Mat_Prop_Mat = log(19)*(MAT_SIZE - (mat_p50 - (mat_site_sigma*unscaled_site_fx[1,MAT_SITE])'))/ mat_p95; } model{ mat_p50 ~ normal(mat_p50_PRIOR_MEAN, mat_p50_PRIOR_SD); mat_p95 ~ normal(mat_p95_PRIOR_MEAN, mat_p95_PRIOR_SD); mat_scaler ~ beta(mat_scaler_PRIOR_A, mat_scaler_PRIOR_B); mat_site_sigma ~ normal(5,mat_site_sigma_PRIOR_SD); OBS_MAT ~ binomial_logit(N_MAT_OBS, Pred_Mat_Prop_Mat); // growth params growth_k_mean ~ normal(growth_k_mean_PRIOR_LOG_MEAN, growth_k_mean_PRIOR_LOG_SD); growth_alpha ~ normal(growth_alpha_PRIOR_MEAN, growth_alpha_PRIOR_SD); to_vector(unscaled_qma_fx) ~ std_normal(); qma_Omega ~ lkj_corr_cholesky(1); qma_scales[1] ~ normal(0, growth_k_qma_sd_PRIOR_SD); qma_scales[2] ~ normal(0, growth_a_qma_sd_PRIOR_SD); to_vector(unscaled_site_fx) ~ std_normal(); site_Omega ~ lkj_corr_cholesky(1); site_scales[1] ~ normal(0, growth_k_site_sd_PRIOR_SD); site_scales[2] ~ normal(0, growth_a_site_sd_PRIOR_SD); Ea ~ normal(0, 5); stratum_fast ~ normal(0, 10); stratum_avg ~ normal(0, 10); growth_a_ind_fx ~ normal(0, growth_a_ind_sd); growth_a_ind_sd ~ normal(0, growth_a_ind_sd_PRIOR_SD); growth_obs_error_sd ~ normal(0, growth_obs_error_sd_PRIOR); //growth_method ~ normal(0, 2); expo ~ normal(10,expo_PRIOR); // growth observation model target += reduce_sum(tags_sum_lupdf, OBS_INCR, grainsize, growth_obs_error_sd, mat_scaler, mat_p95, expo, log_growth_a, log_growth_k, Pred_Tag_Alloc, LENGTH_AT_RELEASE, DAYS); } generated quantities{ vector[2] mg; vector[5] theta; // [1]anabolism [2] catabolsim [3] lmat [4] investment vector[L] pred_mean_growth; vector[L] pred_sd_growth; real pred_means_growth[L,N_QMAS,3]; real pred_sds_growth[L,N_QMAS,3]; real pred_zeros_growth[L,N_QMAS,3]; vector[L] p_mat; vector[N_TAGS] TAG_LL; real TAG_LL_sum; vector[N_MAT] MAT_LL; real MAT_LL_sum; matrix[L,N_SITES] PI; matrix[L,N_SITES] PM; matrix[L,N_SITES] PA; real PIS[N_INT,L,N_QMAS,3]; real PMS[N_INT,L,N_QMAS,3]; real PAS[N_INT,L,N_QMAS,3]; real growth_ind_fx_pred; vector[3] gtypes; real DD[1]; vector[2] LS; DD[1] = 1; LS[1] = 75; LS[2] = 125; gtypes = rep_vector(0,3); gtypes[1] = 0; gtypes[2] = stratum_avg; gtypes[3] = stratum_fast; theta[1] = exp(growth_alpha); theta[2] = exp(growth_k_mean); theta[3] = mat_p50; theta[4] = mat_scaler; theta[5] = mat_p95/mat_scaler^expo; for (i in 1:2){ vector[1] mgt[1] = ode_rk45_tol(growth, LS[i:i], 0.0, DD, 1e-5, 1e-5, 10000000, theta); mg[i] = mgt[1,1]-LS[i]; } for (i in 1:L){ for (s in 1:N_SITES){ vector[5] eta; // [1]anabolism [2] catabolsim [3] lmat [4] investment growth_ind_fx_pred = normal_rng(0,growth_a_ind_sd); eta[1] = exp(growth_alpha + site_fx[s,2] + qma_fx[SITE_QMA[s],2]+ growth_ind_fx_pred + stratum_avg * AVG[s] + stratum_fast * FAST[s]); eta[2] = exp(growth_k_mean + site_fx[s,1]+ qma_fx[SITE_QMA[s],1] + Ea* TEMP[s]); eta[3] = mat_p50 + mat_site_sigma*unscaled_site_fx[s,1]; eta[4] = mat_scaler; eta[5] = mat_p95/mat_scaler^expo; vector[1] mgt[1] = ode_rk45_tol(growth, LENGTHS[i:i], 0.0, DD, 1e-6, 1e-6, 100000, eta); PI[i,s] = mgt[1,1]-LENGTHS[i]; PA[i,s] = inv_logit(log(19)*(LENGTHS[i] - theta[3]/theta[4]) / theta[5]); PM[i,s] = inv_logit(log(19)*(LENGTHS[i] - (mat_p50 + (mat_site_sigma*unscaled_site_fx[s,2])))/ mat_p95); } p_mat[i] = log(19)*(LENGTHS[i] - mat_p50) / mat_p95; pred_mean_growth[i] = mean(PI[i,]); pred_sd_growth[i] = sd(PI[i,]); } for (n in 1:N_INT){ real growth_ind_fx_ipred = normal_rng(0,growth_a_ind_sd); vector[2] unscaled_site_fx_pred; vector[2] site_fx_pred; unscaled_site_fx_pred[1] = normal_rng(0,1); unscaled_site_fx_pred[2] = normal_rng(0,1); site_fx_pred = site_scales .* unscaled_site_fx_pred; for (i in 1:L){ for (s in 1:N_QMAS){ for(t in 1:3){ vector[5] eta; // [1]anabolism [2] catabolsim [3] lmat [4] investment eta[1] = exp(growth_alpha + site_fx_pred[2] + qma_fx[s,2]+ growth_ind_fx_pred + gtypes[t]); eta[2] = exp(growth_k_mean + site_fx_pred[1]+ qma_fx[s,1] + Ea*QTEMP[s]); eta[3] = mat_p50 + mat_site_sigma*unscaled_site_fx_pred[2]; eta[4] = mat_scaler; eta[5] = mat_p95/mat_scaler^expo; vector[1] mgt[1] = ode_rk45_tol(growth, LENGTHS[i:i], 0.0, DD, 1e-6, 1e-6, 100000, eta); PIS[n,i,s,t] = mgt[1,1]-LENGTHS[i]; PAS[n,i,s,t] = inv_logit(log(19)*(LENGTHS[i] - theta[3]/theta[4]) / theta[5]); PMS[n,i,s,t] = inv_logit(log(19)*(LENGTHS[i] - (mat_p50 + (mat_site_sigma*unscaled_site_fx_pred[2])))/ mat_p95); } } } } for (i in 1:L){ for (s in 1:N_QMAS){ for(t in 1:3){ real SS=0; real VV=0; real PP=0; SS+= PIS[1,i,s,t]>0 ? log(PIS[1,i,s,t]) : 0; PP+= PIS[1,i,s,t]>0 ? 0 : 1; for (n in 2:N_INT){ real SSU = PIS[n,i,s,t]>0 ? SS + (log(PIS[n,i,s,t])-SS)/n : SS; VV+= PIS[n,i,s,t]>0 ? (log(PIS[n,i,s,t])-SS)*(log(PIS[n,i,s,t])-SSU) : 0; SS = SSU; PP+= PIS[n,i,s,t]>0 ? 0 : 1; } pred_means_growth[i,s,t] = SS; pred_zeros_growth[i,s,t] = PP/N_INT; pred_sds_growth[i,s,t] = sqrt(VV/(N_INT-1)); } } } // for (i in 1:N_MAT) MAT_LL[i] = binomial_logit_lpmf(OBS_MAT[i] | N_MAT_OBS[i], Pred_Mat_Prop_Mat[i]); //MAT_LL_sum = sum(MAT_LL); //for (i in 1:N_TAGS) TAG_LL[i] = cauchy_lpdf(OBS_INCR[i] | Pred_Tag_Inc[i], growth_obs_error_sd); //TAG_LL_sum = sum(TAG_LL); }