Session aborts when sampling commences with Rstan

Hi, every time I try to run this MARSS model (code below), the model compiles and then R aborts the session just as sampling begins. I have uninstalled and re-installed rstan and rtools, I have tried running the model both with and without a makevars file (and when running with a makevars file I removed the terms that cause windows to melt down), and I have tried running the model after deleting the .rds file. I have ensured that I am using rstan v2.26.22, and rtools v4.3.

Anyone have any other suggestions to get around these session aborts?

// 10 timeseries and 3 states MARSS model

data {
  int<lower=0> obs;       //total number of real observations across all timeseries
  int<lower=0> N;         //number of years in the longest timeseries - max-min year + 1
  int<lower=0> K;         //number of timeseries in the observation portion of the model
  int<lower=0> J;         //number of states
  int<lower=0> P[K];      //number of real data points in each timeseries of observations
  vector[J] p_st;         //mean of first timestep in the timeseries by state grouping
  vector[obs] tsx;        //vector of observations
  int juv_type[obs];      //integers defining 2(smolt) or 1(fry)
  int<lower=0> ii[obs];   //vector of integers describing the indexed locations of each observation relative to the state timeseries

parameters {
  matrix[J,N] x;                //estimated true underlying data (the states)
  vector[K] beta;               //persistent difference between observations and a shared state
  vector[2] beta2;              //categorical variable for effect of smolt vs fry counts in the calc of produ/surv
  real<lower=0> sigma;          //standard deviation of process
  vector[K] sigma_o;            //standard deviation of observation
  cholesky_factor_corr[J] L;    //lower triangular of a positive definite correlation matrix for the process
  cholesky_factor_corr[K] L_o;  //lower triangular of a positive definite correlation matrix for the observation

//calculating the residual error between time-steps
transformed parameters {
  matrix[J,N] mu;
  matrix[K,N] x_m;
  vector[J] sigma_u;
  mu[,1] = p_st;
  for (i in 2:N){
    mu[,i] = x[,i-1];
  x_m[1,] = x[1,];
  x_m[2,] = x[1,];
  x_m[3,] = x[1,];
  x_m[4,] = x[1,];
  x_m[5,] = x[2,];
  x_m[6,] = x[1,];
  x_m[7,] = x[1,];
  x_m[8,] = x[2,];
  x_m[9,] = x[3,];
  x_m[10,] = x[3,];
  sigma_u[1] = sigma;
  sigma_u[2] = sigma;
  sigma_u[3] = sigma;


model {
  int pos;
  pos = 1;
  for(k in 1:K){
    for(i in 1:P[k]){
      segment(tsx, pos, P[k])[i] ~ normal(beta[k] + beta2[segment(juv_type, pos, P[k])[i]] + x_m[k, segment(ii, pos, P[k])[i]], sigma_o[k]); //observation portion of the model (no intercept)
      pos = pos + P[k]; //segment takes a long vector and chops it into segments according to the lengths of each segment
    }                   //here we chop the observations and the positions of the observations relative to the state time-series 
  }                     //beta is a categorical variable that allows observations to be higher or lower than average (i.e., allows for persistent bias)
                        //beta2 is the effect of whether the calculation of productivity/survival uses smolts or fry
  x[,1] ~ multi_normal_cholesky(p_st, diag_pre_multiply(sigma_u, L));
  for(n in 2:N){
    x[,n] ~ multi_normal_cholesky(mu[,n], diag_pre_multiply(sigma_u, L));
  beta ~ normal(0, 2);
  beta2 ~ normal(0, 2);
  sigma ~ cauchy(0.3, 0.05);       //SD of process
  sigma_o ~ cauchy(0.3, 0.05);     //SD of observations
  L ~ lkj_corr_cholesky(2);        //prior for lower triangular of process
  L_o ~ lkj_corr_cholesky(1.5);      //prior for lower triangular of observation

generated quantities {
  matrix[J, N] nuS;
  vector[obs] nuY;
  nuS[,1] = multi_normal_cholesky_rng(p_st, diag_pre_multiply(sigma_u, L));
  for(i in 2:N){
    nuS[,i] = multi_normal_cholesky_rng(mu[,i], diag_pre_multiply(sigma_u, L));

If possible, add also code to simulate data or attach a (subset of) the dataset you work with.

Operating System: Windows 10 Enterprise
Interface Version: R-4.3.0
Compiler/Toolkit: Rtools 4.3

I recently experienced a similar issue with rstan and was able to resolve it by installing a different StanHeaders version. Which StanHeaders version are you using?

Hi, I am using version 2.26.27 of Stan Headers. I just tried going back to version 2.26.26 but the same error occurred. I’ll keep trying moving back through versions and see if one of them works.

Issue resolved! I removed and re-installed rstan and StanHeaders to older versions which produced slightly different failures. Then I removed and re-installed with the latest versions of rstan and StanHeaders and tried a standard rstan test model which did work.

Next I tried running my model, which provided an error rather than immediate R session abort (improvement over before). I discovered that the code line “pos = pos + P[k];” needed to be placed outside the “for(i in 1:P[k]){}” for loop, but inside the “for(i in 1:P[k]){}” for loop for the segment function to work properly.

1 Like