Błąd nlme

Question

Dla IGF danych nlme bibliotece, Dostaję komunikat o błędzie:

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) 

Ale wszystko jest w porządku z tym kodem

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 

Jak IGF jest groupedData, więc oba kody są identyczne. Jestem zdezorientowany, dlaczego pierwszy kod generuje błąd. Dziękuję za poświęcony czas i pomoc.

Answer 1

Jeśli narysujesz dane, możesz zauważyć, że nie ma to wpływu na age, więc wydaje się dziwne, że próbujesz dopasować losowy efekt age pomimo tego. Nic dziwnego, że się nie zbliża.

library(nlme) 
library(ggplot2) 

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

myślę, co chcesz zrobić, to wzór losowy efekt Lot na przecięcia. Możemy spróbować tym age jako stały efekt, ale zobaczymy, że nie jest istotne i może być wyrzucony:

> 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

znajdę drugiego, starszy odpowiedź tutaj niezadowalające. Rozróżniam przypadki, w których statystycznie wiek nie ma żadnego wpływu i na odwrót napotykamy błąd obliczeniowy. Osobiście popełniłem błędy zawodowe, łącząc te dwa przypadki. R zasygnalizował to ostatnie i chciałbym się zastanowić, dlaczego tak jest.

Model określony przez OP to model wzrostu z losowymi nachyleniami i przechwytami. Uwzględniono duży punkt przecięcia, ale nie jest to stok wielkiej epoki. Jednym z niesmacznych ograniczeń nałożonych przez dopasowanie losowego nachylenia bez dodawania jego "wielkiego" terminu jest wymuszenie, że nachylenie losowe ma 0, co jest bardzo trudne do zoptymalizowania. Modele marginalne wskazują, że wiek nie ma statystycznie istotnej różnej wartości od 0 w modelu. Ponadto dodanie wieku jako efektu stałego nie rozwiązuje problemu.

> 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) 

Tutaj błąd jest oczywisty. Może być kuszące ustawienie liczby iteracji. lmeControl ma wiele iteracyjnych oszacowań. Ale nawet to nie działa:

> 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) 

Nie jest to rzeczą precyzyjną, optymalizator działa poza wyznaczonymi granicami.

Pomiędzy dwoma modelami, które zaproponowałeś, muszą występować kluczowe różnice, a także sposób na zdiagnozowanie znalezionego błędu. Jednym prostym podejściem jest określenie "gęsiego" dopasowania do problematycznego modelu:

> 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 

Pierwszy termin to REML (myślę). Druga do czwarta kategoria to parametry obiektu o nazwie lmeSt klasy lmeStructInt, lmeStruct i modelStruct. Jeśli użyjesz debuggera Rstudio do sprawdzenia atrybutów tego obiektu (zlinczowania problemu), zobaczysz, że tutaj wybucha składnik losowych efektów.coef(lmeSt) po 50 powtórzeń produkuje reStruct.Lot1 reStruct.Lot2 reStruct.Lot3 6.316295 4.809975 4638.570586

jak widać powyżej i produkuje

> coef(lmeSt, unconstrained = FALSE) 

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

która jest taka sama jak

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

Więc jasne jest kowariancja efektów losowych jest coś, co wybucha tutaj ten konkretny optymalizator. Procedury PORT w nlminb zostały skrytykowane za ich nieinformacyjne błędy. Tekst z David Gay (Bell Labs) znajduje się tutaj http://ms.mcmaster.ca/~bolker/misc/port.pdf Dokumentacja PORT sugeruje, że nasz błąd 7 wynika z użycia 1 miliarda iter max "x może mieć zbyt wiele wolnych elementów (patrz §5."). Zamiast naprawiać algorytm, wypada zapytać, czy istnieją przybliżone wyniki, które powinny generować podobne wyniki. Jest, na przykład, łatwo dopasować obiekt lmList wymyślić losowego przechwycić i losowych zbocza wariancji:

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

choć idealnie byłyby one ważone przez swoich wagach precyzyjnych:

Aby użyć nlme pakiet Pragnę zauważyć, że swobodne optymalizacji przy użyciu BFGS nie wywołuje takiego błędu i daje podobne wyniki:

> 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.9976 (Intr) 
age   0.005647296 -0.698 
Residual 0.820819785  

Number of Observations: 237 
Number of Groups: 10 

alternatywą składniowej deklaracja takiego modelu można zrobić z t on znacznie łatwiejsze lme4 pakiet:

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

co daje:

> 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 

atrybutem lmer i jego Optimizer to, że efektów losowych korelacje, które są bardzo zbliżone do 1, 0 lub -1 po prostu ustawić do tych wartości, ponieważ znacznie upraszcza optymalizację (i statystyczną efektywność estymacji).

Podsumowując, to nie , a nie sugeruje, że wiek nie ma wpływu, jak wspomniano wcześniej, i ten argument może być obsługiwany przez wyniki liczbowe.