Ajustando um modelo misto especial em R - alternativas para optim ()

Eu gostaria de fazer algo em R que o SAS possa fazer usando o proc misto do SAS (existe alguma maneira de se fazer isso no STATA), ou seja, ajustar o chamado modelo bivariado de Reitsma et al (2005). Este modelo é um modelo misto especial em que a variação depende do estudo (veja abaixo). Pesquisando e conversando com algumas pessoas familiarizadas com o modelo não produzimos uma abordagem direta e rápida ao mesmo tempo (ou seja, uma boa função de ajuste de modelo de alto nível). No entanto, tenho certeza de que há algo rápido no R que se pode construir.

Em poucas palavras, nos deparamos com a seguinte situação: Dados pares de proporções em , gostaríamos de ajustar uma normal bivariada aos pares transformados em logit. Como as proporções vêm de uma tabela 2x2 (isto é, dados binomiais), cada proporção observada transformada em logit possui uma estimativa de variação que deve ser incluída no processo de ajuste, digamos . Então, alguém gostaria de ajustar uma normal bivariada aos pares, onde a matriz de covariância depende da observação, ou seja,$(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)$

onde S é a matriz diagonal com e depende inteiramente dos dados, mas varia de observação para observação. mu e Sigma são os mesmos para todas as observações.$(s_1, s_2)$

No momento, estou usando uma chamada para optim()(usando BFGS) para estimar os cinco parâmetros ( , e três parâmetros para ). No entanto, isso é dolorosamente lento e especialmente inadequado para simulação. Também um dos meus objetivos é introduzir coeficientes de regressão para mu posteriormente, aumentando o número de parâmetros.$\mu_1$$\mu_2$$\Sigma$

Tentei acelerar o ajuste fornecendo valores iniciais e também pensei em calcular gradientes para os cinco parâmetros. Como a probabilidade se torna bastante complexa devido à adição de , senti o risco de introduzir erros dessa maneira era muito grande e não tentei ainda, nem vi uma maneira de verificar meus cálculos.$S$

O cálculo dos gradientes normalmente vale a pena? Como você os verifica?

Além disso optim(), conheço outro otimizador , ou seja, nlm()e também conheço a visualização Tarefa CRAN: Otimização. Quais valem a pena tentar?

Que tipo de truques existem para acelerar, optim()além de reduzir a precisão?

Ficaria muito grato por qualquer dica.

Talvez essa não seja a solução que você antecipou, mas acho que você poderia ajustar esse modelo com brms (como uma interface para Stan), com o que ele chama de 'modelo distributivo' ... veja https: //cran.r- project.org/web/packages/brms/vignettes/brms_distreg.html

Veja também a vinheta de visão geral aqui, que mostra como ajustar o modelo binomial. https://cran.r-project.org/web/packages/brms/vignettes/brms_overview.pdf