更新R后如何纠正错误收敛mgcv::gamm

Question

当我尝试使用 corARMA 函数为相关误差结构拟合模型时，使用 mgcv::gamm 在 GAMM 上出现错误收敛。响应数据是降水量，并且是 Gamma 分布的。预测变量是年份和年份 ('julian')。我正在将模型拟合到发生降水的观测子集（即，我只对非 0 降水日进行建模）。你会看到数据中有很多 NA，这是因为我只对夏季的降水模式感兴趣。一位同事建议我用 NA 填充其他日期而不是删除它们（有人告诉我，如果我省略 9 月至 4 月的日期，模型可能会尝试以循环模式将 8 月底与 5 月初联系起来）。

相关误差结构旨在描述时间自相关，遵循 Gavin Simpson 的优秀博客 post 中概述的 advice/process。当我使用 R 版本 3.6.1 时，此代码很容易适合模型，但在更新到 4.1.2 后，我收到错误的收敛错误和“系数矩阵不可逆”错误。我试过更改一些控件（主要是迭代次数），但似乎无济于事。我不确定在这种情况下如何调整模型收敛的控制（以及为什么 R 更新会导致此问题），而且我对“系数矩阵不可逆”错误一无所知。非常感谢任何建议或想法。

Link to download the data（行太多直接post）

没有自相关误差结构的建模工作正常：

require("mgcv")
cv.data <- read.csv("cv.data.csv", header = T)
#model precipitation amount for days on which it did rain:
gam.precip.non0<-gamm(PRCP ~ s(julian) + s(year, k = 23), data =  cv.data[cv.data$prcp_binary!=0,], family = Gamma, method = "REML")
   summary(gam.precip.non0$gam)

拟合几个自相关函数进行比较：

    ctrl<-list(niterEM = 0, optimMethod = "L-BFGS-B", maxIter = 100, msMaxIter = 100)
gam.non01<-gamm(PRCP ~ s(julian) + s(year, k = 23), 
                data = cv.data[cv.data$prcp_binary!=0,], 
                family = Gamma, method = "REML", control = ctrl, 
                correlation = corARMA(form = ~1|year, p = 1))
#false convergence

gam.non02<-gamm(PRCP ~ s(julian) + s(year, k = 23), 
                data = cv.data[cv.data$prcp_binary!=0,], 
                family = Gamma, method = "REML", control = ctrl, 
                correlation = corARMA(form = ~1|year, p = 2))
#other error: "coefficient matrix not invertible"

gam.non03<-gamm(PRCP ~ s(julian) + s(year, k = 23), 
                data = cv.data[cv.data$prcp_binary!=0,], 
                family = Gamma, method = "REML", control = ctrl, 
                correlation = corARMA(form = ~1|year, p = 3))
#other error: "coefficient matrix not invertible"

Answer 1

我不确定这是否是您想要的，但我怀疑问题的发生是因为 year 的结构具有重复值和间隙。

如果我们使用连续时间变量，它（从技术上）有效：

require("mgcv")
cv.data <- read.csv("cv.data.csv", header = T)
#model precipitation amount for days on which it did rain:
gam.precip.non0<-gamm(PRCP ~ s(julian) + s(year, k = 23), data =  cv.data[cv.data$prcp_binary!=0,], family = Gamma, method = "REML")
summary(gam.precip.non0$gam)

plot(PRCP ~ X, data=cv.data)

## check if X is time, approximately
#cv.data$time <- cv.data$year * 365 + cv.data$julian
#plot(PRCP ~  time, data=cv.data)

## add msVerbose to show trace
ctrl <- list(niterEM = 0, optimMethod = "L-BFGS-B", maxIter = 100, msMaxIter = 100, msVerbose=TRUE)

## use time or X for autocorelation
gam.non01 <- gamm(PRCP ~ s(julian) + s(year, k = 23), 
                data = cv.data[cv.data$prcp_binary!=0,], 
                family = Gamma, method = "REML", control = ctrl, 
                correlation = corARMA(form = ~1|X, p = 1))

summary(gam.non01$gam)
vis.gam(gam.non01$gam)

我怀疑这只是一半，但仔细重新考虑预期的相关结构可能会有所帮助。

更新R后如何纠正错误收敛mgcv::gamm

How to correct false convergence mgcv::gamm after updating R

r

mgcv

convergence

autocorrelation