Errore con NLME

Question

Per IGF dati da nlme biblioteca, sto ottenendo questo messaggio di errore:

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) 

ma tutto è bene con questo codice

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 

Come IGF è groupedData, in modo da entrambi i codici sono identici. Sono confuso perché il primo codice produce errori. Grazie per il tuo tempo e aiuto.

Answer 1

Se si tracciano i dati, è possibile vedere che non vi è alcun effetto di age, quindi sembra strano provare a montare un effetto casuale di age nonostante ciò. Non c'è da stupirsi che non converga.

library(nlme) 
library(ggplot2) 

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

Penso che ciò che si vuole fare è il modello un effetto casuale di Lot sul intercetta. Possiamo provare compreso age come effetto fisso, ma vedremo che non è significativo e può essere buttato fuori:

> 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

trovo l'altro, la risposta più vecchio qui insoddisfacente. Distinguo tra casi in cui, statisticamente, l'età non ha alcun impatto e al contrario incontriamo un errore computazionale. Personalmente, ho commesso errori di carriera confondendo questi due casi. R ha segnalato quest'ultimo e mi piacerebbe immergermi nel perché è così.

Il modello che OP ha specificato è un modello di crescita, con pendenze e intercettazioni casuali. Una grande intercettazione è inclusa ma non una pendenza grand age. Un vincolo sgradevole che viene imposto adattando una pendenza casuale senza l'aggiunta del suo termine "grande" è che si sta costringendo la pendenza casuale ad avere la media 0, che è molto difficile da ottimizzare. I modelli marginali indicano che l'età non ha un valore diverso statisticamente significativo da 0 nel modello. Inoltre, l'aggiunta dell'età come effetto fisso non risolve il problema.

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

Qui l'errore è ovvio. Potrebbe essere allettante impostare il numero di iterazioni. lmeControl ha molti stimatori iterativi. Ma anche questo non funziona:

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

Quindi non è una cosa di precisione, l'ottimizzatore sta finendo i limiti.

Ci devono essere differenze chiave tra i due modelli che hai proposto di montare e un modo per diagnosticare l'errore che hai trovato. Un approccio semplice è specificare una misura "verbose" per il modello problematico:

> 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 

Il primo termine è la REML (credo). I termini dal secondo al quarto sono i parametri di un oggetto chiamato lmeSt della classe lmeStructInt, lmeStruct e modelStruct. Se usi il debugger di Rstudio per ispezionare gli attributi di questo oggetto (il lynchpin del problema), vedrai che è il componente degli effetti casuali che esplode qui.coef(lmeSt) dopo 50 iterazioni produce reStruct.Lot1 reStruct.Lot2 reStruct.Lot3 6.316295 4.809975 4638.570586

come visto sopra e produce

> coef(lmeSt, unconstrained = FALSE) 

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

che è lo stesso del

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

Quindi è chiaro che la covarianza degli effetti casuali, è qualcosa che è esplode qui per questo particolare ottimizzatore. Le routine PORT in nlminb sono state criticate per i loro errori non informativi. Il testo di David Gay (Bell Labs) è qui http://ms.mcmaster.ca/~bolker/misc/port.pdf La documentazione PORT suggerisce che il nostro errore 7 di utilizzare un miliardo di iter max "x potrebbe avere troppi componenti liberi. Piuttosto che correggere l'algoritmo, è necessario chiedere se ci sono risultati approssimativi che dovrebbero generare risultati simili. È, per esempio, facile da montare un oggetto lmList a venire con l'intercetta casuale e varianza pendenza casuale:

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

sebbene idealmente questi sarebbero ponderati dai rispettivi pesi di precisione:

Per utilizzare la nlme pacchetto Prendo atto che l'ottimizzazione non vincolata usando BFGS non produce un errore del genere e dà risultati simili:

> 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 

una dichiarazione sintattica alternativa di un tale modello può essere fatto con t ha molto più facile pacchetto lme4:

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

che produce:

> 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 

Un attributo di lmer e la sua ottimizzatore è che le correlazioni effetti casuali che sono molto vicino a 1, 0, -1 o vengono semplicemente inseriti a quei valori poiché semplifica notevolmente l'ottimizzazione (e l'efficienza statistica della stima).

Prendendo insieme, questo non suggerisce che l'età non ha un effetto, come è stato detto in precedenza, e questo argomento può essere supportato dai risultati numerici.