Hello all,

I am implementing the Besag York Mollié (BYM) spatial model as specified by Mitzi Morris in the iCAR case-study. https://mc-stan.org/users/documentation/case-studies/icar_stan.html

There is plenty of background at that link, but in general it is a Poisson GLM with both an ordinary random effects, as well as a spatial random effect via an iCAR prior. The code for the BYM2 model is included below.

My question is once a model with covariates is fit to an area (ex. voting districts of NYC), how can this model be predicted to a new area (ex. voting districts of Boston given the parameters distributions of the NYC model)? For non-spatial models, I have done this either as generated quantities or with posterior samples, but I am unsure how to do it when there is a spatial weights matrix involved.

If I try to do this in the generated quantities block, I can easily supply `N_test`

, `y_test`

, `E_test`

, and `x_test`

data to represent the predicted geography. However, the convolution prior is a vector if length `N`

; the training data. Since I cannot sample a model with parameters vectors of both length `N`

and `N_test`

, I cannot predict to a new area based on the `phi`

, `rho`

, or `theta`

parameters of the initial model. Perhaps the answer is that I need to re-fit the entire ordinary random and spatial random effects components with the spatial weight matrix of the new area (e.g. Boston) and use the `betas`

from the initial model fit (e.g. NYC)?

Any thoughts or help would be appreciated. Apologies if this is incomplete, unclear, trivial.

Thank you,

Matt

```
data {
int<lower=0> N;
int<lower=0> N_edges;
int<lower=1, upper=N> node1[N_edges]; // node1[i] adjacent to node2[i]
int<lower=1, upper=N> node2[N_edges]; // and node1[i] < node2[i]
int<lower=0> y[N]; // count outcomes
vector<lower=0>[N] E; // exposure
int<lower=1> K; // num covariates
matrix[N, K] x; // design matrix
real<lower=0> scaling_factor; // scales the variance of the spatial effects
}
transformed data {
vector[N] log_E = log(E);
}
parameters {
real beta0; // intercept
vector[K] betas; // covariates
real<lower=0> sigma; // overall standard deviation
real<lower=0, upper=1> rho; // proportion unstructured vs. spatially structured variance
vector[N] theta; // heterogeneous effects
vector[N] phi; // spatial effects
}
transformed parameters {
vector[N] convolved_re;
// variance of each component should be approximately equal to 1
convolved_re = sqrt(1 - rho) * theta + sqrt(rho / scaling_factor) * phi;
}
model {
y ~ poisson_log(log_E + beta0 + x * betas + convolved_re * sigma); // co-variates
// This is the prior for phi! (up to proportionality)
target += -0.5 * dot_self(phi[node1] - phi[node2]);
beta0 ~ normal(0.0, 2.5);
betas ~ normal(0.0, 2.5);
theta ~ normal(0.0, 1.0);
sigma ~ normal(0,5);
rho ~ beta(0.5, 0.5);
// soft sum-to-zero constraint on phi)
sum(phi) ~ normal(0, 0.001 * N); // equivalent to mean(phi) ~ normal(0,0.001)
}
generated quantities {
real logit_rho = log(rho / (1.0 - rho));
vector[N] eta = log_E + beta0 + x * betas + convolved_re * sigma; // co-variates
vector[N] mu = exp(eta);
}
```