拟合LGCP模型遇计算奇异错误，求问题排查与解决方案

问题重现 #

拟合对数高斯Cox过程（LGCP）模型时，触发以下计算奇异错误及警告：

Inhomogeneous Cox point process model
Fitted to point pattern dataset ‘Y’
Fitted by maximum second order composite likelihood
        rmax = 153.75
        weight function: Indicator(distance <= 76.875)

Error in solve.default(M) : 
  system is computationally singular: reciprocal condition number = 4.3093e-22
Error in solve.default(M) : 
  system is computationally singular: reciprocal condition number = 4.3093e-22
In addition: Warning message:
Cannot compute variance: Fisher information matrix is singular 
Log intensity:  ~DF + DP + (x + y + I(x^2) + I(x * y) + I(y^2))

Fitted trend coefficients:
  (Intercept)            DF            DP             x             y 
-1.025219e+03 -3.115902e-02 -2.003711e-02 -1.016590e-01  4.949720e-01 
       I(x^2)      I(x * y)        I(y^2) 
-3.557955e-05  2.898502e-05 -6.005674e-05 

Cox model: log-Gaussian Cox process
        Covariance model: exponential
Fitted covariance parameters:
      var     scale 
 2.024728 18.913952 
Fitted mean of log of random intensity: [pixel image]
Warning message:
Cannot compute variance: Fisher information matrix is singular

已执行Y <- rescale(X,1000)缩放数据，尝试移除协变量后问题仍存在，拟合采用《空间点模式》推荐的mincon、palm、clik复合似然方法。

排查方向与解决方案 #

1. 解决协变量共线性问题 #

• 二次项（I(x^2)、I(x*y)、I(y^2)）与一次项x、y易产生严重共线性，即使做了数据缩放也可能无法缓解。可通过car包的vif()函数计算方差膨胀因子（VIF>10则共线性显著），或直接计算协变量矩阵的条件数验证。
• 处理方案：
  • 简化趋势模型：先移除交互项I(x*y)，仅保留x、y及各自二次项；若仍有问题，进一步简化为仅线性项~x+y，验证是否还会触发奇异错误。
  • 坐标中心化：将x、y转换为均值为0的中心化变量（如x_centered <- x - mean(x)），再构建二次项，能大幅降低共线性程度。

2. 调整复合似然参数设置 #

• 信息矩阵奇异可能源于权重函数或rmax的设置：
  • 调整rmax值：当前rmax=153.75，可尝试缩小或扩大范围，比如设置为数据最近邻距离的中位数2倍，或直接测试rmax=100、rmax=200等取值。
  • 更换平滑权重函数：替代默认的指示函数权重，在调用lgcp()时指定weight="biweight"或weight="triweight"，避免权重突变导致的矩阵计算问题。

3. 检查数据空间分布与质量 #

• 极端聚类、大面积空白区域、重复点或异常坐标都可能引发矩阵奇异：
  • 用plot(Y)直观查看空间分布，若存在严重不均匀聚集，可考虑划分子区域分别拟合，或加入断层分布、地壳应力等更贴合的空间协变量解释异质性。
  • 清理数据：用unique(Y)去除重复点，用range(Y$x)、range(Y$y)检查坐标范围是否合理，排除异常值。

4. 更换复合似然方法 #

• 尝试《空间点模式》中提到的其他复合似然实现：
  • 用method="palm"替代默认的method="clik"，或调整mincon参数（如mincon=0.1）限制信息矩阵的最小条件数。示例调用：

    fit <- lgcp(Y, trend=~DF+DP+x+y+I(x^2)+I(y^2), method="palm", rmax=100)
    

验证步骤 #

1. 先拟合最简模型：lgcp(Y, trend=~1)，确认是否还会出现奇异错误。
2. 逐步添加协变量（先线性项，再二次项，最后交互项），每一步检查模型是否能正常输出方差信息。
3. 对比不同rmax和权重函数下的结果，选择信息矩阵非奇异的参数组合。

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