模拟非齐次泊松点过程并检索协变量系数
Simulate inhomogenous poisson point process and retrieve the covariate coefficients
这个问题对你们中的一些人来说可能微不足道,但就在这里。我实际上是在尝试用截距和单个变量模拟一个简单的非齐次点过程。我的目标只是使用包 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))
}
现在,模拟似乎工作正常并且 return 连贯的情节。但是,当我尝试使用 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,这是完全错误的。
有没有错误或者我只是用错误的角度解释结果?
如果有人能给出答案,我将不胜感激!
本
此问题与 spatstat
包有关。
在您的示例中,强度函数的形式为 a + b*Z
,其中 Z
是协变量函数,a
和 b
是数字。
在model-fitting函数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 章中有进一步解释。
这个问题对你们中的一些人来说可能微不足道,但就在这里。我实际上是在尝试用截距和单个变量模拟一个简单的非齐次点过程。我的目标只是使用包 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))
}
现在,模拟似乎工作正常并且 return 连贯的情节。但是,当我尝试使用 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,这是完全错误的。
有没有错误或者我只是用错误的角度解释结果?
如果有人能给出答案,我将不胜感激!
本
此问题与 spatstat
包有关。
在您的示例中,强度函数的形式为 a + b*Z
,其中 Z
是协变量函数,a
和 b
是数字。
在model-fitting函数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 章中有进一步解释。