Jak mogę przyspieszyć obliczanie ustalonych efektów w GLMM?

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 a błędy nie są ididalne. Statystyka, którą chcemy uruchomić to . Dlatego każda replikacja bootstrap wymaga dopasowania LMM i GLMM. Mój kod R to (używając lme4)

M = α_{0} + α_{1} X + ϵ_{1}

$M=\alpha_0+\alpha_1X+\epsilon_1$

Y^{*} = β_{0} + β_{1} X + β_{2} M + ϵ_{2}

$Y^*=\beta_0+\beta_1X+\beta_2M+\epsilon_2$

^{*}

$^*$

Y

$Y$

α_{1} β_{2}

$\alpha_1\beta_2$

    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 $\alpha_1$ jeśli po prostu dopasuję go jako model liniowy, dzięki czemu zaoszczędzisz trochę czasu, ale ta sama sztuczka nie działa dla $\beta_2$ .

Czy muszę po prostu kupić szybszy komputer? :)

r mixed-model

— BR
źródło

@BR co to jest szyjka butelki tutaj? Zasadniczo, co zajmuje czas Rprof.

— suncoolsu

Jednym ze sposobów jest po prostu zignorowanie „mieszanej” części GLMM. Wystarczy dopasować zwykły GLM, szacunki niewiele się zmienią, ale prawdopodobnie zmienią się ich standardowe błędy

— probabilislogic

@probabilityislogic. Oprócz twojej uwagi, myślę również, czy odpowiedź różni się znacznie, zależy od wielkości grupy i indywidualnych zachowań w grupie. Jak twierdzą Gelman i Hill: wyniki modelu mieszanych efektów byłyby między pulą a brakiem puli. (Oczywiście dotyczy to bayesowskich modeli hierarchicznych, ale modele mieszane to klasyczny sposób robienia tego samego.)

— suncoolsu

@probabilityislogic: Działa to w przypadku LMM, ale wydaje się, że zawodzi w przypadku GLMM (co oznacza, że uruchomiłem modele z dodatkowym M na tych samych danych i bez nich, co zakończyło się znacząco innymi wynikami). Chyba że, oczywiście, wystąpił błąd w implementacji glmer.

— BR

@suncoolsu: wąskim gardłem jest oszacowanie GLMM, które może potrwać kilka sekund (szczególnie z kilkoma efektami losowymi). Ale zrób to 1000 * 1000 razy, a to 280 godzin obliczeń. Montaż GLM zajmuje około 1/100 czasu.

— BR

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 =opcji glmer.

Możesz również rozważyć, czy tolerancje dla deklarowania konwergencji są bardziej rygorystyczne niż potrzebujesz. Nie jestem jednak pewien, jak je zmienić w lme4dokumentacji.

— jeden przystanek
źródło

Przed zakupem nowego komputera rozważ też dwie inne możliwości.

Przetwarzanie równoległe - ładowanie jest łatwe do uruchomienia równoległego. Jeśli twój komputer jest dość nowy, prawdopodobnie masz cztery rdzenie. Spójrz na bibliotekę wielordzeniową w R.
Przetwarzanie w chmurze jest również możliwe i niedrogie. Mam kolegów, którzy używali chmury Amazon do uruchamiania skryptów R. Stwierdzili, że było to dość opłacalne.

— csgillespie
źródło

Dziękuję za odpowiedź. Jakoś przeoczyłem fakt, że mam dwa rdzenie (mój komputer nie jest bardzo nowy). Dawno powinienem był spojrzeć na multicore.

— BR

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 . $g(.)$ $Y=1$ $Y^*$ $g(z)=log \Big(\frac{z}{1-z}\Big)$

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\rightarrow Y\sim Bernoulli(p)$ $p\sim Beta(Y_{obs}+\frac{1}{2},1-Y_{obs}+\frac{1}{2})$ $Y's$ $1$ $n_i$ $Y_i$

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). $p_{sim}$ $Y_{sim}=g(p_{sim})$ $Y_{sim}$ $N\times 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 rodzaju $gmler()$ $lmer()$ $Y$ $lmer()$ $b$

a = \dots

$a=\dots$

b = 0

$b=0$

d o s = 1, \dots, S

$do \ s=1,\dots,S$

b_{e s t} = l m e r (Y_{s} \dots)

$b_{est}=lmer(Y_s\dots)$

b = b + \frac{1}{s} (b_{e s t} - b)

$b=b+\frac{1}{s}(b_{est}-b)$

e n d

$end$

r e t u r n (a * b)

$return(a*b)$

Daj mi znać, jeśli będę musiał coś wyjaśnić

— prawdopodobieństwo prawdopodobieństwa
źródło

Dzięki za odpowiedź, zajmie mi to trochę czasu (i mam już plany na sobotnią noc). Jest na tyle inny, że nie jest dla mnie jasne, czy daje taką samą odpowiedź jak podejście GLMM, ale muszę o tym więcej pomyśleć.

— BR