mgcv包"incorrect number of dimensions"错误排查求助

环境与错误信息 #

• 运行环境：Linux Ubuntu 22.04，R 4.3.1，mgcv包
• 报错内容：

Error in backsolve(L, forwardsolve(t(L), (D * A)[piv, ]))[ipiv, , drop = FALSE] : 
  incorrect number of dimensions

模型代码 #

model_ndvi <- gam(
  list(
    NDVI ~ cutblock + s(long, lat, bs = 'ds', k = 100, by = cutblock) + 
      s(year, bs = 'cr', k = 7, by = cutblock) + s(doy, bs = 'cc', k = 8, by = cutblock) + 
      s(cb_age, bs = "cr", k = 3) + ti(year, doy, bs = 'cc') + ti(cb_age, year, bs = 'cr'),
    ~ cutblock + s(long, lat, bs = 'ds', k = 50, by = cutblock) + 
      s(year, bs = 'cr', k = 5, by = cutblock) + s(doy, bs = 'cc', k = 4, by = cutblock)+ 
      s(cb_age, bs = "cr", k = 3) + ti(year, doy, bs = 'cc') + ti(cb_age, year, bs = 'cr')
  ),
  family = betals(), data = babyd, method = 'REML',
  knots = list(doy = c(0.5, 366.5)), control = gam.control(nthreads = 30), trace = TRUE
)

数据情况 #

• 原数据集大小87GB，测试使用90MB子集（共1,474,926行）
• 数据片段（head()输出）：

long      lat        NDVI    cutblock harvestyr  date    dec_date year doy cb_age
1 -127.2368 51.60493 -0.0275    FALSE        NA 2000-02-18 2000.131 2000  49     NA
2 -127.2368 51.60712 -0.0053    FALSE        NA 2000-02-18 2000.131 2000  49     NA
3 -127.2368 51.60932  0.1173    FALSE        NA 2000-02-18 2000.131 2000  49     NA
4 -127.2368 51.61151  0.2855    FALSE        NA 2000-02-18 2000.131 2000  49     NA
5 -127.2368 51.61371  0.2417    FALSE        NA 2000-02-18 2000.131 2000  49     NA
6 -127.2368 51.61591  0.4240    FALSE        NA 2000-02-18 2000.131 2000  49     NA

• 数据结构（str()输出）：

'data.frame':   1474926 obs. of  10 variables:
 $ long     : num  -127 -127 -127 -127 -127 ...
 $ lat      : num  51.6 54 53.4 52.9 52.1 ...
 $ NDVI     : num  -0.0275 0.4184 -0.0103 -0.0292 0.5711 ...
 $ cutblock : Factor w/ 2 levels "FALSE","TRUE": 1 1 1 1 1 1 1 1 1 1 ...
 $ harvestyr: num  NA NA NA NA NA NA NA NA NA NA ...
 $ date     : Date, format: "2000-02-18" "2000-02-18" "2000-02-18" "2000-02-18" ...
 $ dec_date : num  2000 2000 2000 2000 2000 ...
 $ year     : int  2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 ...
 $ doy      : int  49 49 49 49 49 49 49 49 49 49 ...
 $ cb_age   : int  NA NA NA NA NA NA NA NA NA NA ...

• 关键变量说明：cutblock为二分类因子，cb_age和harvestyr存在大量缺失值，NDVI包含负值。

问题背景 #

该模型用于分析森林皆伐区域（cutblock=TRUE）与非皆伐区域（cutblock=FALSE）的NDVI长期趋势，探究皆伐区域恢复至全生产力的时间点，使用betals()同时建模NDVI的均值和方差变化。

已做排查 #

已检查所有变量的唯一值数量，确认满足模型拟合要求。

解决方案建议 #

1. 移除无信息变量：从数据片段和结构看，cb_age全为NA，模型中包含s(cb_age, ...)和ti(cb_age, year, ...)会导致设计矩阵退化，先移除这些项，测试模型能否正常运行，再逐步添加有效变量。
2. 修复响应变量范围：标准beta分布要求响应值在(0,1)区间，但你的NDVI存在负值，betals()可能无法处理。先将NDVI缩放至(0,1)：

   babyd$NDVI_scaled <- (babyd$NDVI - min(babyd$NDVI, na.rm=TRUE)) / 
     (max(babyd$NDVI, na.rm=TRUE) - min(babyd$NDVI, na.rm=TRUE))
   

   然后用NDVI_scaled作为响应变量重新拟合模型。
3. 简化模型逐步验证：先运行最简版betals模型，比如：

   model_simple <- gam(list(NDVI_scaled ~ cutblock + s(year, bs='cr'), ~ cutblock),
                       family=betals(), data=babyd, method='REML')
   

   确认模型能正常拟合后，再逐步添加平滑项、交互项和by参数。
4. 降低平滑项复杂度：部分平滑项的k值设置过高（如ds平滑k=100），对于子集数据可能导致矩阵维度问题，先降低k值（如k=30）测试。
5. 禁用多线程：暂时移除control = gam.control(nthreads = 30)，多线程在复杂模型下可能引发内存或矩阵操作异常，单线程运行排查问题。
6. 更新mgcv包：运行update.packages("mgcv")更新到最新版本，旧版本可能存在betals()或多线程相关的bug。

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