#This version allows the offspring env to already effect mode in the offspring generation. data { int n; // number of data points real growth_rate[n]; // growth rate per generation int offspring_env_val[n]; // the treatment, binary 0 is norm, 1 is anoxia int maternal_env_val[n]; // the maternal treatment, binary 0 is norm, 1 is anoxia int maternal_cue_val[n]; // Was a cue given to the mother? binary 0 is no, 1 is yes int location_val[n]; // Where was the experiment done, 0 is Lisbon, 1 is Paris int generation[n]; real lba; real lbn; real upa; real upn; } parameters { real del_loc; real del_cue; real lambda_A; #rate of carryover real lambda_N; real carry_eff; # the magnitude of the carryover real s0; real mode_A; real mode_N; real sd_N; real sd_A; real err[n]; # individual error effects per observation real sigma_err; #variance parameter for error term } transformed parameters{ real rate[n]; real shape[n]; real mode[n]; real sd_t[n]; real s[n]; #The current env carryover effect s[1] = (s0 +err[1]) ; for(i in 2:n){ s[i] = generation[i]-3 <0.1 ? (s0 + err[i]) :( (s[i-1] * (1-lambda_N) + carry_eff * lambda_N ) * (1 - offspring_env_val[i]) + (s[i-1] * (1-lambda_A) - carry_eff * lambda_A ) * (offspring_env_val[i]) + err[i]) ; #s0 if gen is 2, otherwise carryover in an open-ended way } for(i in 1:n){ mode[i] = exp( log( mode_N * (1 - offspring_env_val[i]) + mode_A * ( offspring_env_val[i]) )+ del_cue*(maternal_cue_val[i]) - del_cue*(1-maternal_cue_val[i]) + del_loc * location_val[i]- del_loc * (1-location_val[i])+ s[i] ); sd_t[i] = sd_N * (1-offspring_env_val[i]) + sd_A * (offspring_env_val[i]); rate[i]= (mode[i] + sqrt( mode[i]^2 + 4*sd_t[i]^2 ) )/( 2 * sd_t[i]^2 ); shape[i]= 1 + mode[i] * rate[i]; } } model { #priors del_loc ~normal(0,1); sd_N ~ cauchy(0,2.5); sd_A ~ cauchy(0,2.5); carry_eff ~ normal(0,1); s0 ~ normal(0,1); err ~ normal(0,sigma_err); # growth_rate ~ gamma(shape,rate); # version that keeps full log likelihood for comparison.... target += gamma_lpdf( growth_rate | shape, rate); # growth_rate ~ gamma(shape,rate); } generated quantities{ real log_lik[n]; real simulated_data[n]; real simulated_err[n]; real simulated_data2[n]; #shadow parameters for the simulation real sim_rate[n]; real sim_shape[n]; real sim_mode[n]; real sim_sd_t[n]; real sim_s[n]; #The current env carryover effect #even deeper shadow parameters to show the effects of multiple generations of anoxia real multi_mode_A[100]; #track mode real multi_mode_N[100]; #track mode real multi_mode_S[100]; #track mode for (i in 1:n) { log_lik[i] = gamma_lpdf(growth_rate[i] | shape[i], rate[i] ) ; simulated_data[i] = gamma_rng(shape[i],rate[i]); simulated_err[i] = normal_rng(0,sigma_err); } sim_s[1] = (s0 +simulated_err[1]) ; for(i in 2:n){ sim_s[i] = generation[i]-3 <0.1 ? (s0 + simulated_err[i]) :( (sim_s[i-1] * (1-lambda_N) + carry_eff * lambda_N ) * (1 - offspring_env_val[i]) + (sim_s[i-1] * (1-lambda_A) - carry_eff * lambda_A ) * (offspring_env_val[i]) + simulated_err[i]) ; #s0 if gen is 2, otherwise carryover in an open-ended way } for(i in 1:n){ sim_mode[i] = exp( log( mode_N * (1 - offspring_env_val[i]) + mode_A * ( offspring_env_val[i]) )+ del_loc * location_val[i]- del_loc * (1-location_val[i])+ sim_s[i] ); sim_sd_t[i] = sd_N * (1-offspring_env_val[i]) + sd_A * (offspring_env_val[i]); sim_rate[i]= (sim_mode[i] + sqrt( sim_mode[i]^2 + 4*sim_sd_t[i]^2 ) )/( 2 * sim_sd_t[i]^2 ); sim_shape[i]= 1 + sim_mode[i] * sim_rate[i]; } for (i in 1:n) { simulated_data2[i] = gamma_rng(sim_shape[i],sim_rate[i]); } #Forward simulation of the way that the carryover decays with time sim_s[1] = (0 ) ; for(i in 2:10){ sim_s[i] = (sim_s[i-1] * (1-lambda_A) - carry_eff * lambda_A ) ; #do some generations of anoxia } for(i in 11:20){ sim_s[i] = (sim_s[i-1] * (1-lambda_N) + carry_eff * lambda_N ) ; #do some generations of normoxia } for(i in 21:25){ sim_s[i] = (sim_s[i-1] * (1-lambda_A) - carry_eff * lambda_A ) ; #do some generations of anoxia } for(i in 26:30){ sim_s[i] = (sim_s[i-1] * (1-lambda_N) + carry_eff * lambda_N ) ; #do some generations of normoxia } for(i in 31:33){ sim_s[i] = (sim_s[i-1] * (1-lambda_A) - carry_eff * lambda_A ) ; #do some generations of anoxia } for(i in 34:36){ sim_s[i] = (sim_s[i-1] * (1-lambda_N) + carry_eff * lambda_N ) ; #do some generations of normoxia } for(i in 37:38){ sim_s[i] = (sim_s[i-1] * (1-lambda_A) - carry_eff * lambda_A ) ; #do some generations of anoxia } for(i in 39:40){ sim_s[i] = (sim_s[i-1] * (1-lambda_N) + carry_eff * lambda_N ) ; #do some generations of normoxia } for(i in 41:41){ sim_s[i] = (sim_s[i-1] * (1-lambda_A) - carry_eff * lambda_A ) ; #do some generations of anoxia } for(i in 42:42){ sim_s[i] = (sim_s[i-1] * (1-lambda_N) + carry_eff * lambda_N ) ; #do some generations of normoxia } for(i in 43:43){ sim_s[i] = (sim_s[i-1] * (1-lambda_A) - carry_eff * lambda_A ) ; #do some generations of anoxia } for(i in 44:44){ sim_s[i] = (sim_s[i-1] * (1-lambda_N) + carry_eff * lambda_N ) ; #do some generations of normoxia } for(i in 45:45){ sim_s[i] = (sim_s[i-1] * (1-lambda_A) - carry_eff * lambda_A ) ; #do some generations of anoxia } for(i in 46:46){ sim_s[i] = (sim_s[i-1] * (1-lambda_N) + carry_eff * lambda_N ) ; #do some generations of normoxia } for(i in 47:47){ sim_s[i] = (sim_s[i-1] * (1-lambda_A) - carry_eff * lambda_A ) ; #do some generations of anoxia } for(i in 48:48){ sim_s[i] = (sim_s[i-1] * (1-lambda_N) + carry_eff * lambda_N ) ; #do some generations of normoxia } for(i in 49:49){ sim_s[i] = (sim_s[i-1] * (1-lambda_A) - carry_eff * lambda_A ) ; #do some generations of anoxia } for(i in 50:50){ sim_s[i] = (sim_s[i-1] * (1-lambda_N) + carry_eff * lambda_N ) ; #do some generations of normoxia } for(i in 51:75){ sim_s[i] = (sim_s[i-1] * (1-lambda_A) - carry_eff * lambda_A ) ; #do some generations of anoxia } for(i in 76:100){ sim_s[i] = (sim_s[i-1] * (1-lambda_N) + carry_eff * lambda_N ) ; #do some generations of normoxia } for(i in 1:100){ multi_mode_A[i] = exp( log( mode_A )+ sim_s[i] ); multi_mode_N[i] = exp( log( mode_N )+ sim_s[i] ); multi_mode_S[i] = exp( sim_s[i] ); } }