Ошибка с 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

Если вы нанесете данные на график, вы увидите, что эффект 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 

Answer 2

Я нахожу другой, более старый ответ здесь неудовлетворительным.Я различаю случаи, когда, по статистике, возраст не оказывает влияния, и наоборот, мы сталкиваемся с вычислительной ошибкой.Лично я допустил ошибки в карьере, соединив эти два случая.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 были подвергнуты критике за их неинформативные ошибки.Текст от Дэвида Гея (Bell Labs) находится здесь http://ms.mcmaster.ca/~bolker/misc/port.pdf Документация PORT предлагает нашу ошибку 7 от использования максимального значения в 1 миллиард итеров "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 просто устанавливаются на эти значения, поскольку это существенно упрощает оптимизацию (и статистическую эффективность оценки).

В совокупности это не предполагает, что возраст не имеетэффект, как было сказано ранее, и этот аргумент может быть подтвержден числовыми результатами.