Montage d'un modèle mixte spécial dans R - alternatives à optim ()

Je voudrais faire quelque chose en R que SAS peut faire en utilisant le proc mixte de SAS (il y a bien un moyen de le faire en STATA es), à savoir adapter le soi-disant modèle bivarié de Reitsma et al (2005). Ce modèle est un modèle mixte spécial où la variance dépend de l'étude (voir ci-dessous). La recherche sur Google et la discussion avec des personnes familières avec le modèle n'ont pas donné lieu à une approche simple et rapide en même temps (c'est-à-dire une belle fonction d'ajustement de modèle de haut niveau). Je suis néanmoins sûr qu'il y a quelque chose de rapide dans R sur lequel on peut s'appuyer.

En un mot, on est confronté à la situation suivante: étant donné les paires de proportions dans on voudrait ajuster un bivarié normal aux paires transformées en logit. Étant donné que les proportions proviennent d'un tableau 2x2 (c'est-à-dire des données binomiales), chaque proportion observée transformée logit a une estimation de la variance qui doit être incluse dans le processus d'ajustement, disons . On voudrait donc ajuster une normale bivariée aux paires, où la matrice de covariance dépend de l'observation, ie$(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)$

où S est la matrice diagonale avec et dépend entièrement des données mais varie d'une observation à l'autre. mu et Sigma sont les mêmes pour toutes les observations.$(s_1, s_2)$

En ce moment, j'utilise un appel à optim()(en utilisant BFGS) pour estimer les cinq paramètres ( , et trois paramètres pour ). Néanmoins, ceci est douloureusement lent, et particulièrement impropre à la simulation. Un de mes objectifs est également d'introduire des coefficients de régression pour mu plus tard, en augmentant le nombre de paramètres.$\mu_1$$\mu_2$$\Sigma$

J'ai essayé d'accélérer l'ajustement en fournissant des valeurs de départ et j'ai également pensé à calculer des gradients pour les cinq paramètres. Étant donné que la probabilité devient assez complexe en raison de l'ajout de , j'ai senti que le risque d'introduire des erreurs de cette manière était trop grand et je ne l'ai pas encore essayé, et je n'ai pas vu de moyen de vérifier mes calculs.$S$

Le calcul des gradients vaut-il généralement la peine? Comment les vérifiez-vous?

Je connais d'autres optimiseurs en plus optim(), c'est-à nlm()- dire et je connais également la vue Tâche CRAN: Optimisation. Lesquelles valent la peine d'être essayées?

Quels types d'astuces sont disponibles pour accélérer en optim()plus de réduire la précision?

Je serais très reconnaissant pour tout indice.

Ce n'est peut-être pas la solution que vous attendiez, mais je pense que vous pourriez équiper ce modèle de brms (c'est-à-dire en tant qu'interface avec Stan), avec ce qu'il appelle un «modèle distributionnel» ... voir https: //cran.r- project.org/web/packages/brms/vignettes/brms_distreg.html

Voir également la vignette de présentation ici qui montre comment ajuster le modèle binomial. https://cran.r-project.org/web/packages/brms/vignettes/brms_overview.pdf