重复测量分析问题：lme模型出现NA/NaN/Inf错误

问题与报错 #

使用nlme包的lme()函数拟合线性混合效应模型分析重复测量数据时，完整数据集上模型运行正常，但按性别过滤数据集后出现如下错误：

mod_lin <- lme(variable ~ seg * group + fumador + diabetes + hipercol + hta + educ,
               random = ~1 | id2,
               correlation = corCAR1(form = ~seg | id2),
               control = lmeControl(opt = "optim"),
               data = xxx,
               method = 'REML',
               na.action = na.exclude)

Error in logLik.reStruct(object, conLin): 
  NA/NaN/Inf in foreign function call (arg 3)

排查步骤 #

• 确认过滤后的数据集无缺失值
• 检查预测变量间共线性：通过构造模型矩阵并使用caret::findLinearCombos()检测，未发现线性组合

  X <- model.matrix(~ seg * group + fumador + diabetes + hipercol + hta + educ, data = xxx)
  caret::findLinearCombos(X)
  

数据集结构与统计信息 #

过滤后数据集包含51683行记录，结构与统计信息如下：

str(xxx)
'data.frame':   51683 obs. of  10 variables:
 $ id      : chr  "01_110097" "01_110097" "01_110097" "01_110097" ...
 $ variable: num  24.5 25.1 25.8 27.2 23.9 ...  # 因变量（BMI）
 $ group   : Factor w/ 4 levels "1","2","3","4": 2 2 2 2 3 3 3 3 1 1 ...
 $ seg     : num  47.2 49.2 51.2 52.5 55.6 ...  # 年龄（连续变量）
 $ fumador : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 2 2 ...
 $ diabetes: Factor w/ 2 levels "1","2": 2 2 2 2 2 2 2 2 2 2 ...
 $ hipercol: Factor w/ 2 levels "1","2": 2 2 2 2 2 2 2 2 2 2 ...
 $ hta     : Factor w/ 2 levels "1","2": 2 2 2 2 2 2 2 2 2 2 ...
 $ educ    : Factor w/ 3 levels "0","1","3": 3 3 3 3 1 1 1 1 1 1 ...
 $ id2     : Factor w/ 6414 levels "01_110097","01_112224",..: 1 1 1 1 2 2 2 2 3 3 ...

# 各分类变量频数统计
sapply(xxx[, c("variable", "seg", "group", "fumador", "diabetes", "hipercol", "hta")], function(x) table(x))

$group
    1     2     3     4 
11599 12184 12294 15606 

$fumador
    1     2     3 
14224 14688 22771 

$diabetes
    1     2 
10645 41038 

$hipercol
    1     2 
20178 31505 

$hta
    1     2 
22447 29236 

# 部分参与者的重复测量记录数示例
 codigo_unique  n
     01_110097  4
     01_112224  4
     01_112923  4
     01_113024 13
     01_114243  3
     01_114244 14

解决方法 #

通过修改control参数，更换优化器为nlminb并调整迭代参数，解决了模型收敛问题：

mod_lin <- lme(variable ~ seg * group + fumador + diabetes + hipercol + hta + educ, 
               random = ~1 | id2, 
               correlation = corCAR1(form = ~seg | id2), 
               control = lmeControl(opt = "nlminb", 
                                    msMaxIter = 200, 
                                    msVerbose = TRUE,
                                    msTol = 1e-5), 
               data = xxx,
               method = 'REML',
               na.action = na.exclude)

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