[R-sig-Geo] Simulation of inhomogeneous point process

[Rolf Turner]

On 26/03/15 16:15, ASANTOS wrote:
> Dear Members,
>
> - I have ant nests position in UTM in:
>
> #Ant nests position-----------------------
> a.x<-c(371349.2,371351.4,371277.1,371200.9,371139.9,371435.7,371344,371263.5,
>
> 371105.2,371023.2,371262.4,371273.2,371231.8,370959.5,370961.7,371020.3,371048.4,
>
> 371039.9,371086.8,371124.5,371034.7,371040,371025.8,370979.7,370954.3,370966.4,
>
> 370950.7,370914.1,370828.4,371294.4,371446.4,371468,371294.5,371276.1,371268.5,
>
> 371213.2,371055.2,371073.1,371333.6,371347.7,371242.4,371247,371081.9,371054.7,
>
> 370866.6,370885.9,370890,370916.2,370956.4,371125.5,371138.9,371094.8,371048.3,
>
> 371015.8,370990.2,370926.1,370915.8,370902.2,370882.6,371531.1)
> a.y<-c(8247031,8247030,8247100,8247187,8247085,8246710,8246664,8246813,8247012,
>
> 8246970,8246750,8246696,8246644,8246803,8246802,8246748,8246711,8246696,8246624,
>
> 8246509,8246508,8246500,8246539,8246537,8246546,8246601,8246629,8246652,8246782,
>
> 8247238,8246805,8246742,8246811,8246829,8246836,8246903,8247033,8247014,8246730,
>
> 8246660,8246575,8246550,8246806,8246802,8246865,8246828,8246794,8246794,8246754,
>
> 8246474,8246452,8246460,8246533,8246615,8246615,8246661,8246692,8246710,8246738,
>
> 8246867)
> ant.pos<-cbind(a.x,a.y)
>
>
> - This nest are distributed in a stand area:
>
> #Stand contour ----------------------------------
> b.x<-c(371299.9,371266.4,371205.6,371111.8,371047.6,371018.2,371014,371009.3,
>
> 370983.1,370919.7,370853.6,370785.6,370748.8,370711.8,370687.8,370696.4,370785.9,
>
> 370885.5,371035.8,371148.1,371205.2,371231.7,371236.5,371240.3,371285.8,371326.5,
>
> 371397.2,371417.1,371432.9,371445,371455.7,371466.4,371476.6,371502.6,371536,371550,
>
> 371546.8,371528.3,371470,371393.3,371299.9)
> b.y<-c(8246589,8246560,8246508,8246428,8246373,8246349,8246348,8246352,8246385,
>
> 8246465,8246551,8246638,8246685,8246732,8246764,8246771,8246846,8246932,8247062,
>
> 8247160,8247209,8247230,8247224,8247221,8247160,8247107,8247016,8246991,8246967,
>
> 8246939,8246914,8246892,8246875,8246846,8246821,8246809,8246802,8246785,8246735,
>
> 8246669,8246589)
> bord<-cbind(b.x,b.y)
>
>
> - I suspect that this ant nests have inhomogeneous Poisson distribution,
> than I make a CSR for inhomogeneous Poisson process using spastat package:
>
> require(spatstat)
> require(gpclib)
> require(rgeos)
> #Spatstat object
> ---------------------------------------------------------------------
> a <- as(bord, "gpc.poly")
> w <- as.owin(a)
> formi05.ppp<-ppp(x=coordinates(ant.pos)[,1],y=coordinates(ant.pos)[,2],window=w)
>
> plot(formi05.ppp,main="")
>
> ## CSR inhomogeneous
> test-----------------------------------------------------
> env.formi<-envelope(formi05.ppp,nsim=999,fun=Kinhom)
> par(mfrow=c(1,2))
> plot(env.formi,lwd=list(3,1,1,1),xlim=c(0,200), main="Stand 4")
> #
>
>
> - Ok, I have a  inhomogeneous Poisson process. If my process is
> homogeneous Poisson is very simple generate a random point pattern for
> ant nests using rpoispp() function, but now based in my results, what is
> the better approach for find in observed ant nests position the
> retention probability for lambda(x,y)/lmax, when my purposes is create a
> simulation of inhomogeneous point process?

Considerable confusion is indicated by your code and your question.
Note that there is no need for rgeos or gpclib.

(1) You can create your point pattern much more simply by:

W <- owin(poly=list(x=rev(b.x),y=rev(b.y))
X <- ppp(a.x,a.y,window=W)

(2) But ... there is a warning indicating that two points lie outside 
the window; some investigation indicates that these are points 30 and 
60.  You should make sure you understand why these points are present in 
the data.  Something is wrong; it would appear that your specified 
observation window was *not* the actual observation window.

(3) Try again:
a.x <- a.x[-c(30,60)]
a.y <- a.y[-c(30,60)]
X   <- ppp(a.x,a.y,window=W)

(4) Now is there evidence that the pattern is inhomogeneous?  You 
*don't* start with Kinhom()!!! That is for analysing a pattern when you 
have concluded that the pattern is inhomogeneous and you want to decide 
whether there is any interaction between the points.

Testing for inhomogeneity is not all that easy, especially when/if you 
have no particular alternative hypothesis in mind.  The function 
quadrat.test() can be applied, although this test is not particularly 
powerful:

quadrat.test(X,method="MonteCarlo") # The default method produces a
                                     # warning about expected counts
                                     # being too small.

I got a p-value of 0.121 --- no real evidence against a homogeneous 
Poisson process.

Now investigate dependence on the Cartesian coordinates:

cdf.test(X,"x")

gives a p-value of 0.3214, but

cdf.test(X,"y")

gives a p-value of 0.01545 --- here is a *bit* of evidence against 
constant intensity.  However if we try

fit0 <- ppm(X ~ 1)
fit1 <- ppm(X ~ y)
anova(fit0,fit1,test="Chi")

we only get a p-value of 0.05472 --- so "not quite significant" by this 
test.

Finally we might do

E <- envelope(X,savefuns=TRUE)
dclf.test(E)

which gave me p-value of 0.15.  This test is not really appropriate for 
inhomogeneity.  It is essentially testing for the presence of 
interaction *assuming* homogeneity.  However the fact that the null
hypothesis is not rejected says that the data are *consistent* with the 
assumption of a homogeneous Poisson process.

Thus the evidence against a homogeneous Poisson process is meagre at best.

(5) I am not really sure what you are asking in your question at the end 
of the post.  I *think* that you just want to simulate inhomogeneous 
Poisson patterns that "imitate" the behaviour of your observed pattern. 
  In order for that to make much sense you really
should have some model in mind for the intensity of your observed pattern.

One way to proceed without such a model would be to apply smoothing 
methods to get a non-parametric estimate of the intensity of your 
observed pattern:

est.int <- density(X)
Y <- rpoispp(est.int)

This produces a simulated realisation of an inhomogeneous Poisson 
process with intensity equal to the non-parametric estimate of the 
intensity of X.

Another way would be to use the model fitted in terms of the Cartesian
coordinate "y" --- given that there was a little evidence that this
coordinate has an influence on the intensity of X.

Y <- rmh(fit1)

Whether either of these two approaches makes any sense at all depends on 
what you are trying to accomplish --- which is unclear.

cheers,

Rolf Turner

P. S.  Don't send posts like this to r-sig-geo-request at r-project.org; 
that address is for matters pertaining to the administration of the 
r-sig-geo list, e.g. subscribing or unsubscribing.

R. T.

-- 
Rolf Turner
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276
Home phone: +64-9-480-4619