含随机截距、连续AR1相关结构与异方差的线性混合模型（LMM）收敛问题排查求助

大家好，我在用nlme拟合重复测量模型时遇到了收敛问题，想请各位帮忙梳理下排查思路，先说说我的情况：

数据概况 #

• 我的数据是长格式，涵盖104名患者的304次测量记录：其中68名患者只有2次测量，32名有5次，3名4次，1名3次，1名6次。
• 测量间隔差异很大：有74次测量的间隔约0.5年，而部分仅做2次测量的患者间隔长达5年，整体平均随访时长大概3.5年。

已成功拟合的模型 #

我已经顺利拟合了随机截距+斜率模型，代码如下：

library(nlme)
ris <- lme(fixed=outcome~ 1 + fu_time + age*fu_time + sex*fu_time + smoking*fu_time + obesity*fu_time + diab*fu_time + hypt*fu_time + hyperchol*fu_time + ckd*fu_time, random=~1 + fu_time|patid, data=data, na.action="na.omit", method="ML")

尝试的模型及收敛问题 #

为了考虑测量间隔的差异，同时保留方差随时间递增的特性，我尝试构建一个包含以下要素的模型：

1. 随机截距
2. 连续AR1相关结构
3. 异方差（随fu_time变化）

对应的代码如下，但这个模型始终无法收敛：

ri_cAR1_hetvar <- update(ris, random=~1|patid, correlation=corCAR1(form=~1|patid), weight=varIdent(form=~1|fu_time), control=lmeControl(opt="optim")) # does not converge

调试输出 #

按照Ben Bolker之前的建议，我运行debug(nlme:::logLik.reStruct)后执行相关命令，得到了以下输出：

运行str(object)的结果 #

List of 1
 $ patid: 'pdLogChol' num 0.0752
 ..- attr(*, "formula")=Class 'formula' language ~1
 .. .. ..- attr(*, ".Environment")=<environment: 0x000000001ec31628>
 ..- attr(*, "Dimnames")=List of 2
 .. ..$ : chr "(Intercept)"
 .. ..$ : chr "(Intercept)"
 - attr(*, "settings")= int [1:4] 0 1 0 4
 - attr(*, "class")= chr "reStruct"
 - attr(*, "plen")= Named int 1
 ..- attr(*, "names")= chr "patid"

运行str(conLin)的结果 #

List of 5
 $ Xy : num [1:314, 1:22] 1 0.816 0.816 0.816 0.816 ...
 ..- attr(*, "dimnames")=List of 2
 .. ..$ : chr [1:314] "1" "2" "3" "4" ...
 .. ..$ : chr [1:22] "(Intercept)" "(Intercept)" "fu_time" "age" ...
 $ dims :List of 15
 ..$ N : int 314
 ..$ ZXrows : int 314
 ..$ ZXcols : num 22
 ..$ Q : int 1
 ..$ StrRows: num 125
 ..$ qvec : Named num [1:3] 1 0 0
 .. ..- attr(*, "names")= chr [1:3] "patid" "" ""
 ..$ ngrps : Named int [1:3] 104 1 1
 .. ..- attr(*, "names")= chr [1:3] "patid" "X" "y"
 ..$ DmOff : Named num [1:3] 0 1 401
 .. ..- attr(*, "names")= chr [1:3] "" "patid" ""
 ..$ ncol : Named num [1:3] 1 20 1
 .. ..- attr(*, "names")= chr [1:3] "patid" "" ""
 ..$ nrot : Named num [1:3] 21 1 0
 .. ..- attr(*, "names")= chr [1:3] "" "" ""
 ..$ ZXoff :List of 3
 .. ..$ patid: num [1:104] 0 6 8 13 15 18 23 28 30 32 ...
 .. ..$ X : Named num 314
 .. .. ..- attr(*, "names")= chr "patid"
 .. ..$ y : Named num 6594
 .. .. ..- attr(*, "names")= chr ""
 ..$ ZXlen :List of 3
 .. ..$ patid: num [1:104] 6 2 5 2 3 5 5 2 2 2 ...
 .. ..$ X : num 314
 .. ..$ y : num 314
 ..$ SToff :List of 3
 .. ..$ patid: num [1:104] 0 1 2 3 4 5 6 7 8 9 ...
 .. ..$ X : Named num 229
 .. .. ..- attr(*, "names")= chr "patid"
 .. ..$ y : Named num 2749
 .. .. ..- attr(*, "names")= chr ""
 ..$ DecOff :List of 3
 .. ..$ PXE_nr: num [1:104] 0 1 2 3 4 5 6 7 8 9 ...
 .. ..$ X : Named num 125
 .. .. ..- attr(*, "names")= chr "patid"
 .. ..$ y : Named num 2625
 .. .. ..- attr(*, "names")= chr ""
 ..$ DecLen :List of 3
 .. ..$ PXE_nr: num [1:104] 1 1 1 1 1 1 1 1 1 1 ...
 .. ..$ X : num 125
 .. ..$ y : num 125
 $ logLik :Class 'logLik' : 4.3 (df=178)
 $ sigma : num 0
 $ auxSigma: num 0

运行dput(object)的结果 #

structure(list(patid= structure(-1.29387589179439, formula = ~1, Dimnames = list( "(Intercept)", "(Intercept)"), class = c("pdLogChol", "pdSymm", "pdMat")), settings = c(0L, 1L, 0L, 4L), class = "reStruct", plen = c(patid= 1L))

最新发现 #

更新：我排查后发现，单独添加连续AR1相关结构的模型是可以正常收敛的，问题出在加上weight=varIdent(form=~1|fu_time)指定随fu_time变化的异方差之后。我现在有点怀疑是参数估计量的问题：原本以为异方差只需要估计phi参数，但突然想到我的fu_time是连续变量，有大量不同的取值水平，是不是意味着模型要为每个独特的时间点估计一个参数？这样参数数量会不会太多，导致模型无法收敛？

我知道没有可复现的示例会让大家排查起来更麻烦，但实在没办法提供更多数据，而且这个模型在我提取的小数据子集上是能正常收敛的。真心恳请各位给我一些排查这个收敛问题的思路和建议！

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