Using pMEM for Spatial Modelling with Predictive Moran’s Eigenvector Maps

Guillaume Guénard

Pierre Legendre

2024-09-26

Introduction

In the present paper, we show how to model the values of variables in space using a representation of their spatial variation. These representations are provided by predictive Moran’s eigenveector maps (pMEM), an extension of Moran’s eigenvector maps (Dray, Legendre, and Peres-Neto 2006), which take their name from Moran’s autocorrelation index (Moran 1948, 1950). For pMEM, we use a simplified version of MEM whereby the connectivity is decided only by the distance between the sampling locations. Also, package pMEM implements functionalities, such as supplementary distance weighting functions, that were not defined by Dray, Legendre, and Peres-Neto (2006) in their original definition of the method.

Data: Atlantic salmon parr distribution

The exemplary data set features observations (performed by snorkelling) of the presence and numbers of Atlantic salmon parr (juvenile) in \(76\) \(20\,\mathrm{m}\) river segments forming a \(1\,520\,\mathrm{m}\) transect of the St. Marguerite River, Québec, Canada. Sampled together with these observations were the channel depth (\(D\)) the water velocity (\(V\)), and mean substrate grain size (\(D_{50}\)), which were taken at the thalweg in the middle of each section (i.e., \(10\,\mathrm{m}\) after the beginning and \(10\,\mathrm{m}\) before the end of each section). Observations were performed while moving in the upstream direction in a zigzag manner in order not to disturb the fish prior to their observation. This data package is loaded in the R environment as follows:

data("salmon", package = "pMEM")

Calculations

Les us first load the necessary packages for this example:

library("pMEM")         ## To calculate pMEM
## Le chargement a nécessité le package : sf
## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
library("magrittr")     ## For its very handy pipe operateur (%>%)
library("glmnet")       ## To calculate elastic net regression
## Le chargement a nécessité le package : Matrix
## Loaded glmnet 4.1-8

Packages magrittr and glmnet are suggested package of pMEM, and are therefore not automatically installed along with the latter. You may have to install them manually in order to execute this example code. Also, we need to draw indices for the training and the testing sets, which we are obtain as follows:

set.seed(1234567890)                     ## For the drawing to be repeatable.
sort(sample(nrow(salmon),25)) -> test    ## Drawing the testing set (N = 25).
(1:nrow(salmon))[-test] -> train         ## The training set is the remainder.

Here, we set a random seed to a specific value to make subsequent analyses repeatable. Feel free to skip this line, in which case your results will differ somewhat from the ones that follows.

Function genSEF

A pMEM is generated using function genSEF() as follows:

genSEF(
  x = salmon$Position[train],   ## The set of locations.
  m = genDistMetric(),          ## The distance metric function.
  f = genDWF("linear", 500)     ## The distance weighting function.
) -> sef0

sef0                            ## Show the resulting object.
## A SEMap-class object
## --------------------
## Number of sites: 51
## Directional: no
## Number of components: 50
## Eigenvalues: 12.36926,9.94984,5.29604,...,0.01169,0.00863
## --------------------

For this call, argument x is the set of locations coordinates; they are given to the distance metric function that is given to genSEF through its argument m. In this example, it is the returned value of function genDistMetric. In this simple case, genDistMetric is called without an argument and returns a two argument function that calculate the Euclidean distance between the elements (or rows) of the vectors or matrices given to these two arguments. Argument f is also given a one argument spatial weighting function; which take the distances and returning weights. The later is generated by a function called genDWF that itself has three arguments: fun which specify the type of weighting function, dmax which specify the threshold distance (or scale parameters when fun is one of “exponential”, “Gaussian”, or “hole_effet”), and shape for any additional shape parameters (currently used when fun is one of “power” or “hyperbolic”). The shape of the resulting spatial eigenfunctions can be shown as follows:

## A regular transect of points 1 m apart:
salmon$Position %>%
  {seq(min(.) - 20, max(.) + 20, 1)} -> xx

## Custom plotting function:
plotSEF <- function(sef, xTrain, xTest, xx, wh, ...) {
  plot(x = xx, y = predict(sef, xx, wh), type = "l", ...)
  points(x = xTrain, y = as.matrix(sef, wh), pch=21, bg="black")
  points(x = xTest, y = predict(sef, xTest)[,wh], pch=21, bg="red")
  invisible(NULL)
}

## Storing the graphical parameters:
p <- par(no.readonly = TRUE)

## Changing the graphical parameters:
par(mfrow=c(3,2), mar=c(4.6,4.6,3,1))

## Generate a six-inset plot:
for(fun in c("power","hyperbolic","spherical","exponential","Gaussian",
             "hole_effect"))
  genSEF(
    x = salmon$Position[train],
    m = genDistMetric(),
    f = genDWF(fun, 250, 0.75)
  ) %>%
    plotSEF(salmon$Position[train], salmon$Position[test], xx, 1,
            xlab="Location (m)", ylab="pMEM score", main=fun, lwd=2)

Fig. 1. Examples of pMEM spatial eigenfunctions of order 1 with descriptor (back markers) and prediction scores (red markers). The black continuous line is calculated for 1-m intervals to show the continuity of the spatial eigenfunctions.

## Restoring the graphical parameters:
par(p)

Function getMinMSE

Another function that we need to introduce before proceeding any further is called getMinMSE() and is a utility function for fitting simple linear models involving only SEFs and a single normally-distributed response variables. It needs a training and a testing set and search for set of SEF that, when fitted in the training set allows one to best predict the values of the testing set. Also, it is implemented in C++ through the Rcpp package and thus runs very fast. Function getMinMSE() is called as follows:

getMinMSE(
  U = as.matrix(sef0),
  y = salmon$Depth[train],
  Up = predict(sef0, salmon$Position[test]),
  yy = salmon$Depth[test],
  complete = FALSE
)
## $betasq
## [1] 0.0009423604
## 
## $mse
## [1] 0.008758679

From the documentation, arguments of that function are:

U
A matrix of spatial eigenvectors to be used as training data.
y
A numeric vector containing a single response variable to be used as training labels.
Up
A numeric matrix of spatial eigenvector scores to be used as testing data.
yy
A numeric vector containing a single response variable to be used as testing labels.
complete
A boolean specifying whether to return the complete data of the selection procedure (complete=TRUE; the default) or only the resulting mean square error and beta threshold (complete=FALSE).

Called with argument complete=FALSE, the function returns only the (out of the sample) mean square error (MSE) of the model and the minimum standardized regression coefficient of the model used to obtain it.

Building models for channel depth, current velocity, and substrate grain size

In the following section, we build three models that are purely spatial, involving only pMEM spatial eigenfunctions. To simplify the scripts, we will begin by creating lists for storing the results as follows:

sefTrain <- list()  ## For storing the "SEMap" objects.
mseRes <- list()    ## For storing results from function getMinMSE().
sel <- list()       ## For storing selected pMEM eigenfunctions.
lm <- list()        ## For storing the linear models.
prd <- list()       ## For storing the predictions.

To make predictions in a spatially-explicit manner, we first need to figure out which distance weighting function (DWF) to use and what parameter values to use with it. That decision can be carried out in many ways, one of which consists in performing a global search using the previously-defined objective function (objf). For this example, the global search is itself carried out using the simulated annealing implemented in R package stat. Since we have to repeat the same analyses three times for the different parr habitat descriptors, we defined a function performing the required calculations as follows:

estimateSEF <- function(x, xx, y, yy, lower, upper) {
  
  res <- list(optim = list())  ## A list to contain the results.
  
  ## This loop tries the seven DWF one by one, estimating 'dmax' (and, when
  ## necessary, 'shape') using simulated annealing.
  for(fun in c("linear","power","hyperbolic","spherical","exponential",
               "Gaussian","hole_effect")) {
    optim(
      par = c(0,if(fun %in% c("power","hyperbolic")) 0), fn = objf,
      method = "SANN", m = genDistMetric(), fun = fun,
      x = x, xx = xx, y = y, yy = yy,
      lb = c(lower[1],if(fun %in% c("power","hyperbolic")) lower[2]),
      ub = c(upper[1],if(fun %in% c("power","hyperbolic")) upper[2])
      
    ) -> res$optim[[fun]]
  }
  
  ## Extract the minimum values from the list of optimization:
  unlist(
    lapply(
      res$optim,
      function(x) x$value
    )
  ) -> res$bestval
  
  ## Find which DWF had the minimum objective criterion value:
  names(
    which.min(
      res$bestval
    )
  ) -> res$fun
  
  ## Back-transform the parameter values:
  res %>%
    {.$optim[[.$fun]]$par} %>%
    {(upper - lower) * (1 + exp(-.))^(-1) + lower} -> res$par
  
  ## Calculate the SEF using the optimized DWF parameters:
  res %>%
    {genSEF(
      x = x,
      m = genDistMetric(),
      f = genDWF(.$fun, .$par[1L], if(length(.$par) > 1) .$par[1L])
    )} -> res$sef
  
  ## Return the result list:
  res
}

Channel depth

The channel depth is an important descriptor of juvenile Atlantic salmon (parr) habitat. For instance, these fish use riffles to feed on drifting preys; it is while feeding that they are the most readily observable by snorkelers. They may also use pools as a refuge from predators or, perhaps, head for riffles with a close by pool in order to benefit from favourable hydrological conditions as well as readily available refuge. In salmon rivers, pools and riffles alternate successively, in a more or less regular manners; it might thus by possible to model part of the variation in channel depth using pMEM.

Let us estimate the optimal SEF for modelling channel depth as follows:

estimateSEF(
  x = salmon$Position[train],
  xx = salmon$Position[test],
  y = salmon[["Depth"]][train],
  yy = salmon[["Depth"]][test],
  lower = c(20,0.25),
  upper = c(1000,1.75)
) -> sefTrain[["Depth"]]

The best DWF found has been linear, with a \(d_{max}\) of \(65\,\mathrm{m}\). The minMSE model estimated for the channel depth using these parameters is obtained as follows:

## Calculate the channel depth model:
sefTrain[["Depth"]]$sef %>%
  {getMinMSE(
    U = as.matrix(.),
    y = salmon[["Depth"]][train],
    Up = predict(., salmon$Position[test]),
    yy = salmon[["Depth"]][test]
  )} -> mseRes[["Depth"]]

## Extract the selected SEF:
mseRes[["Depth"]] %>% {sort(.$ord[1:.$wh])} -> sel[["Depth"]]

which has a coefficient of prediction of 0.7003. Now, we can calculate a linear (regression) model from the selected spatial eigenfunctions:

## Calculate a linear model from the selected SEF:
lm(
  formula = y~.,
  data = cbind(
    y = salmon[["Depth"]][train],
    as.data.frame(sefTrain[["Depth"]]$sef, wh=sel[["Depth"]])
  )
) -> lm[["Depth"]]

## Calculate the predictions:
predict(
  lm[["Depth"]],
  newdata = as.data.frame(
    predict(
      object = sefTrain[["Depth"]]$sef,
      newdata = xx,
      wh = sel[["Depth"]]
    )
  )
) -> prd[["Depth"]]

The predictions of the channel depth can be displayed as follows:

plot(x=xx, y=prd[["Depth"]], type="l",
     ylim=range(salmon[["Depth"]], prd[["Depth"]]), las=1,
     ylab="Channel depth (m)", xlab="Location along the transect (m)")
points(x = salmon$Position[train], y = salmon[["Depth"]][train], pch=21,
       bg="black")
points(x = salmon$Position[test], y = salmon[["Depth"]][test], pch=21, bg="red")

Fig. 2. Spatially-explicit predictions of channel depth (black solid line) with observations used training (black markers) and testing (red markers) the model.

Current velocity

The second parr habitat descriptor is current velocity. Parrs fed on drifting prey animals (mostly insects) and thus rely on water current to bring about these preys. Whereas a fast current will bring preys at a high rate, is also entails less time reach for them and having to work harder in order to do so. As for channel depth, sections of slow and fast current alternate in a more or less consistent manner which might, to some extent, be modelled in a spatially-explicit manner. Current velocity is modelled as for the channel depth as follows:

## Estimate the most adequate predictive Moran's eigenvector map:
estimateSEF(
  x = salmon$Position[train],
  xx = salmon$Position[test],
  y = salmon[["Velocity"]][train],
  yy = salmon[["Velocity"]][test],
  lower = c(20,0.25),
  upper = c(1000,1.75)
) -> sefTrain[["Velocity"]]
## Calculate the current velocity model:
sefTrain[["Velocity"]]$sef %>%
  {getMinMSE(
    U = as.matrix(.),
    y = salmon[["Velocity"]][train],
    Up = predict(., salmon$Position[test]),
    yy = salmon[["Velocity"]][test]
  )} -> mseRes[["Velocity"]]

## Extract the selected SEF:
mseRes[["Velocity"]] %>% {sort(.$ord[1:.$wh])} -> sel[["Velocity"]]

## Calculate a linear model from the selected SEF:
lm(
  formula = y~.,
  data = cbind(
    y = salmon[["Velocity"]][train],
    as.data.frame(sefTrain[["Velocity"]]$sef, wh=sel[["Velocity"]])
  )
) -> lm[["Velocity"]]

## Calculate the predictions:
predict(
  lm[["Velocity"]],
  newdata = as.data.frame(
    predict(
      object = sefTrain[["Velocity"]]$sef,
      newdata = xx,
      wh = sel[["Velocity"]]
    )
  )
) -> prd[["Velocity"]]

This time, the best DWF found has been linear, with a \(d_{max}\) of \(742\,\mathrm{m}\), a coefficient of predictions of 0.4914, and the predictions appear as follows:

plot(x=xx, y=prd[["Velocity"]], type="l",
     ylim=range(salmon[["Velocity"]], prd[["Velocity"]]), las=1,
     ylab="Velocity (m/s)", xlab="Location along the transect (m)")
points(x = salmon$Position[train], y = salmon[["Velocity"]][train], pch=21,
       bg="black")
points(x = salmon$Position[test], y = salmon[["Velocity"]][test], pch=21,
       bg="red")

Fig. 3. Spatially-explicit predictions of current velocity (black solid line) with observations used training (black markers) and testing (red markers) the model.

Substrate grain size

Atlantic salmon parr is bottom dwelling fish that use the hydrodynamic conditions in the vicinity of the substrate in various manner. While feeding, for instance, it typically choose a cobble against which it lays its pectoral fins to help in holding a steady position in the stream. Also, it uses zones of back-flow close to the rough bottom surface in order to swim back to its ambush position. Substrate composition may thus be a dependable descriptor of the parr habitat. Substrate mean grain size is modelled as for the two previous descriptors as follows:

## Estimate the most adequate predictive Moran's eigenvector map:
estimateSEF(
  x = salmon$Position[train],
  xx = salmon$Position[test],
  y = salmon[["Substrate"]][train],
  yy = salmon[["Substrate"]][test],
  lower = c(20,0.25),
  upper = c(1000,1.75)
) -> sefTrain[["Substrate"]]
## Calculate the mean substrate grain size model:
sefTrain[["Substrate"]]$sef %>%
  {getMinMSE(
    U = as.matrix(.),
    y = salmon[["Substrate"]][train],
    Up = predict(., salmon$Position[test]),
    yy = salmon[["Substrate"]][test]
  )} -> mseRes[["Substrate"]]

## Extract the selected SEF:
mseRes[["Substrate"]] %>% {sort(.$ord[1:.$wh])} -> sel[["Substrate"]]

## Calculate a linear model from the selected SEF:
lm(
  formula = y~.,
  data = cbind(
    y = salmon[["Substrate"]][train],
    as.data.frame(sefTrain[["Substrate"]]$sef, wh=sel[["Substrate"]])
  )
) -> lm[["Substrate"]]

## Calculate the predictions:
predict(
  lm[["Substrate"]],
  newdata = as.data.frame(
    predict(
      object = sefTrain[["Substrate"]]$sef,
      newdata = xx,
      wh = sel[["Substrate"]]
    )
  )
) -> prd[["Substrate"]]

For mean substrate grain size, the best DWF found has been power, with a \(d_{max}\) of \(881\,\mathrm{m}\), a \(\alpha\) of \(1.67\) coefficient of predictions of 0.5726, and the predictions appear as follows:

Fig. 4. Spatially-explicit predictions of mean substrate grain size (black solid line) with observations used training (black markers) and testing (red markers) the model.

Atlantic salmon parr abundance

The Atlantic parr abundances are count data, suggesting it is the result of a Poisson process. A Poisson GLM is a straightforward way to model a Poisson process. Here, we will take this situation as an opportunity to exemplify another way whereby pMEM could be used for modelling, this time using an “elasticnet”-regularized generalized linear model. For that purpose, we need another objective function using function glmnet for model estimation, and having arguments w and ww to to pass auxiliary descriptors (training and testing sets, respectively), to be used alongside the SEF in modelling parr, abundance as follows:

objf2 <- function(par, m, fun, x, xx, y, yy, w, ww, lb, ub) {
  par <- (ub - lb) * (1 + exp(-par))^(-1) + lb
  if(fun %in% c("power","hyperbolic")) {
    sef <- genSEF(x, m, genDWF(fun, range = par[3L], shape = par[4L]))
  } else
    sef <- genSEF(x, m, genDWF(fun, range = par[3L]))
  glm1 <- glmnet(x = cbind(w, as.matrix(sef)), y = y, family = "poisson",
                 alpha = par[1L], lambda = par[2L])
  pp <- predict(glm1, newx = cbind(ww, predict(sef, xx)), type="response")
  -2*sum(dpois(yy, pp, log = TRUE))
}

That objective function has three or four parameters (depending on the presence of the shape parameter), rather than 1 or 2 for objf:

  1. the \(\alpha\) parameter of the elastic net regression that define the amount of \(L_1\) with respect to \(L_2\) norm used for regularization,

  2. the \(\lambda\) parameter, which is the overall amount of regularization,

  3. the \(d_{max}\) parameter of the DWF, and

  4. the shape parameter of the DWF, when necessary.

The value returned is the out of the sample deviance value rather than the out of the sample mean square error that was returned by the first objective function.

The auxiliary descriptors are the channel depth, current velocity, and mean substrate grain size. There may be non-linear relationships between these descriptors and parr abndance, because parr may prefer sites with intermediate values of these descriptors over extremes. To allow the parr abundance model the opportunity to exploit that possibility, we calculated orthogonal polynomials from the training data using R function poly as follows:

## Implement a list of orthogonal polynomial objects:
plist <- list()
plist[["Depth"]] <- poly(salmon[train,"Depth"],2)
plist[["Velocity"]] <- poly(salmon[train,"Velocity"],2)
plist[["Substrate"]] <- poly(salmon[train,"Substrate"],2)

## The matrix of auxiliary descriptor for the training set:
cbind(
  as.matrix(plist[["Depth"]]),
  as.matrix(plist[["Velocity"]]),
  as.matrix(plist[["Substrate"]])
) -> w

## Generate suitable column names:
c("Depth^1","Depth^2",
  "Velocity^1","Velocity^2",
  "Substrate^1","Substrate^2") -> colnames(w)

## The matrix of auxiliary descriptor for the testing set:
cbind(
  predict(plist[["Depth"]], newdata=salmon[test,"Depth"]),
  predict(plist[["Velocity"]], newdata=salmon[test,"Velocity"]),
  predict(plist[["Substrate"]], newdata=salmon[test,"Substrate"])
) -> ww

## Copying the column names:
colnames(ww) <- colnames(w)

The objective function is executed as follows:

objf2(
  par = c(0, 0, 0, 0),
  m = genDistMetric(),
  fun = "Gaussian",
  x = salmon$Position[train],
  xx = salmon$Position[test],
  y = salmon[["Abundance"]][train],
  yy = salmon[["Abundance"]][test],
  w = w,
  ww = ww,
  lb = c(0,0,20,0.25),
  ub=c(1,1,1000,1.75)
) -> res2

res2
## [1] 135.7557

yielding a parr abundance model with an out of the sample deviance of \(135.8\).

estimateSEF2 <- function(x, xx, y, yy, w, ww, lower, upper) {
  
  res <- list(optim = list())  ## A list to contain the results.
  
  ## This loop tries the seven DWF one by one, estimating 'dmax' (and, when
  ## necessary, 'shape') using simulated annealing.
  for(fun in c("linear","power","hyperbolic","spherical","exponential",
               "Gaussian","hole_effect")) {
    optim(
      par = c(0,0,0,if(fun %in% c("power","hyperbolic")) 0), fn = objf2,
      method = "SANN", m = genDistMetric(), fun = fun,
      x = x, xx = xx, y = y, yy = yy, w = w, ww = ww,
      lb = c(lower[1:3],if(fun %in% c("power","hyperbolic")) lower[4]),
      ub = c(upper[1:3],if(fun %in% c("power","hyperbolic")) upper[4])
    ) -> res$optim[[fun]]
  }
  
  ## Extract the minimum values from the list of optimization:
  unlist(
    lapply(
      res$optim,
      function(x) x$value
    )
  ) -> res$bestval
  
  ## Find which DWF had the minimum objective criterion value:
  names(
    which.min(
      res$bestval
    )
  ) -> res$fun
  
  ## Back-transform the parameter values:
  res %>%
    {.$optim[[.$fun]]$par} %>%
    {(upper[1:length(.)] - lower[1:length(.)]) * (1 + exp(-.))^(-1) +
        lower[1:length(.)]} -> res$par
  
  ## Calculate the SEF using the optimized DWF parameters:
  res %>%
    {genSEF(
      x = x,
      m = genDistMetric(),
      f = genDWF(.$fun, .$par[3], if(length(.$par) > 3) .$par[4])
    )} -> res$sef
  
  ## Return the result list:
  res
}

We can now estimate the optimal SEF for modelling parr abundance as follows:

estimateSEF2(
  x = salmon$Position[train],
  xx = salmon$Position[test],
  y = salmon[["Abundance"]][train],
  yy = salmon[["Abundance"]][test],
  w = w,
  ww = ww,
  lower = c(0,0,20,0.25),
  upper = c(1,1,1000,1.75)
) -> sefTrain[["Abundance"]]

The best DWF found has been hole_effect, with an \(\alpha\) of \(0.993\), a \(\lambda\) of \(0.1054\), and a \(d_{max}\) of \(91\,\mathrm{m}\).

The elasticnet model is estimated as follows:

cbind(w, as.matrix(sefTrain[["Abundance"]]$sef)) %>%
  glmnet(
    y = salmon$Abundance[train], family = "poisson",
    alpha = sefTrain[["Abundance"]]$par[1L],
    lambda = sefTrain[["Abundance"]]$par[2L]) -> lm[["Abundance"]]

## Model coefficients:
coef(lm[["Abundance"]])
## 31 x 1 sparse Matrix of class "dgCMatrix"
##                      s0
## (Intercept)  0.48775321
## Depth^1      .         
## Depth^2      .         
## Velocity^1   .         
## Velocity^2   .         
## Substrate^1  0.77749010
## Substrate^2  .         
## pMEM_1       0.56231401
## pMEM_2      -0.42979149
## pMEM_3       2.60281237
## pMEM_4       0.59190038
## pMEM_5      -1.33793350
## pMEM_6       0.69311351
## pMEM_7       .         
## pMEM_8       1.12570665
## pMEM_9      -0.32720176
## pMEM_10      .         
## pMEM_11      .         
## pMEM_12      0.14375648
## pMEM_13     -0.24449250
## pMEM_14     -0.38446542
## pMEM_15     -0.70611753
## pMEM_16     -2.64039180
## pMEM_17      .         
## pMEM_18      0.58608499
## pMEM_19     -0.27411596
## pMEM_20     -0.13680399
## pMEM_21      0.75884282
## pMEM_22      2.06077847
## pMEM_23     -0.08334478
## pMEM_24     -0.38344153

the continuous predictions are obtained as follows:

lm[["Abundance"]] %>%
  predict(
    cbind(
      predict(plist[["Depth"]], prd[["Depth"]]),
      predict(plist[["Velocity"]], prd[["Velocity"]]),
      predict(plist[["Substrate"]], prd[["Substrate"]]),
      predict(sefTrain[["Abundance"]]$sef, xx)
    ),
    type="response"
  ) ->  prd[["Abundance"]]

and the results are displayed as follows:

plot(x=xx, y=prd[["Abundance"]], type="l",
     ylim=range(salmon[["Abundance"]], prd[["Abundance"]]), las=1,
     ylab="Parr abundance (fish)",
     xlab="Location along the transect (m)")
points(x = salmon$Position[train], y = salmon[["Abundance"]][train], pch=21,
       bg="black")
points(x = salmon$Position[test], y = salmon[["Abundance"]][test], pch=21,
       bg="red")

References

Dray, S., P. Legendre, and P. R. Peres-Neto. 2006. “Spatial Modelling: A Comprehensive Framework for Principal Coordinate Analysis of Neighbour Matrices (PCNM).” Ecological Modelling, 483–93. https://doi.org/10.1016/j.ecolmodel.2006.02.015.
Moran, P. A. P. 1948. “The Interpretation of Statistical Maps.” Journal of the Royal Statistical Society: Series B (Methodological) 10: 243–51. https://doi.org/10.1111/j.2517-6161.1948.tb00012.x.
———. 1950. “Notes on Continuous Stochastic Phenomena.” Biometrika 37: 17–23. https://doi.org/10.1093/biomet/37.1-2.17.