Initialization issues with Rstan 2.21.2 on MacOS Catalina

Modeling
, ,
1

If this is a duplicate, please apologize in advance. But I don’t seem to find other posts with the same problem here, although it seems that others have encountered similar situations.

Summary of the problem
After compiling stan models, I get an error message at the beginning of sampling from the sampler.
[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”
error occurred during calling the sampler; sampling not done
here are whatever error messages were returned

I tried running with cores=1，error showed that：
error occurred during calling the sampler; sampling not done

There is some information attached, which may help to understand：
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: binomial_lpmf: Successes variable is 176, but must be in the interval [0, 55] (in ‘model7a72a5b9a43_af543013862bc307cac5e99e31ad53e6’ at line 111)

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.”
error occurred during calling the sampler; sampling not done

code_to_run_your_model(if_applicable)
library(plyr)
library(StanHeaders)
library(ggplot2)
library(rstan)
library(loo)

abmodelnb_f = "
data{
int N;
int Nt;
int Ns;

int TP[N];
int Dis[N];
int TN[N];
int NDis[N];
int Study[N];
int Test[N];

}
parameters{
matrix[2, Nt] logitmu;
vector[Ns] nu[2];
matrix[Ns, Nt] delta[2];
vector<lower=0>[Nt] tau[2]; //*
vector<lower=0>[2] sigmab;
real<lower=-1, upper=1> rho;
}
transformed parameters{
matrix[Ns, 2] p_i[Nt];
matrix[2, Nt] MU;
matrix[2, Nt] RR;
matrix[2, Nt] OR;
vector[Nt] DOR;
vector[Nt] S;
matrix[Nt, Nt] A;
matrix[Nt, Nt] B;
matrix[Nt, Nt] C;

vector<lower=0>[Nt] tausq[2];
vector<lower=0>[2] sigmabsq;

matrix[Nt, Nt] sigmasq[2];
matrix[Nt, Nt] rhow[2];


for (i in 1:Ns){
    for (j in 1:2){
        for (k in 1:Nt)
            p_i[k][i,j] = inv_logit(logitmu[j,k] +  nu[j][i] + delta[j][i,k]);
    }
}

for (j in 1:2){
    for (k in 1:Nt){
        MU[j,k] = mean(col(p_i[k], j));
    }
    tausq[j] = (tau[j]).*(tau[j]);
}

for (j in 1:2){
    for (k in 1:Nt){
        RR[j, k] = MU[j, k]/MU[j, 1]; 
        OR[j, k] = (MU[j, k]*(1 - MU[j, 1]))/(MU[j, 1]*(1 - MU[j, k]));
     }
}

for (l in 1:Nt){
    DOR[l] = (MU[1, l]*MU[2, l])/((1 - MU[1, l])*(1 - MU[2, l]));

    for(m in 1:Nt){
        A[l, m] = if_else((MU[1, l] > MU[1, m]) && (MU[2, l] > MU[2, m]), 1, 0);
        B[l, m] = if_else((MU[1, l] < MU[1, m]) && (MU[2, l] < MU[2, m]), 1, 0);
        C[l, m] = if_else((MU[1, l] == MU[1, m]) && (MU[2, l] == MU[2, m]), 1, 0);
    }

    S[l] = (2*sum(row(A, l)) + sum(row(C, l)))/(2*sum(row(B, l)) + sum(row(C, l)));
}

sigmabsq = (sigmab).*(sigmab);

for (j in 1:2){
    for (k in 1:Nt){
        for (l in 1:Nt){
            sigmasq[j][k,l] = (sigmabsq[j] + tausq[j][k])*((sigmabsq[j] + tausq[j][l]));
            rhow[j][k,l] = sigmabsq[j]/sqrt(sigmasq[j][k,l]);
        }
    }
}

}
model{
//Priors
for (j in 1:2){
logitmu[j] ~ normal(0, 5);
tau[j] ~ cauchy(0, 2.5);
}

sigmab ~ cauchy(0, 2.5);
rho ~ uniform(-1, 1);

nu[2] ~ normal(0, sigmab[2]);


for (i in 1:Ns){

    nu[1][i] ~ normal((sigmab[1]/sigmab[2])*rho*nu[2][i], sqrt(sigmabsq[1]*(1 - (rho*rho))));

    for (j in 1:2){
        for (k in 1:Nt)
            delta[j][i,k] ~ normal(0, tau[j][k]);
    }
}

for (n in 1:N){
    TP[n] ~ binomial(Dis[n], p_i[Test[n]][Study[n], 1]);
    TN[n] ~ binomial(NDis[n], p_i[Test[n]][Study[n], 2]);
}

}
generated quantities{
vector[2*N] loglik;

for (n in 1:N)
    loglik[n] = binomial_lpmf(TN[n]| NDis[n], p_i[Test[n]][Study[n], 1]);

for (n in (N+1):(2*N))
    loglik[n] = binomial_lpmf(TN[n-N]| NDis[n-N], p_i[Test[n-N]][Study[n-N], 2]);

}
"
N <- nrow(df1)
Ns <- max(df1$Study)
Nt <- max(df1$Test)

datalist <- list(
N = N,
Ns = Ns,
Nt = Nt,
TP = df1$TP,
Dis = df1$Dis,
TN = df1$TN,
NDis = df1$NDis,
Test = df1$Test,
Study = df1$Study)
model1.0 <- stan(model_code = abmodelnb_f,
data=datalist,
chains = 3,
iter = 4000,
warmup = 2000,
thin = 5,
seed = 1,
cores=3,
verbose=FALSE)

RStan Version
2.21.2

R Version
4.02

Operating System:
MacOS Catalina 10.15.5

2

Your number of successes is greater than the number of trials when you call the binomial likelihood.

3

Thank you, I will try to change the value.