library(spatstat)
#> Loading required package: spatstat.data
#> Loading required package: spatstat.univar
#> spatstat.univar 3.1-3
#> Loading required package: spatstat.geom
#> spatstat.geom 3.4-1
#> Loading required package: spatstat.random
#> spatstat.random 3.4-1
#> Loading required package: spatstat.explore
#> Loading required package: nlme
#> spatstat.explore 3.4-3
#> Loading required package: spatstat.model
#> Loading required package: rpart
#> spatstat.model 3.3-6
#> Loading required package: spatstat.linnet
#> spatstat.linnet 3.2-5
#>
#> spatstat 3.3-2
#> For an introduction to spatstat, type 'beginner'
# First set of points
set.seed(1)
X1 <- rpoispp(lambda = \(x, y) exp(4 + x))
npt1 <- npoints(X1)
ppm(X1 ~ x)
#> Nonstationary Poisson process
#> Fitted to point pattern dataset 'X1'
#>
#> Log intensity: ~x
#>
#> Fitted trend coefficients:
#> (Intercept) x
#> 4.0207596 0.8124458
#>
#> Estimate S.E. CI95.lo CI95.hi Ztest Zval
#> (Intercept) 4.0207596 0.2408283 3.5487447 4.492774 *** 16.695542
#> x 0.8124458 0.3798072 0.0680374 1.556854 * 2.139101
# Second set of points
X2 <- rpoispp(lambda = \(x, y) exp(4 + x))
npt2 <- npoints(X2)
ppm(X2 ~ x)
#> Nonstationary Poisson process
#> Fitted to point pattern dataset 'X2'
#>
#> Log intensity: ~x
#>
#> Fitted trend coefficients:
#> (Intercept) x
#> 4.2655978 0.6626184
#>
#> Estimate S.E. CI95.lo CI95.hi Ztest Zval
#> (Intercept) 4.2655978 0.2174979 3.83930972 4.691886 *** 19.612131
#> x 0.6626184 0.3486133 -0.02065114 1.345888 1.900726
# Superimpose
X <- superimpose(X1 = X1, X2 = X2)
qs <- quadscheme(unmark(X))
ppm(qs, trend = ~ 1 + x)
#> Nonstationary Poisson process
#> Fitted to point pattern dataset 'qs'
#>
#> Log intensity: ~1 + x
#>
#> Fitted trend coefficients:
#> (Intercept) x
#> 4.8433365 0.7311765
#>
#> Estimate S.E. CI95.lo CI95.hi Ztest Zval
#> (Intercept) 4.8433365 0.1614216 4.526956 5.159717 *** 30.004266
#> x 0.7311765 0.2568024 0.227853 1.234500 ** 2.847234
# Now I would like to recover the intensity of X1 with quadscheme, modifying the
# weights associated to the quadrature scheme. Based on my understanding, I
# would need to
# 1) Give approx 0 weight to the points in X which belong to X2 and recompute
# the weights of the points in X1
qs1 <- quadscheme(X1)
isd <- is.data(qs)
qs$w[isd] <- qs$w[isd] * c(rep(1, npoints(X1)), rep(1e-10, npoints(X2)))
# The following replaces the weights associated to the points in X which
# originated from X1 removing the influence of X2
isX1 <- c(rep(TRUE, npoints(X1)), rep(FALSE, npoints(qs$data) + npoints(qs$dummy) - npoints(X1)))
qs$w[isX1] <- qs1$w[seq_len(npoints(X1))]
# 2) Remove the influence of X2 from the weights assigned to the dummy points.
isd1 <- is.data(qs1)
qs$w[!isd] <- qs1$w[!isd1]
qs$param <- NULL
# Refit the model.
ppm(qs, trend = ~ 1 + x)
#> Nonstationary Poisson process
#> Fitted to point pattern dataset 'qs'
#>
#> Log intensity: ~1 + x
#>
#> Fitted trend coefficients:
#> (Intercept) x
#> 4.8434199 0.7310179
#>
#> Estimate S.E. CI95.lo CI95.hi Ztest Zval
#> (Intercept) 4.8434199 0.1614071 4.5270679 5.159772 *** 30.007480
#> x 0.7310179 0.2567734 0.2277514 1.234284 ** 2.846938
# The real parameters should be (Intercept) = 4 and x = 1. CI are not ok.Created on 2025-05-30 with reprex v2.1.1.9000
Updates after the mails from Adrian and Ege:
Created on 2025-06-03 with reprex v2.1.1.9000