Intrinsic Gaussian CAR model using NIMBLE

CAR_nimble.Rmd

---
title: "R Notebook"
output: html_notebook
---

# Intrinsic Gaussian CAR model using NIMBLE

```{r setup}
library(readr)
library(nimble)
library(coda)
library(stringr)
library(ggplot2)
```

## Data

Distribution Data of Quercus glauca from Iwanami Data Science Vol. 1.

```{r}
data_source <- url("https://raw.githubusercontent.com/iwanami-datascience/vol1/master/ito/Qglauca.csv")
qgl <- read_csv(data_source, col_types = "ddi")
```

### Plot

```{r}
ggplot(qgl) +
  geom_tile(aes(x = X - 0.5, y = Y - 0.5, fill = N)) +
  labs(x = "X", y = "Y") +
  scale_y_reverse(breaks = c(0, 5, 10)) +
  coord_fixed()
```

### Geometry and adjacent list

```{r}
## Adjacent matrix
n.x <- max(qgl$X)
n.y <- max(qgl$Y)
n.sites <- n.x * n.y

num <- c(2, rep(3, n.y - 2), 2,
         rep(c(3, rep(4, n.y - 2), 3), n.x - 2),
         2, rep(3, n.y - 2), 2)

adj <- rep(0, sum(num))
i <- 1
adj[i:(i + 1)] <- c(2, n.y + 1)
i <- i + 2
for (y in 2:(n.y - 1)) {
  adj[i:(i + 2)] <- c(y - 1, y + n.y, y + 1)
  i <- i + 3
}
adj[i:(i + 1)] <- c(n.y - 1, n.y * 2)
i <- i + 2
for (x in 1:(n.x - 2)) {
  adj[i:(i + 2)] <- c((x - 1) * n.y + 1, x * n.y + 2, (x + 1) * n.y + 1)
  i <- i + 3
  for (y in 2:(n.y - 1)) {
    adj[i:(i + 3)] <- c(x * n.y + y - 1, (x - 1) * n.y + y,
                        (x + 1) * n.y + y, x * n.y + y + 1)
    i <- i + 4
  }
  adj[i:(i + 2)] <- c(x * n.y, (x + 1) * n.y - 1, (x + 2) * n.y)
  i <- i + 3
}
adj[i:(i + 1)] <- c((n.x - 2) * n.y + 1, (n.x - 1) * n.y + 2)
i <- i + 2
for (y in 2:(n.y - 1)) {
  adj[i:(i + 2)] <- c((n.x - 1) * n.y + y - 1, (n.x - 2) * n.y + y, 
                      (n.x - 1) * n.y + y + 1)
  i <- i + 3
}
adj[i:(i + 1)] <- c((n.x - 1) * n.y, n.x * n.y - 1)
```

check the adjacent list

```{r}
p1 <- c(0, cumsum(num))[seq_along(num)] + 1
p2 <- cumsum(num)
for (i in seq_along(num))
  print(adj[p1[i]:p2[i]])
```



## NIMBLE model

```{r}
mc <- nimbleCode({
  for (n in 1:N) {
    Y[n] ~ dpois(lambda[n]);
    log(lambda[n]) <- S[n];
  }
  S[1:N] ~ dcar_normal(adj = ADJ[], num = NUM[], tau = tau);
  sigma ~ dunif(0, 10);
  tau <- 1 / (sigma * sigma);
})
```

## Fitting

```{r fit}
nimble_const <- list(N = nrow(qgl))
nimble_data <- list(Y = qgl$N,
                    ADJ = adj,
                    NUM = num)
nimble_init <- function() {
  list(sigma = runif(1, 0, 2),
       S = runif(nrow(qgl), -1, 1))}
nimble_fit <- nimbleMCMC(code = mc,
                         constants = nimble_const,
                         data = nimble_data,
                         inits = nimble_init,
                         niter = 20000,
                         nburnin = 10000,
                         thin = 10,
                         nchains = 3,
                         monitors = c("sigma", "lambda", "S"),
                         setSeed = 1:3,
                         samplesAsCodaMCMC = TRUE)
```

## Results

```{r results}
plot(nimble_fit[, c(1:2, 401)])
```

R-hat

```{r}
gelman.diag(nimble_fit[, "sigma"])

# maximum in S[]
max(gelman.diag(nimble_fit[, str_starts(colnames(nimble_fit$chain1), "S\\[")])$psrf)
```

## Plot

```{r}
pos_lambda <- str_starts(colnames(nimble_fit$chain1), "lambda")
qgl2 <- cbind(qgl, lambda = summary(nimble_fit)$statistics[pos_lambda, "Mean"])
ggplot(qgl2) +
  geom_tile(aes(x = X - 0.5, y = Y - 0.5, fill = lambda)) +
  labs(x = "X", y = "Y") +
  scale_y_reverse(breaks = c(0, 5, 10)) +
  coord_fixed()
```