[R-sig-Geo] problem simulating random point patterns

[Roger Bivand]

On Tue, 9 Oct 2007, Dale Steele wrote:

> I'm trying to reproduce a randomization test for complete spatial
> randomness (described in Manly, 2007 Randomization, Bootstrap and
> Monte Carlo Methods in Biology). As a first step, I am generating
> random point processes in a rectangle.  I can reproduce results in the
> book when I generate data as follows:
>
> # (1)
> nsim <- 999
> xmin <- 0; xmax <- 40
> ymin <- 0; ymax <- 50
> npoints <- 22
> xpoly <- c(0,40,40,0)
> ypoly <- c(0,0,50,50)
> poly <- as.points(xpoly, ypoly)
>
> # generate an array of nsim random spatial point patterns, each with npoints
> rpoints <- array(NA, dim=c(npoints,2,nsim))
>  for (i in 1:nsim) {
>  rpoints[ , ,i] <- cbind(runif(npoints,min=xmin,
> max=xmax),runif(npoints,min=xmin,max=xmax))
> }
>
> However, when I use either the csr() function from splancs or
> runifpoint from spatstat, the randomization distributions of my
> statistic (mean kth nearest neighbor distance) are similar, but do not
> reproduce the results of simulations using (1) to simulate the data.
>
> I would expect that (1), (2) and (3) would be identical for a
> rectangular sampling area.  Would like to extend to irregular
> polygons.  Am I missing something simple here?  Thanks!  --Dale

One of the features of open source software is that you can read the code. 
Just entering the function name with no brackets at the prompt shows the 
R code of the function.

Tracing back from csr() in splancs, you get to ranpts() then gen(), where 
you can see that runif() is being used slightly differently, to 
accommodate irregular polygons; if you has used set.seed(), you would 
probably have seen detail differences. pip() at the end of gen() can also 
give different answers to the in/out question for points falling exactly 
on the boundary than you might expect.

Similarly, looking at runifpoint() in spatstat takes you to runifrect(), 
which again makes heuristic adjustments to the number of points being 
sampled to suit a range of windows, with some sampled points being 
discarded if they are superfluous.

Finally, the spsample methods in the sp package also provide a solution, 
based on Brian Ripley's 1981 Spatial Statistics, for which you can read 
the R code. The sample.Spatial() method also uses runif() internally, in 
your case:

W <- new("Spatial", bbox=matrix(c(0,0,40,50), nrow=2, dimnames=list(NULL,
   c("min", "max"))))
coordinates(spsample(W, npoints, type="random"))

retrieves the coordinates. So you can trace through what all of these 
variants are doing in R code, and setting set.seed() to a fixed number, 
you can provide all the methods with the same stream of random numbers. 
Using debug() to step through an example run, you'll be able to find out 
exactly why what you expect isn't happening.

Please recall that splancs and spatstat are based on decades of 
experience. For corroboration, you can also read the underlying C code 
under Psim() in spatial, written by Brian Ripley (and supporting his book 
with Bill Venables: Modern Applied Statistics with S) - hint: it is using 
the C API to runif() as you see in spatial/src/pps.c.

Hope this helps,

Roger

>
> # (1) generate an array of nsim random spatial point patterns within a
> specific polygon area
> library(splancs)
> rpoints <- array(NA, dim=c(npoints,2,nsim))
>  for (i in 1:nsim) {
>    rpoints[, ,i] <- csr(poly, npoints)
>  }
>
> ###  (2) same as above using spatstat runifpoint
> library(spatstat)
> rpoints <- array(NA, dim=c(npoints,2,nsim))
>  for (i in 1:nsim) {
>    rpoints[, ,i] <-  runifpoint(npoints,
> win=owin(c(xmin,xmax),c(ymin,ymax)), giveup=1000)
>  }
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, 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