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

Question

这个问题对你们中的一些人来说可能是微不足道的，但在这里。我实际上是在尝试用一个截距和一个变量来模拟一个简单的非同质点过程。我的目标只是使用包 spatstat 和 rstan 恢复这两个系数。

这是模拟的代码：

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 的强度的log 被假定为的线性函数协变量Z。因此强度函数的形式为exp(a + b * Z)，其中a 和b 是在调用ppm 时估计的参数。

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

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