I’m working on a model that I want to try slight variations of, and I don’t want to recode the entire thing for each variation. Thus, I wish to include a Boolean variable (in below code denoted includeKeffectonalpha) that determines whether some parameters are to be included in the model or not. That is, I want to use if [if(includeKeffectonalpha){do this}] and if/else statements [if(includeKeffectonalpha){do this}else{do the other thing}]. I have tried this in several ways, none of which seems to work both when includeKeffectonalpha=0 and when includeKeffectonalpha=1. The issue seems to be declaration.
Strategy one: Declare all variable as usual, independent on whether they will be used in the model. This appears to works fine when includeKeffectonalpha=1, but gives the following errors when includeKeffectonalpha=0:
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: model26dc5f757c4a_Hunt_stan_code_krets_lognstd_if_namespace::write_array: talpha[k0__] is nan, but must be greater than or equal to 0 (in ‘model26dc5f757c4a_Hunt_stan_code_krets_lognstd_if’ at line 44)
Chain 1:
Chain 1: Initialization between (-2, 2) failed after 100 attempts.
Chain 1: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”
[1] “error occurred during calling the sampler; sampling not done”
Strategy two: use the trick proposed here: https://dev.to/martinmodrak/optional-parametersdata-in-stan-4o33, where parameters are declared as having dimensions dependent on the Boolean input variable, in this case making them zero sized if includeKeffectonalpha=0. Independent on includeKeffectonalpha=0 or includeKeffectonalpha=1, I however get the following error.
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
No matches for:
real[] * real
Available argument signatures for operator*:
real * real
vector * real
row vector * real
matrix * real
row vector * vector
vector * row vector
matrix * vector
row vector * matrix
matrix * matrix
real * vector
real * row vector
real * matrix
Expression is ill formed
error in ‘model26dc44c91d4b_Hunt_stan_code_krets_lognstd_if’ at line 75, column 44
73:
74: if (includeKeffectonalpha){
75: talpha[il] = taoSalpha * taoLalpha[il]; //
^
76: }
Error in stanc(file = file, model_code = model_code, model_name = model_name, :
failed to parse Stan model ‘Hunt_stan_code_krets_lognstd_if’ due to the above error.
Any ideas on how to solve this? Are there other strategies that works better? The code as per my strategy two and strategy one commented out is pasted below.
data {
int<lower=1> NrOfLans;
int<lower=1> NrOfKretsar;
int<lower=1> Reports;
int<lower=0> K[Reports];
int<lower=0> A[Reports];
int<lower=1> KretsNrSeqID[Reports];
int<lower=1> LansSeqID[NrOfLans];
int<lower=1> LanSeqID[Reports];
int<lower=1> KretsarSeqID[NrOfKretsar];
int<lower=1> KretsListLanSeqID[NrOfKretsar];
int<lower=0, upper=1> includeKeffectonalpha;//Boolean declaration
}
parameters {
real logLm[NrOfLans];
real logKm[NrOfKretsar];
real logLalpha[NrOfLans];
real logSm;
real logSalpha;
real<lower=0> taoLm[NrOfLans];
real<lower=0> Tm;
real<lower=0> Talpha;
real<lower=0> sigmataoLm;
real<lower=0> taoSm;
real logKalpha[includeKeffectonalpha ? NrOfKretsar : 0];// size of array is zero if includeKeffectonalpha is false (0)
real<lower=0> taoLalpha[includeKeffectonalpha ? NrOfLans : 0];
real<lower=0> sigmataoLalpha[includeKeffectonalpha ? 1 : 0];
real<lower=0> taoSalpha[includeKeffectonalpha ? 1 : 0];
// real logKalpha[NrOfKretsar];//
// real<lower=0> taoSalpha;
// real<lower=0> sigmataoLalpha;
// real<lower=0> taoLalpha[NrOfLans];
}
transformed parameters {
real<lower=0> alpha[NrOfKretsar];
real<lower=0> m[NrOfKretsar];
real<lower=0> tm[NrOfLans];
// real<lower=0> talpha[NrOfLans];
//
real<lower=0> talpha[includeKeffectonalpha ? NrOfLans : 0];
for (ik in 1:NrOfKretsar){
m[ik] = exp( logSm + logLm[KretsListLanSeqID[ik]] + logKm[ik] );
if (includeKeffectonalpha) {
alpha[ik] = exp(logSalpha + logLalpha[KretsListLanSeqID[ik]] + logKalpha[ik] );
} else {
alpha[ik] = exp(logSalpha + logLalpha[KretsListLanSeqID[ik]]);
}
}
for (il in 1:NrOfLans){
tm[il] = taoSm * taoLm[il]; //tm is the precission of between circuit variability of logK for each county
if (includeKeffectonalpha){
talpha[il] = taoSalpha * taoLalpha[il]; //
}
}
}
model{
for (r in 1:Reports){
// Number of killed animals in report r is negative binomial distributed,
// but defined from the mean (m[KretsNrSeqID[r]] * A[r]) and shape (c[Lan[r]]) of the gamma distribution of the equivalent gamma-poisson mixture
K[r] ~ neg_binomial_2(m[KretsNrSeqID[r]] * A[r] ,alpha[KretsNrSeqID[r]]);
}
for(ik in 1:NrOfKretsar){
logKm[ik] ~ normal( 0 , tm[KretsListLanSeqID[ik]] ); // logKm[Kretsar[ik]] is the circuit effect on logm
if (includeKeffectonalpha){
logKalpha[ik] ~ normal( 0 , talpha[KretsListLanSeqID[ik]] ); // logKalpha[Kretsar[ik]] is the circuit effect on logalpha
}
}
for (il in 1:NrOfLans){
logLm[il] ~ normal( 0 , Tm );// Tm is the std of between county variability for logm
logLalpha[il] ~ normal( 0 , Talpha );// Talpha is the std of between county variability for logalpha
taoLm[il] ~ lognormal( 0 , sigmataoLm); // taoLm[il] is the std of variability of logm between circuits in county il
if (includeKeffectonalpha){
taoLalpha[il] ~ lognormal( 0 , sigmataoLalpha); //taoLalpha [il] is the std of variability of logalpha between circuits in county il and is log-normally distributed with std sigmataoLalpha
}
}
logSm ~ normal(0,100);
logSalpha ~ normal(0,100);
Tm ~ lognormal(0,5);
Talpha ~ lognormal(0,5);
sigmataoLm ~ lognormal(0,1);
taoSm ~ lognormal(0,2);
if (includeKeffectonalpha){
sigmataoLalpha ~ lognormal(0,1);
taoSalpha ~ lognormal(0,2);
}
}
generated quantities {
vector[Reports] log_lik;
for (r in 1:Reports){
log_lik[r] = neg_binomial_2_lpmf(K[r] | m[KretsNrSeqID[r]] * A[r], alpha[KretsNrSeqID[r]]);
}
}