Error: accessing element out of range

Modeling
1

Hello,

I have run into the following progress error message when running RStan for a sparse ICAR model:

SAMPLING FOR MODEL 'sparse_icar_nocovariates' NOW (CHAIN 1).
Unrecoverable error evaluating the log probability at the initial value.
Exception: Exception: []: accessing element out of range. index -2147483648 out of range; expecting index to be between 1 and 56; index position = 1phit_W  (in 'model20d86748338d_sparse_icar_nocovariates' at line 25)
  (in 'model20d86748338d_sparse_icar_nocovariates' at line 79)

This is repeated for every chain.

My searching found others had asked about this kind of error previously, but the solutions did not seem relevant to my situation (normally their index number was within the index range).

My R code is below (using Scotland lip cancer data for reproducibility), followed by my Stan model (which is directly from http://mc-stan.org/users/documentation/case-studies/mbjoseph-CARStan.html (under the postcript at the end). The only change made to the Stan code was to remove the covariate X details).

I am using R version 3.3.3 and rstan 2.17.3.

Any suggestions will be appreciated - thank you!

##R code
library(rstan)  
library(CARBayesdata)
library(CARBayes)
library(shapefiles)
library(spdep)

options(mc.cores = parallel::detectCores())  

#Create adjacency matrix
data(lipdata)
data(lipshp)
data(lipdbf)

lipdbf$dbf <- lipdbf$dbf[ ,c(2,1)]
data.combined <- combine.data.shapefile(data=lipdata, shp=lipshp, dbf=lipdbf)
W.nb <- poly2nb(data.combined, row.names = rownames(lipdata))
W <- nb2mat(W.nb, style="B")

# Define MCMC parameters 
niter <- 1E4   # definitely overkill, but good for comparison
nchains <- 4

full_d <- list(n = 56,                 # number of observations
               p = 1,                  # number of beta0 (here, just the intercept)
               y = lipdata$observed,               # observed number of cases
               log_offset = log(lipdata$expected), # log(expected) num. cases
               W_n = 264,                 #SumNumNeigh
               W = W)                         # adjacency matrix

full_fit <- stan('sparse_icar_nocovariates.stan', data = full_d, 
                 iter = niter, chains = nchains, verbose = FALSE)

##Stan model

functions {
  /**
  * Return the log probability of a proper intrinsic autoregressive (IAR) prior 
  * with a sparse representation for the adjacency matrix
  *
  * @param phi Vector containing the parameters with a IAR prior
  * @param tau Precision parameter for the IAR prior (real)
  * @param W_sparse Sparse representation of adjacency matrix (int array)
  * @param n Length of phi (int)
  * @param W_n Number of adjacent pairs (int)
  * @param D_sparse Number of neighbors for each location (vector)
  * @param lambda Eigenvalues of D^{-1/2}*W*D^{-1/2} (vector)
  *
  * @return Log probability density of IAR prior up to additive constant
  */
  real sparse_iar_lpdf(vector phi, real tau,
    int[,] W_sparse, vector D_sparse, vector lambda, int n, int W_n) {
      row_vector[n] phit_D; // phi' * D
      row_vector[n] phit_W; // phi' * W
      vector[n] ldet_terms;
    
      phit_D = (phi .* D_sparse)';
      phit_W = rep_row_vector(0, n);
      for (i in 1:W_n) {
        phit_W[W_sparse[i, 1]] = phit_W[W_sparse[i, 1]] + phi[W_sparse[i, 2]];
        phit_W[W_sparse[i, 2]] = phit_W[W_sparse[i, 2]] + phi[W_sparse[i, 1]];
      }
    
      return 0.5 * ((n-1) * log(tau)
                    - tau * (phit_D * phi - (phit_W * phi)));
  }
}
data {
  int<lower = 1> n;
  int<lower = 1> p;
  int<lower = 0> y[n];
  vector[n] log_offset;
  matrix<lower = 0, upper = 1>[n, n] W; // adjacency matrix
  int W_n;                // number of adjacent region pairs
}
transformed data {
  int W_sparse[W_n, 2];   // adjacency pairs
  vector[n] D_sparse;     // diagonal of D (number of neigbors for each site)
  vector[n] lambda;       // eigenvalues of invsqrtD * W * invsqrtD
  
  { // generate sparse representation for W
  int counter;
  counter = 1;
  // loop over upper triangular part of W to identify neighbor pairs
    for (i in 1:(n - 1)) {
      for (j in (i + 1):n) {
        if (W[i, j] == 1) {
          W_sparse[counter, 1] = i;
          W_sparse[counter, 2] = j;
          counter = counter + 1;
        }
      }
    }
  }
  for (i in 1:n) D_sparse[i] = sum(W[i]);
  {
    vector[n] invsqrtD;  
    for (i in 1:n) {
      invsqrtD[i] = 1 / sqrt(D_sparse[i]);
    }
    lambda = eigenvalues_sym(quad_form(W, diag_matrix(invsqrtD)));
  }
}
parameters {
  vector[p] beta0;
  vector[n] phi_unscaled;
  real<lower = 0> tau;
}
transformed parameters {
  vector[n] phi; // brute force centering
  phi = phi_unscaled - mean(phi_unscaled);
}
model {
  phi_unscaled ~ sparse_iar(tau, W_sparse, D_sparse, lambda, n, W_n);
  beta0 ~ normal(0, 2.5);
  tau ~ gamma(2, 2);
  y ~ poisson_log(beta0 + phi + log_offset);
}

[edit: escaped code]

2

That means on line 79 it called something on line 25 and you tried to index phit_W with an unintialized integer. That would imply that some elements of W_sparse are uninitialized.