Clarification on CAR Pairwise Difference Derivation

Modeling
#1

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?

#2

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.

#3

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
#4

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.