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.
