nlme出错

Question

对于IGF来自nlme库的数据,我收到此错误消息:

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（我认为）。至第四方面的第二是参数称为对象lmeSt类的lmeStructInt，lmeStruct和modelStruct。如果您使用 Rstudio 的调试器来检查此对象的属性（问题的关键），您会看到它是随机效应组件在这里爆炸。coef(lmeSt)经过 50 次迭代产生 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

所以很明显，随机效应的协方差对于这个特定的优化器来说是爆炸性的。中的 PORT 例程nlminb因其无信息错误而受到批评。David Gay（贝尔实验室）的文本在这里http://ms.mcmaster.ca/~bolker/misc/port.pdf PORT 文档建议我们使用 10 亿迭代最大的错误 7 “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