icar_sparse.stan

functions {
  /**
   * Do a row-wise sum over a matrix
   */
  vector rowwise_sum(matrix mat) {
    int mat_rows = rows(mat);
    vector[mat_rows] D;
    for (i in 1:mat_rows) {
      D[i] = sum(mat[i]);
    } 
    return D;
  }
  /**
  * Return the log probability of a proper conditional autoregressive (CAR) prior
  * with a sparse representation for the adjacency matrix
  *
  * @param phi Vector containing the parameters with a CAR prior
  * @param tau Precision parameter for the CAR prior (real)
  * @param alpha Dependence (usually spatial) parameter for the CAR 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 CAR prior up to additive constant
  */
  real sparse_car_lpdf(vector phi, real tau, real alpha,
    sparse_matrix data W, vector D_sparse, vector lambda, int n, int W_n) {
      // row_vector[n] phit_W = rep_row_vector(0, n); // phi' * W
      // Is this right???
      row_vector[n] phit_W = phi' * W; 
      row_vector[n] phit_D = (phi .* D_sparse)'; // phi' * D
      vector[n] ldet_terms = log1m(alpha * lambda);
      /*
      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 * log(tau)
                    + sum(ldet_terms)
                    - tau * (phit_D * phi - alpha * (phit_W * phi)));
  }
}
data {
  int<lower = 1> n;
  int<lower = 1> p;
  matrix[n, p] X;
  int<lower = 0> y[n];
  vector[n] log_offset;
  int W_n;                // number of adjacent region pairs
  int W_sparse[W_n, 2];   // adjacency pairs
  sparse_matrix[n, n, W_sparse] W; // adjacency matrix
}
transformed data {
  // diagonal of D (number of neigbors for each site)
  // less efficient to put in func but easier for me to read
  vector[n] D_sparse = rowwise_sum(W);     
  // eigenvalues of (1/sqrt(rowsum(W))) * W * (1/sqrt(rowsum(W)))
  vector[n] lambda = eigenvalues_sym(quad_form(W, diag_matrix(1.0 ./ sqrt(D_sparse))));
}
parameters {
  vector[p] beta;
  vector[n] phi;
  real<lower = 0> tau;
  real<lower = 0, upper = 1> alpha;
}
model {
  phi ~ sparse_car(tau, alpha, W_sparse, D_sparse, lambda, n, W_n);
  beta ~ normal(0, 1);
  tau ~ gamma(2, 2);
  y ~ poisson_log(X * beta + phi + log_offset);
}