Ajuste de un modelo mixto especial en R: alternativas a optim ()

Me gustaría hacer algo en R que SAS pueda hacer usando el proceso mixto de SAS (hay alguna manera de hacerlo en STATA es bien), es decir, ajustar el llamado modelo Bivariate de Reitsma et al (2005). Este modelo es un modelo mixto especial donde la varianza depende del estudio (ver más abajo). Buscar en Google y hablar con algunas personas familiarizadas con el modelo no produjo un enfoque directo que sea rápido al mismo tiempo (es decir, una buena función de ajuste del modelo de alto nivel). Sin embargo, estoy seguro de que hay algo rápido en R que se puede construir.

En pocas palabras, uno se enfrenta a la siguiente situación: pares de proporciones dados $(p_1,p_2)$ en $[0,1]^2$A uno le gustaría ajustar una normal bivariada a los pares transformados logit. Dado que las proporciones provienen de una tabla de 2x2 (es decir, datos binomiales), cada proporción observada transformada logit tiene una estimación de la varianza que se debe incluir en el proceso de ajuste, digamos$(s_1, s_2)$. Entonces, uno quisiera ajustar una bivariada normal a los pares, donde la matriz de covarianza$\Sigma$ depende de la observación, es decir

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

donde S es la matriz diagonal con $(s_1, s_2)$y depende completamente de los datos, pero varía de una observación a otra. Sin embargo, mu y Sigma son iguales para toda observación.

En este momento estoy usando una llamada a optim()(usando BFGS) para estimar los cinco parámetros ($\mu_1$, $\mu_2$y tres parámetros para $\Sigma$) Sin embargo, esto es dolorosamente lento y especialmente inadecuado para la simulación. También uno de mis objetivos es introducir coeficientes de regresión para mu más tarde, aumentando el número de parámetros.

Intenté acelerar el ajuste al proporcionar valores iniciales y también pensé en calcular gradientes para los cinco parámetros. Dado que la probabilidad se vuelve bastante compleja debido a la adición de$S$, Sentí que el riesgo de introducir errores de esta manera era demasiado grande y aún no lo intenté, ni vi una forma de verificar mis cálculos.

¿Vale la pena el cálculo de los gradientes? ¿Cómo los revisas?

Además optim(), conozco otro optimizador , es decir, nlm()y también conozco la vista de tareas de CRAN: Optimización. ¿Cuáles valen la pena probar?

¿Qué tipo de trucos hay para acelerar optim()además de reducir la precisión?

Estaría muy agradecido por cualquier pista.

Tal vez esta no sea la solución que esperaba, pero creo que podría ajustar este modelo con brms (es decir, como una interfaz para Stan), con lo que llama un "modelo de distribución" ... consulte https: //cran.r- project.org/web/packages/brms/vignettes/brms_distreg.html

Vea también la viñeta de resumen aquí que muestra cómo ajustar el modelo binomial. https://cran.r-project.org/web/packages/brms/vignettes/brms_overview.pdf