rLGCP模拟及点过程模型拟合报错、结果异常问题咨询

我尝试使用rLGCP函数生成点，假设这些点在观测窗口内的分布受两个协变量ras1和ras2控制，因此需要先计算log-lambda。

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

步骤1：创建研究域D与两个栅格 #

D <- c(300, 300)  # 边长为300的正方形研究域D
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) # 栅格的范围
# 第一个栅格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)

# 第二个栅格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()

步骤2：将栅格转换为im对象 #

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)

# 绘制im对象
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)

步骤3：计算log-Lambda #

# 标准化
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)

# 计算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) 

虽然可以继续运行，但后续得到的结果不符合预期：

# 稀疏化
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')

# 点模式
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)

# 提取每个采样点的协变量值
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)))

### 选择缺失点模式
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)
# 构建LGCP数据集
LGCP.Data.Set <- rbind(presence.data, absence.data)

#' 模型
#' 我们使用非平稳公式
covar.formula <- as.formula(paste("~", paste(names(LGCP.Data.Set[,3:4]), collapse = "+")))

# quadrature方案
Q.lgcp <- quadscheme(ppp.lgcp, ppp.abs.points, method = 'grid')
plot(Q.lgcp)

警告信息：
In countingweights(id, areas) :
部分面积为正的瓦片不包含任何正交点：相对误差 = 94.2%

# 非齐次泊松过程模型
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')

操作中的问题主要集中在以下4个环节：

• 首先是spatstat.options(npixel=c(D[1],D[2]))的设置不合理：将全局像素分辨率设为300×300后，rLGCP模拟时需要生成对应尺寸的高斯随机场网格，300×300共9万个点，加上协方差矩阵计算的内存开销，自然会出现内存分配报错。该精度远高于你栅格10m分辨率、300m范围的实际数据精度，把这个参数改成spatstat.options(npixel=c(100,100))即可。
• 其次是log.lambda的定义错误：你标准化后的log lambda均值几乎为0，意味着大部分区域的强度exp(log.lambda)≈1，且后续直接把log.lambda>0的二值掩膜作为mu传入rLGCP，完全丢失了协变量的梯度信息，模拟出来的点分布和你预设的协变量控制逻辑完全无关，后续拟合自然不可能得到符合预期的结果。正确的做法是给log lambda加一个常数偏移量，调整整体点密度到合理范围，比如log.lambda <- 2 + norm.im.ras1 + 2*norm.im.ras2，保证整体平均强度在exp(2)≈7.39每单位面积，控制模拟点总数在合理区间，也能减少内存压力。
• 第三是 quadrature 方案构造错误：你自定义的存在点和缺失点的窗口范围不一致，且缺失点直接用栅格转点的网格点，和quadscheme的grid方法逻辑冲突，导致大部分网格瓦片没有匹配的积分点，才会出现94.2%的误差警告。你不需要手动构造缺失点，直接用Q.lgcp <- quadscheme(ppp.lgcp, method = 'grid', nd = 50)即可，spatstat会自动生成匹配窗口的积分点。
• 第四是ppm拟合时协变量传入错误：你直接传入数据框形式的协变量，ppm要求协变量是和窗口匹配的im对象或者命名列表，你已经构造了covar.list，直接传入covariates = covar.list即可，否则模型无法匹配每个积分点对应的协变量值，自然无法收敛。

修正后的核心代码参考 #

# 修正全局像素设置
spatstat.options(npixel=c(100,100))

# 修正log.lambda计算
log.lambda <- 2 + norm.im.ras1 + 2*norm.im.ras2

# 正常模拟LGCP
gen.lgcp <- rLGCP("matern", mu=log.lambda, var=0.5, scale=0.05, nu=1)

# 构造积分方案不需要手动生成缺失点
Q.lgcp <- quadscheme(gen.lgcp, method = 'grid', nd = 50)

# 拟合模型传入正确的协变量列表
fit.ipp <- ppm(Q.lgcp, trend = ~ Radiation.sim + Precipitation.sim, covariates = covar.list)

# 预测
prediction.ipp <- predict(fit.ipp, type = 'intensity')

内容的提问来源于stack exchange，提问作者Luc