---
title: "Models for areal data"
subtitle: "Lecture 18"
author: "Dr. Colin Rundel"
footer: "Sta 344/644 - Fall 2023"
format:
  revealjs:
    theme: slides.scss
    transition: fade
    slide-number: true
    self-contained: true
    html-math-method: mathjax
    include-in-header:
      - text: |
          <script>
          MathJax = {
            options: {
              menuOptions: {
                settings: {
                  assistiveMml: false
                }
              }
            }
          };
          </script>
execute:
  echo: true
  warning: false
  collapse: true
---

```{r setup}
#| include: false
library(dukestm)
library(sf)
library(raster)

library(dplyr)
library(ggplot2)
library(patchwork)

knitr::opts_chunk$set(
  fig.align = "center"
)

set.seed(20221109)

ggplot2::theme_set(ggplot2::theme_bw())
```

# areal / lattice data

## Example - NC SIDS

```{r echo=FALSE}
nc = st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE) |> 
  select(-(AREA:CNTY_ID), -(FIPS:CRESS_ID)) |>
  st_transform(st_crs("+proj=utm +zone=17 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"))

plot(nc[,"SID79"], axes=FALSE, graticule=st_crs(4326), las=1)
plot(st_centroid(st_geometry(nc)), pch=16, add=TRUE)
```


## Adjacency Matrix

```{r}
1*st_touches(nc[1:12,], sparse=FALSE)
```


## Normalized spatial weight matrix

::: {.small}
```{r}
dukestm::normalize_weights( st_touches(nc[1:12,], sparse=FALSE) )
```
:::


## EDA - Moran's I

If we have observations at $n$ spatial locations $(s_1, \ldots s_n)$

$$ I = \frac{n}{\sum_{i=1}^n \sum_{j=1}^n w_{ij}} \frac{\sum_{i=1}^n \sum_{j=1}^n w_{ij} \big(y(s_i)-\bar{y}\big)\big(y(s_j)-\bar{y}\big)}{\sum_{i=1}^n \big(y(s_i) - \bar{y}\big)^2} $$ 

where $\boldsymbol{w}$ is a normalized spatial weight matrix.

. . .

Some properties of Moran's I when there is no spatial autocorrelation / dependence:

* $E(I) = -1 / (n-1)$

* $Var(I) = (\text{Something ugly but closed form})- E(I)^2$

* $\underset{n\to\infty}{\lim} \frac{I - E(I)}{\sqrt{Var(I)}} \sim N(0,1)$ (via the CLT)


## NC SIDS & Moran's I

Lets start by using a normalized spatial weight matrix for $\boldsymbol{w}$ (basedd on shared county borders).


```{r warning=FALSE, message=FALSE}
morans_I = function(y, w) {
  # w here can be either the adjacency matrix or the normalized weight matrix
  w = normalize_weights(w)
  n = length(y)
  num = sum(w * (y-mean(y)) %*% t(y-mean(y)))  
  denom = sum( (y-mean(y))^2 )
  (n/sum(w)) * (num/denom)
}

A = 1*st_touches(nc, sparse=FALSE)
morans_I(y = nc$SID74, A)
```

## 

```{r}
ape::Moran.I(nc$SID74, weight = A)
```



## EDA - Geary's C

Like Moran's I, if we have observations at $n$ spatial locations $(s_1, \ldots s_n)$

$$ C = \frac{n-1}{2\sum_{i=1}^n \sum_{j=1}^n w_{ij}} \frac{\sum_{i=1}^n \sum_{j=1}^n w_{ij} \big(y(s_i)-y(s_j)\big)^2}{\sum_{i=1}^n \big(y(s_i) - \bar{y}\big)^2} $$ 

where $\boldsymbol{w}$ is a normalized spatial weights matrix.


. . .

\vspace{7mm}

Some properties of Geary's C:

* $0 < C < 2$
    * If $C \approx 1$ then no spatial autocorrelation
    * If $C > 1$ then negative spatial autocorrelation
    * If $C < 1$ then positive spatial autocorrelation

* Geary's C is inversely related to Moran's I


## NC SIDS & Geary's C

Again using an adjacency matrix for $\boldsymbol{w}$ (shared county borders).

```{r message=FALSE, warning=FALSE}
gearys_C = function(y, w) {
  # w here can be either the adjacency matrix or the normalized weight matrix
  w = normalize_weights(w)
  
  n = length(y)
  y_i = y %*% t(rep(1,n))
  y_j = t(y_i)
  ((n-1)/(2*sum(w))) * (sum(w * (y_i-y_j)^2) / sum( (y - mean(y))^2 ))
}

A = 1*st_touches(nc, sparse=FALSE)
gearys_C(y = nc$SID74, w = A)
```


## Spatial Correlogram

:::: {.columns .medium}
::: {.column width='50%'}
Rather than using the touches predicate to determine adjacency, we can alternatively use a distance based approach.

In this case we will define the elements of the adjacency matrix as

$$
\{\boldsymbol A\}_{ij} = \mathbb1_{\{\text{distance between } s_i \text{ and } s_j \text{ is less than } d\}}
$$

We can then construct a correlogram by varying $d$ and then calculating Moran's I and Geary's C for each $\boldsymbol A$.
:::

::: {.column width='50%'}
```{r}
nc_pt = st_centroid(nc)
plot(nc_pt[,"SID74"], pch=16)
```
:::
::::



##

```{r echo=FALSE}
d = nc_pt |> st_distance() |> unclass()
breaks = seq(0, max(d), length.out = 31)
d_cut = cut(d, breaks)

adj_mats = purrr::map(
  levels(d_cut), 
  function(l) {
    m = matrix(d_cut == l, ncol=ncol(d), nrow=nrow(d))
    diag(m) = 0
    m
  }
)

d = tibble(
  dist   = breaks[-1],
  `Moran's I` = purrr::map_dbl(adj_mats, morans_I, y = nc$SID74),
  `Geary's C` = purrr::map_dbl(adj_mats, gearys_C, y = nc$SID74)
)

d |>
  tidyr::gather(var, value, -dist) |>
  mutate(var = forcats::as_factor(var)) |>
  ggplot(aes(x=dist/1000, y=value, color=var)) + 
    geom_line() + 
    facet_wrap(~var, scales = "free_y") +
    labs(x = "dist (km)") +
    guides(color="none")
```

## Inverse correlation of Moran's I and Geary's C

```{r echo=FALSE}
ggplot(d, aes(x=`Moran's I`, y=`Geary's C`)) +
  geom_point() + 
  geom_smooth(method="lm", se=FALSE, fullrange=TRUE)
```


# Autoregressive Models

## AR Models - Time

Lets return to the simplest case, an $AR(1)$ process

$$ y_t = \delta + \phi \, y_{t-1} + w_t $$

where $w_t \sim N(0,\sigma^2_w)$ and $|\phi| < 1$, then

:::: {.columns}
::: {.column width='50%'}
$$
\begin{aligned}
E(y_t) &= \frac{\delta}{1-\phi} \\
Var(y_t) &= \frac{\sigma^2}{1-\phi}
\end{aligned}
$$
:::

::: {.column width='50%'}
$$
\begin{aligned}
\gamma(h) &=  \phi^h \frac{\sigma^2}{1-\phi} \\
\rho(h) &= \phi^{h} 
\end{aligned}
$$
:::
::::




## AR Models - Time - Joint Distribution

We also saw that this $AR(1)$ model can be represented using a multivariate normal distribution,

$$
\begin{pmatrix}
y_1 \\ y_2 \\ \vdots \\ y_n
\end{pmatrix}
\sim N \begin{pmatrix}
\frac{\delta}{1-\phi} \begin{pmatrix}1\\ 1\\ \vdots\\ 1\end{pmatrix},~
\frac{\sigma^2}{1-\phi}
\begin{pmatrix}
1      & \phi   & \cdots & \phi^{n-1} \\
\phi   & 1      & \cdots & \phi^{n-2} \\
\vdots & \vdots & \ddots & \vdots     \\
\phi^{n-1} & \phi^{n-2}  & \cdots & 1 \\
\end{pmatrix}
\end{pmatrix}
$$

. . .



In writing down the likelihood we also saw that an $AR(1)$ is 1st order Markovian,

$$ \begin{aligned}
f(y_1, \ldots, y_n) 
  &= f(y_1) \, f(y_2 \mid y_1) \,  f(y_3 \mid y_2,y_1) \,\cdots\, f(y_n|y_{n-1},y_{n-2},\ldots,y_1) \\
  &= f(y_1) \, f(y_2 \mid y_1) \,  f(y_3 \mid y_2) \,\cdots\, f(y_n \mid y_{n-1})
\end{aligned} $$


## Alternative Definitions for $y_t$

$$ y_t = \delta + \phi \, y_{t-1} + w_t $$

::: {.center}
vs.
:::

$$ y_t\mid y_{t-1} \sim N(\delta + \phi \, y_{t-1},~\sigma^2) $$

. . .

In the case of time, both of these definitions result in the same multivariate distribution for $\boldsymbol{y}$ given on the previous slide.


## AR in Space

```{r echo=FALSE, fig.width=15, fig.height=3, out.width="100%"}
sq = st_polygon(list(matrix(c(0,0,0,1,1,1,1,0,0,0),ncol=2, byrow=TRUE)))

sqs = purrr::map(1:10, ~ sq + .*c(1,0)) |> st_sfc()

plot(sqs)
plot(st_centroid(sqs), add=TRUE, pch=16)
text( (sqs+c(0,-0.25)) |> st_centroid() |> st_coordinates(),labels=paste0("s",1:10), adj=c(0.5,0.5))
```

. . .

Even in the simplest spatial case there is no clear / unique ordering,

::: {.small}
$$
\begin{aligned}
f\big(y(s_1), \ldots, y(s_{10})\big) 
  &= f\big(y(s_1)\big) \, f\big(y(s_2) \mid y(s_1)\big) \, \cdots \, f\big(y(s_{10}\mid y(s_{9}),y(s_{8}),\ldots,y(s_1)\big)  \\
  &= f\big(y(s_{10})\big) \, f\big(y(s_9) \mid y(s_{10})\big) \, \cdots \, f\big(y(s_{1}\mid y(s_{2}),y(s_{3}),\ldots,y(s_{10})\big)  \\
  &= ~?
\end{aligned} 
$$
:::

. . .

Instead we need to think about things in terms of their neighbors / neighborhoods. We define $N(s_i)$ to be the set of neighbors of location $s_i$.

* If we define the neighborhood based on "touching" then $N(s_3) = \{s_2, s_4\}$

* If we include 2nd order neighbors, then $N(s_3) = \{s_1,s_2,s_3,s_4\}$


## Defining the Spatial AR model

Here we will consider a simple average of neighboring observations, just like with the temporal AR model we have two options in terms of defining the autoregressive process, 

. . .

* Simultaneous Autogressve (SAR)

$$ y(s) = \delta + \phi \frac{1}{|N(s)|}\sum_{s' \in N(s)} y(s') + N(0,\sigma^2) $$

. . .

* Conditional Autoregressive (CAR)

$$ y(s) \mid y(-s) \sim N\left(\delta + \phi \frac{1}{|N(s)|}\sum_{s' \in N(s)} y(s'),~ \sigma^2 \right) $$


# SAR

## Simultaneous Autogressve (SAR)

Using
$$ y(s) = \phi \frac{1}{|N(s)|}\sum_{s' \in N(s)} y(s') + N(0,\sigma^2) $$
we want to find the distribution of $\boldsymbol{y} = \Big(y(s_1),\, y(s_2),\,\ldots,\,y(s_n)\Big)^t$.

##

First we can define a weight matrix $\boldsymbol{W}$ where
$$ 
\{\boldsymbol{W}\}_{ij} = \begin{cases}
1/|N(s_i)| & \text{if $j \in N(s_i)$} \\
0        & \text{otherwise}
\end{cases}
$$
then we can write $\boldsymbol{y}$ as follows,
$$ \boldsymbol{y} = \phi \, \boldsymbol{W} \, \boldsymbol{y} + \boldsymbol{\epsilon} $$
where
$$ \boldsymbol{\epsilon} \sim N(0,\sigma^2 \, \boldsymbol{I}) $$

::: {.aside}
Note $\boldsymbol W$ here is the *normalized* weight matrix we mentioned previously.
:::


## A toy example

:::: {.columns}
::: {.column width='66%'}
![](imgs/triangle_adj.png){fig-align="center" width="90%"}
:::

::: {.column width='33%' .fragment}
<br/><br/>
$$ 
\begin{aligned}
\boldsymbol{A} &= \begin{pmatrix}
0 & 1 & 1 \\
1 & 0 & 1 \\
1 & 1 & 0 \\
\end{pmatrix} \\
\\
\boldsymbol{W} &= \begin{pmatrix}
0 & 1/2 & 1/2 \\
1/2 & 0 & 1/2 \\
1/2 & 1/2 & 0 \\
\end{pmatrix}
\end{aligned}
$$
:::
::::




## Back to SAR

$$ \boldsymbol{y} = \phi \, \boldsymbol{W} \, \boldsymbol{y} + \boldsymbol{\epsilon} $$

<!--
$$ \begin{aligned} 
\boldsymbol{y} &= \boldsymbol{\delta} + \phi \, \boldsymbol{W} \, \boldsymbol{y} + \boldsymbol{\epsilon} \\
\boldsymbol{y} - \phi \, \boldsymbol{W} \, \boldsymbol{y} &= \boldsymbol{\delta} + \boldsymbol{\epsilon} \\
(I-\phi \, \boldsymbol{W}) \, \boldsymbol{y} &= \boldsymbol{\delta} + \boldsymbol{\epsilon} \\
\boldsymbol{y} &= (I-\phi \, \boldsymbol{W})^{-1} \boldsymbol{\delta} + (I-\phi \, \boldsymbol{W})^{-1} \boldsymbol{\epsilon} \\
\end{aligned}$$

$$\begin{aligned}
E(\boldsymbol{y}) &= (I-\phi \, \boldsymbol{W})^{-1} \boldsymbol{\delta} \\
Var(\boldsymbol{y}) 
  &= \left((I-\phi \, \boldsymbol{W})^{-1}\right) \sigma^2 I \left((I-\phi \, \boldsymbol{W})^{-1}\right)^{t} \\
  &= \sigma^2 \left((I-\phi \, \boldsymbol{W})^{-1}\right) \left((I-\phi \, \boldsymbol{W})^{-1}\right)^{t} \\
\end{aligned}$$

$$ \boldsymbol{y} \sim N((I-\phi \, \boldsymbol{W})^{-1} \boldsymbol{\delta},~\sigma^2 \left((I-\phi \, \boldsymbol{W})^{-1}\right) \left((I-\phi \, \boldsymbol{W})^{-1}\right)^{t})$$
-->


# CAR


## Conditional Autogressve (CAR)

This is a bit trickier, in the case of the temporal AR process we actually went from joint distribution $\to$ conditional distributions (which we were then able to simplify).

. . .

Since we don't have a natural ordering we can't get away with this (at least not easily).

. . .

Going the other way, conditional distributions $\to$ joint distribution is difficult because it is possible to specify conditional distributions that lead to an improper joint distribution.


## Brooks' Lemma

For sets of observations $\boldsymbol{x}$ and $\boldsymbol{y}$ where $p(x) > 0~~\forall ~ x\in\boldsymbol{x}$ and $p(y) > 0~~\forall ~ y\in\boldsymbol{y}$ then

$$
\begin{aligned}
\frac{p(\boldsymbol{y})}{p(\boldsymbol{x})} 
  &= \prod_{i=1}^n \frac{p(y_i \mid  y_1,\ldots,y_{i-1},x_{i+1},\ldots,x_n)}{p(x_i \mid y_1,\ldots,y_{i-1},x_{i+1},\ldots,x_n)} \\
  &= \prod_{i=1}^n \frac{p(y_i \mid x_1,\ldots,x_{i-1},y_{i+1},\ldots,y_n)}{p(x_i \mid x_1,\ldots,x_{i-1},y_{i+1},\ldots,y_n)} \\
\end{aligned}
$$


## Derivation for $n=2$

Let $\boldsymbol{y} = (y_1,y_2)$ and $\boldsymbol{x} = (x_1,x_2)$ then we can derive Brook's Lemma for this case,

::: {.small}
$$
\begin{aligned}
p (y_1,y_2)  
    & \class{fragment}{= p(y_1 \mid y_2) p(y_2)} \\[3px]
    & \class{fragment}{= p(y_1 \mid y_2) \frac{p(y_2  \mid  x_1)}{p(x_1 \mid y_2)} p(x_1)} \\[3px]
    & \class{fragment}{= \frac{p(y_1 \mid y_2)}{p(x_1 \mid y_2)} p(y_2  \mid  x_1) \, p(x_1)} \\[3px]
    & \class{fragment}{= \frac{p(y_1 \mid y_2)}{p(x_1 \mid y_2)} p(y_2  \mid  x_1) \, p(x_1)   \left(\frac{p(x_2  \mid  x_1)}{p(x_2  \mid  x_1)}\right)} \\[3px]
    & \class{fragment}{= \frac{p(y_1 \mid y_2)}{p(x_1 \mid y_2)} \frac{p(y_2  \mid  x_1)}{p(x_2  \mid  x_1)} \, p(x_1,x_2)}
\end{aligned}
$$
:::

##

$$
\begin{aligned}
p (y_1,y_2)
    & = \frac{p(y_1 \mid y_2)}{p(x_1 \mid y_2)} \frac{p(y_2  \mid  x_1)}{p(x_2  \mid  x_1)} \, p(x_1,x_2) \\
  \\
  \frac{p (y_1,y_2) }{p(x_1,x_2)} 
    & = \frac{p(y_1 \mid y_2)}{p(x_1 \mid y_2)} \frac{p(y_2  \mid  x_1)}{p(x_2  \mid  x_1)}
\end{aligned}
$$

. . .

<br/>

From which we can generalize for $n=3$,

$$ 
\begin{aligned}
\frac{p(y_1,y_2,y_3)}{p(x_1,x_2,x_3)}
  = \frac{p(y_1 \mid y_2,y_3)}{p(x_1 \mid y_2,y_3)} \frac{p(y_2 \mid x_1,y_3)}{p(x_2 \mid x_1,y_3} \frac{p(y_3 \mid x_1,x_2)}{p(x_3 \mid x_1,x_2)}
\end{aligned} 
$$
and so on.


## Usefulness?

Lets repeat that last example but consider the case where $\boldsymbol{y} = (y_1,y_2)$ but now we let $\boldsymbol{x} = (y_1=0,y_2=0)$

$$
\begin{aligned}
\frac{p (y_1,y_2) }{p(x_1,x_2)} 
    &= \frac{p (y_1,y_2) }{p(y_1=0,y_2=0)}
\end{aligned}
$$

. . .

$$
\begin{aligned}
  p(y_1,y_2) &= \frac{p(y_1 \mid y_2)}{p(y_1=0 \mid y_2)} \frac{p(y_2  \mid  y_1=0)}{p(y_2=0 \mid y_1=0)} ~ p(y_1=0,y_2=0)
\end{aligned}
$$

. . .

$$
\begin{aligned}
  p(y_1,y_2) 
    &\propto \frac{p(y_1 \mid y_2) ~ p(y_2  \mid  y_1=0) }{ p(y_1=0 \mid y_2)} \\
%    &\propto \frac{p(y_2 \mid y_1) ~ p(y_1 \mid y_2=0) }{ p(y_2=0 \mid y_1)}
\end{aligned}
$$

## As applied to a **simple** CAR model

```{r}
#| echo: false
#| out-width: "50%"
sq = st_polygon(list(matrix(c(0,0,0,1,1,1,1,0,0,0),ncol=2, byrow=TRUE)))

sqs = purrr::map(1:2, ~ sq + .*c(1,0)) |> st_sfc()

par(mar=c(0.2,0.2,0.2,0.2))
plot(sqs)
plot(st_centroid(sqs), add=TRUE, pch=16)
text( (sqs+c(0,-0.25)) |> st_centroid() |> st_coordinates(),labels=paste0("s",seq_along(sqs)), adj=c(0.5,0.5))
```

$$
\begin{aligned}
y(s_1)\mid y(s_2) \sim N(\phi W_{12}\, y(s_2),\, \sigma^2) \\
y(s_2)\mid y(s_1) \sim N(\phi W_{21}\, y(s_1),\, \sigma^2)
\end{aligned}
$$

##

::: {.medium}

$$
\begin{aligned}
p\big(y(s_1),y(s_2)\big) 
    &\propto \frac{p\big(y(s_1) \mid y(s_2)\big) ~   p\big(y(s_2) \mid y(s_1)=0\big)}{p\big(y(s_1)=0  \mid y(s_2)\big)}\\
%    
    &\class{fragment}{{} \propto 
      \frac{
        \exp\left(-\frac{1}{2\sigma^2}\left(y(s_1)-\phi \, W_{12} \,   y(s_2)\right)^2\right)
        \exp\left(-\frac{1}{2\sigma^2}\left(y(s_2)-\phi \, W_{21} \, 0\right)^2\right) 
      }{
        \exp\left(-\frac{1}{2\sigma^2}\left(0-\phi W_{12} y(s_2)\right)^2 \right)
      }}\\
%    
    &\class{fragment}{{} \propto \exp\left(-\frac{1}{2\sigma^2}\left(\left(y(s_1)-\phi \, W_{12} \,   y(s_2)\right)^2 + y(s_2)^2- (\phi W_{21} y(s_2))^2\right)\right)} \\
%    
    &\class{fragment}{{} \propto \exp\left(-\frac{1}{2\sigma^2}\left(y(s_1)^2-\phi \, W_{12} \,   y(s_1)\,y(s_2) -\phi \, W_{21} \,   y(s_1)\,y(s_2) + y(s_2)^2\right)\right)} \\
%    
    &\class{fragment}{{} \propto \exp\left(-\frac{1}{2\sigma^2} (\boldsymbol{y}-0)
      \begin{pmatrix} 
      1 & -\phi W_{12} \\
      -\phi W_{21} & 1
      \end{pmatrix}
      (\boldsymbol{y}-0)^{t}
    \right)}
\end{aligned}
$$
:::


## Implications for $\boldsymbol{y}$

$$ \boldsymbol{\mu} = 0 $$

. . .

$$
\begin{aligned}
\boldsymbol{\Sigma}^{-1} &= \frac{1}{\sigma^2}
  \begin{pmatrix} 
    1 & -\phi W_{12} \\
    -\phi W_{21} & 1
  \end{pmatrix} \\
  &= \frac{1}{\sigma^2}(\boldsymbol{I} - \phi \, \boldsymbol{W})
\end{aligned}
$$

. . .


$$
\boldsymbol{\Sigma} = \sigma^2(\boldsymbol{I} - \phi \, \boldsymbol{W})^{-1}
$$

. . .

we can then conclude that for $\boldsymbol{y} = (y(s_1),~y(s_2))^t$,

$$
\boldsymbol{y} \sim N\left(
\boldsymbol{0}, ~
\sigma^2 (\boldsymbol{I} - \phi \, \boldsymbol{W})^{-1}
\right)
$$

which generalizes for all mean $0$ CAR models.


## General Proof

Let $\boldsymbol{y} = (y(s_1),\ldots,y(s_n))$ and $\boldsymbol{0} = (y(s_1) = 0, \ldots, y(s_n)=0)$ then by Brook's lemma,


::: {.small}
$$
\begin{aligned}
\frac{p(\boldsymbol{y})}{p(\boldsymbol{0})} 
  &= \prod_{i=1}^n \frac{p(y_i \mid y_1,\ldots,y_{i-1},0_{i+1},\ldots,0_{n})}{p(0_i \mid y_1,\ldots,y_{i-1},0_{i+1},\ldots,0_{n})}
  = \prod_{i=1}^n 
    \frac{
      \exp\left(-\frac{1}{2\sigma^2} \left(y_i - \phi \sum_{j<i} W_{ij} \, y_j - \phi \sum_{j>i} 0_j \right)^2 \right)
    }{
      \exp\left(-\frac{1}{2\sigma^2} \left(0_i - \phi \sum_{j<i} W_{ij} \, y_j - \phi \sum_{j>i} 0_j \right)^2 \right)
    } \\
  &= \exp\left(-\frac{1}{2\sigma^2} \sum_{i=1}^n \left(y_i - \phi \sum_{j<i} W_{ij} \, y_j\right)^2 + \frac{1}{2\sigma^2} \sum_{i=1}^n \left( \phi \sum_{j<i} W_{ij} \, y_j \right)^2 \right) \\
  &= \exp\left(-\frac{1}{2\sigma^2} \sum_{i=1}^n y_i^2 - 2 \phi y_i \,\sum_{j<i} W_{ij} \, y_j \right) \\
  &= \exp\left(-\frac{1}{2\sigma^2} \sum_{i=1}^n y_i^2 - \phi \sum_{i=1}^n \sum_{j=1}^n y_i \, W_{ij} \, y_j \right) \quad \mathit{\big(\text{if } W_{ij} = W_{ji}\big)} \\
  &= \exp\left(-\frac{1}{2\sigma^2} (\boldsymbol{y}-0)^t (\boldsymbol{I} - \phi \boldsymbol{W}) (\boldsymbol{y}-0)  \right)
\end{aligned}
$$
:::


## CAR vs SAR

* Simultaneous Autogressve (SAR)

$$ y(s) = \phi \sum_{s'} W_{s\;s'} \, y(s') + \epsilon $$

$$ \boldsymbol{y} \sim N(0,~\sigma^2 \, ((\boldsymbol{I}-\phi \boldsymbol{W})^{-1}) ((\boldsymbol{I}-\phi \boldsymbol{W})^{-1})^t )$$

* Conditional Autoregressive (CAR)

$$ y(s) \mid y(-s) \sim N\left(\sum_{s'} W_{s\;s'} \, y(s'),~ \sigma^2 \right) $$

$$ \boldsymbol{y} \sim N(0,~\sigma^2 \, (\boldsymbol{I}-\phi \boldsymbol{W})^{-1})$$

## Generalizations

* Adopting different weight matrices $(\boldsymbol{W})$
  
    * Between SAR and CAR model we move to a generic weight matrix definition (beyond average of nearest neighbors)
    
    * In time we varied $p$ in the $AR(p)$ model, in space we adjust the weight matrix.
    
    * In general having a symmetric W is helpful, but not required

. . .
    
* More complex Variance (beyond $\sigma^2 \, I$)
  
    * $\sigma^2$ can be a vector (differences between areal locations)
    
    * i.e. since areal data tends to be aggregated - adjust variance based on sample size
    
    * i.e. scale based on the number of neighbors
    
## Some specific generalizations

Generally speaking we want to work with a scaled / normalized version of the weight matrix,
$$ W_{ij}/W_{i\boldsymbol{\cdot}}  $$

When $W$ is derived from an adjacency matrix, $\boldsymbol{A}$, we can express this as 
$$ \boldsymbol{W} = \boldsymbol{D}^{-1} \boldsymbol{A} $$
where $\boldsymbol{D}^{-1} = \text{diag}(1/|N(s_i)|)$. 


We can also allow $\sigma^2$ to vary between locations, we can define this as $\boldsymbol{D}_{\sigma^2} = \text{diag}(\sigma^2_i)$ and most often we use
$$ \boldsymbol{D}_{\sigma^2}^{-1} = \text{diag}\left(\frac{\sigma^2}{|N(s_i)|}\right) =  \sigma^2 \boldsymbol{D}^{-1}.  $$

## Revised SAR Model

* Formula Model

$$
y(s_i) = X_{i\cdot}\beta + \phi \sum_{j=1}^n D^{-1}_{jj} \, A_{ij} \, \big(y(s_j) - X_{j\cdot}\beta\big) + \epsilon_i
$$
$$ 
\boldsymbol{\epsilon} \sim N(\boldsymbol{0},\,\boldsymbol{D}_{\sigma^2}^{-1}) =  N(\boldsymbol{0},\, \sigma^2 \boldsymbol{D}^{-1}) 
$$

. . .

* Joint Model

$$
\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \phi \boldsymbol{D}^{-1} \boldsymbol{A} ~\big(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\big) + \boldsymbol{\epsilon}
$$

. . .

$$
\boldsymbol{y} \sim N\left(\boldsymbol{X}\boldsymbol{\beta}, (\boldsymbol{I} - \phi \boldsymbol{D}^{-1} \boldsymbol{A})^{-1} \sigma^2 \boldsymbol{D}^{-1} \big((\boldsymbol{I} - \phi \boldsymbol{D}^{-1} \boldsymbol{A})^{-1}\big)^t \right)
$$



## Revised CAR Model

* Conditional Model

:::  {.medium}
$$ y(s_i) \mid y_{-s_i} \sim N\left(X_{i\cdot}\beta + \phi\sum_{j=1}^n \frac{W_{ij}}{D_{ii}} ~ \big(y(s_j)-X_{j\cdot}\beta\big),~ \sigma^2 D^{-1}_{ii} \right) $$
:::

. . .

* Joint Model

::: {.medium}
$$
\boldsymbol{y} \sim N(\boldsymbol{X}\boldsymbol{\beta},~\Sigma_{CAR})
$$
:::

. . .

::: {.medium}
$$
\begin{aligned}
\Sigma_{CAR}
  &= (\boldsymbol{D}_{\sigma} \, (\boldsymbol{I}-\phi \boldsymbol{D}^{-1}\boldsymbol{A}))^{-1} \\
  &= (1/\sigma^2 \boldsymbol{D} \, (\boldsymbol{I}-\phi \boldsymbol{D}^{-1}\boldsymbol{A}))^{-1} \\
  &= (1/\sigma^2 (\boldsymbol{D}-\phi \boldsymbol{A}))^{-1} \\
  &= \sigma^2(\boldsymbol{D}-\phi \boldsymbol{A})^{-1}
\end{aligned}
$$
:::




## Toy CAR Example

```{r echo=FALSE}
x = c(0,1,2)
y = c(0,1,0)

plot(x, y, pch=16, type="b", axes=FALSE,xlab="",ylab="", xlim=c(-0.2,2.2),ylim=c(-0.2,1.5))
text(x+c(-0.1,0,0.1),y+c(0,0.2,0), labels = c("s1","s2","s3"))
```

. . .

$$
\boldsymbol{A} = \begin{pmatrix}
0 & 1 & 0 \\
1 & 0 & 1 \\
0 & 1 & 0 
\end{pmatrix}
\qquad\qquad
\boldsymbol{D} = \begin{pmatrix}
1 & 0 & 0 \\
0 & 2 & 0 \\
0 & 0 & 1 
\end{pmatrix}
$$

. . .

$$
\boldsymbol{\Sigma} = \sigma^2 \, (\boldsymbol{D} - \phi \, \boldsymbol{W})^{-1} = \sigma^2~\begin{pmatrix}
1 & -\phi & 0 \\
-\phi & 2 & -\phi \\
0 & -\phi & 1 
\end{pmatrix}^{-1}
$$

## When does $\Sigma$ exist?

::: {.small}
```{r error=TRUE}
check_sigma = function(phi) {
  Sigma_inv = matrix(c(1,-phi,0,-phi,2,-phi,0,-phi,1), ncol=3, byrow=TRUE) 
  solve(Sigma_inv)
}

check_sigma(phi=0)

check_sigma(phi=0.5)

check_sigma(phi=-0.6)
```
:::

##

::: {.small}
```{r error=TRUE}
check_sigma(phi=1)

check_sigma(phi=-1)

check_sigma(phi=1.2)

check_sigma(phi=-1.2)
```
:::

## When is $\Sigma$ positive semidefinite?

::: {.small}
```{r error=TRUE}
check_sigma_pd = function(phi) {
  Sigma_inv = matrix(c(1,-phi,0,-phi,2,-phi,0,-phi,1), ncol=3, byrow=TRUE) 
  chol(solve(Sigma_inv))
}

check_sigma_pd(phi=0)

check_sigma_pd(phi=0.5)

check_sigma_pd(phi=-0.6)
```
:::

##

::: {.small}
```{r error=TRUE}
check_sigma_pd(phi=1)

check_sigma_pd(phi=-1)

check_sigma_pd(phi=1.2)

check_sigma_pd(phi=-1.2)
```
:::


## Conclusions

Generally speaking just like the AR(1) model for time series we require that $|\phi| < 1$ for the CAR model to be proper.

These results for $\phi$ also apply in the context where $\sigma^2_i$ is constant across locations, i.e. 
$$
\boldsymbol{\Sigma} = \big(\sigma^2 \, (\boldsymbol{I}-\phi \boldsymbol{D}^{-1}\boldsymbol{A})\big)^{-1}
$$

. . .

As a side note, the special case where $\phi=1$ is known as an intrinsic autoregressive (IAR) model and they are popular as an *improper* prior for spatial random effects. An additional sum constraint is necessary for identifiability 
$$\sum_{i=1}^n y(s_i) = 0$$