Co zrobić z korelacją efektów losowych równą 1 lub -1?

Nierzadko zdarza się, gdy mamy do czynienia ze złożonymi maksymalnymi modelami mieszanymi (szacowanie wszystkich możliwych efektów losowych dla danych i modelu) jest idealna (+1 lub -1) lub prawie idealna korelacja między niektórymi efektami losowymi. Na potrzeby dyskusji przyjrzyjmy się poniższemu modelowi i podsumowaniu modelu

Model: Y ~ X*Cond + (X*Cond|subj)

# Y = logit variable  
# X = continuous variable  
# Condition = values A and B, dummy coded; the design is repeated 
#             so all participants go through both Conditions  
# subject = random effects for different subjects  

Random effects:
 Groups  Name             Variance Std.Dev. Corr             
 subject (Intercept)      0.85052  0.9222                    
         X                0.08427  0.2903   -1.00            
         CondB            0.54367  0.7373   -0.37  0.37      
         X:CondB          0.14812  0.3849    0.26 -0.26 -0.56
Number of obs: 39401, groups:  subject, 219

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)       2.49686    0.06909   36.14  < 2e-16 ***
X                -1.03854    0.03812  -27.24  < 2e-16 ***
CondB            -0.19707    0.06382   -3.09  0.00202 ** 
X:CondB           0.22809    0.05356    4.26 2.06e-05 ***

Przypuszczalną przyczyną tych idealnych korelacji jest to, że stworzyliśmy model, który jest zbyt złożony dla danych, które mamy. Powszechną radą, która jest podawana w takich sytuacjach jest (np. Matuschek i in., 2017; paper ), aby naprawić nadparametryzowane współczynniki na 0, ponieważ takie zdegenerowane modele mają tendencję do obniżania mocy. Jeśli zauważymy wyraźną zmianę w ustalonych efektach w modelu zredukowanym, powinniśmy to zaakceptować; jeśli nie ma zmian, nie ma problemu z zaakceptowaniem oryginalnego.

Załóżmy jednak, że nie interesują nas nie tylko ustalone efekty kontrolowane dla RE (efekty losowe), ale także struktura RE. W danym przypadku teoretycznie rozsądnie byłoby założyć, że Intercepti nachylenie Xma niezerową korelację ujemną. Oto kilka pytań:

1. Co robić w takich sytuacjach? Czy powinniśmy zgłosić idealną korelację i powiedzieć, że nasze dane nie są „wystarczająco dobre”, aby oszacować „prawdziwą” korelację? Czy też powinniśmy zgłosić model korelacji 0? A może powinniśmy spróbować ustawić inną korelację na 0 w nadziei, że ta „ważna” nie będzie już idealna? Nie wydaje mi się, żeby tu były w 100% poprawne odpowiedzi, w większości chciałbym usłyszeć wasze opinie.

2. Jak napisać kod, który naprawia korelację 2 konkretnych losowych efektów do 0, bez wpływu na korelacje między innymi parametrami?

Pojedyncze macierze kowariancji efektu losowego

Uzyskanie oszacowania korelacji efektu losowego +1 lub -1 oznacza, że ​​algorytm optymalizacyjny uderzył „granicę”: korelacje nie mogą być wyższe niż +1 ani niższe niż -1. Nawet jeśli nie ma wyraźnych błędów lub ostrzeżeń o konwergencji, może to potencjalnie wskazywać na pewne problemy z konwergencją, ponieważ nie oczekujemy, że prawdziwe zależności leżą na granicy. Jak powiedziałeś, zwykle oznacza to, że nie ma wystarczających danych, aby wiarygodnie oszacować wszystkie parametry. Matuschek i in. 2017 twierdzą, że w tej sytuacji moc może zostać zagrożona.

Innym sposobem na przekroczenie granicy jest oszacowanie wariancji na 0: Dlaczego otrzymuję zerową wariancję losowego efektu w moim modelu mieszanym, pomimo pewnych różnic w danych?

Obie sytuacje mogą być postrzegane jako uzyskiwanie zdegenerowanej macierzy kowariancji efektów losowych (w twoim przykładzie wyjściowa macierz kowariancji wynosi ); zerowa wariancja lub idealna korelacja oznacza, że ​​macierz kowariancji nie ma pełnej rangi, a [przynajmniej] jedna z jej wartości własnych wynosi zero. Ta obserwacja natychmiast sugeruje, że istnieją inne , bardziej złożone sposoby uzyskania zdegenerowanej macierzy kowariancji: można mieć macierz kowariancji bez zer i zerowych korelacji, ale mimo to z niedoborem rang (liczba pojedyncza). Bates i in. Parsimonious Mixed Models 2015$4\times 4$$4\times 4$(nieopublikowany przedruk) zaleca użycie analizy głównego składnika (PCA), aby sprawdzić, czy uzyskana macierz kowariancji jest pojedyncza. Jeśli tak, sugerują potraktowanie tej sytuacji w taki sam sposób, jak powyższe sytuacje pojedyncze.

Co więc zrobić?

Jeśli nie ma wystarczających danych do wiarygodnego oszacowania wszystkich parametrów modelu, powinniśmy rozważyć uproszczenie modelu. Biorąc przykładowy model, X*Cond + (X*Cond|subj)istnieją różne możliwe sposoby jego uproszczenia:

1. Usuń jeden z efektów losowych, zwykle korelację najwyższego rzędu:

   X*Cond + (X+Cond|subj)

2. Pozbądź się wszystkich parametrów korelacji:

   X*Cond + (X*Cond||subj)

   Aktualizacja: jak zauważa @Henrik, ||składnia usunie korelacje tylko wtedy, gdy wszystkie zmienne po lewej stronie są liczbowe. Jeśli w grę Condwchodzą zmienne jakościowe (takie jak ), należy raczej użyć jego wygodnego afexpakietu (lub uciążliwych obejść ręcznych). Zobacz jego odpowiedź, aby uzyskać więcej informacji.

3. Pozbądź się niektórych parametrów korelacji, dzieląc ten termin na kilka, np .:

   X*Cond + (X+Cond|subj) + (0+X:Cond|subj)

4. Ogranicz macierz kowariancji w określony sposób, np. Ustawiając jedną konkretną korelację (tę, która osiągnęła granicę) na zero, jak sugerujesz. Nie ma wbudowanego sposobu lme4na osiągnięcie tego. Zobacz odpowiedź @ BenBolkera na temat SO, aby zobaczyć, jak to osiągnąć poprzez inteligentne hakowanie.

Wbrew temu, co powiedziałeś, nie sądzę, że Matuschek i in. 2017 szczególnie polecam # 4. Istota Matuschek i in. 2017 oraz Bates i in. Wydaje się, że rok 2015 zaczyna się od maksymalnego modelu a la Barr i in. 2013, a następnie zmniejsza złożoność, aż macierz kowariancji osiągnie pełną rangę. (Co więcej, często zalecają jeszcze bardziej zmniejszenie złożoności, aby zwiększyć moc.) Aktualizacja: W przeciwieństwie do tego, Barr i in. zaleca się zmniejszenie złożoności WYŁĄCZNIE, jeśli model nie jest zbieżny; są gotowi tolerować pojedyncze macierze kowariancji. Zobacz odpowiedź @ Henrika.

Jeśli ktoś zgadza się z Bates / Matuschek, to uważam, że dobrze jest wypróbować różne sposoby zmniejszenia złożoności, aby znaleźć ten, który wykona zadanie, wyrządzając jednocześnie „najmniej obrażeń”. Patrząc na moją listę powyżej, pierwotna macierz kowariancji ma 10 parametrów; # 1 ma 6 parametrów, # 2 ma 4 parametry, # 3 ma 7 parametrów. Który model pozbyje się idealnych korelacji, nie można powiedzieć bez ich dopasowania.

Ale co jeśli jesteś zainteresowany tym parametrem?

Powyższa dyskusja traktuje macierz kowariancji efektu losowego jako parametr uciążliwy. Zadajesz interesujące pytanie, co zrobić, jeśli jesteś szczególnie zainteresowany parametrem korelacji, który musisz „zrezygnować”, aby uzyskać sensowne rozwiązanie pełnej rangi.

Zauważ, że ustawienie parametru korelacji na zero niekoniecznie da BLUP ( ranef), które są nieskorelowane; w rzeczywistości mogą one nawet nie mieć tak dużego wpływu ( demonstracja znajduje się w odpowiedzi na @ Placidia ). Tak więc jedną opcją byłoby przyjrzenie się korelacjom BLUP i zgłoszenie tego.

Inną, być może mniej atrakcyjną opcją byłoby zastosowanie traktowania subjectjako stałego efektu Y~X*cond*subj, uzyskanie szacunków dla każdego pacjenta i obliczenie korelacji między nimi. Jest to równoważne z uruchomieniem osobnych Y~X*condregresji dla każdego podmiotu osobno i uzyskanie od nich oszacowań korelacji.

Zobacz także sekcję dotyczącą pojedynczych modeli w często zadawanych pytaniach dotyczących modelu mieszanego Bena Bolkera:

Bardzo często w przypadku zbyt wielu modeli mieszanych powstają pojedyncze dopasowania. Technicznie, osobliwość oznacza, że ​​niektóre parametry (rozkład wariancji Choleskiego-kowariancji) odpowiadające elementom diagonalnym współczynnika Choleskiego są dokładnie zerowe, co stanowi krawędź możliwej przestrzeni, lub równoważnie, że macierz wariancji-kowariancji ma jakieś zero wartości własne (tj. jest dodatnim półfinałem raczej niż dodatnim określonym) lub (prawie równoważnie), że niektóre wariancje są szacowane na zero lub niektóre korelacje są szacowane na +/- 1.$\theta$

Zgadzam się ze wszystkim, co zostało powiedziane w odpowiedzi ameby, która stanowi doskonałe podsumowanie bieżącej dyskusji na ten temat. Spróbuję dodać kilka dodatkowych punktów i w innym przypadku odnieść się do materiałów informacyjnych mojego niedawnego kursu modelu mieszanego, który również podsumowuje te punkty.

Tłumienie parametrów korelacji (opcje 2 i 3 w odpowiedzi ameby) za pomocą ||działa tylko dla współzmiennych liczbowych lmerdla czynników, a nie dla czynników. Jest to szczegółowo omówione w kodzie Reinholda Kliegl .

Jednak mój afexpakiet zapewnia funkcję tłumienia korelacji również wśród czynników, jeśli argument expand_re = TRUEw wywołaniu do mixed()(patrz także funkcja lmer_alt()). Zasadniczo robi to, wdrażając podejście omówione przez Reinholda Kliegl (tj. Przenosząc czynniki na zmienne liczbowe i określając na nich strukturę efektów losowych).

Prosty przykład:

library("afex")
data("Machines", package = "MEMSS") # same data as in Kliegl code

# with correlation:
summary(lmer(score ~ Machine + (Machine  | Worker), data=Machines))
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr       
#  Worker   (Intercept) 16.6405  4.0793              
#           MachineB    34.5467  5.8776    0.48      
#           MachineC    13.6150  3.6899   -0.37  0.30
#  Residual              0.9246  0.9616              
# Number of obs: 54, groups:  Worker, 6

## crazy results:
summary(lmer(score ~ Machine + (Machine  || Worker), data=Machines))
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr     
#  Worker   (Intercept)  0.2576  0.5076            
#  Worker.1 MachineA    16.3829  4.0476            
#           MachineB    74.1381  8.6103   0.80     
#           MachineC    19.0099  4.3600   0.62 0.77
#  Residual              0.9246  0.9616            
# Number of obs: 54, groups:  Worker, 6

## as expected:
summary(lmer_alt(score ~ Machine + (Machine  || Worker), data=Machines))
# Random effects:
#  Groups   Name         Variance Std.Dev.
#  Worker   (Intercept)  16.600   4.0743  
#  Worker.1 re1.MachineB 34.684   5.8894  
#  Worker.2 re1.MachineC 13.301   3.6471  
#  Residual               0.926   0.9623  
# Number of obs: 54, groups:  Worker, 6

Dla tych, którzy nie wiedzą afex, główną funkcją modeli mieszanych jest zapewnienie wartości p dla ustalonych efektów, np .:

(m1 <- mixed(score ~ Machine + (Machine  || Worker), data=Machines, expand_re = TRUE))
# Mixed Model Anova Table (Type 3 tests, KR-method)
# 
# Model: score ~ Machine + (Machine || Worker)
# Data: Machines
#    Effect      df        F p.value
# 1 Machine 2, 5.98 20.96 **    .002
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

summary(m1)  
# [...]
# Random effects:
#  Groups   Name         Variance Std.Dev.
#  Worker   (Intercept)  27.4947  5.2435  
#  Worker.1 re1.Machine1  6.6794  2.5845  
#  Worker.2 re1.Machine2 13.8015  3.7150  
#  Residual               0.9265  0.9626  
# Number of obs: 54, groups:  Worker, 6
# [...]

Dale Barr z Barr i in. (2013) praca jest bardziej ostrożna w zalecaniu ograniczenia struktury efektów losowych niż w odpowiedzi Amoeby. W ostatniej wymianie na Twitterze napisał:

• „redukcja modelu wprowadza nieznane ryzyko antykonserwatywności i należy to robić ostrożnie, jeśli w ogóle”. i
• „Moim głównym zmartwieniem jest to, że ludzie rozumieją ryzyko związane z redukcją modelu i że minimalizacja tego ryzyka wymaga bardziej konserwatywnego podejścia niż jest to powszechnie przyjęte (np. Każde nachylenie testowane przy .05)”.

Dlatego zaleca się ostrożność.

Jako jeden z recenzentów mogę również wyjaśnić, dlaczego my, Bates i in. (2015) papier pozostał niepublikowany. Ja i pozostali dwaj recenzenci (którzy podpisali, ale pozostaną tutaj bezimienni) krytykowaliśmy podejście PCA (wydaje się, że nie jest to zasadą podstawową i nie ma dowodów na to, że jest on lepszy pod względem siły). Ponadto uważam, że wszyscy trzej skrytykowali, że artykuł nie skupił się na tym, jak określić strukturę efektów losowych, ale również próbował wprowadzić GAMM. Tak więc artykuł Batesa i in. (2015) przekształcił się w Matuschek i in. (2017), który porusza kwestię struktury efektów losowych za pomocą symulacji oraz Baayen i in. (2017) dokument wprowadzający GAMM.

Moja pełna recenzja Batesa i in. szkic można znaleźć tutaj . IIRC, inne recenzje miały podobne podobne główne punkty.

Ja również miałem ten problem, gdy korzystałem z oszacowania maksymalnego prawdopodobieństwa - tylko używam algorytmu Goldstein IGLS zaimplementowanego za pomocą oprogramowania MLwiN, a nie LME4 w R. Jednak w każdym przypadku problem został rozwiązany, gdy przełączyłem się na oszacowanie MCMC przy użyciu tego samego oprogramowanie. Miałem nawet korelację powyżej 3, która rozwiązała się, gdy zmieniłem oszacowanie. Za pomocą IGLS korelacja jest obliczana po oszacowaniu jako kowariancja podzielona przez iloczyn pierwiastka kwadratowego iloczynu powiązanych wariancji - i to nie uwzględnia niepewności w każdym z szacunków składowych.

Oprogramowanie IGLS nie „wie”, że kowariancja implikuje korelację i po prostu oblicza oszacowania stałej, liniowej, kwadratowej funkcji wariancji itp. Natomiast podejście MCMC opiera się na założeniu próbek z wielowymiarowego rozkładu normalnego, który odpowiada wariancjom i kowariancjom o dobrych właściwościach i pełnym propagowaniu błędów, tak że niepewność przy szacowaniu kowariancji jest uwzględniana przy szacowaniu wariancji i wzajemnie.

MLwin jest łańcuchem szacunkowym MCMC z oszacowaniami IGLS i nieujemną macierzą kowariancji wariancji z wyraźną wariancją może być konieczna zmiana przez zmianę kowariancji na zero przed rozpoczęciem próbkowania.

Dla sprawdzonego przykładu patrz

Opracowanie wielopoziomowych modeli do analizy kontekstu, heterogeniczności i zmian za pomocą MLwiN 3, tom 1 (zaktualizowany wrzesień 2017 r.); Tom 2 znajduje się również w RGate

https://www.researchgate.net/publication/320197425_Vol1Training_manualRevisedSept2017

Dodatek do rozdziału 10