# Testing species distribution model residuals for spatial autocorrelation?

I build a species distribution model with data from location A and tested it on location B. The picture shows a rastermap of location B and what I guess should be the model residuals.

• Green - are the areas where the model was right (value 0)
• Blue - are the areas where the model predicted species occurrence but there was no observation (value -1)
• Black - are the areas where the species was observed but the model predicted no occurrence (value 1)

I think the clustering is obvious but I also want to measure the spatial auto correlation and give it a value. In Arcgis and R I found the morans I functions. However both need point data as input.

Should I transform the Raster to points (one point each cel) and use this as basis for the calculation?

What should I do with the areas where the model was right (value 0)?

should I delete them before I perform the test?

This question is related to Calculating residuals of species distribution model and test for spatial auto correlation at Cross Validated.

Edit:

I asked the question because I always read that you have to take care of spatial autocorrelation in the analysis of species distributional data when building a SDM:

"Species distributional or trait data [...] often display spatial autocorrelation, i.e. locations close to each other exhibit more similar values than those further apart. If this pattern remains present in the residuals of a statistical model based on such data, one of the key assumptions of standard statistical analyses, that residuals are independent and identically distributed (i.i.d), is violated. "(Dormann, 2007).

However I cannot figure out whether it is an issue in my case. I did an systematic, full sampling of species occurrence/absence and environmental variables in Area A. This data (not georeferenced) I used to build a deterministic model to predict species distribution based on habitat preference. I just built habitat preference curves for each predictor based on species occurrence and habitat availability at each sample point and combined the habitat preferences for each predictor to a combined preference by multiplication. @Jeffrey Evans: this could also be interpreted as occurrence probability.

Then I applied the model to the independent area B, where I have a grid/raster of the environmental variables. I also mapped the species occurrence in B to test the model against independent data. I tested the occurrence probability against the observed occurrences. My AUC values are good (around 0.8).

Now I wanted to see whether there is some kind of spatial pattern in the residuals of the model. Iam not even sure whether I should call it residuals because i didn't build a statistical model and do the test against independent data. Maybe better call it error, the error between observed and predicted occurrence. My understanding is that when there is an error which shows spatial autocorrelation then my model omits an important predictor or process that explains the spatial distribution

My main worry is whether I can I trust the AUC value and assess the predictive performance as good?

Because I could not figure out on which data I should do the morans I test I askesd the initial question.

• What model did you use? Autocorrelation is something to address during the model fitting process and not post hoc. Violations of iid and multivariate normality are really only a concern in linear, frequentist statistics. These methods do not perform well with complex high dimensional space which is why the sdm community has moved towards nonparametric methods such as random forests. – Jeffrey Evans Jan 19 '18 at 19:56
• I understand that post hoc I will not be able to account for spatial auto correlation. However I applied a systematic sampling and made sure that the model building and test data are from two different locations. I built two models. The first is based on the geometric mean of the species preference for each environmental variable. The second is a random forest model. Are you saying that with these two approaches I do not have to consider spatial auto correlation in the model building data? Does it has to do with the trees being a randomly selected subspace of the data? thanks! – Konstantin_o Jan 22 '18 at 16:25
• From a habitat utilization, niche volume standpoint I would seriously question your methodology and would most certainly call it into question in review and, I review many niche related papers for several journals. An additive or multiplicative set of habitat use curves does not account for the hypervolume nor nonlinearity inherent in this type of data. Additionally, since the spatial variation is not sampled, this model cannot be extrapolated. – Jeffrey Evans Jan 22 '18 at 16:26
• A geometric mean would be highly suspect because of expected skewness in the data and how it aggregates in a binomial process. A random forests model is nonparametric and does not assume parametric norms as such, is not notably effected by autocorrelation in the y or x side of the equation. Based on your explanation of the data, it is quite unclear how you would specify a valid RF model. You really need to talk with a statistician, versed in this type of modeling, because there are some serious holes in your background statistical understanding that the literature is just not filling in. – Jeffrey Evans Jan 22 '18 at 16:31
• I sure realize now that my approach is somehow antiquated. However the model does a decent job in predicting the species occurrence. The AUC is 0.8 tested against data from a different (adjacent) location. The RF I build on 960 sample points with the information of species occurrence (response) and 3 environmental variable (predictor). Why would the model not be valid? – Konstantin_o Jan 22 '18 at 16:41

Given an inevitable nonparametric model structure, a residual in the traditional sense is not possible. You could approximate the response, Pearson's or deviance (likelihood ratio chi-squared) residual error using standard methods found in the logistic regression literature, eg., `response = [y - y-hat]`, `Pearson's = [(y - y-hat) / sqrt(y-hat)]` where; y = observed, y-hat = estimated probability. However, you do not know the curvilinear nature of the model so, cannot account for it in evaluating the residual error.

As @Senshi pointed out, "don't do statistics for statistics sake". In this case start with a testable hypothesis, evaluate appropriate test statistics and finally, understand application and the underlying assumptions of the selected statistic.

An autocorrelation statistic is not supported in this case. You functionally have a nominal outcome and statistics such as Moran's-I assume a continuous distribution. You are quite limited in model evaluation when using the binomial response and should really be estimating a continuous outcome (ie., probability, log likelihood, density). The only applicable autocorrelation statistic for binomial process would be a Joins-Count. You could try some of the point process statistics (eg., K-hat, Besag's-L, F-hat) but a density dependence would be in question thus, complicating the analysis. Besides, for testing any relevant hypothesis you would have to perform a cross statistic.

If you estimate the probability, then you could evaluate the autocorrelation of the resulting estimate and get an idea of the global and localized spatial structure. Honestly, the local autocorrelation would be much more informative here than the global.

I would ask, what are you trying to portray in evaluating the error of your model? I would assume that understanding the uncertainty (estimate or spatial) would be a relevant outcome here. You could do this through a Monte Carlo simulation or simple randomization of your model. In practice, you would apply simulation conditions to your data, rerun your model and make an estimate, to the observations or ideally, as a spatial prediction. This is done many times (eg., 999 at a p of 0.05) and a test statistic/moment, such as variance, is used to evaluate the uncertainty in the model or spatial estimate. It would be necessary to estimate the probability or log likelihood, and not the binomial response, as the variance of a series of [0,1] is nonsensical. Areas of high uncertainty would exhibit high variances. If you were using a randomization approach then a test statistic such as RMSE would be used (comparing the model estimate to each model randomization) to derive an error distribution.

I have also evaluated uncertainty of spatial estimates, for single model realizations, using focal nonmetric Multidimensional Scaling (nMDS). This indicates the "dissimilarity" with a specified spatial distance and can be interpreted as localized uncertainty. The advantage to nMDS is that you can look at dissimilarities in model estimates or parameter space. However, if only evaluating the spatial predictions, much the same can be done using a local autocorrelation (eg., LISA) to indicate the actual spatial structure ie., juxtaposition of high-low, high-high, low-high and low-low values. Please keep in mind that high variance produce in a Monte Carlo could be perceived as estimate error or model instability but, local evaluation of uncertainty via methods like local autocorrelation can be "real" in terms of the identification of spatial outliers, representing constrained niche conditions. Uncertainty is exactly that, and not always representing error per se.

• Much more detailled and precise answer than mine. Thanks for explaining the applicability of different methods. Discussing statistical methods always is enlightening. – Senshi Jan 19 '18 at 7:41
• thank a lot @ Jeffrey Evans @Senshi and everybody else. I added some details to my question. now it changed from a technical question to a more methodological one. I hope that is ok. – Konstantin_o Jan 19 '18 at 19:30

Yes. You already have point data, it's just presented in the form of a raster matrix, and not the vector-based points that your available Morans I algorithm expects. Morans I works fine on raster data as well, in principle, so it really is just an implementation issue.

You can safely convert forth and back between raster cells and points.

Regarding on what attribute to test for: That depends on what your goal is. If you use Morans I on all the cells (-1 to 1), you'll just learn how "clustered" your data is geographically. If you test for each value individually (e.g. just use cells with 0 / green), you can test the randomness of your matches.

Remember one of the most important questions in scientific work: "Is this relevant"? Reading your question, it seems to me you are a bit unsure of what you're doing and what conclusions you can even draw from Morans I.

Don't do statistics for statistic's sake, use it to reinforce or falsify(!) your assumptions and interpretations.

• Thanks for the comments. Your are totally right, Iam a bit unsure about what Iam doing. My first doubt is whether I can calculate the residuals of my species distribution model like I did and whether I have to care for the areas where the model predicted right (green, residual = 0) when analyzing the spatial auto correlation. – Konstantin_o Jan 18 '18 at 15:28