ICAR Spatial Model with Hierarchical Data

Modeling
, ,
1

Hi, I am working on building a spatial model using an ICAR component ala this case study. However, all the tutorials on this topic I could find use data with only one observation per areal unit. I’m looking to build a hierarchical model where I have multiple observations per areal unit, while incorporating the ICAR priors (or maybe BYM2).

I am mostly confused on how to build the adjacency matrix, and how to properly index the data for such a model.

Here is my first attempt:

data {
  int<lower=0> N;
  int<lower=1> K; // K counties
  array[N] int<lower=1, upper=K> county;
  int<lower=0> N_edges;
  array[N_edges] int<lower=1, upper=N> node1; 
  array[N_edges] int<lower=1, upper=N> node2; 
  vector[N] y;
}
parameters {
  vector[N] phi;
  array[K] real alpha_county;
  real<lower=0> sigma;
}
model {
  target += -0.5 * dot_self(phi[node1] - phi[node2]);
  
  sum(phi) ~ normal(0, 0.01 * N);

  to_vector(alpha_county) ~ normal(0, 5);
  sigma ~ exponential(1);
  
  for (n in 1:N)
    y[n] ~ normal(alpha_county[county[n]], sigma);
}

Using simulated data from the geostan package:

library(geostan)
library(dplyr)
library(rstan)

data("georgia")

georgia_long = georgia %>%
  select(GEOID, income, college) %>% 
  # expand the data to include multiple obs per county
  tidyr::crossing(x = 1:10) %>% 
  mutate(x = rnorm(n(), mean = college),
         .by = GEOID) %>% 
  # not all counties will have the same no. of observations
  slice_sample(prop = .75) %>% 
  sf::st_as_sf()

C <- shape2mat(georgia_long, "B")
nodes <- prep_icar_data(C)

icar_stan = stan_model("icar.stan")
icar_fit = sampling(
  icar_stan, 
  data = list(
    N = nrow(georgia_long),
    K = length(unique(georgia_long$GEOID)),
    county = as.numeric(as.factor(georgia_long$GEOID)),
    N_edges = nodes$n_edges,
    node1 = nodes$node1,
    node2 = nodes$node2,
    y = georgia_long$income
  ), 
  chains = 4,
  cores = 4,
  iter = 1000
)

Should the adjacency matrix be constructed using the “long” data like I have done it here? Or should I use a county-level adjacency matrix and index the nodes by county?

Thanks!

2

I don’t have time to carefully engage with your Stan code, but in case it’s clarifying, you can think about your model as follows.

Step 1: Build an ICAR model for the county averages as if you have only one observation per county. This will look exactly like the case studies with one observation per county, except that instead of declaring the outcome as data, you will declare the outcome as a parameter.

Step 2: Build a data model that gives the likelihood of your observed data conditional on the county average from step 1.

3 Likes
3

Thanks, that definitely helps clarify things!

I am still not comfortable enough with Stan syntax to know how to translate the ICAR prior target += -0.5 * dot_self(phi[node1] - phi[node2]); into something that is operating at the county-level as opposed to observation-level.

1 Like
4

Forget about the observation level entirely for a moment. Make up some data, if you like, that have a single observation per county, and follow the case studies to get a model for these. Suppose you call those data mu. Now, remove your declaration of mu as data, and instead declare mu as a parameter. Now you’re done with the ICAR part, and you can put target += -0.5 * dot_self(phi[node1] - phi[node2]); out of mind. All you need to do now is write down a likelihood for the observed data in each county conditional on mu.

5

Here’s some example code you can use. Its not using the ICAR but I think you might get the idea. I haven’t tried to run this code, its just adapted from something else I was working on. Other Stan code for ‘varying intercepts’/‘random effects’ might help too, you’re just adding autocorrelation to them.

Or should I use a county-level adjacency matrix and index the nodes by county?

Yeah, the trick is to create an index that connects each individual to a geographic area.

So my code is kind of like this:

# create cross-walk from geoid to index position
geo_index = data.frame(GEOID =  geo_df$GEOID,  geo_index = 1:nrow(geo_df))

# join to your individual-level observations 
dat <- merge(dat, geo_index, by = 'GEOID')

# create your data list for stan, respecting order
A <- shape2mat(geo_df)
dl <- prep_car_data(A)
dl$y <- dat$y
dl$geo_index = dat$geo_index
dl$N = nrow(dat)
dl$N_geo = nrow(geo_df)
...

You use the index like this in Stan:

  transformed parameters {
  vector[N] mu;
  for (i in 1:N) mu[i] = phi[geo_index[i]] + alpha;
  mu += X * beta;
}

Here’s the example:


// hiererchical model: individuals nested in geographies

functions {
  // spatial model
   //  https://connordonegan.github.io/geostan/articles/custom-spatial-models.html
  #include functions.stan
}

data {
  int<lower=1> N;
  int<lower=0> K;
  vector[N] y;
  matrix[N, K] X;
  
  // geo identifier
  array[N] int geo_index;    

 // size of geog
  int N_geo;

  // CAR parts
  int nAx_w;
  int nC;
  vector[nAx_w] Ax_w;
  array[nAx_w] int Ax_v;
  array[N_geo + 1] int Ax_u;
  array[nC] int Cidx;
  vector[N_geo] Delta_inv;
  real log_det_Delta_inv;
  vector[N_geo] lambda;

}

parameters {
  // observation-level regression
  vector[K] beta;
  real alpha;
  real<lower=0> sigma;  

  // CAR 
  vector[N_geo] phi;  
  real<lower=0> sigma_phi;
  real<lower=1/min(lambda), upper=1/max(lambda)> rho_phi;  
}

transformed parameters {
  vector[N] mu;
  for (i in 1:N) mu[i] = phi[geo_index[i]] + alpha;
  mu += X * beta;
}

  //
  // alternative is to drop alpha from the equation for mu and put it in phi:
  // vector[N_geo] mu_phi = rep_vector(alpha, N_geo);
  // you can have covariates at this level also, like: matrix[N_geo, K_geo] X_geo;
  //

model {
  // phi centered on zero:
  vector[N_geo] mu_phi = rep_vector(0, N_geo);

  // likelihood
  target += normal_lpdf(y | mu, sigma);
  
  // regression priors
  target += normal_lpdf(beta0 | 0, 25);
  target += normal_lpdf(sigma | 0, 5);
  target += normal_lpdf(beta | 0, 5);
  
   // CAR: phi ~ Normal(mu_phi, Sigma)
  target += wcar_normal_lpdf(phi |
                 mu_phi, sigma_geo, rho_geo, 
                 Ax_w, Ax_v, Ax_u, 
                 Delta_inv, 
                 log_det_Delta_inv,
                 lambda, 
                 N_geo);    
  
   // CAR model priors [implicit uniform on rho in its support]
  target += student_t_lpdf(sigma_geo | 10, 0, 5);
}

generated quantities {
  vector[N] log_lik;
    for (i in 1:N_obs) log_lik[i] = normal_lpdf(y[i] | mu[i], sigma);    
}
6

Thanks for your help @jsocolar and @cmcd! I think I now have something that works. Following @jsocolar’s advice I attempted to model the BYM2 part (comprised of the spatial component and the unstructured component) on its own, and then condition on those parameters in the observation-level model. Rather than calling those parameters mu, I call them convolved_re to match the case study. But I hope they reflect the intuition you were trying to provide!

data {
  int<lower=1> N; // no. observations
  int<lower=1> K; // no. counties
  array[N] int<lower=1, upper=K> county;
  int<lower=1> N_edges; // no. edges created from county-level adjacency matrix
  array[N_edges] int<lower=1, upper=K> node1;
  array[N_edges] int<lower=1, upper=K> node2; 
  
  vector[N] x; // observation-level predictor
  vector[N] y; // outcome
  
  real<lower=0> scaling_factor;
}
parameters {
  real<lower=0, upper=1> rho;
  vector[K] phi; // spatial component
  vector[K] theta; // unstructured component
  
  real beta; // fixed effect
  real alpha; // global intercept
  
  vector[K] county_alpha;
  real<lower=0> county_sigma;
  
  real<lower=0> sigma; // residual sd
}
transformed parameters {
  vector[K] convolved_re;
  convolved_re = sqrt(1 - rho) * theta + sqrt(rho / scaling_factor) * phi;
}
model {
  rho ~ beta(0.5, 0.5);
  target += -0.5 * dot_self(phi[node1] - phi[node2]);
  sum(phi) ~ normal(0, 0.01 * K);
  theta ~ normal(0, 2);
  
  beta ~ normal(0, 2);
  alpha ~ normal(0, 2);
  sigma ~ normal(0, 1); 

  county_sigma ~ normal(0, 1); // hyperprior for hierarchical term
  
  for (k in 1:K) {
    county_alpha[k] ~ normal(convolved_re[k], county_sigma);
  }
  
  y ~ normal(alpha + county_alpha[county] + x * beta, sigma);
}