Absent prior draws?

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.

  1. 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?

  2. 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]);}
	}}

}

Just in general, I think you are carrying over some JAGS assumptions that don’t hold with Stan. There is an appendix to the Stan User Manual that covers a lot of these.

(1) In line 103, yy[1,2,2] = d_multB*yy[1,2,1];, the yy[1,2,1] value has not been defined yet and thus yy[1,2,2] is assigned NaN (not a number). This is also the first question in the FAQ. I’m not sure what you were trying to do here. Maybe the first point in the Bezier curve was supposed to be taken from the data or a parameter?

(2) In general, the bounds on one parameter can be (a function of) another parameter that is declared earlier. However, in your case the bounds seem to vary in which case, that won’t work because bounds in Stan can only be scalars. Thus, in your case what you should be doing is declaring a primitive parameter with a lower bound of 0 and an upper bound of 1, and then shift and scale it to the correct interval (i.e. if z = a + b * u where u is standard uniform, then z is uniform between a and b + a. It would look something like

parameters {
  real<lower=0, upper=1> std_u[4,nages,ngenders];
  ...
}
transformed parameters {
  real<lower=0, upper=2> yy[4,nages,ngenders];
  // map elements of std_u to those of yy
}

If you are using the uniform density in Stan, you are doing it wrong and missing what makes the MCMC algorithm Stan uses fundamentally different from a Gibbs sampler. In Stan, you define a posterior kernel (in log-units), which characterizes the probability that each leapfrog step is chosen as the next posterior realization. Since the uniform density is equal to a constant (within its support), you are just adding the logarithm of that constant to the posterior kernel function, which cannot have any effect on which leapfrog proposals are realized. Instead, you need to be using transformations like those above to ensure that the elements of yy fall in the appropriate intervals.

Thank you so much Ben, this explanation has been very useful and has allowed me to reshape my model.

I have another question now, which I can’t wrap my head around:
I have this matrix of parameters on line 91 that should be probabilities (and go into the binomial likelihood, on line 216), so defined (hopefully correctly) as being between 0 and 1.
Turns out that they are not, from what I see from preliminary runs and prints of the chain.
And they are not, also because another parameter (on line 90) turns out to be >1, which it also shouldn’t (being, again, a percent).
The relations between the two are on line 160 and 161, and I can’t understand why this happens.

On top of that, a log probability equals to log(0) shows up, but that’s probably related.

Again apologies for any naive questions, thank you so much for the help you’ve given me and for any future help.
Best, EM

data: https://www.dropbox.com/s/8vmqeuyt1jlhi8d/datastan.rda?dl=0
possible initial values for chain: https://www.dropbox.com/s/xr3qw8og1lsatv6/initstan.rda?dl=0

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{
// indices (defined as real) to avoid integer division in Bezier definition
	real<lower=0, upper=n_single_years+pred> yp[n_single_years+pred];
	real<lower=0> yp1; 
	// assignments
	for(w in 1:(n_single_years+pred)) yp[w]= n_single_years+pred-w;
	yp1 = (n_single_years+pred-1);
}

parameters{
// numbers
	real<lower=0, upper=1> 		psy;								//multiplier for prob sympt
	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; 	//multiplier for Bezier males vs. females
//vectors	
	vector<lower=dmultN_parms[1], upper=dmultN_parms[2]>[ngenders] d_multN; //multiplier for duration trt asympt
// matrices
	real 					a[4,nages,ngenders]; 							//coeff cubic polyn for prevalence
	real<lower=0, upper=1> 	base_bzp[4,nages,ngenders]; 					// bases for unif distrib on Bezier points
	real<lower=0, upper=1> 	base_beta[(n_single_years+pred),nages,ngenders];// bases for unif distrib on reporting
	} 

transformed parameters {
// numbers
	real<lower=0, upper=1> 		m;										 		//prob.symptomatic
// vectors	
	vector<lower=0>[ngenders] 	dnew; 									 		//dur.treated, asymptomatic
	vector<lower=0>[ngenders] 	du;										 		//dur.untreated
	vector<lower=0>[ngenders] 	dt;										 		//dur.treated, symptomatic				
// matrices
	real<lower=0, upper=2> 		yy[4,nages,ngenders]; 						 	//support points for Bezier curves (transformed from base_bzp)
	real<lower=0, upper=2> 		P[4,nages,ngenders];							//transformed support points for Bezier (transformed from yy)
	real<lower=0> 				dur[(n_single_years+pred),nages,ngenders]; 		//overall duration
	real<lower=0, upper=1> 		pi_w[(n_single_years+pred),nages,ngenders];	 	//prevalence_w (transf of the coefficients a (cubic polyn))
	real<lower=0, upper=1> 		pi_j[nyrs,nages,ngenders];					 	//prevalence_j (prevalence on the paired years)
	real 						bezcurve[(n_single_years+pred),nages,ngenders]; //uncapped Bezier (from Ps)
	real<lower=0, upper=1> 		bez[(n_single_years+pred),nages,ngenders]; 	 	//capped (at 0 and 1) Bezier curve  (from uncapped Bezier)	
	real<lower=0, upper=1> 		gam[(n_single_years+pred),nages,ngenders]; 	 	//prob.detection (fraction treated)
	real<lower=0, upper=0.95> 	beta[(n_single_years+pred),nages,ngenders];		//prob.case is correctly reported in recording system (transf from base_beta)
	real<lower=0, upper=2> 		eta[(n_single_years+pred),nages,ngenders]; 	 	//incidence (transf from prevalence (prevalence/duration))
	real<lower=0, upper=1> 		theta[(n_single_years+pred),nages,ngenders];  	//rate case reports (from eta, beta, gamma)

// ASSIGMENTS //
// 1. map elements of base_bzp into those of yy
// these allow me to only put uniform(0,1) priors, and the transformations
// take care of all the relationships-constraints
	yy[1,2,1] = base_bzp[1,2,1];
	yy[1,2,2] = base_bzp[1,2,1]*(dmultB_parms[1]+(dmultB_parms[2]-dmultB_parms[1])*base_bzp[1,2,2]);
for(g in 1:ngenders){
	yy[3,2,g] = yy[1,2,g]+yy[1,2,g]*base_bzp[3,2,g];
	yy[2,2,g] = yy[1,2,g]+(yy[3,2,g]-yy[1,2,g])*base_bzp[2,2,g];
	yy[4,2,g] = yy[1,2,g]+yy[1,2,g]*base_bzp[4,2,g];
	//
	yy[1,1,g] = yy[1,2,g]*base_bzp[1,1,g];
	yy[3,1,g] = yy[1,1,g]*(1+pow(base_bzp[3,1,g],2));
	yy[2,1,g] = yy[1,1,g]+(yy[3,1,g]-yy[1,1,g])*base_bzp[2,1,g];
	yy[4,1,g] = yy[1,1,g]*(1+pow(base_bzp[4,1,g],2));
	//
	yy[1,3,g] = yy[1,1,g]*base_bzp[1,3,g];
	yy[3,3,g] = yy[1,3,g]*(1+pow(base_bzp[3,3,g],3));
	yy[2,3,g] = yy[1,3,g]+(yy[3,3,g]-yy[1,3,g])*base_bzp[2,3,g];
	yy[4,3,g] = yy[1,3,g]*(1+pow(base_bzp[4,3,g],3));

// 2. map elements of base_beta into beta
for(k in 1:nages){				   beta[1,k,g] = base_beta[1,k,g];
for(w in 2:(n_single_years+pred)){ beta[w,k,g] = beta[(w-1),k,g]-less_beta+(0.95-beta[(w-1),k,g]+less_beta)*base_beta[(w-1),k,g]; } //closes w
	
// 3. calculate Ps from yys (Bezier support points)
	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];

// 4. Bezier curve (to be truncated at 0 and 1 below)
for(w in 1:(n_single_years+pred)){					// opens w again
	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];

	// 4.1 truncation
	bez[w,k,g]  = fmax(0, fmin(1, bezcurve[w,k,g]) ); 
}}} //closes w, closes k, closes g opened above

// 5. prob symptomatic (from literature)
	m = m_int+m_slope*psy;

for(g in 1:ngenders){for(k in 1:nages){for(w in 1:(n_single_years+pred)){
// 6. detection
	gam[w,k,g]  = m*vs+(1-m)*bez[w,k,g];
}}}

// 7. durations(untreated and treated, symptomatic, from literature)
	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;

for(g in 1:ngenders){
// 8. duration treated, asymptomatic (as a transformation of dt, du, dmultN)
	dnew[g] = dt[g]+(d_multN[g])*(du[g]-dt[g]);

for(k in 1:nages){for(w in 1:(n_single_years+pred)){
// 9. 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];

// 10. prevalence (logit(prevalence)=cubic.polynomial => prevalence=inv_logit(cubic.polynomial))
	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) );
	for(j in 1:nyrs){ pi_j[j,k,g]=(pi_w[(2*j-1),k,g]+pi_w[(2*j),k,g])/2; }
	
// 11. incidence (prevalence/duration)
	eta[w,k,g]  = pi_w[w,k,g]/dur[w,k,g];

// 12. rate of case reports (true incidence*p(incident case is diagnosed)*p(diagnosed case captured in reporting system))
// 	   OR p(incident)*p(diagnosed|incident)*p(reported|incident and diagnosed)
	theta[w,k,g]= eta[w,k,g]*beta[w,k,g]*gam[w,k,g];
}}} // closes w, closes k, closes g

   print("eta: ",eta[,2,1]);
   print("theta: ",theta[,2,1]);

}

model {
// PRIORS //
	//Bezier points
	//females [,,1] first one is beta (for age group 2), the others uniform
	base_bzp[1,2,1] ~ beta(pbetaBz[1], pbetaBz[2]);
	base_bzp[1,1,1]	~ uniform(0,1);
	base_bzp[1,3,1]	~ uniform(0,1);
	//males [,,2]
	for(k in 1:nages){base_bzp[1,k,2]~ uniform(0,1);}
		
for(s in 2:4){for(k in 1:nages){for(g in 1:ngenders){ base_bzp[s,k,g] ~ uniform(0,1);}}}

	// Reporting
for(k in 1:nages){
	base_beta[1,k,1] ~ beta(pbetaF[1],pbetaF[2]);
	base_beta[1,k,2] ~ beta(pbetaM[1],pbetaM[2]);
for(g in 1:ngenders){ for(w in 2:(n_single_years+pred)){ base_beta[w,k,g] ~ uniform(0,1); }}
}	
	
	psy ~ beta(pm[1],pm[2]);
	vs 	~ beta(pvs[1], pvs[2]);

	// cubic polyn 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))
}}

// 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]);}
	}}

}

One thing you can (and should) do with Stan is define your own functions in the functions block that do anything non-trivial. You can then export those functions to R and verify that they are doing what they are supposed to do. In your case, it could start off like

functions {
  real[, , ] make_yy(real [, ,] base_bzp, vector dmultB_parms, 
                                  int nages, int ngenders) {
    real yy[4, nages, ngenders];
    yy[1,2,1] = base_bzp[1,2,1];
    yy[1,2,2] = base_bzp[1,2,1] * (dmultB_parms[1]+
                                  (dmultB_parms[2]-dmultB_parms[1])*base_bzp[1,2,2]);
    for (g in 1:ngenders){
	yy[3,2,g] = yy[1,2,g]+yy[1,2,g]*base_bzp[3,2,g];
	yy[2,2,g] = yy[1,2,g]+(yy[3,2,g]-yy[1,2,g])*base_bzp[2,2,g];
	yy[4,2,g] = yy[1,2,g]+yy[1,2,g]*base_bzp[4,2,g];
	yy[1,1,g] = yy[1,2,g]*base_bzp[1,1,g];
	yy[3,1,g] = yy[1,1,g]*(1+pow(base_bzp[3,1,g],2));
	yy[2,1,g] = yy[1,1,g]+(yy[3,1,g]-yy[1,1,g])*base_bzp[2,1,g];
	yy[4,1,g] = yy[1,1,g]*(1+pow(base_bzp[4,1,g],2));
	yy[1,3,g] = yy[1,1,g]*base_bzp[1,3,g];
	yy[3,3,g] = yy[1,3,g]*(1+pow(base_bzp[3,3,g],3));
	yy[2,3,g] = yy[1,3,g]+(yy[3,3,g]-yy[1,3,g])*base_bzp[2,3,g];
	yy[4,3,g] = yy[1,3,g]*(1+pow(base_bzp[4,3,g],3));
    }
    return yy;
}
data {
  ...
}

but you need a separate function for each moderately complicated transformation. Then, in R, do,

library(rstan)
expose_stan_functions("complicated_model.stan")

and now there will be a make_yy function in the R environment that you can call with inputs that are known to be valid.

Thank you so much Ben, your help is always much appreciated.

I have tried what you suggest, and it actually kind of “worked”, as you can verify from the code attached below.
There are some issues I still see though: I am not able to actually verify the functions in R (i.e. as “real” R functions) because I get the

> Error in make_yy(base_bzp, data$dmultB_parms, 3, 2) : 
>   VECTOR_ELT() can only be applied to a 'list', not a 'double'

This also happens when I try to transform, say, the “base_bzp” array into a list with a

lapply

on every dimension of the array.
Any idea of what’s going on?

Also, I still get some log(0) error from the chain running.

Thank you very much, you’re helping me a lot.

data: https://www.dropbox.com/s/gq3nmxz4cd1y2da/datastan.rda?dl=0
initial chain values: https://www.dropbox.com/s/vwanl598n5bek92/initstan.rda?dl=0

code:

functions {

real[,,] make_yy(real[,,] base_bzp, vector dmultB_parms, int nages, int ngenders){
  real yy[4,nages,ngenders];
 	yy[1,2,1] = base_bzp[1,2,1];
	yy[1,2,2] = base_bzp[1,2,1]*(dmultB_parms[1]+(dmultB_parms[2]-dmultB_parms[1])*base_bzp[1,2,2]);
  for(g in 1:ngenders){
	yy[3,2,g] = yy[1,2,g]+yy[1,2,g]*base_bzp[3,2,g];
	yy[2,2,g] = yy[1,2,g]+(yy[3,2,g]-yy[1,2,g])*base_bzp[2,2,g];
	yy[4,2,g] = yy[1,2,g]+yy[1,2,g]*base_bzp[4,2,g];
	//
	yy[1,1,g] = yy[1,2,g]*base_bzp[1,1,g];
	yy[3,1,g] = yy[1,1,g]*(1+pow(base_bzp[3,1,g],2));
	yy[2,1,g] = yy[1,1,g]+(yy[3,1,g]-yy[1,1,g])*base_bzp[2,1,g];
	yy[4,1,g] = yy[1,1,g]*(1+pow(base_bzp[4,1,g],2));
	//
	yy[1,3,g] = yy[1,1,g]*base_bzp[1,3,g];
	yy[3,3,g] = yy[1,3,g]*(1+pow(base_bzp[3,3,g],3));
	yy[2,3,g] = yy[1,3,g]+(yy[3,3,g]-yy[1,3,g])*base_bzp[2,3,g];
	yy[4,3,g] = yy[1,3,g]*(1+pow(base_bzp[4,3,g],3)); }
  return yy;}
  
real[,,] make_beta(real[,,] base_beta, int ngenders, int nages, int nsyp, real less_beta){
  real beta[nsyp,nages,ngenders];
  for(g in 1:ngenders){ for(k in 1:nages){ beta[1,k,g] = base_beta[1,k,g];
  for(w in 2:nsyp){ beta[w,k,g] = beta[(w-1),k,g]-less_beta+(0.95-beta[(w-1),k,g]+less_beta)*base_beta[(w-1),k,g]; }   }}
  return beta;}

real[,,] make_bezier_curve(real[,,] P, real[] yp, real ypo, int nsyp, int nages, int ngenders){
  real bez[nsyp,nages,ngenders];
  for(w in 1:nsyp){for(k in 1:nages){for (g in 1:ngenders){
	bez[w,k,g] = fmax(0, fmin(1, pow( yp[w]/ypo, 3)*P[1,k,g]+ 3*pow(yp[w]/ypo, 2)*((w-1)/ypo)*P[2,k,g]+
	3*(yp[w]/ypo)*pow((w-1)/ypo,2)*P[3,k,g]+ pow((w-1)/ypo,3)*P[4,k,g] )); }}}
  return bez;
  }

}

data{
// prevalence, standard dev., # case rep., populations
	real<lower=0, upper=1> 	P_hat[8,3,2];         //prevalence
	real<lower=0> 			    s_over_rootn[8,3,2];  //std.error from prevalence data
	int<lower=0> 			      R_by_year[18,3,2];    //case reports
	int<lower=0> 		    	  N_by_year[18,3,2];    //populations
// n.years(pairs), n.years(singles), n.age groups, n.genders
	int<lower=1> 			      nyrs;                 // #year pairs
	int<lower=1> 			      n_single_years;       // #single years
	int<lower=1, upper=3> 	nages;                // #age groups (3)
	int<lower=1, upper=2> 	ngenders;             // #genders (2)

// 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, upper=0.25> less_beta;
}

transformed data{
// indices (defined as real) to avoid integer division in Bezier definition
	real<lower=0, upper=n_single_years+pred> yp[n_single_years+pred];
	real<lower=0> yp1; 
// assignments
	for(w in 1:(n_single_years+pred)) yp[w]= n_single_years+pred-w;
	yp1 = (n_single_years+pred-1);
}

parameters{
// numbers
	real<lower=0, upper=1> psy;				//multiplier for prob sympt, gets a prior distrib beta
	real<lower=0, upper=1> vs;				//prob.treatment given symptomatic, gets a prior distrib beta	
	real<lower=0, upper=1> d_unt_1;	  //duration untreated, F, gets a prior distrib beta
	real<lower=0, upper=1> d_unt_2;   //duration untreated, M, gets a prior distrib beta
	real<lower=0, upper=1> d_ts_1;    //duration trt, sympt F, gets a prior distrib beta
	real<lower=0, upper=1> d_ts_2;    //duration trt, sympt M, gets a prior distrib beta
	real<lower=dmultB_parms[1], upper=dmultB_parms[2]> 	d_multB; 	//multiplier for Bezier males vs. females, gets uniform distrib. betw. 0.1 & 0.8
//vectors	
	vector<lower=dmultN_parms[1], upper=dmultN_parms[2]>[ngenders] d_multN; //multiplier for duration trt asympt, gets uniform distrib. betw. 0 & 1
// matrices
	real 					a[4,nages,ngenders]; 							                        //coeff cubic polyn for prevalence
	real<lower=0, upper=1> 	base_bzp[4,nages,ngenders]; 					          // bases for unif distrib on Bezier points
	//these are the ones who get a U(0,1) distribution that will be transformed to obtain the "real" distribution
	//of the support points of the beziér curve
	real<lower=0, upper=1> 	base_beta[(n_single_years+pred),nages,ngenders];// bases for unif distrib on reporting
	//same idea applied to the completeness of reporting
	} 

transformed parameters {
// numbers
	real<lower=m_int, upper=m_int+m_slope> 					m;	        //prob.symptomatic (transf.from psy & data)
// vectors	
	vector<lower=du_int, upper=du_int+du_slope>[ngenders] 	du;	//dur.untreated (transf.from d_unt_1, d_unt_2 & data)
	vector<lower=dt_int, upper=dt_int+dt_slope>[ngenders] 	dt;	//dur.treated, symptomatic (transf.from d_ts_1, d_ts_2 & data)				
	vector<lower=dt_int, upper=du_int+du_slope>[ngenders] dnew; //dur.treated, asymptomatic (transf.from du,dt,d_multN)	
// matrices
	real<lower=0, upper=2> 		yy[4,nages,ngenders];                         //support points for Bezier curves (transformed from base_bzp)
	real           				    P[4,nages,ngenders];	                        //transformed support points for Bezier (transformed from yy)
	real<lower=0>	 			      dur[(n_single_years+pred),nages,ngenders]; 		//overall duration 	
	real<lower=0>			        pi_w[(n_single_years+pred),nages,ngenders];	 	//prevalence_w (transf of the coefficients a (cubic polyn))
	real<lower=0> 				    pi_j[nyrs,nages,ngenders];					        	//prevalence_j (prevalence on the paired years)
	real 						          bezcurve[(n_single_years+pred),nages,ngenders];//uncapped Bezier (from Ps)
	real<lower=0, upper=1> 		bez[(n_single_years+pred),nages,ngenders]; 	 	//capped (at 0 and 1) Bezier curve  (from uncapped Bezier)	
	real<lower=0>          		gam[(n_single_years+pred),nages,ngenders]; 	 	//prob.detection (fraction treated, from m, vs, bez)
	real<lower=0, upper=0.95> beta[(n_single_years+pred),nages,ngenders];		//prob.case is correctly reported in recording system (transf from base_beta)
	real<lower=0>				      eta[(n_single_years+pred),nages,ngenders]; 	 	//incidence (transf from prevalence (prev/dur))
	real<lower=0, upper=1>		theta[(n_single_years+pred),nages,ngenders];  //rate case reports (from eta, beta, gamma)

// ASSIGMENTS //
// 1. map elements of base_bzp into those of yy
// these allow me to only put uniform(0,1) priors, and the transformations
// take care of all the relationships-constraints

//	yy[1,2,1] = base_bzp[1,2,1];
//	yy[1,2,2] = base_bzp[1,2,1]*(dmultB_parms[1]+(dmultB_parms[2]-dmultB_parms[1])*base_bzp[1,2,2]);
//  for(g in 1:ngenders){
//	yy[3,2,g] = yy[1,2,g]+yy[1,2,g]*base_bzp[3,2,g];
//	yy[2,2,g] = yy[1,2,g]+(yy[3,2,g]-yy[1,2,g])*base_bzp[2,2,g];
//	yy[4,2,g] = yy[1,2,g]+yy[1,2,g]*base_bzp[4,2,g];
	//
//	yy[1,1,g] = yy[1,2,g]*base_bzp[1,1,g];
//	yy[3,1,g] = yy[1,1,g]*(1+pow(base_bzp[3,1,g],2));
//	yy[2,1,g] = yy[1,1,g]+(yy[3,1,g]-yy[1,1,g])*base_bzp[2,1,g];
//	yy[4,1,g] = yy[1,1,g]*(1+pow(base_bzp[4,1,g],2));
	//
//	yy[1,3,g] = yy[1,1,g]*base_bzp[1,3,g];
//	yy[3,3,g] = yy[1,3,g]*(1+pow(base_bzp[3,3,g],3));
//	yy[2,3,g] = yy[1,3,g]+(yy[3,3,g]-yy[1,3,g])*base_bzp[2,3,g];
//	yy[4,3,g] = yy[1,3,g]*(1+pow(base_bzp[4,3,g],3));

 yy = make_yy(base_bzp, dmultB_parms, nages, ngenders);

// 2. map elements of base_beta into beta
//  for(k in 1:nages){ beta[1,k,g] = base_beta[1,k,g];
//  for(w in 2:(n_single_years+pred)){ beta[w,k,g] = beta[(w-1),k,g]-less_beta+(0.95-beta[(w-1),k,g]+less_beta)*base_beta[(w-1),k,g]; } //closes w
	
	beta = make_beta(base_beta, ngenders, nages, (n_single_years+pred), less_beta);

for(g in 1:ngenders){for(k in 1:nages){	
// 3. calculate Ps from yys (Beziér support points)
	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];

// 4. Beziér curve (to be truncated at 0 and 1 below)
// the curve is extended to the whole definition set (i.e.#years+#prediction years) to allow for smooth predictions
// for(w in 1:(n_single_years+pred)){					// opens w again
//	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];
// 4.1 truncation at 0 and 1
//	bez[w,k,g]  = fmax(0, fmin(1, bezcurve[w,k,g]) ); } //closes w
  }}  //closes k, closes g opened above

  bez = make_bezier_curve(P, yp, yp1, (n_single_years+pred), nages, ngenders);

// 5. prob symptomatic (from literature)
	m = m_int+m_slope*psy;

for(g in 1:ngenders){for(k in 1:nages){for(w in 1:(n_single_years+pred)){
// 6. detection
	gam[w,k,g]  = m*vs+(1-m)*bez[w,k,g];
}}}
  //print("gam: ", gam);

// 7. durations(untreated and treated, symptomatic, from literature)
	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;

for(g in 1:ngenders){
// 8. duration treated, asymptomatic (as a transformation of dt, du, dmultN)
	dnew[g] = dt[g]+(d_multN[g])*(du[g]-dt[g]);
	
for(k in 1:nages){for(w in 1:(n_single_years+pred)){
// 9. 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];
	
// 10. prevalence (logit(prevalence)=cubic.polynomial => prevalence=inv_logit(cubic.polynomial))
	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) );
	for(j in 1:nyrs){ pi_j[j,k,g]=(pi_w[(2*j-1),k,g]+pi_w[(2*j),k,g])/2; }

// 11. incidence (prevalence/duration)
	eta[w,k,g]  = pi_w[w,k,g]/dur[w,k,g];

// 12. rate of case reports (true incidence*p(incident case is diagnosed)*p(diagnosed case captured in reporting system))
// 	   OR p(incident)*p(diagnosed|incident)*p(reported|incident and diagnosed)
	theta[w,k,g]= eta[w,k,g]*beta[w,k,g]*gam[w,k,g];
}}} // closes w, closes k, closes g
  //print("eta: ", eta);
  //print("theta: ", theta);
  
}

model {
// PRIORS //
	//Bezier points
	//females [,,1] first one is beta (for age group 2), the others uniform, already defined in the transf.parameters block
	base_bzp[1,2,1] ~ beta(pbetaBz[1], pbetaBz[2]);
	
	// should be redundant as already defined in the transf.parameters block
	//base_bzp[1,1,1]	~ uniform(0,1);
	//base_bzp[1,3,1]	~ uniform(0,1);
	//males [,,2]
	//for(k in 1:nages){base_bzp[1,k,2]~ uniform(0,1);}
	//for(s in 2:4){for(k in 1:nages){for(g in 1:ngenders){ base_bzp[s,k,g] ~ uniform(0,1);}}}

	// Reporting
for(k in 1:nages){
	base_beta[1,k,1] ~ beta(pbetaF[1],pbetaF[2]);
	base_beta[1,k,2] ~ beta(pbetaM[1],pbetaM[2]);
	// as above, already defined
	// for(g in 1:ngenders){ for(w in 2:(n_single_years+pred)){ base_beta[w,k,g] ~ uniform(0,1); }}
}	
	
	psy ~ beta(pm[1],pm[2]);
	vs 	~ beta(pvs[1], pvs[2]);

	// cubic polyn 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))
}}

// 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
	// should be redundant as already specified in the parameters to have
	// limits between the "real" uniform limits here below
	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]);}
	}}

}

I have no idea what’s going on this case, but it’s clearly calling a function designed for lists with a double, so you just need to figure out where. Just look at the types that make_yy expects and the types of the things you’re passing in.