模拟非齐次泊松点过程并检索协变量系数

Question

这个问题对于你们中的一些人来说可能是微不足道的，但事实就是如此。我实际上试图用截距和单个变量来模拟一个简单的非齐次点过程。我的目标只是使用 spatstat 和 rstan 包恢复这 2 个系数。

这是模拟的代码：

library(spatstat)
library(sf)
library(sp)
library(maptools)
library(raster)
library(fields)
library(rstan)
library(tidyverse)

# Generate species distribution
genDat_pp <- function(b1, b2, dim, plotdat = TRUE){

# Define the window of interest
win <- owin(c(0,dim[1]), c(0,dim[2]))

# set number of pixels to simulate an environmental covariate
spatstat.options(npixel=c(dim[1],dim[2]))

y0 <- seq(win$yrange[1], win$yrange[2],
        length=spatstat.options()$npixel[2])
x0 <- seq(win$xrange[1], win$xrange[2],
        length=spatstat.options()$npixel[1])
multiplier <- 1/dim[2]

# Make the environmental covariate
gridcov <- outer(x0,y0, function (x,y) multiplier*y + 0*x)

# Set the coefficients
beta0 <- b1
beta1 <- b2

# Simulate the point pattern
pp <- rpoispp(im(beta0 + beta1*gridcov, xcol=x0, yrow=y0))
qcounts <- quadratcount(pp, ny=dim[1], nx=dim[2])
dens <- density(pp)
Lambda <- as.vector(t(qcounts))

if(plotdat == TRUE){
  par(mfrow=c(1,2), mar=c(2,2,1,1), mgp=c(1,0.5,0))
  plot(im(gridcov), main = 'Covariate')
  plot(dens, main = 'Intensity')
}
# Return a list with which I can play with
return(list(Lambda = Lambda, pp = pp, gridcov = gridcov))
}

现在，模拟似乎可以工作并返回连贯的图。然而，当我尝试使用 spatstat“分析”数据时，我无法恢复截距或 beta 系数：

b1 <- 2
b2 <- 5
dim <- c(20,20)

# Generate data
pp <- genDat_pp(b1, b2, dim)

# Fit a poisson point process model in spatstat
cov <- im(pp$gridcov)
fit <- ppm(pp$pp ~ 1 + cov)
summary(fit)

在这个例子中，对象“fit”告诉我截距是 2.65（相对来说还不错），而 beta 系数是 2.68，这是完全错误的。

是否有任何错误或者我只是以错误的角度解释结果？

我真的很感激任何人能给出答案！

本

Answer 1

这个问题涉及到spatstat包。

在您的示例中，强度函数的形式a + b*Z为Z协变量函数， 和a是b数字。

在模型拟合函数 中ppm，模型公式描述了强度函数的对数。模型公式X ~ 1 + Z或等效公式X ~ Z（1 在这种情况下是多余的）意味着点模式强度的对数X被假定为协变量 的线性函数Z。因此强度函数的形式为exp(a + b * Z)其中a和是调用b时将被估计的参数。ppm

您拟合的模型lambda = exp(a+bZ)与您模拟的模型 不一致lambda = a+bZ，因此结果存在差异。

这在spatstat 书的帮助文件ppm和第 9 章中有进一步解释。