时空计数数据GAM建模四大技术问题咨询

1. 含时空张量积交互项的GAM拟合失败/不收敛的解决方法 #

• 简化张量积复杂度：先降低时空维度的平滑度参数k，比如设置te(time, x, y, k=c(4,3,3))（时间维度k=4，空间x/y各k=3），收敛后再逐步调高k值。也可以用ti()（张量积交互项，不含主效应）替代te()，先拟合单独的时间/空间主效应，再加入交互项，降低模型复杂度。
• 调整估计方法与迭代参数：将默认的GCV估计换成REML（method="REML"），REML对复杂模型的稳定性更好；同时增加迭代次数，比如control=gam.control(maxit=1000, trace=TRUE)，trace参数可以查看迭代过程，定位收敛瓶颈。
• 用简单模型初始化：先拟合不含交互项的GAM（比如仅含时间和空间单独平滑项），将其拟合结果作为初始值传入复杂模型，比如gam(..., start=fitted(simple_model))。
• 检查数据质量：排查是否存在极端离群值、时空维度的缺失/异常分组（比如某时间点仅少数网格有数据），这类数据会干扰拟合稳定性。

2. 负二项GAM的Pearson离散参数0.85的解读 #

• 负二项模型的离散参数接近1时，本质上趋近于泊松模型。0.85属于轻微欠离散，但负二项模型对这种程度的偏离鲁棒性较强，若模型的偏差解释率（Deviance explained）、预测效果表现良好，可认为模型合适。
• 若担心欠离散影响推断，可改用准泊松模型（quasi-Poisson），它允许离散参数自由调整，无需假设特定的分布形式，能适配过离散或欠离散场景。
• 同时需检查离散参数的计算方式：确保是用Pearson残差平方和除以残差自由度（n-p，p为模型参数个数），避免因计算方式错误导致误判。

3. ziP零膨胀模型提取Pearson残差/检验过离散的方法 #

手动计算Pearson残差 #

如果使用mgcv拟合ziP模型（family=ziP()），可按以下步骤手动计算：

# 提取模型拟合值（零膨胀后的总体均值）
mu_total <- fitted(zip_model)
# 提取计数部分的均值（泊松部分的μ）
mu_count <- exp(predict(zip_model, type = "link")[, 1])
# 提取零膨胀概率p
p_zero <- plogis(predict(zip_model, type = "link")[, 2])
# 计算ziP的方差：Var(Y) = μ_total + (p_zero/(1-p_zero)) * mu_count^2
var_y <- mu_total + (p_zero / (1 - p_zero)) * mu_count^2
# 计算Pearson残差
pearson_resid <- (y_obs - mu_total) / sqrt(var_y)
# 计算离散参数：残差平方和/残差自由度
dispersion <- sum(pearson_resid^2) / (length(y_obs) - length(coef(zip_model)))

用模拟残差检验（更可靠） #

推荐使用DHARMa包生成标准化模拟残差，适配几乎所有模型类型：

library(DHARMa)
# 生成模拟残差
sim_resids <- simulateResiduals(model = zip_model, n = 1000)
# 检验离散性
testDispersion(sim_resids)
# 可视化残差
plot(sim_resids)

4. 网格质心vs网格ID随机效应，GAM中体现空间相似性的方法 #

• 两种方式的差异：
  • 质心坐标平滑项（s(x,y)）：适合捕捉连续的空间趋势，但1000km低分辨率网格下，空间趋势可能不明显，平滑效果有限。
  • 网格ID普通随机效应（s(grid_id, bs="re")）：仅能捕捉网格间的异质性，但无法体现相邻网格的空间相似性（随机效应独立）。
• 体现空间相似性的最优方案：
  1. 马尔可夫随机场（MRF）：构建网格的邻接矩阵，用bs="mrf"指定空间结构化随机效应，强制相邻网格的效应相关：

     # 先构建邻接矩阵adj_matrix（需根据网格拓扑生成）
     gam_model <- gam(y ~ s(grid_id, bs="mrf", xt=list(adj=adj_matrix)) + s(time, bs="cr"), 
                      family=ziP(), data=df, method="REML")
     

  2. 高斯过程（GP）：基于质心坐标的空间距离建模相关性，用bs="gp"：

     gam_model <- gam(y ~ s(x, y, bs="gp", k=50) + s(time, bs="cr"), 
                      family=ziP(), data=df, method="REML")
     

  3. 时空结合的话，可将结构化空间效应与时间平滑项做张量积，比如te(x,y,time, bs=c("gp", "cr"))，同时捕捉时空交互与空间相似性。

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