[R-sig-Geo] Simulating spatially autocorrelated data

[Roger Bivand]

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