Ostrzeżenie „Model nie zbiegło się” w lmer ()

Za pomocą następującego zestawu danych chciałem sprawdzić, czy odpowiedź (efekt) zmienia się w odniesieniu do witryn, sezonu, czasu trwania i ich interakcji. Niektóre internetowe fora poświęcone statystykom sugerowały, żebym kontynuował liniowe modele z efektami mieszanymi, ale problem polega na tym, że ponieważ replikacje są losowe w każdej stacji, nie mam szans na pobranie próbki z dokładnie tego samego miejsca w kolejnych sezonach (na przykład repl-1 z s1 po monsunu może nie być taki sam jak monsun). W przeciwieństwie do badań klinicznych (z projektem wewnątrz przedmiotu), w których mierzysz ten sam przedmiot wielokrotnie w różnych porach roku. Jednak biorąc pod uwagę losowość witryn i sezonu, uruchomiłem następujące polecenia i otrzymałem komunikat ostrzegawczy:

Warning messages:
1: In checkConv(attr(opt, "derivs"), optpar,ctrl=controlpar,ctrl=controlcheckConv, 
: unable to evaluate scaled gradient
2: In checkConv(attr(opt, "derivs"), optpar,ctrl=controlpar,ctrl=controlcheckConv, 
: Model failed to converge: degenerate Hessian with 1 negative eigenvalues

Czy ktoś może mi pomóc rozwiązać problem? Kody podano poniżej:

library(lme4)
read.table(textConnection("duration season  sites   effect
                          4d    mon s1  7305.91
                          4d    mon s2  856.297
                          4d    mon s3  649.93
                          4d    mon s1  10121.62
                          4d    mon s2  5137.85
                          4d    mon s3  3059.89
                          4d    mon s1  5384.3
                          4d    mon s2  5014.66
                          4d    mon s3  3378.15
                          4d    post    s1  6475.53
                          4d    post    s2  2923.15
                          4d    post    s3  554.05
                          4d    post    s1  7590.8
                          4d    post    s2  3888.01
                          4d    post    s3  600.07
                          4d    post    s1  6717.63
                          4d    post    s2  1542.93
                          4d    post    s3  1001.4
                          4d    pre s1  9290.84
                          4d    pre s2  2199.05
                          4d    pre s3  1149.99
                          4d    pre s1  5864.29
                          4d    pre s2  4847.92
                          4d    pre s3  4172.71
                          4d    pre s1  8419.88
                          4d    pre s2  685.18
                          4d    pre s3  4133.15
                          7d    mon s1  11129.86
                          7d    mon s2  1492.36
                          7d    mon s3  1375
                          7d    mon s1  10927.16
                          7d    mon s2  8131.14
                          7d    mon s3  9610.08
                          7d    mon s1  13732.55
                          7d    mon s2  13314.01
                          7d    mon s3  4075.65
                          7d    post    s1  11770.79
                          7d    post    s2  4254.88
                          7d    post    s3  753.2
                          7d    post    s1  11324.95
                          7d    post    s2  5133.76
                          7d    post    s3  2156.2
                          7d    post    s1  12103.76
                          7d    post    s2  3143.72
                          7d    post    s3  2603.23
                          7d    pre s1  13928.88
                          7d    pre s2  3208.28
                          7d    pre s3  8015.04
                          7d    pre s1  11851.47
                          7d    pre s2  6815.31
                          7d    pre s3  8478.77
                          7d    pre s1  13600.48
                          7d    pre s2  1219.46
                          7d    pre s3  6987.5
                          "),header=T)->dat1


m1 = lmer(effect ~ duration + (1+duration|sites) +(1+duration|season),
          data=dat1, REML=FALSE)

„Rozwiązanie” napotkanego problemu polegającego na nieotrzymywaniu ostrzeżeń o nieudanej konwergencji jest raczej proste: nie używasz domyślnego optymalizatora BOBYQA, ale zamiast tego decydujesz się na procedurę optymalizacji Neldera-Meada używaną domyślnie we wcześniejszych 1.0.xwersjach. Lub instalujesz pakiet optimx, abyś mógł bezpośrednio wykonać procedurę L-BFGS-B lub nlminb(taką samą jak w lme4wersjach wcześniejszych niż wer. 1). Na przykład:

m1 = lmer(effect~duration+(1+duration|sites)+(1+duration|season), 
          data = dat1, REML = FALSE, 
          control = lmerControl(optimizer ="Nelder_Mead")
library(optimx)
m1 = lmer(effect~duration+(1+duration|sites)+(1+duration|season), 
          data = dat1, REML = FALSE, 
          control = lmerControl(
                           optimizer ='optimx', optCtrl=list(method='L-BFGS-B')))
m1 = lmer(effect~duration+(1+duration|sites)+(1+duration|season), 
          data = dat1, REML = FALSE, 
          control = lmerControl(
                           optimizer ='optimx', optCtrl=list(method='nlminb')))

wszystko działa dobrze (bez ostrzeżeń). Interesujące pytania to:

1. dlaczego masz te ostrzeżenia na początek i
2. dlaczego przy użyciu REML = TRUEnie masz ostrzeżeń.

W skrócie: 1. otrzymałeś te ostrzeżenia, ponieważ zdefiniowałeś durationzarówno jako stały efekt, jak i losowe nachylenie dla tego czynnika, sitesa także season. W modelu skutecznie zabrakło stopni swobody w celu oszacowania korelacji między nachyleniami a zdefiniowanymi przez siebie punktami przecięcia. Jeśli użyłeś nieco prostszego modelu, takiego jak:

m1 = lmer(effect~duration+ (1+duration|sites) + (0+duration|season) + (1|season),
          data=dat1, REML = FALSE)

nie wystąpiłyby problemy z konwergencją. Ten model skutecznie oszacowałby nieskorelowane losowe przechwyty i losowe nachylenia dla każdego z nich season.

REML = FALSE$X$$y = X\beta + Z\gamma + \epsilon$$K$$KX = 0$$y$$Ky$$Z$$KZ$$Z$

Ostatnia uwaga jest taka, że ​​nie jestem pewien, czy seasonna początek warto zastosować jako efekt losowy. W końcu jest tylko tyle pór roku, więc równie dobrze możesz traktować je jako ustalone efekty.

Pytanie jest raczej statystyczne niż techniczne. Właściwie użyłem modelu efektu losowego zamiast modelu efektu stałego. Myślę, że żaden z czynników nie powinien być traktowany jako czynnik losowy, ponieważ potrzebujemy co najmniej 5 lub 6 poziomów lub powtórzeń, aby traktować czynnik jako efekt losowy (patrz tutaj Jaka jest minimalna zalecana liczba grup dla współczynnika efektów losowych? ).

Powyższy zestaw danych zawiera tylko trzykrotne próbki / miejsce / sezon, co jest niewystarczające dla modelu efektu losowego. W zestawie danych czas trwania, 4-dniowy i 7-dniowy, należą do dwóch oddzielnych równoległych eksperymentów przeprowadzanych w tym samym czasie. Więc plując zestawem danych według czasu trwania (4-dniowy i 7-dniowy) i wykonując 2-kierunkową anovę dla każdego czasu trwania z sezonem i miejscami, ponieważ czynniki byłyby wystarczające do modelowania efektu (zmienna odpowiedzi) tutaj. Model powinien wyglądać następująco:

lm(day_4_effect~sites*season, data=dat1)

lm(day_7_effect~sites*season, data=dat1)

Podziękowania dla Bodo Winter ( http://www.bodowinter.com/tutorial/bw_LME_tutorial2.pdf ) i @ usεr11852, którzy pomogli mi rozwiązać problem.