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.

  • This really depends on what you are after in identifying "spatial outliers". Are you worried about nonstationarity in a geostatistical model, them effecting you inference, wanting to draw inference on 2nd order process, etc... It is not enough to want to do something, there should be a rational behind it and, in this case a statistical justification. If you focused your question, providing context to your analysis goals/model, what you have already tried, perhaps with a reproducible example, you will get a much more relevant answer. Commented Feb 6, 2018 at 15:38
  • 1
    @JeffreyEvans thank for your comments. I just simply want to detect some points that are very different from their neighbors to be considered an spatial outlier. I don't understand why an answer to my question depends on what I want to do next. In my case, an echosounder may introduce bad measures of depths and I want to detect them but the goals and objectives can be diverse.
    – Guz
    Commented Feb 6, 2018 at 15:50
  • 2
    Because, the type of statistic really matters. For example a lisa would indicate the juxtaposition of high and low values within the neighborhood whereas a local Geary's-C would be comparable to the locally decomposed semivariance and would represent a different type of outlier effect that would be more relevant to violation of nonstatiionarity in geostatistical models. The lisa would be relevant in drawing inference regarding 2nd order spatial variation or could be used as weights in a conditional autoregressive model. Commented Feb 6, 2018 at 17:27
  • 1
    Im not familiar with the papers you cited but I would start with a local Moran analysis. Then it’s a subjective decision where to cut off autocorrelation to be considered an outlier. Maybe there are some literature values for that.
    – Dominik
    Commented Feb 6, 2018 at 20:22

1 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.


  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))), 
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.

  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") 

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


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) )
  • 1
    Thanks a lot for your very detailed and complete answer! Also, the first part is a very intuitive approach.
    – Guz
    Commented Feb 13, 2018 at 17:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.