Simulating spatially autocorrelated data
On Thu, 1 Sep 2011, Downey, Patrick wrote:
Hello all,
I'm trying to simulate a spatially autocorrelated random variable, and I
cannot figure out what the problem is. All I want is a simple spatial lag
model where
Y = rho*W*Y + e
Where e is a vector of iid normal random variables, rho is the
autocorrelation, W is a row-normalized distance matrix (a spatial weights
matrix), and Y is the random variable.
I thought the following program should do it, but it's not working. At the
end of the program, I calculate Moran's I, and it is not even close to
rejecting the null hypothesis of no spatial autocorrelation, even when rho
is very high (for example, below, rho is 0.95). Can someone please identify
what the problem is and offer some guidance on how to fix it?
PS - I apologize in advance, but I am not familiar with R's spatial
packages. I've done very little spatial analysis in R, so if there's a
package that can already do this, please recommend.
BEGIN PROGRAM:
install.packages("fields");library(fields)
install.packages("ape");library(ape)
N <- 200
rho <- 0.95
x.coord <- runif(N,0,100)
y.coord <- runif(N,0,100)
points <- cbind(x.coord,y.coord)
e <- rnorm(N,0,1)
dist.nonnorm <- rdist(points,points) # Matrix of Euclidean distances
dist <- dist.nonnorm/rowSums(dist.nonnorm) # Row normalizing the distance
matrix
diag(dist) <- 0 # Ensuring that the main diagonal is exactly 0
I think that you are using the distances as weights, not inverse distances, which seems more sensible.
I <- diag(N) # Identity matrix (not Moran's I) inv <- solve(I-rho.lag*dist) # Inverting (I - rho*W) y <- as.vector(inv %*% e) # Generating data that is supposed to be spatially autocorrelated Moran.I(y,dist) # Does not reject null hypothesis of no spatial autocorrelation
As Terry Griffin says, you can use spdep for this: library(spdep) rho <- 0.95 N <- 200 x.coord <- runif(N,0,100) y.coord <- runif(N,0,100) points <- cbind(x.coord,y.coord) e <- rnorm(N,0,1) dnb <- dnearneigh(points, 0, 150) dsts <- nbdists(dnb, points) idw <- lapply(dsts, function(x) 1/x) lw <- nb2listw(dnb, glist=idw, style="W") inv <- invIrW(lw, rho) y <- inv %*% e moran.test(y, lw) to reproduce your analysis with IDW, here without: lw <- nb2listw(dnb, glist=dsts, style="W") inv <- invIrW(lw, rho) y <- inv %*% e moran.test(y, lw) # no autocorrelation and here with a less inclusive distance threshold: dnb <- dnearneigh(points, 0, 15) dsts <- nbdists(dnb, points) idw <- lapply(dsts, function(x) 1/x) lw <- nb2listw(dnb, glist=idw, style="W") inv <- invIrW(lw, rho) y <- inv %*% e moran.test(y, lw) the larger the distance threshold, the less well the spatial process is captured, alternatively use idw <- lapply(dsts, function(x) 1/(x^2)), for example, to attenuate the weights more sharply. Hope this clarifies, Roger
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand Department of Economics, NHH Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: Roger.Bivand at nhh.no