Error while comparing two raster layer in R

Question

I have got a geographic point distribution with N prehistoric settlements with a size of x and want to compare the distribution with an economic simulation, which produces M settlements (Usally N > M). I've done about 110 runs with different parameters and now I'm looking for a best-fit model. Therefore I calculate the Kernel Density Estimation for the real settlemts (x_dens), the simulation results (y_dens) and for the M-biggest real settlements (xM_dens) and want to compare the cell values of two kernel density estimations in R.

As a first step, i am using the function density.ppp from the Package "spatstat".

x_dens <- density(x_ppp, sd, eps=rw, edge=TRUE, at="pixels")
x_ras <- raster(x_dens, crs="+proj=utm +zone=35 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0")

y_dens <- density(y_ppp, sd, eps=rw, edge=TRUE, at="pixels")
y_ras <- raster(y_dens, crs="+proj=utm +zone=35 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0")

xM_dens <- density(xM_ppp, sd, eps=rw, edge=TRUE, at="pixels")
xM_ras <- raster(xM_dens, crs="+proj=utm +zone=35 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0")

The next step includes the calculation of the cell differences for x with xM, to see which would be minimum error for my M simulated settlements:

xxM_comp <- sqrt(((x_ras - xM_ras)^2 / x_ras))
xxM_comp [!is.finite(xxM_comp )] <- 0
xxM_comp_result <- sum(xxM_comp @data@values)

As a next step, i'm doing the same for my Y-simulated settlements:

xy_comp <- sqrt(((x_ras - y_ras)^2 / x_ras))
xy_comp [!is.finite(xy_comp )] <- 0
xy_comp_result <- sum(xy_comp @data@values)

To evaluate the best-fit model, i caluclate:

results <- xy_comp_results/xxM_comp_result

xy_comp shall result in a new raster layer, which shows the differences in both point pattern, but the resulting cell values are enormously high at some points and as a consequence xy_comp_result has a illogically high value. It is due to the fact, that some pixel values of x_ras much to small compared with y_ras. Do you have got any suggestions, how to solve this problem. I've tried to remove outliers, but I did not get proper results.

The methodology is described in D. Stelder, Where do Cities Form? A Geographical Agglomeration Model for Europa, Journal of Regional Science 45 (4), 2005, 657-679 (http://www.regroningen.nl/stelder/doc/JRS_nov2005_c.pdf) Starting with page 669.

Answer 1

I am curious as to why you are not approaching this problem using a point pattern analysis? It is apparent that you are after a multiscale comparision but, it is not clear as to what end or what type of supported inference would be made. The type of standardization that your are attempting is hinting that a PPA would be a more supported methodology.

Depending on your hypothesis, I would highly recommend taking a look at the Geits-Ord or Ripley's-K statistic(s). If you have covariates there are correlogram methods implementing the partial Mantel test, which if permuted, support exploratory analysis. Another option for a marked point process would be the family of scan statistics, which can be specified using assumed point process distributions (ie., Poisson, Gaussian, Binomial) and as a multiscale model. All of these methods provide a spatially lagged (multiscale) evaluation of the spatial process and avoid many of the issues associated with kernel density estimates, which is never really considered a supported inferential method. For a robust inference you could formalize a point process model, in a hierarchical Bayesian framework (MCMC), using the estimated intensity function.

Now, if I have misunderstood and you, in fact, want to integrate the volume of two kernel density estimates you could follow one of the methodologies presented in Hurlbert (1978) that provides mathematical definitions for integrating niche volume overlap. The equations are relevant to any volume integration and the paper is available on JSTOR.

Hurlbert, S.H., (1978) The Measurement of Niche Overlap and Some Relatives. Ecology 59(1):67-77