Specifying Bayesian Spatially lagged dependent (SAR) variable regression in winBUGS or Stan

Question

I have been looking for code example on how to specify a SAR spatially lagged dependent variable regression in winBUGS (geoBUGS) or, if possible, in Stan.

I mean, something of the format:

Y = beta_0 + beta_1*Y*W + ... + e

Where W is the weight matrix identifying neighborhood of geographical areas.

Notice that I am not looking for a CAR model, i.e. when the errors are spatially auto correlated, but specifically for a SAR lag model.

Can anyone help with some code example or any source suggestion?

Answer 1

I am assuming that you would like an AR(II) process term to account for serial autocorrelation.

Here is example WinBUGS code for running a Simultaneous Autoregressive model with an AR(II) term that could be adapted for your problem.

model;
{
# Defines likelihood
for(i in 1:I){ D[i, 2:(I+2)] ~ dmulti(C[i,], D[i, 1]); }

for(i in 1:(I-1)){
    lphi[i] <- log(phi[i])
    logit(phi[i]) <- beta + e[i]
    for(j in (i+1):I){
            C[i, j] <- lambda[j]*exp(sum(lphi[i:(j-1)]))
    }
    for (j in 1:i){
        C[i+1, j] <- 0
    }
}
for(i in 1:I){
    C[i, i] <- lambda[i]
    C[i, I+1] <- 1 - sum(C[i, 1:I])
}

# Defines epsilon    
e[1] ~ dnorm(mu[1], tau1)
e[2] ~ dnorm(mu[2], tau2)
mu[1] <- 0
mu[2] <- (rho[1]/(1-rho[2]))*e[1]
tau1 <- ((1+rho[2])/(1-rho[2]))*((1-rho[2])*(1-rho[2]) - rho[1]*rho[1])*tau
tau2 <- tau*(1 - rho[2]*rho[2])
for(i in 3:(I-1)){
    e[i] ~ dnorm(mu[i],tau)
    mu[i] <- rho[1]*e[i-1] + rho[2]*e[i-2]
}
sigma <- 1/sqrt(tau)

#  Prior distribution
for(i in 1:2){beta[i] ~ dnorm(0, 0.01)}
beta ~ dnorm(0, 0.01)
for(i in 1:I){lambda[i] ~ dunif(0, 1)}
tau ~ dgamma(0.001, 0.001)

## Prior for rho is approx (uniform on the AR(II) stationary triangle)
rho[1] ~ dunif(l, u)
u <- abs(1 - rho[2])
l <- -u
rho[2] ~ dunif(-1, 1)    
}

Answer 2

Here is a specification in RJAGS, which you can probably port over to winbugs or stan pretty easily. It is taken from these lecture notes by Colin Rundel. The accompanying slides can be found here. Note that the model is incorrect on the lecture pdf, (and is crossed out in the annotated pdf), but is different in the rmd file.

```{r echo=FALSE}
sar_model = "model{
  C <- (I - phi * W_tilde)
  Sigma_inv <- tau * t(C) %*% D %*% C
  mu <- beta0 + beta1*x

  y ~ dmnorm(mu, Sigma_inv)
  y_pred ~ dmnorm(mu, Sigma_inv)

  beta0 ~ dnorm(0, 1/100)
  beta1 ~ dnorm(0, 1/100)

  tau <- 1 / sigma2
  sigma2 ~ dnorm(0, 1/100) T(0,)
  phi ~ dunif(-0.99, 0.99)
}"
cat(sar_model,"\n")
```

```{r}
D_inv = diag(1/diag(D))
W_tilde = D_inv %*% W
I = diag(1, ncol=length(y), nrow=length(y))
```