Hello,
I am very new to rstan after having used rjags for quite a while. I find myself dealing with a fairly complicated model which I attach here below, with its own dataset.
I have essentially two questions.
-
It looks like that despite having set all the possible boundaries to the parameters and the transformed parameters, I don’t seem to be able to get a prior draw for the quantity yy[1,2,1] (on line 103), which should be distributed as a beta with parameters taken from the data (line 179).
Can anybody explain me why this happens and how to solve the issue? -
You see I have many variables distributed as uniform random variables with extremes taken in turn from other variables (eg code line 181): how do I set the limits for these? Is it possible to write them like I have in my code?
Thank you very much in advance for all the help, and sorry for the naive questions
Best,
EM
data: Dropbox - File Deleted - Simplify your life
model:
data{ // prevalence, standard dev., case rep., populations real<lower=0> P_hat[8,3,2]; real<lower=0> s_over_rootn[8,3,2]; int<lower=0> R_by_year[18,3,2]; int<lower=0> N_by_year[18,3,2]; // n.years(pairs), n.years(singles), n.age groups, n.genders int<lower=0> nyrs; int<lower=0> n_single_years; int<lower=0> nages; int<lower=0> ngenders; // parameters for beta distributions of: // 1st point bezier, beta, prob.detection, treated given symptomatic and durations vector<lower=0>[2] pbetaBz; vector<lower=0>[2] pbetaF; vector<lower=0>[2] pbetaM; vector<lower=0>[2] pm; vector<lower=0>[2] pvs; vector<lower=0>[2] pdt; vector<lower=0>[2] pdu; // to build duration priors (from confidence intervals in literature) real<lower=0> du_int; real<lower=0, upper=1> du_slope; real<lower=0, upper=1> dt_int; real<lower=0, upper=1> dt_slope; real<lower=0, upper=1> m_int; real<lower=0, upper=1> m_slope; // parameters for priors coefficients cubic polyn matrix[3,2] a1_mean; real<lower=0> a1_sigma; matrix[3,2] a2_mean; real<lower=0> a2_sigma; matrix[3,2] a3_mean; real<lower=0> a3_sigma; matrix[3,2] a4_mean; real<lower=0> a4_sigma; // final durations vector[2] dmultN_parms; vector[2] dmultB_parms; // n.predicted years, beta.decrease int<lower=0> pred; real<lower=0> less_beta; } transformed data{ real<lower=0, upper=n_single_years+pred> yp[n_single_years+pred]; real<lower=0> yp1; // definitions for(w in 1:(n_single_years+pred)) yp[w]= n_single_years+pred-w; yp1 = (n_single_years+pred-1); } parameters{ // matrices real a[4,nages,ngenders]; //coeff.cubic polyn.for prevalence // vectors vector<lower=1, upper=2>[ngenders] yy3fac; vector<lower=1, upper=2>[ngenders] factor_y3; vector<lower=1, upper=2>[ngenders] yy3fac_35; vector<lower=1, upper=2>[ngenders] yy4fac; vector<lower=1, upper=2>[ngenders] factor_y4; vector<lower=1, upper=2>[ngenders] yy4fac_35; vector<lower=dmultN_parms[1], upper=dmultN_parms[2]>[ngenders] d_multN; // numbers real<lower=0, upper=1> psy; real<lower=0, upper=1> vs; real<lower=0, upper=1> d_unt_1; real<lower=0, upper=1> d_unt_2; real<lower=0, upper=1> d_ts_1; real<lower=0, upper=1> d_ts_2; real<lower=dmultB_parms[1], upper=dmultB_parms[2]> d_multB; } transformed parameters { real<lower=0, upper=2> yy[4,nages,ngenders]; //support points for Bezier curve real P[4,nages,ngenders]; real bezcurve[(n_single_years+pred),nages,ngenders]; real<lower=0, upper=1> bez[(n_single_years+pred),nages,ngenders]; //values of the Bezier curve real<lower=0, upper=1> gam[(n_single_years+pred),nages,ngenders]; //prob.detection (fraction treated) real<lower=0, upper=1> theta[(n_single_years+pred),nages,ngenders]; //rate case reports real<lower=0, upper=1> beta[(n_single_years+pred),nages,ngenders]; //prob.case is correctly reported in recording system real<lower=0> dur[(n_single_years+pred),nages,ngenders]; //overall duration real<lower=0, upper=1> eta[(n_single_years+pred),nages,ngenders]; //incidence real<lower=0, upper=1> pi_w[(n_single_years+pred),nages,ngenders]; //prevalence_w real<lower=0, upper=1> pi_j[nyrs,nages,ngenders]; //prevalence_j real<lower=0, upper=1> beta_uF[(n_single_years+pred),nages]; real<lower=0, upper=1> beta_uM[(n_single_years+pred),nages]; real<lower=0, upper=1> m; //prob.symptomatic vector<lower=0>[ngenders] dnew; //dur.treated, asymptomatic vector<lower=0>[ngenders] du; //dur.untreated vector<lower=0>[ngenders] dt; //dur.treated, symptomatic // ASSIGMENTS // // 1. transform Bezier support points to comply with constraints // females & males together yy[1,2,2] = d_multB*yy[1,2,1]; for(g in 1:ngenders){ yy[3,1,g] = yy[1,1,g]*yy3fac[g]; yy[3,2,g] = yy[1,2,g]*factor_y3[g]; yy[3,3,g] = yy[1,3,g]*yy3fac_35[g]; yy[4,1,g] = yy[1,1,g]*yy4fac[g]; yy[4,2,g] = yy[1,2,g]*factor_y4[g]; yy[4,3,g] = yy[1,3,g]*yy4fac_35[g]; // 2. calculate Ps from yys (Bezier support points) for(k in 1:nages){ P[1,k,g] = yy[1,k,g]; P[2,k,g] = (-5*yy[1,k,g]+18*yy[2,k,g]- 9*yy[3,k,g]+2*yy[4,k,g])/6; P[3,k,g] = ( 2*yy[1,k,g]- 9*yy[2,k,g]+18*yy[3,k,g]-5*yy[4,k,g])/6; P[4,k,g] = yy[4,k,g]; // 3. Bezier curve (to be truncated at 0 and 1) for(w in 1:(n_single_years+pred)){ bezcurve[w,k,g] = pow( yp[w]/yp1, 3)*P[1,k,g]+ 3*pow(yp[w]/yp1, 2)*((w-1)/yp1)*P[2,k,g]+ 3*(yp[w]/yp1)*pow((w-1)/yp1,2)*P[3,k,g]+ pow((w-1)/yp1,3)*P[4,k,g]; bez[w,k,g] = fmax(0, fmin(1, bezcurve[w,k,g]) ); // 4. detection gam[w,k,g] = m*vs+(1-m)*bez[w,k,g]; // 5. overall duration dur[w,k,g] = m*vs*dt[g] + m*(1-vs)*du[g] + (1-m)*bez[w,k,g]*dnew[g] + (1-m)*(1-bez[w,k,g])*du[g]; // 6. prevalence pi_w[w,k,g] = inv_logit( a[1,k,g]+a[2,k,g]*w+a[3,k,g]*pow(w,2)+a[4,k,g]*pow(w,3) ); // 7. incidence eta[w,k,g] = pi_w[w,k,g]/dur[w,k,g]; // 8. rate of case reports theta[w,k,g]= eta[w,k,g]*beta[w,k,g]*gam[w,k,g]; } // closes w for(j in 1:nyrs){ pi_j[j,k,g]=(pi_w[(2*j-1),k,g]+pi_w[(2*j),k,g])/2; } } // closes k // 9. duration treated, asymptomatic (as a transformation of dt, du, dmultN) dnew[g] = dt[g]+(d_multN[g])*(du[g]-dt[g]); } //closes g //10. prob.case correctly recorded in reporting system (random walk) for(k in 1:nages){ beta_uF[1,k] = 0; beta_uM[1,k] = 0; for(w in 2:(n_single_years+pred)){ beta[w,k,1] = (0.95+less_beta-beta[w-1,k,1])*beta_uF[w,k]+(beta[w-1,k,1]-less_beta); beta[w,k,2] = (0.95+less_beta-beta[w-1,k,2])*beta_uM[w,k]+(beta[w-1,k,2]-less_beta); }} //11. durations(untreated and treated, symptomatic) du[1] = du_int+du_slope*d_unt_1; du[2] = du_int+du_slope*d_unt_2; dt[1] = dt_int+dt_slope*d_ts_1; dt[2] = dt_int+dt_slope*d_ts_2; //12. prob symptomatic m = m_int+m_slope*psy; print(yy[1,2,1]); print(pbetaBz); } model { //// PRIORS //// ///Bezier points females /// /// age group 2 yy[1,2,1] ~ beta(pbetaBz[1], pbetaBz[2]); factor_y3[1]~ uniform(1,2); yy[2,2,1] ~ uniform(yy[1,2,1], yy[3,2,1]); factor_y4[1]~ uniform(1,2); /// age group 1 yy[1,1,1] ~ uniform(0*yy[1,2,1], 1*yy[1,2,1]); yy3fac[1] ~ uniform(1, factor_y3[1]); yy4fac[1] ~ uniform(1, factor_y4[1]); yy[2,1,1] ~ uniform(yy[1,1,1], yy[3,1,1]); /// age group 3 yy[1,3,1] ~ uniform(0*yy[1,1,1], 1*yy[1,1,1]); yy3fac_35[1]~ uniform(1, yy3fac[1]); yy4fac_35[1]~ uniform(1, yy4fac[1]); yy[2,3,1] ~ uniform(yy[1,3,1], yy[3,3,1]); ///Bezier points males /// age group 2 factor_y3[2]~ uniform(1,2); yy[2,2,2] ~ uniform(yy[1,2,2], yy[3,2,2]); factor_y4[2]~ uniform(1,2); /// age group 1 yy[1,1,2] ~ uniform(0*yy[1,2,2], 1*yy[1,2,2]); yy3fac[2] ~ uniform(1, factor_y3[2]); yy4fac[2] ~ uniform(1, factor_y4[2]); yy[2,1,2] ~ uniform(yy[1,1,2], yy[3,1,2]); /// age group 3 yy[1,3,2] ~ uniform(0*yy[1,1,2], 1*yy[1,1,2]); yy3fac_35[2]~ uniform(1, yy3fac[2]); yy4fac_35[2]~ uniform(1, yy4fac[2]); yy[2,3,2] ~ uniform(yy[1,3,2], yy[3,3,2]); /////////////////////////////////////////////////// psy ~ beta(pm[1],pm[2]); vs ~ beta(pvs[1], pvs[2]); /////////////////////////////////////////////////// /// cubic coefficients for(g in 1:ngenders){for(k in 1:nages){ // age-group and gender-specific a[1,k,g]~ normal(a1_mean[k,g], a1_sigma); //dnorm(coef[1,k,g],1) #dnorm(0,1) #dnorm(-2,1)# a[2,k,g]~ normal(a2_mean[k,g], a2_sigma); //dnorm(coef[2,k,g],sqrt(0.05)) #dnorm(0,1) #dnorm(-0.05, sqrt(0.05)) a[3,k,g]~ normal(a3_mean[k,g], a3_sigma); //dnorm(coef[3,k,g],sqrt(0.01)) #dnorm(0,1) #dnorm(0, sqrt(0.01)) a[4,k,g]~ normal(a4_mean[k,g], a4_sigma); //dnorm(coef[4,k,g],sqrt(0.0001)) #dnorm(0,1) #dnorm(0, sqrt(0.0001)) }} /////////////////////////////////////////////////// /// reporting for(k in 1:nages){ beta[1,k,1] ~ beta(pbetaF[1],pbetaF[2]); beta[1,k,2] ~ beta(pbetaM[1],pbetaM[2]); for(w in 2:(n_single_years+pred)){ beta_uF[w,k] ~ uniform(0,1); beta_uM[w,k] ~ uniform(0,1); }} /////////////////////////////////////////////////// /// durations /// untreated d_unt_1 ~ beta(pdu[1],pdu[2]); d_unt_2 ~ beta(pdu[1],pdu[2]); /// treated symptomatic d_ts_1 ~ beta(pdt[1],pdt[2]); d_ts_2 ~ beta(pdt[1],pdt[2]); /// multipliers for duration treated, asymptomatic d_multB ~ uniform(dmultB_parms[1],dmultB_parms[2]); d_multN[1] ~ uniform(dmultN_parms[1],dmultN_parms[2]); d_multN[2] ~ uniform(dmultN_parms[1],dmultN_parms[2]); /// likelihoods for(g in 1:ngenders){ for(k in 1:nages){ for(j in 1:nyrs){ P_hat[j,k,g] ~ normal( pi_j[j,k,g], s_over_rootn[j,k,g] );} for(w in 1:n_single_years){ R_by_year[w,k,g] ~ binomial(N_by_year[w,k,g], theta[w,k,g]);} }} }