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