Wie kann ich die Berechnung der festen Effekte in einem GLMM beschleunigen?

Ich mache eine Simulationsstudie, die Bootstrapping-Schätzungen erfordert, die aus einem verallgemeinerten linearen gemischten Modell erhalten wurden (tatsächlich das Produkt von zwei Schätzungen für feste Effekte, eine von einem GLMM und eine von einem LMM). Um die Studie gut durchzuführen, wären etwa 1000 Simulationen mit jeweils 1000 oder 1500 Bootstrap-Replikationen erforderlich. Dies nimmt auf meinem Computer viel Zeit in Anspruch (viele Tage).

How can I speed up the computation of these fixed effects?

Um genauer zu sein, ich habe Probanden, die wiederholt auf drei Arten gemessen werden, was zu Variablen X, M und Y führt, wobei X und M stetig und Y binär sind. Wir haben zwei Regressionsgleichungen wobei Y die zugrunde liegende latente kontinuierliche Variable für und die Fehler nicht iid sind. Die Statistik, die wir booten möchten, lautet . Daher erfordert jede Bootstrap-Replikation die Anpassung eines LMM und eines GLMM. Mein R-Code ist (mit lme4)Y ∗ = β 0 + β 1 X + β 2 M + ϵ 2 ∗ Y α 1 β

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

Mir ist klar, dass ich die gleiche Schätzung für erhalte, wenn ich es nur als lineares Modell anpasse, was einige Zeit spart, aber der gleiche Trick funktioniert nicht für .β $\alpha_1$$\beta_2$

Muss ich nur einen schnelleren Computer kaufen? :) :)

Es sollte hilfreich sein, Startwerte anzugeben, obwohl es schwierig ist zu wissen, wie viel. Während Sie Simulation und Bootstrapping durchführen, sollten Sie die "wahren" Werte oder die nicht gebooteten Schätzungen oder beides kennen. Versuchen Sie, diese in der start =Option von zu verwenden glmer.

Sie können auch prüfen, ob die Toleranzen für die Konvergenzerklärung strenger sind als erforderlich. Ich bin mir jedoch nicht sicher, wie ich sie aus der lme4Dokumentation ändern soll .

Zwei weitere Möglichkeiten sollten vor dem Kauf eines neuen Computers in Betracht gezogen werden.

1. Paralleles Rechnen - Bootstrapping kann einfach parallel ausgeführt werden. Wenn Ihr Computer relativ neu ist, haben Sie wahrscheinlich vier Kerne. Schauen Sie sich die Multicore- Bibliothek in R an.
2. Cloud Computing ist ebenfalls eine Möglichkeit und relativ billig. Ich habe Kollegen, die die Amazon Cloud zum Ausführen von R-Skripten verwendet haben. Sie fanden, dass es ziemlich kostengünstig war.

Es könnte möglicherweise ein schnellerer Computer sein. Aber hier ist ein Trick, der funktionieren kann.

Generieren Sie eine Simulation von , jedoch nur unter der Bedingung von , und führen Sie dann OLS oder LMM für die simulierten -Werte aus.$Y^*$$Y$$Y^*$

Angenommen, Ihre Link-Funktion ist . Dies sagt aus, wie Sie von der Wahrscheinlichkeit von zum -Wert gelangen, und ist höchstwahrscheinlich die logistische Funktion .$g(.)$$Y=1$$Y^*$$g(z)=log \Big(\frac{z}{1-z}\Big)$

Wenn Sie also eine Bernouli-Stichprobenverteilung für annehmen und dann die Jeffreys vor der Wahrscheinlichkeit verwenden, erhalten Sie einen Beta-Posterior für . Daraus zu simulieren sollte wie Beleuchtung sein, und wenn dies nicht der Fall ist, benötigen Sie einen schnelleren Computer. Darüber hinaus sind die Stichproben unabhängig, sodass Sie keine "Konvergenz" -Diagnose wie in MCMC überprüfen müssen und wahrscheinlich nicht so viele Stichproben benötigen - 100 funktionieren möglicherweise gut für Ihren Fall. Wenn Sie Binomial- , ersetzen Sie einfach die im obigen Posterior durch , die Anzahl der Versuche des Binomials für jedes .$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$

Sie haben also eine Reihe von simulierten Werten . Anschließend wenden Sie die Verknüpfungsfunktion auf jeden dieser Werte an, um . Passen Sie ein LMM an , was wahrscheinlich schneller ist als das GLMM-Programm. Sie können die ursprünglichen Binärwerte grundsätzlich ignorieren (aber nicht löschen!) Und einfach mit der "Simulationsmatrix" arbeiten ( , wobei die Stichprobengröße und die Anzahl der Simulationen ist).$p_{sim}$$Y_{sim}=g(p_{sim})$$Y_{sim}$$N\times S$$N$$S$

In Ihrem Programm würde ich also die Funktion Funktion und durch eine einzige Simultation ersetzen. Sie würden dann eine Art Schleife erstellen, die die Funktion auf jede Simulation , und dann die Durchschnitt als Schätzung von . So etwas wie $gmler()$$lmer()$$Y$$lmer()$$b$

Lassen Sie mich wissen, wenn ich etwas klarer erklären muss