Spatial outliers detection in R?

Question

What is the best way to detect spatial outliers in R? With spatial outlier detection I mean:

Discover the objects whose non-spatial attribute values are significantly different from the values of their spatial neighbors.

Maybe performing univariate Z algorithm, Moran scatterplot or Scatterplot approach like this paper describes: On Detecting Spatial Outliers?

Or performing Metric Learning like this paper describes: Contextual Spatial Outlier Detection with Metric Learning?

Other approaches are welcome! I'm looking for a working implementation using R code.

Answer 1

We can run through some straight-forward approaches starting with looking at the modified z-score on the variable of interest (a-spatial) then move to calculating the local z-score and variance within a specified distance and finally calculate the local Moran's-I statistic. Our final example is evaluating the autocorrelation and heteroskedasticity in the residuals in an OLS regression.

First, let's add libraries and data.

library(sp)
library(spdep)
library(RANN)
library(spatialEco)

data(meuse)
  coordinates(meuse) <- ~x+y

Here we can calculate global outliers using modified z-score. None of the z values in this variable exceed the common threshold (ie., z>9) but, by plotting the z values you can see the spatial distribution of z.

( meuse$Zscore <- spatialEco::outliers(meuse$cadmium) )  
  spplot(meuse, "Zscore", col.regions=cm.colors(10))

Here we calculate local (1000m radius) variation using a distance-based neighbor variance and the local z-score. This is simply a statistic of each point including all neighbors within the specified distance. You could apply a large variety of statistics using this approach. In this example, we do identify one significant local outlier using the local modified z-score.

dnn <- RANN::nn2(coordinates(meuse), searchtype="radius", 
                 radius = 1000)$nn.idx
  var.cadmium <- rep(NA,nrow(meuse))
  z.cadmium <- rep(NA,nrow(meuse))  
    for(i in 1:nrow(dnn)){
      dnn.idx <- dnn[i,] 
      var.cadmium[i] <- var( meuse[dnn.idx[dnn.idx != 0],]$cadmium, na.rm=TRUE)
      z.cadmium[i] <- outliers(meuse[dnn.idx[dnn.idx != 0],]$cadmium)[1]
    }
z.cadmium[!is.finite(z.cadmium)] <- 0 

meuse$var.cadmium <- var.cadmium
spplot(meuse, "var.cadmium", col.regions=cm.colors(10))

meuse$z.cadmium <- z.cadmium
spplot(meuse, "z.cadmium", col.regions=cm.colors(10))

Here is where we calculate the local autocorrelation using the Local Moran's-I or LISA statistic. First, we need to calculate a minimum search distance to ensure that we do not have any empty sets (null Wij neighbor matrix) and then use the min distance to build the Wij spatial weights matrix.

all.linked <- max(unlist(nbdists(knn2nb(knearneigh(coordinates(meuse))), 
                  coordinates(meuse))))
nb <- dnearneigh(meuse, 0, all.linked)

Now we can use the Wij matrix to calculate the Local Moran's-I and create a point dataset, corresponding to meuse, to explore results. We also print and plot the number of observations that are statistically significant.

 mI <- localmoran(meuse@data[,"cadmium"], nb2listw(nb, style="W"))
 LocalI <- meuse
   LocalI@data <- data.frame(ID=rownames(LocalI@data), as.data.frame(mI))
     names(LocalI@data)[6] <- "Pr"
spplot(LocalI, "Z.Ii", xlab="Local Morans-I", col.regions=topo.colors(30))   

cat(nrow( LocalI[LocalI@data[,"Pr"] < 0.05 ,]), "obs of", 
    nrow(LocalI), "are significant at p=0.05","\n")

plot(LocalI, pch=19)
  points(LocalI[which(LocalI$Pr <= 0.05),], pch=19,col="red")

The high and low have no distinct I values so, one must create a vector to distinguish significant and high hotspots. Red represents hot spots or spatial outliers.

LocalI@data <- data.frame(LocalI@data, HotSpots=ifelse( mI[,5] <= 0.05 & mI[,4] >= mean(mI[,4]), 1, 0) )
  LocalI@data$HotSpots <- as.factor(LocalI@data$HotSpots)

spplot(LocalI, "HotSpots", xlab="Local Moran’s-I Hot Spots", col.regions=c("blue","red") )

There are, of course, many ways to address spatial outliers. You could use the Durbin-Watson or Lagrange multiplier test statistics to evaluate the autocorrelation in the residual error of a linear model and perhaps specify a spatial lag model.

Here is an example were we fit an OLS model to the Boston housing data.

library(spdep)
data(boston)
  b <- boston.c
    b$LOGMEDV <- log(b$CMEDV)
  coordinates(b) <- boston.utm
blm <- lm(LOGMEDV ~ RM + LSTAT + CRIM + ZN + CHAS + DIS, data=b@data)
  b$lmresid <- residuals(blm)
    spplot(b, "lmresid")

Now we can test the residual error for spatial dependence using the Morna's-I and the Lagrange multiplier diagnostics for spatial dependence. Results indicate that a spatial lag model is in order. We can also use a Breusch-Pagan test for evaluating heteroskedasticity. These tests will indicate the leverage (influence) that an observation(s) has on the autocorrelation structure in a regression model.

dnn <- spdep::dnearneigh(coordinates(b), d1=0, d2=3.973, row.names=row.names(b))  
  Wij <- spdep::nb2listw(dnn)

spdep::lm.morantest(blm, Wij)  
spdep::lm.LMtests(blm, Wij, test="all") 
lmtest::bptest(blm)

A multivariate case of local spatial outlier effect can be approximated using a partial correlation function. This is available in the ncf package.

library(ncf)

x <- expand.grid(1:20, 1:5)[,1]
y <- expand.grid(1:20, 1:5)[,2]
z <- cbind(rmvn.spa(x=x, y=y, p=2, method="exp"),
           rmvn.spa(x=x, y=y, p=2, method="exp"))

( part.corr <- ncf::lisa.nc(x=x, y=y, z=z, neigh=3) )
  plot(part.corr)