Dear Rolf,

Thank you for insightful answer. Following your explanation I found out that the key problem was that I didn't specify the simulation in Linhom envelope. For some reason I thought that lambda would be overriding the envelope creating point process. So, your example 4 was really helpful for that matter.

However in order to produce a proper-looking graph with an envelope, I had to use reciplambda instead of the regular one. I figured this out after seeing that while envelope was calculated and displayed correctly, the observed Linhom(r) line did not appear (envelope's "obs" values are either 0's or NaN's). Using reciplambda did the trick though, but the values on the y-axis are 1e+07 times smaller than they should. Creating random patterns using rpoispp() does not indicate at all that the intensity surface should be used reciprocally: the number of points and their spatial distribution make sense. Also if using randomly created points there is no need for reciplambda, although it works both ways showing slightly different results. Another worrying issue, but I'll try to look into it.


Thank you and all the best,
Allar


On 23/07/2013 11:07, Rolf Turner wrote:

Dear Allar,

It is impossible for me to tell from the information that you have provided what is going wrong. Basically it seems to me that the code you used will work "in general". It is strange that you are getting *empty* plots. A toy example (see below) indicates that "empty plots"
should not result.

I can see why there might be problems with the envelope, but the plot shouldn't be
*empty*.  At the very least you should get a plot of Linhom(sites).

Toy example:

W <- disc(0.5,c(0.5,0.5),npoly=6)
A  <- as.im(function(x,y){0.1*cos(2*pi*(x^2+y^2))^2},W=W)
i    <- summary(A)$integral
B   <- eval.im(100*A/i)
set.seed(42)
X  <- rpoispp(B)[W]
E1 <- envelope(X,Linhom,lambda=A,nsim=19)
E2 <- envelope(X,Linhom,lambda=B,nsim=19)
E3 <- envelope(X,Linhom,lambda=A,simulate=expression(rpoispp(A)),nsim=19)
plot(E1)
plot(E2)
plot(E3)

The envelopes E1 and E2 are silly looking, being inappropriate, but they are
certainly not empty.

The envelope E3 yields a non-empty plot, although there is no envelope as such. This is due to the fact that the "simulate" expression yields nothing but empty
patterns (from which no L function estimate can be calculated).

I suspect that what you actually did was something like the E3 calculation rather
than what you told us you did.

What you probably should have done is something like:

E4 <- envelope(X,Linhom,lambda=B,simulate=expression(rpoispp(B)),nsim=19)

which produces a sensible envelope.

What Linhom() does for you is to allow you to assess whether "sites" arises from an inhomogeneous Poisson process or if there is evidence of interaction in the data. To apply Linhom() you must know (or have an idea of) the underlying intensity of "sites". To use an intensity such as "model" which you apparently *know* would generate fewer points than there are in "sites" makes no sense. Use an intensity which you actually believe to be the underlying intensity for "sites" (or which you at least believe to be a
*plausible* candidate for this intensity).

    cheers,

        Rolf Turner

P. S. The error that you got from applying Linhom() directly to "sites" with lambda=model might have arisen because "model" has some zero values. I *think* this will be a problem only if it has zero values at one or more points of "sites". Haven't experimented with
this.

        R. T.

On 23/07/13 13:23, Allar Haav wrote:
Hi all,

I am currently trying to use spatstat's inhomogeneous cluster analysis methods (mostly Linhom), but ran into a problem. Namely, I have a pixel image (type "im") with values ranging from 0 to 1 indicating point probabilities. When creating and plotting random points with it, this intensity surface works as it should. I can also use it in regular envelope functions without issues. However, any attempt to use the intensity surface with Linhom function results in empty plot. I suspected that this might have something to do with the integral of the surface, therefore made a crude transformation so that it would be approximately the number of points under study (see below), but of course it didn't change anything. Generally I still feel that the crux lies in the intensity surface's values, but I can't quite put my finger on it.

A simplified code I used is as follows:

area <- as.owin(area_polygon)
sites <- as.ppp(as(readShapeSpatial('sites.shp')))
sites$window <- area
model <- as.im(read.asciigrid('model.asc'))
sites.n <- sites$n
integral <- summary(model)$integral

# This creates a nice plot with simulation envelopes
plot(envelope(sites, Lest, nsim=19, simulate=expression(rpoispp(eval.im(model*(sites.n/integral ))))))

# These produce only empty plots
plot(envelope(sites, Linhom, nsim=19, lambda=model ))
plot(envelope(sites, Linhom, nsim=19, lambda=eval.im(model *(sites.n/integral ))))

# A further check without envelope
Linhom(sites, lambda=model)
# This resulted in:
"Error in Kwtsum(dIJ, bI, wIJ, b, w = reciplambda, breaks) : Weights in K-function were infinite or NA"

That last error message left me a bit perplexed. Indeed, I have NA values in the dataset as the study area is not perfectly rectangular. But this should be a non-issue anyway. Frankly I'm running out of ideas. I've been thinking of alternative approaches, e.g. creating random points from the intensity surface and fitting a model based on it using "ppm", perhaps using some covariates. Yet this seems actually an awful idea considering the probable loss of accuracy and the fact that basically the model already exists.

As a side note, a bit of googling on the issue revealed an older post about inhomogeneous process analysis (http://r-sig-geo.2731867.n2.nabble.com/memory-allocation-problem-with-envelope-function-in-spatstat-td5029218.html). There a command equivalent to my first envelope plot was used: Lin=envelope(SPECIES.ppp,Linhom,nsim=5,simulate=expression(rpoispp(lambda)), ! + correction="border") That made me think that I might be overthinking the problem and I've already found the solution with first try. However I am still skeptical as the difference between Lest and Linhom does not lie merely in simulated patterns - weights are also to be considered. If so, why was this solution used?

Any suggestions or ideas? All sort of feedback would be appreciated.

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Reply via email to