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