nlme出错

Question

对于来自nlme库的IGF数据，我得到了这个错误消息：

lme(conc ~ 1, data=IGF, random=~age|Lot)
Error in lme.formula(conc ~ 1, data = IGF, random = ~age | Lot) : 
  nlminb problem, convergence error code = 1
  message = iteration limit reached without convergence (10)

但是这段代码一切都很好

lme(conc ~ age, data=IGF)
Linear mixed-effects model fit by REML
  Data: IGF 
  Log-restricted-likelihood: -297.1831
  Fixed: conc ~ age 
 (Intercept)          age 
 5.374974367 -0.002535021 

Random effects:
 Formula: ~age | Lot
 Structure: General positive-definite
            StdDev      Corr  
(Intercept) 0.082512196 (Intr)
age         0.008092173 -1    
Residual    0.820627711       

Number of Observations: 237
Number of Groups: 10 

因为IGF是groupedData，所以两个代码是相同的。我搞不懂为什么第一个代码会产生错误。感谢您的时间和帮助。

Answer 1

我发现另一个更老的答案并不令人满意。就我个人而言，我曾将这两个案例混为一谈，从而犯下了职业错误。R已经发出了后者的信号，我想深入探讨为什么会这样。

OP指定的模型是一个增长模型，具有随机斜率和截距。包括一个宏大的截距，但不包括一个宏大的年龄坡度。一个令人不快的约束是，通过拟合随机斜率而不添加其“宏大”术语，您强制随机斜率具有0均值，这是非常难以优化的。边际模型表明，年龄与模型中的0没有统计上显着的差异值。此外，将年龄作为固定效果添加并不能解决这个问题。

> lme(conc~ age, random=~age|Lot, data=IGF)
Error in lme.formula(conc ~ age, random = ~age | Lot, data = IGF) : 
  nlminb problem, convergence error code = 1
  message = iteration limit reached without convergence (10)

这里的错误是显而易见的。设置迭代次数可能很诱人。lmeControl有很多迭代的方法。但即使这样也行不通：

> fit <- lme(conc~ 1, random=~age|Lot, data=IGF, 
control = lmeControl(maxIter = 1e8, msMaxIter = 1e8))

Error in lme.formula(conc ~ 1, random = ~age | Lot, 
data = IGF, control = lmeControl(maxIter = 1e+08,  : 
  nlminb problem, convergence error code = 1
  message = singular convergence (7)

所以这不是精确度的问题，优化器超出了界限。

您提议拟合的两个模型之间必须存在关键差异，并且必须有一种方法来诊断您发现的错误。一种简单的方法是为有问题的模型指定一个“详细”模型：

> lme(conc~ 1, random=~age|Lot, data=IGF, control = lmeControl(msVerbose = TRUE))
  0:     602.96050:  2.63471  4.78706  141.598
  1:     602.85855:  3.09182  4.81754  141.597
  2:     602.85312:  3.12199  4.97587  141.598
  3:     602.83803:  3.23502  4.93514  141.598
   (truncated)
 48:     602.76219:  6.22172  4.81029  4211.89
 49:     602.76217:  6.26814  4.81000  4425.23
 50:     602.76216:  6.31630  4.80997  4638.57
 50:     602.76216:  6.31630  4.80997  4638.57

第一个术语是REML (我认为)。如果您使用Rstudio的调试器来检查这个对象的属性(问题的关键)，您将看到在这里爆炸的是随机效果组件。50次迭代后的coef(lmeSt)生成reStruct.Lot1 reStruct.Lot2 reStruct.Lot3 6.316295 4.809975 4638.570586

如上图所示，并产生

> coef(lmeSt, unconstrained = FALSE)

    reStruct.Lot.var((Intercept)) reStruct.Lot.cov(age,(Intercept)) 
                         306382.7                         2567534.6 
            reStruct.Lot.var(age) 
                       21531399.4 

，它与

Browse[1]> lmeSt$reStruct$Lot
Positive definite matrix structure of class pdLogChol representing
            (Intercept)      age
(Intercept)    306382.7  2567535
age           2567534.6 21531399

所以很明显，随机效应的协方差对于这个特定的优化器来说是爆炸性的。nlminb中的端口例程因其缺乏信息性的错误而受到批评。来自David Gay (贝尔实验室)的文本在这里，http://ms.mcmaster.ca/~bolker/misc/port.pdf端口文档表明我们的错误7来自使用10亿iter最大"x可能有太多的自由组件。见§5。“而不是修复算法，我们理应询问是否有应该产生类似结果的近似结果。例如，很容易将lmList对象拟合为随机截距和随机斜率方差：

> fit <- lmList(conc ~ age | Lot, data=IGF)
> cov(coef(fit))
            (Intercept)          age
(Intercept)  0.13763699 -0.018609973
age         -0.01860997  0.003435819

尽管理想情况下，这些将通过它们各自的精度权重进行加权：

为了使用nlme包，我注意到使用BFGS的无约束优化不会产生这样的错误，并给出类似的结果：

> lme(conc ~ 1, data=IGF, random=~age|Lot, control = lmeControl(opt = 'optim'))
Linear mixed-effects model fit by REML
  Data: IGF 
  Log-restricted-likelihood: -292.9675
  Fixed: conc ~ 1 
(Intercept) 
   5.333577 

Random effects:
 Formula: ~age | Lot
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev      Corr  
(Intercept) 0.032109976 (Intr)
age         0.005647296 -0.698
Residual    0.820819785       

Number of Observations: 237
Number of Groups: 10 

这种模型的另一种语法声明可以使用更简单的lme4包来完成：

library(lme4)
lmer(conc ~ 1 + (age | Lot), data=IGF)

这会产生：

> lmer(conc ~ 1 + (age | Lot), data=IGF)
Linear mixed model fit by REML ['lmerMod']
Formula: conc ~ 1 + (age | Lot)
   Data: IGF
REML criterion at convergence: 585.7987
Random effects:
 Groups   Name        Std.Dev. Corr 
 Lot      (Intercept) 0.056254      
          age         0.006687 -1.00
 Residual             0.820609      
Number of obs: 237, groups:  Lot, 10
Fixed Effects:
(Intercept)  
      5.331 

lmer及其优化器的一个属性是，非常接近1、0或-1的随机效果相关性被简单地设置为这些值，因为它大大简化了优化(以及估计的统计效率)。

综上所述，这并不意味着年龄没有影响，正如前面所说的那样，这一论点可以得到数字结果的支持。

Answer 2

如果绘制数据图，可以看到没有age的影响，所以尽管如此，尝试拟合age的随机效果似乎很奇怪。难怪它不会收敛。

library(nlme)
library(ggplot2)

dev.new(width=6, height=3)
qplot(age, conc, data=IGF) + facet_wrap(~Lot, nrow=2) + geom_smooth(method='lm')

​

​

我认为你想要做的是模拟Lot对截获的随机影响。我们可以尝试将age作为固定效果包含在内，但我们将看到它并不重要，可以将其抛出：

> summary(lme(conc ~ 1 + age, data=IGF, random=~1|Lot))
Linear mixed-effects model fit by REML
 Data: IGF 
       AIC      BIC    logLik
  604.8711 618.7094 -298.4355

Random effects:
 Formula: ~1 | Lot
        (Intercept) Residual
StdDev:  0.07153912 0.829998

Fixed effects: conc ~ 1 + age 
                Value  Std.Error  DF  t-value p-value
(Intercept)  5.354435 0.10619982 226 50.41849  0.0000
age         -0.000817 0.00396984 226 -0.20587  0.8371
 Correlation: 
    (Intr)
age -0.828

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-5.46774548 -0.43073893 -0.01519143  0.30336310  5.28952876 

Number of Observations: 237
Number of Groups: 10