First, consider the sum of univariate log-likelihoods
$$L_1= \sum_{i=1}^n\sum_{j=1}^d\, \log f_1(y_{ij};\nu_{ij},\boldsymbol{\gamma})=\sum_{i=1}^n\sum_{j=1}^d\,
\ell_1(\nu_{ij},\boldsymbol{\gamma},\, y_{ij}),
$$
and then the sum of bivariate log-likelihoods
$$L_2=\sum_{i=1}^{n}\sum_{j<k}
\log{f_2(y_{ij},y_{ik};\nu_{ij},\nu_{ik},\boldsymbol{\gamma},\rho_{jk})}=\sum_{i=1}^{n}\sum_{j<k}\ell_2(\nu_{ij},\nu_{ik},\boldsymbol{\gamma},\rho_{jk};y_{ij},y_{ik}),$$
where \(\ell_2(\cdot)=\log f_2(\cdot)\) and
$$
f_2(y_{ij},y_{ik};\nu_{ij},\nu_{ik},\boldsymbol{\gamma},\rho_{jk})=\int_{\Phi^{-1}[F_1(y_{ij}-1;\nu_{ij},\boldsymbol{\gamma})]}^{\Phi^{-1}[F_1(y_{ij};\nu_{ij},\boldsymbol{\gamma})]}
\int_{\Phi^{-1}[F_1(y_{ik}-1;\nu_{ik}),\boldsymbol{\gamma}]}^{\Phi^{-1}[F_1(y_{ik};\nu_{ik},\boldsymbol{\gamma})]} \phi_2(z_j,z_d;\rho_{jk}) dz_j dz_k;
$$
\(\phi_2(\cdot;\rho)\) denotes the standard bivariate normal density with correlation
\(\rho\).
Let \(\mathbf{a}^\top=(\boldsymbol{\beta}^\top,\boldsymbol{\gamma}^\top)\) be the
column vector of all \(r=p+q\) univariate parameters. Differentiating \(L_1\) with respect to \(\mathbf{a}\) leads to the independent estimating equations or univariate composite score functions:
$$
\mathbf{g}_1=\mathbf{g}_1(\mathbf{a})=
\frac{\partial L_1}{\partial \mathbf{a}}=
\sum_{i=1}^n\mathbf{X}_i^\top\mathbf{s}_i^{(1)}(\mathbf{a})=\bf 0,
$$
Differentiating \(L_2\) with respect to \(\mathbf{R}=\bigl(\rho_{jk},1\leq j<k\leq d\bigr)\)
leads to the bivariate composite score functions (Zhao and Joe, 2005):
$$
\mathbf{g}_2=\frac{\partial L_2}{\partial \mathbf{R}}= \sum_{i=1}^{n}\mathbf{s}_i^{(2)}(\mathbf{a},\mathbf{R})= \sum_{i=1}^n\Bigl(\mathbf{s}_{i,jk}^{(2)}(\mathbf{a},\rho_{jk}),1\leq j<k\leq d\Bigr)=\mathbf{0},
$$
where \(\mathbf{s}_i^{(2)}(\mathbf{a},\mathbf{R})=\frac{\partial\sum_{j<k}\ell_2(\nu_{ij},\nu_{ik},\boldsymbol{\gamma},\rho_{jk};y_{ij},y_{ik})}{\partial \mathbf{R}}\) and \(\mathbf{s}_{i,jk}^{(2)}(\boldsymbol{\gamma},\rho_{jk})=\frac{\partial \ell_2(\nu_{ij},\nu_{ik},\boldsymbol{\gamma},\rho_{jk};y_{ij},y_{ik})}{\partial \rho_{jk}}\).
The CL1 estimates \(\widetilde{\mathbf{a}}\) and \(\widetilde{\mathbf{R}}\) of the discretized MVN model are obtained
by solving the above CL1 estimating functions.
The asymptotic covariance matrix for the estimator that solves them, also known as the inverse Godambe (Godambe, 1991) information matrix, is
$$
\mathbf{V}=(-\mathbf{H}_\mathbf{g})^{-1}\mathbf{J}_\mathbf{g}(-\mathbf{H}_\mathbf{g}^\top)^{-1},
$$
where \(\mathbf{g}=(\mathbf{g}_1,\mathbf{g}_2)^\top\). First set \(\boldsymbol{\theta}=(\mathbf{a},\mathbf{R})^\top\), then
$$
-\mathbf{H}_\mathbf{g}=E\Bigl(\frac{\partial \mathbf{g}}{\partial\boldsymbol{\theta}}\Bigr)=
\left[\begin{array}{cc}
E\Bigl(\frac{\partial \mathbf{g}_1}{\partial\mathbf{a}}\Bigr)
& E\Bigl(\frac{\partial \mathbf{g}_1}{\partial\mathbf{R}}\Bigr)\\
E\Bigl(\frac{\partial \mathbf{g}_2}{\partial\mathbf{a}}\Bigr) &
E\Bigl(\frac{\partial \mathbf{g}_2}{\partial\mathbf{R}}\Bigr)
\end{array}\right]=
\left[\begin{array}{cc}
-\mathbf{H}_{\mathbf{g}_1}&\mathbf{0}\\
-\mathbf{H}_{\mathbf{g}_{2,1}}&-\mathbf{H}_{\mathbf{g}_2}
\end{array}\right],
$$
where \(-\mathbf{H}_{\mathbf{g}_1}=\sum_i^n\mathbf{X}_i^\top\boldsymbol{\Delta}_i^{(1)}
\mathbf{X}_i\),
\(-\mathbf{H}_{\mathbf{g}_{2,1}}=\sum_i^n\boldsymbol{\Delta}_i^{(2,1)}\mathbf{X}_i\),
and
\(-\mathbf{H}_{\mathbf{g}_2}=\sum_i^n\boldsymbol{\Delta}_i^{(2,2)}\).
The covariance matrix \(\mathbf{J}_\mathbf{g}\) of the composite score functions \(\mathbf{g}\) is given as below
$$
\mathbf{J}_\mathbf{g}=\mbox{Cov}(\mathbf{g})=
\left[\begin{array}{cc}
\mbox{Cov}(\mathbf{g}_1) & \mbox{Cov}(\mathbf{g}_1,\mathbf{g}_2)\\
\mbox{Cov}(\mathbf{g}_2,\mathbf{g}_1) & \mbox{Cov}(\mathbf{g}_2)
\end{array}\right]=
\left[\begin{array}{cc}\mathbf{J}_\mathbf{g}^{(1)}& \mathbf{J}_\mathbf{g}^{(1,2)}\\
\mathbf{J}_\mathbf{g}^{(2,1)}& \mathbf{J}_\mathbf{g}^{(2)}\end{array}\right]=
\sum_i\left[\begin{array}{cc}
\mathbf{X}_i^\top\boldsymbol{\Omega}_i^{(1)}
\mathbf{X}_i & \mathbf{X}_i^\top\boldsymbol{\Omega}_i^{(1,2)}\\
\boldsymbol{\Omega}_i^{(2,1)}\mathbf{X}_i& \boldsymbol{\Omega}_i^{(2)}\end{array}\right],
$$
where
$$\left[\begin{array}{cc}
\boldsymbol{\Omega}_i^{(1)}& \boldsymbol{\Omega}_i^{(1,2)}\\
\boldsymbol{\Omega}_i^{(2,1)}& \boldsymbol{\Omega}_i^{(2)}
\end{array}\right]=
\left[\begin{array}{cc}
\mbox{Cov}\Bigl(\mathbf{s}_i^{(1)}(\mathbf{a})\Bigr) &
\mbox{Cov}\Bigl(\mathbf{s}_i^{(1)}(\mathbf{a}),
\mathbf{s}_i^{(2)}(\mathbf{a},\mathbf{R})\Bigr)\\
\mbox{Cov}\Bigl(\mathbf{s}_i^{(2)}(\mathbf{a},\mathbf{R}),
\mathbf{s}_i^{(1)}(\mathbf{a})\Bigr) & \mbox{Cov}\Bigl(\mathbf{s}_i^{(2)}(\mathbf{a},\mathbf{R})\Bigr)
\end{array}\right].
$$
To this end, the composite AIC (Varin and Vidoni, 2005) and BIC (Gao and Song, 2011) criteria have the forms:
$$
\mbox{CL1AIC} = -2L_2 + 2\mbox{tr}\Bigl(\mathbf{J}_\mathbf{g}\mathbf{H}_\mathbf{g}^{-1}\Bigr),
$$
$$
\mbox{CL1BIC} =-2L_2+\log(n)\mbox{tr}\Bigl(\mathbf{J}_\mathbf{g}\mathbf{H}_\mathbf{g}^{-1}\Bigr).
$$