#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 in anoxia real lambda_N; #rate of carryover in normoxia 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); }