1

I have two points pattern (ppp) objects p1 and p2. There are X and Y points in p1 and p2 respectively. I have fitted a ppm model (with location coordinates as independent variables) in p1 and then used it to predict "intensity" for each of the Y points in p2.

Now I want to get the probability for event occurrence at that point/zone in p2. How can I use the predicted intensities for this purpose?

Can I do this using Spatstat? Are there any other alternative.

Lesnar
  • 501
  • 3
  • 16

2 Answers2

0

Your question, objective and current method are not very clear to me. It would be beneficial, if you could provide code and graphics, that explains more clearly what you have done, and what you are trying to obtain. If you cannot share your data you can use e.g. the built-in dataset chorley as an example (or simply simulate artificial data):

library(spatstat)
plot(chorley, cols = c(rgb(0,0,0,1), rgb(.8,0,0,.2)))


X <- split(chorley)
X1 <- X$lung
X2 <- X$larynx
mod <- ppm(X1 ~ polynom(x, y, 2))
inten <- predict(mod)
summary(inten)
#> real-valued pixel image
#> 128 x 128 pixel array (ny, nx)
#> enclosing rectangle: [343.45, 366.45] x [410.41, 431.79] km
#> dimensions of each pixel: 0.18 x 0.1670312 km
#> Image is defined on a subset of the rectangular grid
#> Subset area = 315.291058349571 square km
#> Subset area fraction = 0.641
#> Pixel values (inside window):
#>  range = [0.002812544, 11.11172]
#>  integral = 978.5737
#>  mean = 3.103715
plot(inten)

Predicted intensities at the 58 locations in X2

intenX2 <- predict.ppm(mod, locations = X2)
summary(intenX2)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.1372  4.0025  6.0544  6.1012  8.6977 11.0375

These predicted intensities intenX2[i] say that in a small neighbourhood around each point X2[i] the estimated number of points from X1 is Poisson distributed with mean intenX2[i] times the area of the small neighbourhood. So in fact you have estimated a model where in any small area you have a probability distribution for any number of points happening in that area. If you want the distribution in a bigger region you just have to integrate the intensity over that region.

To get a better answer you have to provide more details about your problem. Created on 2018-12-12 by the reprex package (v0.2.1)

Ege Rubak
  • 4,347
  • 1
  • 10
  • 18
  • Thank you for your response! I have edited the question with code examples. Can you please tell me how to solve this. I am still stuck! – Lesnar Dec 13 '18 at 05:39
0

The intensity is the expected number of points per unit area. In small areas (such as pixels) you can just multiply the intensity by the pixel area to get the probability of presence of a point in the pixel. 

fit <- ppm(p1, .......)
inten <- predict(fit)
pixarea <- with(inten, xstep * ystep)
prob <- inten * pixarea

This rule is accurate provided the prob values are smaller than about 0.4.

In a larger region W, the expected number of points is the integral of the intensity function over that region:

EW <- integrate(inten, domain=W)

The result EW is a numeric value, the expected total number of points in W. To get the probability of at least one point,

P <- 1- exp(-EW)

You can also compute prediction intervals for the number of points, using predict.ppm with argument interval="prediction".

Adrian Baddeley
  • 1,956
  • 1
  • 8
  • 7
  • Thank you for your response! I have edited the question with code examples. Can you please tell me how to solve this. I am still stuck! – Lesnar Dec 13 '18 at 05:40