与 rLGCP 模拟和点过程模型拟合相关的一些问题

Question

我尝试从rLGCP 函数生成两个点。我假设这些点在 Window 中的存在由两个 covaiates ras1 和 ras2 控制。因此我需要计算 log-lambda。

rm(list= ls(all=T))
#Libraries
library(spatstat)
library(raster)
library(maptools)
library(fields)

创建域 D 和两个栅格

D <- c(300, 300)  # Square Domaine D of side 300
Win <- owin(xrange =c(0, D[1]), yrange =c(0,D[2])) 
spatstat.options(npixel=c(D[1],D[2]))

ext <- extent(Win$xrange, Win$yrange) # Extent of the rasters
# First raster ras1
par(mfrow=c(1,1))
ras1 <- raster()
extent(ras1) <- ext
res(ras1) <- 10
names(ras1) <- 'Radiation sim'
crs(ras1) <- "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +datum=WGS84"
values(ras1) <- matrix(c(seq(from =0, to =50, length.out=200), seq(from=50, to=100, length.out = 100), seq(from=100, to=150, length.out = 200), seq(from=150, to=200, length.out = 200), seq(from=200, to=290, length.out = 200)), nrow = 30, ncol = 30)
ras1
plot(ras1, asp=1)

# Second Raster ras2
ras2 <- raster()
extent(ras2) <- ext
res(ras2) <- 10
names(ras2) <- 'Precipitation sim'
crs(ras2) <- "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +datum=WGS84"
values(ras2) <- matrix(c(seq(from =-0, to =200, length.out=500), seq(from=400, to=893, length.out = 20), seq(from=200, to=300, length.out = 300),seq(from=300, to = 400, length.out=80)))
ras2
plot(ras2, asp=1)


Rasters.group <- stack(ras1, ras2)
plot(Rasters.group)
graphics.off()

从光栅到即时。对象

im.ras1 <- as.im.RasterLayer(ras1); summary(im.ras1)
im.ras2 <- as.im.RasterLayer(ras2); summary(im.ras2)

covar.list <- list(Radiation.sim=im.ras1, Precipitation.sim=im.ras2)

# plot .im object
par(mfrow=c(1,2))
image.plot(list(x=im.ras1$xcol, y=im.ras1$yrow, z=t(im.ras1$v)), main= "Radiation sim", asp=1)
image.plot(list(x=im.ras2$xcol, y=im.ras2$yrow, z=t(im.ras2$v)), main= "Precipitation sim", asp=1)

现在我可以计算 log-Lambda

#normalization
norm.im.ras1 <- (im.ras1- summary(im.ras1)$mean)/sd(im.ras1) ; summary(norm.im.ras1)
norm.im.ras2 <- (im.ras2- summary(im.ras2)$mean)/sd(im.ras2) ; summary(norm.im.ras2)

#Compute log-lambda

log.lambda <- norm.im.ras1 + 2*norm.im.ras2
summary(log.lambda)

结果显示非常弱的值

像素值 范围 = [-4.657923, 10.94624] 积分 = -9.678445e-12 均值 = -1.075383e-16

当我尝试从 rLGCP 模拟时

gen.lgcp <- rLGCP("matern", mu=log.lambda, var=0.5, scale=0.05, nu=1)

错误：无法分配大小为 181.9 MB 的向量

我试图绕过

log.lambda0 <- as.im(solutionset(log.lambda>0))
gen.lgcp <- rLGCP("matern", mu=log.lambda0, var=0.5, scale=0.05, nu=1)
summary(gen.lgcp) 

我可以继续前进。但进一步，我没有得到相关的结果

#Thinning
image.plot(list(x=log.lambda$xcol, y=log.lambda$yrow, z=t(log.lambda$v)), main= "log.lambda", asp=1)

samp.lgcp <- rthin(gen.lgcp, P=seq(from=0.02, to=0.2, length.out = gen.lgcp$n));  points(samp.lgcp$x, samp.lgcp$y, type = 'p', cex=0.2, lwd=1, col='white')

#point pattern 
pts.locations <- as.data.frame(cbind(longitude=samp.lgcp$x, latitude=samp.lgcp$y))
ppp.lgcp <- ppp(pts.locations$longitude, pts.locations$latitude, window = owin(xrange=c(min(pts.locations [,1]),max(pts.locations [,1])), yrange = c(min(pts.locations[,2]),max(pts.locations[,2]))))
plot(ppp.lgcp)

#Extract value of each sampled point covariate
cov.value  <- extract(Rasters.group, pts.locations)
cov.value <- as.data.frame(cov.value )
presence.data <- data.frame(pts.locations, cov.value, presence=rep(1, nrow(cov.value)))

### Choosing absence point pattern
abs.region <- crop(Virtual.species.domaine, extent(25.28486 , 162.2897 ,181.7417 , 280.7651 ))
im.abs.region <- as.im.RasterLayer(abs.region)
abs.points <- rasterToPoints(abs.region)
ppp.abs.points <- ppp(abs.points[,1], abs.points[,2], window = owin(xrange = c(min(abs.points[,1]), max(abs.points[,1])), yrange =c(min(abs.points[,2]), max(abs.points[,2]))))
plot(ppp.abs.points)

cov.value.abs <- extract(Rasters.group, abs.points[,1:2])
absence.data <- data.frame(abs.points[,1:2], cov.value.abs, presence=rep(0, nrow(abs.points)))
colnames(absence.data)[1:2] <- c("longitude", "latitude")
head(absence.data)
# Get database for LGCP
LGCP.Data.Set <- rbind(presence.data, absence.data)

#' Model
#' we will use non-stationary formula
covar.formula <- as.formula(paste("~", paste(names(LGCP.Data.Set[,3:4]), collapse = "+")))

#Quadrature scheme
Q.lgcp <- quadscheme(ppp.lgcp, ppp.abs.points, method = 'grid')
plot(Q.lgcp)

警告信息： 在 countweights(id, area) 中： 一些具有正面积的图块不包含任何交点：相对误差 = 94.2%

# Inhomogenous poisson process Model
fit.ipp <- ppm(Q.lgcp, trend = covar.formula, covariates = LGCP.Data.Set[,3:4])
summary(fit.ipp)

警告信息： glm.fit: 算法没有收敛

出了什么问题？

我的目标是评估模型和预测

prediction.ipp <- predict.ppm(fit.ipp, log.lambda, type = 'intensity')

Answer 1

这是一个很长且没有重点的问题，但我会尽力提供帮助。

在构建图像log.lambda 后，您说“结果显示非常弱的值”。你什么意思？图像值被分配为从 0 到 200 的一系列值，然后标准化为均值为零，标准差为 1。这有多“弱”？

然后您使用此图像作为平均对数强度调用rLGCP。 log.lambda 的值范围从大约 -4 到 +10。这意味着所需的强度范围从exp(-4) 到exp(+10)，即大约每平方单位的0.01 到20 000 点。图像尺寸为 30 x 30 单位。因此，必须生成大量随机点，但由于内存限制而失败。 （预期点数为integral(exp(log.lambda))。

然后您将log.lambda 更改为另一个仅采用值 0 和 1 的图像。

下一段代码似乎采用了（“不存在”像素的）光栅图像，并尝试使用“不存在”像素作为虚拟点来构建正交方案。这不是quadscheme 的设计目的（对于quadscheme，虚拟点应该是稀疏的）。

您无需构造正交方案即可使用ppm。你可以做类似的事情

D <- solist(A=im.ras.1, B=im.ras.2)
ppm(ppp.logi ~ A+B , data=D)

如果你真的想构造一个正交方案，我建议你改用函数pixelquad。只需执行pixelquad(ppp.lgcp, im.abs.region) 或类似操作即可。然后使用ppm。

由于数据是由 Cox 流程生成的，因此使用 kppm 而不是 ppm 会更合适。

更多信息请参见the spatstat book。