1.1 Autonormal model

This first example is an autonormal model which applies the CAR model directly to the data (as described above). This can also be written as \[y \sim Normal(\alpha + x \beta, \Sigma).\] Here you have \(n\) observations of a continuously measured outcome \(y\), possibly \(k=1\) or more covariates \(x\), and you’re accounting for spatial autocorrelation in \(y\) using the covariance matrix (following the CAR specification).

We fit a model like this using geostan::stan_car:

library(geostan)
data(georgia)
A <- shape2mat(georgia, "B")
car_list <- prep_car_data(A, style = "WCAR")

# prior distributions (to match the Stan model below)
prior_list <- list(intercept = normal(0, 5),
               beta = normal(0, 5),
           sigma = student_t(10, 0, 5)
           )

# an auto-model
# y = mu + rho * C (y - mu) + error
# mu = alpha + beta * .x 
# y = log(income); x = log(population)
# x will be centered: .x = x - mean(x)

fit <- stan_car(log(income / 10e3) ~ log(population / 10e3),
    data = georgia, car = car_list, prior = prior_list, centerx = TRUE)

The following Stan code implements the above model; if you’re following along, copy the Stan code and save it in a file inside your working directory; I’ll name it “autonormal.stan” and I’ll save the name as an R variable:

autonormal_file <- "autonormal.stan"

Here’s the Stan code:

functions {
#include wcar-lpdf.stan
}

data {
  // data
  int<lower=1> n;
  int<lower=1> k;
  vector[n] y;
  matrix[n, k] x;

    // CAR parts
  int nA_w;
  vector[nA_w] A_w;
  array[nA_w] int A_v;
  array[n + 1] int A_u;
  vector[n] Delta_inv;
  real log_det_Delta_inv;
  vector[n] lambda;
}

parameters {
  // spatial autocorrelation (SA) parameter
  real<lower=1/min(lambda), upper=1/max(lambda)> rho;
  
  // scale parameter
  real<lower=0> tau;
  
  // intercept
  real alpha;
  
  // coefficients
  vector[k] beta;  
}

model {
  vector[n] mu = alpha + x * beta;

  // Likelihood: y ~ Normal(Mu, Sigma)
  target += wcar_normal_lpdf(y |
                 mu, tau, rho, // mean, scale, SA
                 A_w, A_v, A_u, // stuff from prep_car_data
                 Delta_inv, 
                 log_det_Delta_inv,
                 lambda, 
                 n);    
  
  // prior for scale parameter
  target += student_t_lpdf(tau | 10, 0, 5);

  // prior for beta
  target += normal_lpdf(beta | 0, 5);

  // prior for intercept
  target += normal_lpdf(alpha | 0, 5);
}

The input file “wcar-lpdf.stan” contains the following code (save this code inside your working directory and name it “wcar-lpdf.stan”):

/**
 * Log probability density of the conditional autoregressive (CAR) model: WCAR specifications only
 *
 * @param y Process to model
 * @param mu Mean vector
 * @param tau Scale parameter
 * @param rho Spatial dependence parameter
 * @param A_w Sparse representation of the symmetric connectivity matrix, A
 * @param A_v Column indices for values in A_w
 * @param A_u Row starting indices for values in A_u
 * @param D_inv The row sums of A; i.e., the diagonal elements from the inverse of Delta, where M = Delta * tau^2 is a diagonal matrix containing the conditional variances.
 * @param log_det_D_inv Log determinant of Delta inverse.
 * @param lambda Eigenvalues of C (or of the symmetric, scaled matrix Delta^{-1/2}*C*Delta^{1/2}); for the WCAR specification, C is the row-standardized version of W.
 * @param n Length of y
 *
 * @return Log probability density of CAR prior up to additive constant
 */
real wcar_normal_lpdf(vector y, vector mu,
              real tau, real rho,
              vector A_w,
              array[] int A_v,
              array[] int A_u,
              vector D_inv,
              real log_det_D_inv,
              vector lambda,
              int n) {
  vector[n] z = y - mu;
  real ztDz = (z .* D_inv)' * z;
  real ztAz = z' * csr_matrix_times_vector(n, n, A_w, A_v, A_u, z);
  real ldet_ImrhoC = sum(log1m(rho * lambda));  
  return 0.5 * (
        -n * log( 2 * pi() )
        -2 * n * log(tau)
        + log_det_D_inv
        + ldet_ImrhoC
        - (1 / tau^2) * (ztDz - rho * ztAz));
}

From R, you can use the following code to prepare the input data (\(y\), \(x\), etc.). I use prep_car_data to get a list of parts for the CAR model, then I append the outcome data to the same list and pass it all to a Stan model to draw samples.

library(rstan)
library(geostan)
data(georgia)

A <- shape2mat(georgia, "B")
car_list <- prep_car_data(A, style = "WCAR")

# add data
## (centering covariates improves sampling efficiency)
car_list$y <- log(georgia$income / 10e3)
car_list$x <- scale(log(georgia$population / 10e3), center = TRUE, scale = FALSE)
car_list$k <- ncol(car_list$x)

# compile Stan model from file
autonormal_file <- "autonormal.stan"
car_model <- stan_model(autonormal_file)

# sample from model
samples <- sampling(car_model, data = car_list)

1.2 Hierarchical model

Disease mapping is a common use for CAR models, though these models have many applications. In this class of models, the CAR model is assigned as the prior distribution to a parameter vector \(\phi\), which is used to model disease incidence rates across small areas like counties. The disease data consist of counts \(y\) together with the size of population at risk \(p\), or possibly a different denominator such as the expected number of cases \(E\) (which occurs when using indirect age-standardization). These models have the form \[\begin{equation} \begin{aligned} &y \sim Poisson( exp((log(p) + \phi) ) \\ &\phi \sim Normal(\mu, \Sigma) \\ &\mu = \alpha + x \beta. \end{aligned} \end{equation}\]

Here you see that the linear predictor \(\mu\) is embedded with the CAR model.

Here is another way of writing down the very same model: \[\begin{equation} \begin{aligned} &y \sim Poisson(e^\eta) \\ &\eta = log(p) + \alpha + x \beta + \phi \\ &\phi \sim Normal(0, \Sigma) \\ \end{aligned} \end{equation}\]

It is possible to build a Stan modeling using either one of these two formulations. They are substantively equivalent, but you may find that one is far more amenable to Stan’s MCMC algorithm. The former specification (\(\mu\) inside the CAR model) is less commonly presented in papers, but in this author’s experience it is often far more stable computationally than forcing the CAR model to have zero mean. (You may well encounter a case where the opposite is true.)

The purpose of the offset term is sometimes unclear in introductory texts. The offset term \(p\) is from the denominator of the rates \(\frac{y}{p}\). The expectation or mean of the model is \[\begin{equation} \begin{aligned} &\mathbb{E}[y] = e^{log(p) + \mu} = e^{log(p)} \cdot e^{\mu} = p \cdot e^\mu \\ &\mathbb{E}\Bigl[ \frac{y}{p} \Bigr] = e^\mu \end{aligned} \end{equation}\] The smaller is the denominator, the less informative is the rate with respect to the characteristic level of risk bearing upon the population of the given period and place. With small denominators, chance renders the rates uninformative (and unstable over time and space), and the Poisson model accounts for this.

We can fit this type of model using geostan::stan_car:

data(georgia)
A <- shape2mat(georgia, "B")
car_list <- prep_car_data(A, style = "WCAR")

# prior distributions (to match the Stan model below)
prior_list <- list(intercept = normal(0, 5),
               beta = normal(0, 5),
           sigma = student_t(10, 0, 5)
           )


# Poisson model
# y ~ Poisson(pop * exp(mu))
# mu = alpha + beta * x + phi
# phi ~ CAR(0, Sigma)
# y = deaths; x = log(income);
# x will be centered: .x = x - mean(x)

fit <- stan_car(deaths.male ~ offset(log(pop.at.risk.male)) + log(income / 1e3),
    data = georgia, 
    car = car_list, 
    centerx = TRUE,
    family = poisson()
    )

The following Stan code provides a template for a simple disease mapping model using the CAR model as a prior for \(\phi\). As before, you can save the code in your working directory and give it a name like “car_poisson.stan”. I’ll store the file name in my R environment:

car_poisson_file <- "car_poisson.stan"

functions {
#include wcar-lpdf.stan
}

data {
  // data
  int<lower=1> n;
  int<lower=1> k;
  array[n] int<lower=0> y;
  matrix[n, k] x;
  vector[n] const_offset;

    // CAR parts
  int nA_w;
  vector[nA_w] A_w;
  array[nA_w] int A_v;
  array[n + 1] int A_u;
  vector[n] Delta_inv;
  real log_det_Delta_inv;
  vector[n] lambda;
}

parameters {
  // spatial autocorrelation (SA) parameter
  real<lower=1/min(lambda), upper=1/max(lambda)> rho;
  
  // scale parameter
  real<lower=0> tau;
  
  // intercept
  real alpha;
  
  // coefficients
  vector[k] beta;
  
  // SA trend component
  vector[n] phi;
}


model {
  vector[n] y_mu = exp(const_offset + phi);
  vector[n] phi_mu = alpha + x * beta;
  
  // Likelihood: y ~ Poisson( population * exp(phi) )
  target += poisson_lpmf(y | y_mu);
    
  // phi ~ Normal(Mu, Sigma)
  target += wcar_normal_lpdf(phi |
                 phi_mu, tau, rho, // mean, scale, SA
                 A_w, A_v, A_u, // stuff from prep_car_data
                 Delta_inv, 
                 log_det_Delta_inv,
                 lambda, 
                 n);

  // prior for scale parameter
  target += student_t_lpdf(tau | 10, 0, 1);

  // prior for beta
  target += normal_lpdf(beta | 0, 5);

  // prior for intercept
  target += normal_lpdf(alpha | 0, 5);
}

Again, you can use geostan::prep_car_data to easily convert the spatial weights matrix into the list of required inputs for the CAR model:

library(rstan)
library(geostan)
data(georgia)

A <- shape2mat(georgia, "B")
car_list <- prep_car_data(A, style = "WCAR")

# add data
car_list$y <- georgia$deaths.male
car_list$const_offset <- log(georgia$pop.at.risk.male)
car_list$x <- scale(log(georgia$income / 1e3), center = TRUE, scale = FALSE)
car_list$k <- ncol(car_list$x)

# compile Stan model from file
car_poisson_file <- "car_poisson.stan"
car_poisson <- stan_model(car_poisson_file)

# sample from model
samples <- sampling(car_poisson, data = car_list)

1.3 Zero-mean parameterization

It was noted above that the hierarchical CAR model can be described in more than one equivalent way. For one, we can apply the CAR model to the log-rates \(\phi\) just like we applied it to the outcome \(y\):

\[\begin{equation} \begin{aligned} y &\sim Poisson(pop * exp(\phi)) \\ \phi &= \mu + \rho C (\phi - \mu) + \epsilon \end{aligned} \end{equation}\]

where \(\mu = \alpha + x * beta\). That is the model we used above for disease mapping.

The alternative is to force \(\phi\) to have a mean of zero:

\[\begin{equation} \begin{aligned} y &\sim Poisson(pop * exp(\mu + \phi)) \\ \phi &\sim Normal(0, \Sigma). \\ \end{aligned} \end{equation}\]

Adopting former method often results in more stable and efficient MCMC sampling. But there are times when the former leads to poor MCMC sampling, particularly when many of the populations at risk and observed counts are very small. In such cases, you will find that the model is slow to converge (more MCMC samples are needed to avoid low ESS warnings) and you may also receive obscure warnings about the ‘Bayesian fraction of missing information’ and perhaps even warnings of divergent transitions. This is a widely-encountered phenomenon in Bayesian analysis. Here we are going to see how to handle the issue with these hierarchical CAR models.

When this happens, we need to make two changes to our hierarchical Poisson model. The first is to force the mean of \(\phi\) to equal zero instead of \(\mu=\alpha + x * \beta\). The second change pertains to the CAR covariance matrix for \(\phi\): we are going to specify the covariance matrix to have a variance parameter of \(\tau^2 = 1\); we symbolize this as \(\Sigma_1\). Then we bring the (actual) scale parameter back in by multiplication:

\[\begin{equation} \begin{aligned} y &\sim Poisson(pop * exp(\mu + \phi * \tau)) \\ \phi &\sim Normal(0, \Sigma_1) \\ \end{aligned} \end{equation}\]

The covariance matrix is, again, \(\Sigma_1 = (I - \rho C)^{-1}M\), and in this case \(M\) contains the variances \(m_{i,i} = \tau^2 * \delta_i = \delta_i\) (where the weights \(\delta_i\) are equal to one divided by the row sums of \(C\)).

In the geostan documentation, this method is refered to as the zero-mean parameterization (ZMP). For hierarchical CAR models, we can use the ZMP by adding zmp = TRUE:

# zero-mean parameterization of the hierarchical CAR model
car_list <- prep_car_data(shape2mat(georgia, "B", quiet = TRUE))
fit <- stan_car(deaths.male ~ offset(log(pop.at.risk.male)) + log(income / 1e3),
    data = georgia, 
    car = car_list, 
    centerx = TRUE,
    family = poisson(),
    zmp = TRUE
    )

(This zmp option is also available in stan_sar.)

To code this in Stan, we take “car_poisson.stan” and change the top of model block to the following:

model {
  vector[n] y_mu = exp(const_offset + alpha + x * beta + phi * tau);
  vector[n] phi_mu = rep_vector(0, n);
  
  // Likelihood: y ~ Poisson( population * exp(phi) )
  target += poisson_lpmf(y | y_mu);
    
  // phi ~ Normal(Mu, Sigma)
  target += wcar_normal_lpdf(phi |
     phi_mu,
     1, // scale: tau = 1
     rho, 
     A_w, A_v, A_u,
     Delta_inv, 
     log_det_Delta_inv,
     lambda, 
     n);
   // ... (priors, etc.)         
}

Everything else is the same. However, when you extract the MCMC samples of \(\phi\) you have to remember to multiple them by the scale parameter. It can be easier if you complete this step inside the Stan model; to do so, you can append the following generated quantity to the model:

generated quantities {
  vector[n] eta = phi * tau;  
}

NB: these ways of parameterizing a hierarchical model are sometimes referred to as ‘centered’ and ‘non-centered’ parameterizations (these are good search terms to learn more about it). We use the term ZMP here because it seems more intuitive: it is obvious which parameterization of the CAR model has mean of zero; it is not obvious which of these two is the ‘non-centered’ parameterization.