Custom distance matrices for spatial autocorrelation?

[BertvanderVeen]

Discussed in #216

^{Originally posted by dhercar March 29, 2025}
Thanks so much for developing this great package and all the recent updates.

I am currently working with river communities (fish, algae, macroinvertebrates, etc. ) across relatively large spatial scales, and I would like to incorporate spatial autocorrelation into our models. However, geographic distance is not very meaningful for taxa groups that can only move within the river network, such as fish. Would it be possible to make dist and distLV take a custom distance matrix directly, instead of spatial coordinates? This might also be useful for users dealing with spatial barriers.

Many thanks!
Daniel

[wlangera]

I also agree with allowing user-specified distance matrices directly for the dist and distLV arguments. This would:

1. Avoids assumptions about coordinate reference systems and helps prevent misuse of unprojected (longitude/latitude) coordinates. Explain how distances are calculated in gllvm() #231
2. Allows integration of custom distance metrics such as river-network distances, resistance surfaces, or ecological connectivity. Custom distance matrices for spatial autocorrelation? #217
3. Makes it possible to model non-metric spatial structures, such as adjacency or nearest-neighbor weights (e.g. binary or ordinal weights).

Perhaps a vignette could explain how to calculate distance matrices in R:

Examples of custom distance or weight matrices in R

1. Euclidean distance from projected coordinates

coords <- data.frame(
    x = c(100, 110, 120),
    y = c(200, 210, 190)
)

# Calculate Euclidean distance matrix
dist_mat <- as.matrix(dist(coords))
dist_mat
#>          1        2        3
#> 1  0.00000 14.14214 22.36068
#> 2 14.14214  0.00000 22.36068
#> 3 22.36068 22.36068  0.00000

^{Created on 2025-07-28 with reprex v2.1.1}

2. Geographic distance from sf object with lon/lat

library(sf)
#> Linking to GEOS 3.13.1, GDAL 3.10.2, PROJ 9.5.1; sf_use_s2() is TRUE

# Coordinates in lon-lat (WGS84)
coords <- data.frame(
  id = 1:3,
  lon = c(4.895, 4.901, 4.880),
  lat = c(52.370, 52.365, 52.390)
)

# Convert to sf
pts_sf <- st_as_sf(coords, coords = c("lon", "lat"), crs = 4326)

# Use st_distance (returns units)
dist_geo <- st_distance(pts_sf)
dist_geo
#> Units: [m]
#>           [,1]      [,2]     [,3]
#> [1,]    0.0000  689.2457 2445.884
#> [2,]  689.2457    0.0000 3124.051
#> [3,] 2445.8835 3124.0513    0.000

^{Created on 2025-07-28 with reprex v2.1.1}

3. k-nearest neighbor weight matrix with ordinal weights

library(FNN)

# Create a grid of 9 points (3x3)
coords <- expand.grid(x = seq(0, 20, by = 10),
                      y = seq(0, 20, by = 10))

# Get 3 nearest neighbors
k <- 3
knn <- get.knn(coords, k = k)

# Create weight matrix (asymmetric)
n <- nrow(coords)
nn_weights <- matrix(0, nrow = n, ncol = n)

for (i in 1:n) {
  nn_weights[i, knn$nn.index[i, 1]] <- 3  # nearest
  nn_weights[i, knn$nn.index[i, 2]] <- 2  # second nearest
  nn_weights[i, knn$nn.index[i, 3]] <- 1  # third nearest
}

nn_weights
#>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#>  [1,]    0    2    0    3    1    0    0    0    0
#>  [2,]    1    0    3    0    2    0    0    0    0
#>  [3,]    0    3    0    0    1    2    0    0    0
#>  [4,]    2    0    0    0    1    0    3    0    0
#>  [5,]    0    2    0    0    0    3    0    1    0
#>  [6,]    0    0    1    0    2    0    0    0    3
#>  [7,]    0    0    0    2    1    0    0    3    0
#>  [8,]    0    0    0    0    3    0    1    0    2
#>  [9,]    0    0    0    0    1    3    0    2    0

^{Created on 2025-07-28 with reprex v2.1.1}