Anpassung eines speziellen gemischten Modells in R - Alternativen zu optim ()

Ich möchte etwas in R tun, das SAS mit SASs gemischten Prozessen tun kann (es gibt eine Möglichkeit, dies in STATA gut zu tun), nämlich das sogenannte bivariate Modell von Reitsma et al. (2005) anzupassen. Dieses Modell ist ein spezielles gemischtes Modell, bei dem die Varianz von der Studie abhängt (siehe unten). Das Googeln und Sprechen mit einigen Personen, die mit dem Modell vertraut sind, ergab keinen einfachen Ansatz, der gleichzeitig schnell ist (dh eine nette Modellanpassungsfunktion auf hohem Niveau). Ich bin mir trotzdem sicher, dass es in R etwas Schnelles gibt, auf das man aufbauen kann.

Kurz gesagt, man sieht sich folgender Situation gegenüber: Bei gegebenen Proportionspaaren in [ 0 , 1 ] 2 möchte man eine bivariate Normale an die logit-transformierten Paare anpassen. Da die Proportionen aus einer 2x2-Tabelle (dh Binomialdaten) stammen, hat jeder logit-transformierte beobachtete Anteil eine Varianzschätzung, die in den Anpassungsprozess einbezogen werden soll, z. B. ( s 1 , s 2 ) . Man möchte also eine bivariate Normale an die Paare anpassen, wobei die Kovarianzmatrix Σ von der Beobachtung abhängt , d. H.$(p_1,p_2)$$[0,1]^2$$(s_1, s_2)$$\Sigma$

$(\text{logit}(p_1),\text{logit}(p_2)) \sim N((mu_1, mu_2), \Sigma + S)$

Dabei ist S die Diagonalmatrix mit und hängt vollständig von den Daten ab, variiert jedoch von Beobachtung zu Beobachtung. mu und Sigma sind jedoch für alle Beobachtungen gleich.$(s_1, s_2)$

Im Moment benutze ich einen Aufruf von optim()(unter Verwendung von BFGS), um die fünf Parameter ( , und drei Parameter für ) zu schätzen . Dies ist jedoch schmerzhaft langsam und für die Simulation besonders ungeeignet. Eines meiner Ziele ist es auch, später Regressionskoeffizienten für mu einzuführen und die Anzahl der Parameter zu erhöhen.$\mu_1$$\mu_2$$\Sigma$

Ich habe versucht, die Anpassung durch Angabe von Startwerten zu beschleunigen, und ich habe auch darüber nachgedacht, Gradienten für die fünf Parameter zu berechnen. Da die Wahrscheinlichkeit aufgrund der Hinzufügung von ziemlich komplex wird , hielt ich das Risiko, Fehler auf diese Weise einzuführen, für zu groß und versuchte es noch nicht, noch sah ich eine Möglichkeit, meine Berechnungen zu überprüfen.$S$

Lohnt sich die Berechnung der Gradienten in der Regel? Wie überprüfst du sie?

Außerdem sind mir andere Optimierer bekannt optim(), dh nlm()ich kenne auch die CRAN-Aufgabenansicht: Optimierung. Welche sind einen Versuch wert?

Welche Tricks gibt es optim()neben der Verringerung der Genauigkeit, um zu beschleunigen ?

Für Hinweise wäre ich sehr dankbar.

Dies ist vielleicht nicht die Lösung, die Sie erwartet haben, aber ich denke, Sie könnten dieses Modell mit brms (dh als Schnittstelle zu Stan) ausstatten, mit dem, was es als "Verteilungsmodell" bezeichnet ... siehe https: //cran.r- project.org/web/packages/brms/vignettes/brms_distreg.html

Siehe auch die Übersichtsvignette hier, die zeigt, wie das Binomialmodell angepasst wird. https://cran.r-project.org/web/packages/brms/vignettes/brms_overview.pdf