Pairwise likelihood inference for the multivariate ordered probit model

TL;DR

This paper derives a closed-form expression for the pairwise score vector in the multivariate ordered probit model, enhancing likelihood-based inference by improving optimization speed and standard error computation.

Contribution

It introduces a novel closed-form expression for the pairwise score vector, facilitating faster optimization and more accurate inference in multivariate ordered probit models.

Findings

Enables faster gradient-based optimization routines.

Provides a new method to compute standard errors and confidence intervals.

Improves likelihood-based inference accuracy.

Abstract

This paper provides a closed form expression for the pairwise score vector for the multivariate ordered probit model. This result has several implications in likelihood-based inference. It is indeed used both to speed-up gradient based optimization routines for point estimation, and to provide a building block to compute standard errors and confidence intervals by means of the Godambe matrix.

Figures

Equations

$$\displaystyle\Sigma=\left(\begin{array}[]{ccccc}1&\rho_{1,2}&\dots&\dots&\rho_{1,q}\\ &1&\dots&\dots&\rho_{2,q}\\ &&1&\rho_{r,s}&\vdots\\ &&&1&\rho_{q-1,q}\\ &&&&1\end{array}\right)$$

$$y_{ij}=k\text{ if and only if }{z}_{ij}\in (a_{k-1}(j),a_k(j)],$$

$$\text{pr}(Y_{i1}=y_{i1},\dots ,Y_{iq}=y_{iq})={\int }_{a_{y_{i1}-1}(1)}^{a_{y_{i1}}(1)}\cdots {\int }_{a_{y_{iq}-1}(q)}^{a_{y_{iq}}(q)}{\varphi }_{\Sigma }(z_1,\dots ,z_q)d{z}_1\dots d{z}_q,$$

$$\theta =({\rho }_{1,2},\dots ,{\rho }_{q-1,q},a_1(1),\dots ,a_{K-1}(1),a_1(2),\dots ,a_{K-1}(2),a_1(q),\dots ,a_{K-1}(q){)}^T,$$

$$\ell (\theta ;y)$$

$${\ell }^P(\theta ;y)=$$

$$=$$

$$\sqrt{n}(\hat{\theta }-\theta )\dot{\sim }N(0,G(\theta {)}^{-1}).$$

$$H(\theta )=\text{mathds}E\{-\partial u^P(\theta ;y)/\partial \theta \}\text{mbox}andJ(\theta )=\mathrm{Var}\{u^P(\theta ;y)\}$$

$$\displaystyle\hat{J}(\theta)=\frac{1}{n}\sum_{i=1}^{n}u^{P}(\theta;y_{i})u^{P}(\theta;y_{i})^{T}\Big{|}_{\theta=\hat{\theta}},$$

$$\displaystyle\hat{H}(\theta)=-\frac{1}{n}\sum_{i=1}^{n}\frac{\partial u^{P}(\theta;y_{i})}{\partial\theta}\Big{|}_{\theta=\hat{\theta}}.$$

$$\displaystyle\hat{H}(\theta)=\frac{1}{n}\sum_{i=1}^{n}\sum_{r=1}^{q-1}\sum_{s=r+1}^{q}u(\theta;y_{ir},y_{is})u(\theta;y_{ir},y_{is})^{T}\Big{|}_{\theta=\hat{\theta}},$$

$${\Phi }_2(a,b,\rho )={\int }_{-\infty }^a{\int }_{-\infty }^b{\varphi }_{\Sigma (\rho )}(u,v)dudv,$$

$${\int }_b^a{\int }_d^c{\varphi }_{\Sigma (\rho )}(u,v)dudv={\Phi }_2(a,c;\rho )-{\Phi }_2(b,c;\rho )-{\Phi }_2(a,d;\rho )+{\Phi }_2(b,d;\rho ).$$

$${\ell }^P(\theta ;y)$$

$$\displaystyle-\Phi_{2}\Bigr{(}a_{y_{ir}}(r),a_{y_{is}-1}(s),\rho_{r,s}\Bigl{)}+\Phi_{2}\Bigr{(}a_{y_{ir}-1}(r),a_{y_{is}-1}(s),\rho_{r,s}\Bigl{)}\biggr{\}}.$$

$$\frac{\partial {\Phi }_2(x_1,x_2,\rho )}{\partial x_1}$$

$$\frac{\partial {\Phi }_2(x_1,x_2,\rho )}{\partial x_2}$$

$$\frac{\partial {\Phi }_2(x_1,x_2,\rho )}{\partial \rho }$$

$$\frac{\partial {\ell }^P(\theta )}{\partial {\rho }_{r,s}}$$

$$\displaystyle\quad\times\biggr{[}\phi_{\Sigma(\rho_{r,s})}\Bigr{(}a_{y_{ir}}(r),a_{y_{is}}(s)\Bigl{)}-\ \phi_{\Sigma(\rho_{r,s})}\Bigr{(}a_{y_{ir}-1}(r),a_{y_{is}}(s)\Bigl{)}\ +$$

$$\displaystyle\quad-\phi_{\Sigma(\rho_{r,s})}\Bigr{(}a_{y_{ir}}(r),a_{y_{is}-1}(s)\Bigl{)}\ +\ \phi_{\Sigma(\rho_{r,s})}\Bigr{(}a_{y_{ir}-1}(r),a_{y_{is}-1}(s)\Bigl{)}\biggl{]}.$$

$$\mathcal{L}^{P}(\theta)=\prod_{l=1}^{K}\ \prod_{m=1}^{K}\ \prod_{r=1}^{q-1}\prod_{s=r+1}^{q}\mbox{pr}\Bigr{(}Y_{r}=l,Y_{s}=m\Bigl{)}^{n_{\{r,l\}\{s,m\}}},$$

$$n_{\{r,l\}\{s,m\}}=\sum _{i=1}^n{1}_{\{Y_{ir}=l,Y_{is}=m\}}.$$

$$\sum _{l=1}^K\sum _{m=1}^K{n}_{\{r,l\}\{s,m\}}=\frac{q(q-1)}{2}n.$$

$$\begin{split}\ell^{P}(\theta)&=\sum_{l=1}^{K}\sum_{m=1}^{K}\sum_{r=1}^{q-1}\sum_{s=r+1}^{q}n_{\{r,l\}\{s,m\}}\log\Bigl{\{}\mbox{pr}\Bigr{(}Y_{r}=l,Y_{s}=m\Bigl{)}\Bigr{\}}\\ &=\sum_{l=1}^{K}\sum_{m=1}^{K}\sum_{r=1}^{q-1}\sum_{s=r+1}^{q}n_{\{r,l\}\{s,m\}}\ \log\biggl{\{}\int_{a_{l-1}(r)}^{a_{l}(r)}\int_{a_{m-1}(s)}^{a_{m}(s)}\phi_{\Sigma(\rho_{r,s})}(z_{r},z_{s})\ dz_{r}dz_{s}\biggr{\}}\\ &=\sum_{l=1}^{K}\sum_{m=1}^{K}\sum_{r=1}^{q-1}\sum_{s=r+1}^{q}n_{\{r,l\}\{s,m\}}\ \text{log}\ \biggl{\{}\Phi_{2}\Bigr{(}a_{l}(r),a_{m}(s),\rho_{r,s}\Bigl{)}\ -\ \Phi_{2}\Bigr{(}a_{l-1}(r),a_{m}(s),\rho_{r,s}\Bigl{)}\ +\\ &\quad-\ \Phi_{2}\Bigr{(}a_{l}(r),a_{m-1}(s),\rho_{r,s}\Bigl{)}\ +\ \Phi_{2}\Bigr{(}a_{l-1}(r),a_{m-1}(s),\rho_{r,s}\Bigl{)}\biggr{\}}.\end{split}$$

$${\ell }^P(\theta ;y)=A_k(j)+B_k(j)+R_k(j),$$

$$\begin{split}A_{k}(j)=&\sum_{l=1}^{K}\sum_{m=k}^{k+1}\ \sum_{r=1}^{j-1}n_{\{r,l\}\{j,m\}}\ \text{log}\ \biggl{\{}\Phi_{2}\Bigr{(}a_{l}(r),a_{m}(j),\rho_{r,j}\Bigl{)}-\Phi_{2}\Bigr{(}a_{l-1}(r),a_{m}(j),\rho_{r,j}\Bigl{)}+\\ &-\ \Phi_{2}\Bigr{(}a_{l}(r),a_{m-1}(j),\rho_{r,j}\Bigl{)}\ +\Phi_{2}\Bigr{(}a_{l-1}(r),a_{m-1}(j),\rho_{r,j}\Bigl{)}\biggr{\}},\\ B_{k}(j)=&\sum_{l=k}^{k+1}\sum_{m=1}^{K}\sum_{s=j+1}^{q}n_{\{j,l\}\{s,m\}}\ \text{log}\ \biggl{\{}\Phi_{2}\Bigr{(}a_{l}(j),a_{m}(s),\rho_{j,s}\Bigl{)}-\Phi_{2}\Bigr{(}a_{l-1}(j),a_{m}(s),\rho_{j,s}\Bigl{)}+\\ &\quad-\Phi_{2}\Bigr{(}a_{l}(j),a_{m-1}(s),\rho_{j,s}\Bigl{)}+\Phi_{2}\Bigr{(}a_{l-1}(j),a_{m-1}(s),\rho_{j,s}\Bigl{)}\biggr{\}},\end{split}$$

$$\frac{\partial {\ell }^P(\theta )}{\partial a_k(j)}=\frac{\partial A_k(j)}{\partial a_k(j)}+\frac{\partial B_k(j)}{\partial a_k(j)}$$

$$\begin{split}\frac{\partial A_{k}(j)}{\partial a_{k}(j)}&=\phi\Bigr{(}a_{k}(j)\Bigl{)}\sum_{l=1}^{K}\sum_{r=1}^{j-1}\left\{\Phi\Biggr{(}\frac{a_{l}(r)-\rho_{r,j}a_{k}(j)}{\sqrt{1-\rho_{r,j}^{2}}}\Biggl{)}-\Phi\Biggr{(}\frac{a_{l-1}(r)-\rho_{r,j}a_{k}(j)}{\sqrt{1-\rho_{r,j}^{2}}}\Biggl{)}\right\}\times\\ &\quad\times\Biggl{[}\ n_{\{r,l\}\{j,k\}}\biggl{\{}\int_{a_{l-1}(r)}^{a_{l}(r)}\int_{a_{k-1}(j)}^{a_{k}(j)}\phi_{\Sigma(\rho_{r,j})}(z_{r},z_{j})\ dz_{r}dz_{j}\biggr{\}}^{-1}+\\ &\quad\quad\quad-\ n_{\{r,l\}\{j,k+1\}}\biggl{\{}\int_{a_{l-1}(r)}^{a_{l}(r)}\int_{a_{k}(j)}^{a_{k+1}(j)}\phi_{\Sigma(\rho_{r,j})}(z_{r},z_{j})\ dz_{r}dz_{j}\biggr{\}}^{-1}\Biggl{]},\end{split}$$

$$\begin{split}\frac{\partial B_{k}(j)}{\partial a_{k}(j)}&=\phi\Bigr{(}a_{k}(j)\Bigl{)}\sum_{m=1}^{K}\sum_{s=j+1}^{q}\Biggl{[}\Phi\Biggr{(}\frac{a_{m}(s)-\rho_{j,s}a_{k}(j)}{\sqrt{1-\rho_{j,s}^{2}}}\Biggl{)}-\Phi\Biggr{(}\frac{a_{m-1}(s)-\rho_{j,s}a_{k}(j)}{\sqrt{1-\rho_{j,s}^{2}}}\Biggl{)}\Biggr{]}\\ &\quad\times\Biggl{[}n_{\{j,k\}\{s,m\}}\ \biggl{\{}\int_{a_{k-1}(j)}^{a_{k}(j)}\int_{a_{m-1}(s)}^{a_{m}(s)}\phi_{\Sigma(\rho_{j,s})}(z_{j},z_{s})\ dz_{j}dz_{s}\biggr{\}}^{-1}\ +\\ &\quad\quad\quad-\ n_{\{j,k+1\}\{s,m\}}\biggl{\{}\int_{a_{k}(j)}^{a_{k+1}(j)}\int_{a_{m-1}(s)}^{a_{m}(s)}\phi_{\Sigma(\rho_{j,s})}(z_{j},z_{s})\ dz_{j}dz_{s}\biggr{\}}^{-1}\Biggl{]}.\end{split}$$

$$\frac{\partial {\ell }^P(\theta )}{\partial a_k(1)}=\frac{\partial B_k(1)}{\partial a_k(1)},\frac{\partial {\ell }^P(\theta )}{\partial a_k(q)}=\frac{\partial A_k(q)}{\partial a_k(q)}.$$

$${\delta }_k(j)=\log (a_k(j)-a_{k-1}(j)),$$

Full text

**Keywords: ** composite likelihood; Godambe information matrix; pairwise score function

1 Introduction

Multivariate ordinal categorical outcomes (Agresti, 2002) arise in many fields of application such as political science, marketing research, educational assessment, and psychometry among others. A common approach consists in considering that the $q$ observed categorical data are related to continuous latent variables. The relation between the latent continuous variables and the observed categorical variables is usually induced by means of thresholds partitioning the latent sample space into a series of regions corresponding to each ordinal category. Popular distributions for the latent variables are the logistic and the Gaussian distributions, leading to the ordered logit and probit models, respectively.

In this paper, we deal with inferential issues related to the multivariate ordered probit class of model following a likelihood-based approach. Specifically, since the multivariate ordered probit model has clear computational problems related to the calculation of a $q$ dimensional integral for each single likelihood contribution, we exploit an approximate approach based on a simple likelihood belonging to the class of composite likelihoods (Lindsay, 1988; Varin et al., 2011): the pairwise likelihood (Cox and Reid, 2004). The use of the pairwise likelihood for the ordered probit model has been shown to have clear advantages in many situations (De Leon, 2005; Varin and Vidoni, 2006; Kenne Pagui and Canale, 2016; Hirk et al., 2018). In particular, in the context of point estimation, Kenne Pagui and Canale (2016) showed a dramatic improvement in the computational time and only moderate bias if compared to a standard likelihood approach. This comparison, however, was only possible for moderate values of $q$ as the standard likelihood approach is not practically feasible in high dimensions. Hirk et al. (2018) confirm this evidence while showing that a triplewise likelihood does not lead to substantial improvement further endorsing the pairwise approximation. In both of the aforementioned papers, the authors try to go beyond point estimation and provide significance testing on some parameters. There is however a clear difficulty in obtaining inferential results beyond point estimation with the approach at hand: an analytic expression for the covariance matrix of the pairwise likelihood estimator is not available. This is mainly due to the difficulties in calculating the Godambe matrix, an estimator for the asymptotic precision of the maximum pairwise likelihood estimator. Indeed, the latter involves the calculation of the derivative of the pairwise log-likelihood, i.e. the derivative of the logs of bivariate integrals of Gaussian densities with respect to the limits of integration and the latent correlation coefficients. To bypass this problem, possible solutions consist in calculating the the derivative of the pairwise log-likelihood numerically or in avoiding asymptotic arguments and instead using a bootstrap-based approach (Kenne Pagui and Canale, 2016). Unfortunately both approaches are computationally demanding in moderate and high dimensions. While Hirk et al. (2018) mention that the estimates standard errors may be obtained “either analytically or by numerical differentiation” (see page 8 of Hirk et al. (2018)) in their paper they do not provide any analytical result.

The analytical expression of the pairwise score vector would be hence of paramount interest both to speed-up gradient based optimization routines for point estimation, and in order to compute both the sensitivity and variability matrices involved in the Godambe’s matrix computation. This paper fills this gap and presents the analytical expression of the pairwise score vector of the multivariate ordered probit model.

The paper is organized as follows: the next section introduces the multivariate ordered probit model. Section 3 —after a brief review of the pairwise likelihood approach to inference— reports the main result of this paper, i.e. the closed form expression for the pairwise score vector. In Section 3.3 some numerical and computational considerations are reported. Section 4 reports a simulation experiment to assess the performance of the pairwise approach to perform confidence interval estimation by means of an empirical approximation of the Godambe’s matrix. Section 5 concludes the paper discussing some ongoing research.

2 The multivariate ordered probit model

Let $Y_{i}=(Y_{i1},\ldots,Y_{iq})^{{\rm T}}$, $Y_{ij}\in\{1,2,\ldots,K\}$ for $j=1,\dots,q$ be a $q$-dimensional ordinal categorical random vector with joint distribution depending on some unknown parameter $\theta$, with $i=1,\dots,n$ defining a collection of iid random vectors. The ordered probit model assumes that, for each $i=1,\dots,n$, there exists a latent random vector ${Z_{i}}=(Z_{i1},\ldots,Z_{iq})^{{\rm T}}$, with ${Z_{i}}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\Sigma)$ where

[TABLE]

is an unknown $q$-dimensional positive definite correlation matrix. In observing $y_{ij}$, the following relation is assumed

[TABLE]

where for every $j=1,\dots,q$, $\{a_{k}(j)\}_{k=0,\dots,K}$ is a sequence of real numbers such that $-\infty=a_{0}(j)<a_{1}(j)<\dots<a_{K-1}(j)<a_{K}(j)=\infty,$ defining a disjoint partition of $\mathbb{R}$. The real vector $a(j)=\{a_{1}(j),\dots,a_{K-1}(j)\}$ is then a $(K-1)$-dimensional vector of threshold parameters referred to the $j$-th marginal.

Under the ordered probit model assumption, the likelihood contribution of a single observation is then proportional to the joint probability

[TABLE]

where $\phi_{\Sigma}(\cdot)$ denotes the multivariate Gaussian distribution with zero mean and variance $\Sigma$. Let $\theta$ be the joint vector of parameters, i.e.

[TABLE]

then, the log-likelihood function is given by

[TABLE]

where $y$ denotes all the observed sample. The total number of model parameters is given by the sum of the number of thresholds ($(K-1)\times q$) plus number of latent correlation coefficients ($q(q-1)/2$). As already discussed, equation (1) involves a $q$-dimensional integral for each observation and hence it is not easily manageable for moderate and high values of $q$.

3 Pairwise inference for the multivariate probit model

3.1 Pairwise likelihood inference

The pairwise likelihood is a likelihood constructed from bivariate marginals. For our multivariate ordered categorical data, the pairwise log-likelihood is

[TABLE]

where $\Sigma(\rho)$ denotes the $2\times 2$ correlation matrix with off-diagonal entries equal to $\rho$. As already discussed, the pairwise approach is particularly appealing because it substitutes the computational challenges related to a $q$ dimensional numerical integration with simpler bivariate integrals.

Generally, composite likelihood inferential procedures shares many properties of standard likelihood methods. For example, the pairwise score vector, $u^{P}(\theta;y)=\partial\ell^{P}(\theta;y)/\partial\theta,$ is still unbiased being the sum of score vectors based on the likelihood contribution of each pair of observations. Under regularity conditions, the maximum pairwise likelihood estimator $\hat{\theta}$ —obtained either by maximizing $\ell^{P}(\theta;y)$ numerically or by solving the pairwise likelihood equations $u^{P}(\theta;y)=0$— is consistent and has asymptotic Gaussian distribution

[TABLE]

In equation above $G(\theta)=H(\theta)J(\theta)^{-1}H(\theta)$ is the Godambe information matrix where

[TABLE]

are the sensitivity and variability matrices, respectively. The sample estimate of the variability matrix is simply

[TABLE]

while the sample estimate of the sensitivity matrix is given by

[TABLE]

However, the calculation of the matrix of the derivatives of the score can be avoided by exploiting the second Bartlett identity, which is still valid as the pairwise likelihood is made of proper likelihood contributions. This yields the alternative expression

[TABLE]

where $u(\theta;y_{ir},y_{is})$ is a suitable adaptation of the score vector of a bivariate Gaussian likelihood with data $(y_{ir},y_{is})^{T}$. This adaptation is needed to match the dimension of the full vector of parameters $\theta$, i.e. $q(q-1)/2+(K-1)q$ as each bivariate Gaussian likelihood involves $1+2(K-1)$ parameters only. Specifically, if a component of $\theta$ is not present in the specific bivariate Gaussian likelihood, a zero is assigned to the related entry of $u(\theta;y_{ir},y_{is})$.

3.2 The pairwise score vector

In this section we are going to obtain closed-form expressions for each element of the pairwise score vector for the pseudolikelihood (2). To our knowledge, this is the first time that such a result is analytically available. In addition, this result is of paramount interest both to speed-up gradient based optimization routines and, even more, to provide sample estimates for the sensitivity and variability matrices used to obtain the asymptotic variance of the pairwise likelihood estimator. The derivation of the following expressions are carefully described in the Master thesis of Bravo (2017) defended at the University of Turin, Italy.

Using linearity in differentiation and the property of the derivation of a composite function, the crucial point related to the calculation of the pairwise score vector is the computation of the derivative of the argument of the logarithm in (2). To this end, consider the standard bivariate Gaussian cumulative distribution function (CDF)

[TABLE]

and note that the general integral of a bivariate Gaussian density can be decomposed in the sum of four CDF, namely

[TABLE]

Hence the arguments of all the logarithms in (2) can be easily rewritten according to equation (5) thus allowing to rewrite the pairwise log-likelihood function as

[TABLE]

It is clear that the partial derivatives of the score vector, follow from the partial derivatives of bivariate Gaussian CDFs. To obtain the latter, consider the following results (Drum and McCullagh, 1993).

Lemma 1

The partial derivatives of the bivariate standard Normal CDF with respect to the limits of integration and the correlation coefficients are

[TABLE]

The application of Lemma 1 allows us to easily compute the partial derivatives with respect to all correlation coefficients, i.e. for every $r=1,\dots,q-1$ and $s=r+1,\dots,q$, we have

[TABLE]

The latter equation contributes to determine the first entries of the score vector $u^{P}(\theta;y)=\partial\ell^{P}(\theta;y)/\partial\theta$.

In computing the partial derivatives with respect to the thresholds parameters, the application of Lemma 1 is not trivial. First of all, note that in (2) the extremes of integration still depend on the observations $y_{ij}$. Hence, it is useful to rewrite (2), grouping the pairs of observations sharing the same values $y_{ij}$. To this end, we can first rewrite the pairwise likelihood (1) as

[TABLE]

where $n_{\{r,l\}\{s,m\}}$ is the number of observations with $Y_{ir}=l$ and $Y_{is}=m$, i.e.

[TABLE]

Clearly the $n_{\{r,l\}\{s,m\}}$ multiplicities are such that

[TABLE]

The full log-likelihood in (2) can be rewritten as

[TABLE]

We now calculate the partial derivative of the pairwise log-likelihood with respect to $a_{k}(j)$. This parameter is referred to the $k$-th response level—and hence it is involved when $l=k,k+1$ and $m=k,k+1$—and to the $j$-th marginal dimension. To emphasise the contribution of $a_{k}(j)$ to the pairwise log-likelihood above, we can rewrite $\ell^{P}(\theta;y)$ as

[TABLE]

where

[TABLE]

and $R_{k}(j)$ is a residual part that does not depend on $a_{k}(j)$. Note that $A_{k}(j)$ is the likelihood contribution where the dimension $j$ is fixed as second element in the pairs with levels $k$ and $k+1$, while the first element is free to vary between all the possible previous dimensions and combinations of levels. Similarly, $B_{k}(j)$ is the likelihood contribution where the dimension $j$ is fixed as first element in the pairs with levels $k$ and $k+1$ and the second is free to vary between all the possible next combinations of levels.

The partial derivative of the pairwise log-likelihood with respect to $a_{k}(j)$ is hence,

[TABLE]

where

[TABLE]

and

[TABLE]

Note that for $j=1$ and $j=q$, equations above simplify to

[TABLE]

3.3 Numerical and computational considerations

The analytical expression of the pairwise score vector equals its numerical equivalent—computed through the grad function in the R library numDeriv. However, its computational burden is dramatically lower. Figure 1 compares the execution times of a direct calculation of the gradient of the pairwise log-likelihood with the analytical expression of the score for different values of $q$ and fixed sample size and number of levels ($n=50$ and $K=5$). The improvement of our result is dramatic. For instance when $q=12$, the user time to numerically evaluate the score is more that 60 times higher than its analytic counterpart. Similar results are obtained comparing the approaches for growing values of $n$ —and fixed $q=5$ and $K=5$. The computational costs behave linearly in both cases, but at a slower rate for the routine using the analytical expression.

As discussed in Kenne Pagui and Canale (2016) and Varin and Vidoni (2006), the optimization of the pairwise likelihood function, to compute point estimates of the parameters, is performed through the quasi-Newton BFGS algorithm. To assure the ordering of the thresholds $a_{0}(j)<\dots<a_{K}(j)$, within the optimization routine, it is useful to define for $j=1,\dots,q$ and $k=2,\dots,K-1$, the new parameters

[TABLE]

and express the pairwise log-likelihood as a function of these new parameters. In order to speed up the optimization procedure, we can provide the analytical expression of the score for the model under the new parameter

[TABLE]

identified by the bijective map $\psi=g(\theta)$. Here $g(\cdot)$ is a vector-to-vector function mapping each correlation coefficient and each $a_{1}(j)$ with the identity function, and all the remaining threshold parameters with (14). The pairwise score vector for the new parameter is simply

[TABLE]

where the matrix of partial derivatives of the inverse of $g(\cdot)$ with respect to $\psi$ is the following block matrix

[TABLE]

where $\mbox{I}_{q(q-1)/2}$ is the identity matrix of dimension ${q(q-1)/2}$, and each $\Delta_{j}$ is the submatrix of partial derivatives of $g^{-1}(\cdot)$ with respect to the new parameters $\delta_{1}(j),\dots,\delta_{K-1}(j)$, i.e.

[TABLE]

Outside the block diagonal, each entry of the matrix of partial derivatives of the inverse of $g(\cdot)$ with respect to $\psi$ is zero.

Figure 2 compares the execution times of the numerical optimization of the pairwise log-likelihood —performed with the nlminb function of R— providing the analytical expression of the score with the numerical optimization in which the gradient is obtained numerically. Data are simulated with different values of dimensions $q$, with $n=50$ and $K=5$ (panel a) and with different sample sizes $n$, fixing $q=3$ and $K=5$ (panel b). In both situations the improvement is considerable. For instance with $q=9$, the optimization with the analytical expression of the gradient performs in less than 3 minutes versus about 20 minutes for the standard procedure. The use of the algebrical expression obtained in this section dramatically lowers the single pseudolikelihood optimization.

4 Simulation study

In this section we describe a simulation study conducted in order to assess the performance of the proposed approach in terms of point and confidence interval estimation. For confidence interval estimation, we consider the Wald-type confidence interval for all the $q(q-1)/2$ polychoric correlations coefficients. These are defined by the subset of $[-1,1]$ with extremes

[TABLE]

where $z_{1-\alpha/2}$ is the quantile of a standard Gaussian distribution of level $1-\alpha/2$, $\hat{\rho}_{r,s}$ is the pairwise maximum likelihood estimates for ${\rho}_{r,s}$ and $\mbox{se}(\hat{\rho}_{r,s})$ it the appropriate element of the diagonal of the inverse of the Godambe’s matrix.

In fact, the pairwise approach for the ordered probit model has already proved to provide satisfying results in terms of bias and variance for point estimation Kenne Pagui and Canale (2016) Hirk et al. (2018) but reporting details on point estimation provides additional insights on the performance in terms of confidence interval estimation.

We compute point and interval estimations for $R=100$ replicated and independent datasets. Then, for each parameter we compute the mean squared error of point prediction, averaging over the $R$ replicates, the average standard errors $\mbox{se}(\hat{\rho}_{r,s})$ and the empirical coverage of the procedure counting how many times over $R$ the true $\rho_{r,s}$ is inside the obtained confidence interval. The confidence level is fixed to $\alpha=0.05$.

In evaluating the performance of the proposed solution in terms of confidence interval estimation, it must be taken into account that we are considering the empirical estimates of the sensitivity and variability matrices provided in equations (3)–(4). However, such empirical estimates may be not very accurate as discussed by Cattelan and Sartori (2016). A workaround solution proposed by the latter authors consists in estimating the sensitivity and variability matrices via Monte Carlo simulation. This solution, despite providing more reliable estimates of the Godambe’s matrix, comes at an increasing computational cost and its implementation must be considered on a case-by-case basis.

To assess the inferential performance in different scenarios, we simulate data sets of different dimension $q$. Specifically we simulated data with $q=10,15$. As discussed, the number of parameter is $O(q^{2})$ and hence in order to have reasonable estimates of the parameter, we need more data points for increasing dimension $q$. The sample sizes are $n=300,400,500$ and $n=400,600,800$ for $q=10$ and for $q=15$, respectively. The latent polychoric correlation matrices are correlation matrices randomly generated assuming a sparse structure with 30% of zeroes. We fix $K=4$ and define two set of thresholds, $(0,0.5,1)^{T}$ and $(-1,0,1)^{T}$ randomly assigned to each marginal variable.

The results of the simulation study are reported in Figure 3. As expected, and consistently with the findings of Kenne Pagui and Canale (2016) and Hirk et al. (2018) the mean squared errors—reported in panels (a)—are low in absolute value and substantially decrease for increasing sample size in both scenarios. The average magnitude of the standard errors also decreases for increasing sample size as can be noticed in panels (b). This is also an expected behaviour as the precision of the estimators is expected to grow with the sample size. The empirical coverage, depicted in the last panels of Figure 3 is close to its nominal values but shows, in all cases, a moderate over confidence. This is probably due to the already discussed inaccuracy of the sample estimates of the variability and sensitivity matrices which may lead to larger standard errors and hence wider intervals. Despite being conservative and showing a considerable variability for each parameter, globally the empirical coverage converges to its nominal value as the sample size increases, i.e. when the averages in (3)–(4) are calculated with more data points.

5 Discussion

In this article, we studied a pairwise likelihood approach for inference in the multivariate ordered probit model. Specifically, we derived and described the analytical expression of the pairwise score vector. This result is of paramount interest both for point and interval estimation. Indeed the analytical expression led to a dramatic reduction of the computational costs related to the numerical maximization of the pairwise log-likelihood function via standard gradient-based numerical optimization as discussed in Section 3.3.

In addition, the pairwise score allowed us to compute an empirical estimate of the Godambe matrix that we used for standard errors quantification and confidence interval estimations. Despite showing promising performance in the simulation study, the empirical estimates of the sensitivity and variability matrices provided in equations (3)–(4) must be used with some care. Indeed, for small sample sizes, they may be very inaccurate and a Monte Carlo estimator may be preferable Cattelan and Sartori (2016). This solution comes at an increasing computational cost which further motivates the need of the fast and reliable routine to compute the pairwise score vector based on our analytical results.

In the specific class of models that we considered here, all the $q$ marginals have the same number of categories $K$. However, extensions to the cases in which each marginal can have different levels of categories $K_{j}$ are straightforward. We also assumed that the latent variables are iid but in many situations it is reasonable to assume some sort of dependence, e.g. from a set of known covariates (Kenne Pagui and Canale, 2016; Hirk et al., 2018). This aspect is of dramatic interest in many application and is subject to ongoing research.

Acknowledgement

This work was initially developed by Martina Bravo as a MSc thesis in Stochastics and Data Science at the University of Torino, Italy.