Clarification on CAR Pairwise Difference Derivation

In the CAR class of spatial models, I’m a little confused as to how the pairwise difference was derived. In the case study, the section about deriving the pairwise difference from the matrix algebra says,

-\frac{1}{2}\big(\phi^{\top}[D-W]\phi\big) + const = -\frac{1}{2}\sum_{i\sim j} (\phi_i - \phi_j)^2\ + const

But if I simulate the math in R using a small artificial lattice (code below) I get the following result (omitting constants and negatives),

\phi^{\top}[D-W]\phi = \frac{1}{2}\sum_{i\sim j} (\phi_i - \phi_j)^2

which says that the matrix algebra result is half of what the pairwise difference computes.

N <- 25
B <- matrix(rbinom(N, 1, 0.5), nrow = sqrt(N))
B[upper.tri(B)] <- 0
B_lwr <- B
B <- B + t(B)
diag(B) <- 0
isSymmetric(B)

# mean vector
m <- rnorm(sqrt(N), 1, 1)

# find neighbors
E <- which(B == 1, arr.ind = TRUE)

# GMRF matrix multiply
t(m) %*% (diag(rowSums(B)) - B) %*% m
# GMRF matrix multiply simplified to pairwise difference
diff <- m[E[,1]] - m[E[,2]]
0.5 * sum(diff^2)

I think I’m misunderstanding something since both the case study and this paper (bottom of p3) have similar mathematical results. Any idea where I’m going wrong in my thinking?

So the 2x2 example is:

D - W = [[1, -1], [-1, 1]]

The terms on the diagonal are going to be the squared terms, and then the off diagonal will be the mixed terms.

So you get x_1^2 + x_2^2 from the diagonal, and then -2 x_1 x_2, one from the upper off diagonal, and one for the lower.

There’s only one term in the sum for two elements, (x_1 - x_2)^2. So in that case the matching -\frac{1}{2} terms seem fine.

Oh, E is double counting.

> E
      row col
 [1,]   2   1
 [2,]   3   1
 [3,]   4   1
 [4,]   5   1
 [5,]   1   2 // duplicate
 [6,]   5   2
 [7,]   1   3 // duplicate
 [8,]   1   4 // duplicate
 [9,]   1   5 // duplicate
[10,]   2   5 // duplicate

edit: fixed what is double counting

1 Like

Wow *face-palm*. Thanks for checking the indices and finding the duplicates Ben.
Now the math matches if you find the indices on the upper/lower triangular matrix of B provided that the diagonal of B is 0.