使用m=1时bam模型summary无法运行的原因排查（广义可加模型）

我想用bam()的method="fREML"加速分析，代码如下：

require(mgcv)
require(gratia)
require(brms)
require(DHARMa)
require(marginaleffects)

mode11 <- bam(Met1 ~ s(Timepoint, k=7) + s(Timepoint, by=subject, k=7) + s(subject, bs="re"), data=df, method="fREML")

model2 <- bam(Met1 ~ s(Timepoint, k=7) + s(Timepoint, by=subject, k=7, m=1) + s(subject, bs="re"), data=df, method="fREML")

两个模型都能正常运行，但调用summary(model2)时出现以下错误：

Error in h(simpleError(msg, call)) : 
 error in evaluating the argument 'x' in selecting a method for function 'chol': non-conformable arrays

m=1应该是问题根源，但不清楚具体原因——用gam()设置m=1时，summary能正常运行。

补充：未找到同类案例，推测和数据有关，但不知道怎么评估model2的输出来定位问题，以下是数据的dput结果：

structure(list(subject = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L), levels = c("1", 
"3", "4", "5", "6", "7", "8", "10", "11", "12", "16", "17", "19", 
"21", "22", "25", "26", "28", "35", "36", "37", "38", "39", "40", 
"42", "43", "44", "45", "46", "48", "50", "51", "52", "55", "56", 
"57", "58", "61", "65", "71"), class = "factor"), Timepoint = c(1L, 
2L, 3L, 4L, 5L, 6L, 7L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 1L, 2L, 3L, 
4L, 5L, 6L, 7L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 1L, 2L, 3L, 4L, 5L, 
6L, 7L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 1L, 2L, 3L, 4L, 5L, 6L, 7L
), Met1 = c(0.012461309, 0.01794298, 0.003863203, -0.006398583, 
0.001760362, 0.000870293, 0.000830637, -0.014258656, -0.01946732, 
0.00389311, 0.000199752, -0.0122629, 0.007792022, 0.007019363, 
0.006641794, 0.012690831, 0.004881705, -0.004663227, 0.002341043, 
-0.008827367, -0.009172269, 0.038058982, -0.002204742, -0.003488286, 
0.004832194, -0.00137924, -0.016501518, 0.006870124, 0.011860637, 
0.005993393, -0.004560137, -0.007139575, 0.000783146, 0.026054272, 
0.000494855, 0.008204079, -0.011764013, -0.006389326, 0.008017244, 
0.004591451, 0.012983199, 0.010973351, 0.005136924, 0.014092113, 
0.002394753, 0.022453531, 0.012407455, 0.001519177, 0.000155242
)), row.names = c(NA, 49L), class = "data.frame")

问题原因 #

bam()是mgcv中针对大数据优化的模型拟合函数，其fREML（快速REML）的内部矩阵处理逻辑和标准gam()不同。当给by=subject的样条设置m=1时：

1. m=1对应一阶导数惩罚（默认m=2是二阶导数惩罚），会改变样条的惩罚矩阵结构；
2. bam()的fREML在处理这种非默认惩罚结构时，可能生成了维度不匹配的矩阵，导致summary()调用chol()进行Cholesky分解时触发"non-conformable arrays"（矩阵维度不兼容）错误；
3. 额外的s(subject, bs="re")随机截距项可能和by=subject的样条存在共线性，进一步加剧了矩阵维度的冲突——毕竟分组样条已经包含了每个subject的时间趋势，随机截距的冗余可能让模型矩阵变得奇异。

解决方案 #

1. 更换拟合方法：暂时改用method="REML"或method="ML"替代fREML，这两种方法的矩阵处理逻辑更接近gam()，大概率能正常生成summary：

   model2 <- bam(Met1 ~ s(Timepoint, k=7) + s(Timepoint, by=subject, k=7, m=1) + s(subject, bs="re"), data=df, method="REML")
   

2. 调整样条自由度：降低k值（比如设为k=6），减少样条的自由度，避免惩罚矩阵出现维度异常：

   model2 <- bam(Met1 ~ s(Timepoint, k=6) + s(Timepoint, by=subject, k=6, m=1) + s(subject, bs="re"), data=df, method="fREML")
   

3. 移除冗余的随机效应：尝试去掉s(subject, bs="re")，因为by=subject的样条已经能捕捉每个subject的时间趋势，随机截距可能是冗余项：

   model2 <- bam(Met1 ~ s(Timepoint, k=7) + s(Timepoint, by=subject, k=7, m=1), data=df, method="fREML")
   

4. 手动提取模型信息：如果必须保留当前模型结构，可以用gam.vcomp(model2)提取方差成分，用predict(model2)获取拟合值，替代summary()的部分功能。

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