Lattice with autocorrelation and predetermined values

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20090419/6538abbb/attachment.pl>
Dear Simon,

Perhaps you could model your data using a multivariate Normal with
spatial autocorrelation (using a SAR or CAR specification). In this way,
you can fix some of the values and then simulate the rest. You can tune
the autocorrelation coefficient to make sure that the spatial
correlation is large enough to give high values of Moran's I.

Hope this helps,

Virgilio

El dom, 19-04-2009 a las 21:16 +0000, Simon Chamaill? escribi?:
Dear all,

I want to distribute a vector of known values onto a lattice, while still making sure that the Moran autocorrelation is higher than a specific threshold.

I know of several ways of getting a lattice with some pre-determined levels of approximate spatial autocorrelation, but all involve some random number generation of one sort or another, whereas I want cell values to have some pre-determined values.

I assume that I could:
1) generate a lattice
2) randomly distribute my pre-determined values
3) compute Moran's I, then randomly permutate two cell values
5) recompute Moran's I and, and if higher that the one calculated in 3) keep the new lattice
6) back to 3) for permutation, or stop when Moran's I higher than my threshold, or kill the simulation when i've reached a specific number of iterations.

The above is obviously so inefficient that it is likely to be of no use for a lattice of even a moderate size or for a high autocorrelation threshold to be reached. I intend to work on approx. 100*100 lattice, with values equally distributed between 0-100, and autocorrelation threshold of approx. 0.3.

Any suggestions (hum..., even 'don't even think of doing something like this'), would be appreciated.
Thanks,
simon

	[[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Dear Simon,

Perhaps you could model your data using a multivariate Normal with
spatial autocorrelation (using a SAR or CAR specification). In this way,
you can fix some of the values and then simulate the rest. You can tune
the autocorrelation coefficient to make sure that the spatial
correlation is large enough to give high values of Moran's I.
Yes, this is the way to go, but raises practical problems for a 100 by 100 
grid (10') observations. In this case, you would need to use a "W" or "C" 
style weights list, convert to sparse matrix format, and power the matrix 
some k times. The inverse of (I - rho W) will then be the sum of:

rho^i W^i, for i=0,...,k

Using the "B" style and/or some kinds of general weights may lead to a 
slow convergence of rho^k W^k to machine zero.

Note also that rho is not Moran's I, but is scaled depending on the 
structure of W and its extreme eigenvalues.

The original idea of swapping values cell-wise will not work, because the 
values propagate through (I - rho W)^{-1}, so you'd need to choose from n 
samples (permutations) of the input vector. Just sample the vector say 
10000 times, compute Moran's I (possibly simplifying the code for speed), 
and choose one of the outputs satisfying your criteria. You need to think 
about storage and about reproducibility (think set.seed() and friends, to 
record the state of the RNG before each permutation, rather than the 
permutation vector.

For a 10 by 10 raster, this might look like:

library(spdep)
n <- 100
set.seed(1)
input <- rnorm(n)
nb <- cell2nb(10, 10, torus=TRUE)
lw <- nb2listw(nb, style="C")
S0 <- Szero(lw)
seeds <- sample.int(1e7, 1e4)
outI <- sapply(seeds, function(x) {
   set.seed(x)
   perm <- sample.int(n, n)
   moran(input[perm], lw, n, S0)$I
})
hist(outI)

What you notice is that you'd actually need a very large number of samples 
to achieve even moderate autocorrelation, so maybe the SAR simulation is 
the way to go, or maybe an SMA scheme:

n <- 100
set.seed(1)
input <- rnorm(n)
nb <- cell2nb(10, 10, torus=TRUE)
lw <- nb2listw(nb, style="C")
S0 <- Szero(lw)
input <- input + 0.9 * lag(lw, input)
moran(input, lw, n, S0)$I

for avoiding the problem of inverting a large matrix, or approximating 
through the sum of a power series.