of the $F\bigl(\cdot;l(\pi),\de\bigr)$ and $l$ in the one-factor copula mixed model.

Relationship with the $2T$-variate GLMM
---------------------------------------

We show what happens when all the bivariate copulas are BVN and the univariate distribution of the random effects is the $N(\mu,\s)$ distribution. One can easily deduce that the within-study model in (\[withinBinom\]) is the same as in the $2T$-variate GLMM.

Furthermore, when all the bivariate copulas are BVN copulas with correlation parameters $\rhobf=(\rho_{11},\ldots,$ $\rho_{1t},\ldots,\rho_{1T},\rho_{01},\ldots,\rho_{0t},\ldots,\rho_{0T})$, the resulting random effects distribution is the $2T$-variate normal distribution with mean vector $\bmu=\bigl(l(\bfpi_1),l(\bfpi_0)\bigr)$ and variance-covariance matrix $$\Sbf=\begin{pmatrix}
\s_{11}^2&\cdots&\rho_{11}\rho_{1T}\s_{11}\s_{1T}&\rho_{11}\rho_{01}\s_{11}\s_{01}&\cdots&\rho_{11}\rho_{0T}\s_{11}\s_{0T}\\
\vdots &\ddots&\vdots& \vdots&\vdots&\vdots\\
\rho_{1T}\rho_{11}\s_{1T}\s_{11}&\cdots&\s_{1T}^2&\rho_{1T}\rho_{01}\s_{1T}\s_{01}&\cdots&\rho_{1T}\rho_{0T}\s_{1T}\s_{0T}\\
\rho_{01}\rho_{11}\s_{01}\s_{11}&\cdots&\rho_{01}\rho_{1T}\s_{01}\s_{1T}&\s_{01}^2&\cdots&\rho_{1T}\rho_{0T}\s_{1T}\s_{0T}\\
\vdots &\vdots&\vdots& \vdots&\ddots&\vdots\\
\rho_{0T}\rho_{11}\s_{0T}\s_{11}&\cdots&\rho_{0T}\rho_{1T}\s_{0T}\s_{1T}&\rho_{0T}\rho_{01}\s_{0T}\s_{01}&\cdots&\s_{0T}^2\\
\end{pmatrix}$$ Hence, the proposed model has as special case the $2T$-variate GLMM.

\[computation\]Maximum likelihood estimation and computational details
----------------------------------------------------------------------

Estimation of the model parameters $(\bfpi_1,\bfpi_0,\debf_1,\debf_0,\thbf_1,\thbf_0)$ can be approached by the standard maximum likelihood (ML) method, by maximizing the logarithm of the joint likelihood in (\[mixed-cop-likelihood\]). The estimated parameters can be obtained by using a quasi-Newton [@nash90] method applied to the logarithm of the joint likelihood. This numerical method requires only the objective function, i.e., the logarithm of the joint likelihood, while the gradients are computed numerically and the Hessian matrix of the second order derivatives is updated in each iteration. 