I would like to do something in R that SAS can do using SAS's proc mixed (there is some way to do in STATA es well), namely fitting the so called Bivariate model from Reitsma et al (2005). This model is a special mixed model where the variance depends on the study (see below). Googling and talking to some people familiar with the model did not yield a straightforward approach that is fast at the same time (i.e. a nice high level model fitting function). I am nevertheless sure, there is something fast in R that one can built on.

In a nutshell one is faced with the following situation: Given pairs of proportions $(p_1,p_2)$ in $[0,1]^2$ one would like to fit a bivariate normal to the logit-transformed pairs. Since the proportions come from a 2x2 table (i.e. binomial data) each logit transformed observed proportion has a variance estimate that is to be included in the fitting process, say $(s_1, s_2)$. So one would like to fit a bivariate normal to the pairs, where the covariance matrix $\Sigma$ depends on the observation, i.e.

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

where S is the diagonal matrix with $(s_1, s_2)$ and depends entirely on the data but varies from observation to observation. mu and Sigma are the same for all observation though.

Right now I am using a call to optim() (using BFGS) to estimate the five parameters ($\mu_1$, $\mu_2$, and three parameters for $\Sigma$). Nevertheless this is painfully slow, and especially unsuitable for simulation. Also one of my aims is to introduce regression coefficients for mu later, increasing the number of parameters.

I tried speeding up fitting by supplying starting values and I also thought about computing gradients for the five parameters. Since the likelihood becomes quite complex due to the addition of $S$, I felt the risk of introducting errors this way was too big and did not attempt it yet, nor did I see a way to check my calculations.

Is the calculation of the gradients typically worthwhile? How do you check them?

I am aware of other optimizer besides optim(), i.e. nlm() and I also know about the CRAN Task view: Optimization. Which ones a are worth a try?

What kind of tricks are there to speed up optim() besides reducing accuracy?

I would be very grateful for any hints.

  • $\begingroup$ Some discussion here $\endgroup$
    – Henry
    Jun 30, 2011 at 23:55
  • $\begingroup$ The answers to this question may be useful: stats.stackexchange.com/questions/9535/… $\endgroup$ Jul 3, 2011 at 15:08
  • $\begingroup$ @mbq: Thanks for migrating the question. Also thank you for hinting at the other questions. $\endgroup$
    – Philipp
    Jul 5, 2011 at 19:05

1 Answer 1


This isn't perhaps the solution you anticipated, but I think you could fit this model with brms (i.e. as an interface to Stan), with what it calls a 'distributional model'... see https://cran.r-project.org/web/packages/brms/vignettes/brms_distreg.html

See also the overview vignette here which shows how to fit the binomial model. https://cran.r-project.org/web/packages/brms/vignettes/brms_overview.pdf


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.