Robię badanie symulacyjne, które wymaga oszacowań ładowania uzyskanych z uogólnionego liniowego modelu mieszanego (w rzeczywistości iloczyn dwóch oszacowań dla ustalonych efektów, jednego z GLMM i jednego z LMM). Prawidłowe wykonanie badania wymagałoby około 1000 symulacji z 1000 lub 1500 replikami ładowania początkowego za każdym razem. To zajmuje dużo czasu na moim komputerze (wiele dni).
How can I speed up the computation of these fixed effects?
Mówiąc ściślej, mam podmioty, które są mierzone wielokrotnie na trzy sposoby, co prowadzi do zmiennych X, M i Y, gdzie X i M są ciągłe, a Y jest binarne. Mamy dwa równania regresji gdzie Y ^ * jest ukrytą ukrytą zmienną ciągłą dla Y, a błędy nie są ididalne. Statystyka, którą chcemy uruchomić to \ alpha_1 \ beta_2 . Dlatego każda replikacja bootstrap wymaga dopasowania LMM i GLMM. Mój kod R to (używając lme4)
stat=function(dat){
a=fixef(lmer(M~X+(X|person),data=dat))["X"]
b=fixef(glmer(Y~X+M+(X+M|person),data=dat,family="binomial"))["M"]
return(a*b)
}
Zdaję sobie sprawę, że otrzymam takie same oszacowanie dla jeśli po prostu dopasuję go jako model liniowy, dzięki czemu zaoszczędzisz trochę czasu, ale ta sama sztuczka nie działa dla .
Czy muszę po prostu kupić szybszy komputer? :)
źródło
Rprof
.Odpowiedzi:
Powinno to pomóc określić wartości początkowe, choć trudno jest ustalić, ile. Podczas przeprowadzania symulacji i ładowania początkowego powinieneś znać „prawdziwe” wartości, szacunki braku ładowania początkowego lub jedno i drugie. Spróbuj użyć tych w
start =
opcjiglmer
.Możesz również rozważyć, czy tolerancje dla deklarowania konwergencji są bardziej rygorystyczne niż potrzebujesz. Nie jestem jednak pewien, jak je zmienić w
lme4
dokumentacji.źródło
Przed zakupem nowego komputera rozważ też dwie inne możliwości.
źródło
Może to być szybszy komputer. Ale oto jedna sztuczka, która może zadziałać.
Wygeneruj równoczesność , ale tylko warunkowo na , a następnie po prostu wykonaj OLS lub LMM na symulowanych wartościach .Y∗ Y Y∗
Załóżmy, że twoją funkcją linku jest . to mówi, w jaki sposób uzyskujesz prawdopodobieństwo do wartości , i najprawdopodobniej jest to funkcja logistyczna .sol( . ) Y= 1 Y∗ sol( z) = l o g(z1 - z)
Więc jeśli przyjmiesz rozkład próbkowania bernouli dla , a następnie użyjesz jeffreys przed prawdopodobieństwem, otrzymasz wersję beta później dla . Symulacja z tego powinna być jak oświetlenie, a jeśli nie jest, potrzebujesz szybszego komputera. Co więcej, próbki są niezależne, więc nie trzeba sprawdzać żadnej diagnostyki „zbieżności”, takiej jak w MCMC, i prawdopodobnie nie potrzebujesz tak wielu próbek - 100 może działać dobrze w twoim przypadku. Jeśli masz dwumianowe , to po prostu zamień w powyższym tylnej na , liczbę prób dwumianu dla każdego .Y→ Y∼ B e r n o u l l i ( p ) p ∼ B e t a (Yo b s+12), 1 -Yo b s+12)) Y′s 1 nja Yi
Masz więc zestaw symulowanych wartości . Następnie zastosuj funkcję link do każdej z tych wartości, aby uzyskać . Dopasuj LMM do , co jest prawdopodobnie szybsze niż program GLMM. Możesz w zasadzie zignorować oryginalne wartości binarne (ale nie usuwaj ich!) I po prostu pracuj z „macierzą symulacji” ( , gdzie to wielkość próbki, a to liczba symulacji).psim Ysim=g(psim) Ysim N×S N S
Tak więc w twoim programie zastąpiłbym funkcję funkcją , a pojedynczą symultanią, a następnie utworzyłbyś pętlę, która stosuje funkcję do każdej symulacji, a następnie bierze średnia jako szacunkowa wartość . Coś w rodzajugmler() lmer() Y lmer() b
Daj mi znać, jeśli będę musiał coś wyjaśnić
źródło