Model Selection in a Composite Likelihood Framework Based on Density Power Divergence

Abstract

This paper presents a model selection criterion in a composite likelihood framework based on density power divergence measures and in the composite minimum density power divergence estimators, which depends on an tuning parameter $\alpha$. After introducing such a criterion, some asymptotic properties are established. We present a simulation study and two numerical examples in order to point out the robustness properties of the introduced model selection criterion.

1. Introduction

Composite likelihood inference is an important approach to deal with those real situations of large data sets or very complex models, in which classical likelihood methods are computationally difficult, or even, not possible to manage. Composite likelihood methods have been successfully used in many applications concerning, for example, genetics ([1]), generalized linear mixed models ([2]), spatial statistics ([3,4,5]), frailty models ([6]), multivariate survival analysis ([7,8]), etc.

Let us introduce the problem, adopting here the notation by [9]. Let $\{f(·;\mathit{\theta }),\mathit{\theta }\in \mathsf{\Theta }\subseteq {\mathbb{R}}^p,p\ge 1\}$ be a parametric identifiable family of distributions for an observation $\mathit{y}={(y_1,\dots ,y_m)}^T$, a realization of a random m-vector $\mathit{Y}$. In this setting, the composite likelihood function based on K different marginal or conditional distributions has the form

$$\mathcal{CL}(\mathit{\theta },\mathit{y})=\prod _{k=1}^K{\left(f_{A_k}(y_j,j\in A_k;\mathit{\theta })\right)}^{w_k}$$

and the corresponding composite log-density

$$\log \mathcal{CL}(\mathit{\theta },\mathit{y})=\sum _{k=1}^K{w}_k{\ell }_{A_k}(\mathit{\theta },\mathit{y}),$$

(1)

with ${\ell }_{A_k}(\mathit{\theta },\mathit{y})=log{f}_{A_k}(y_j,j\in A_k;\mathit{\theta })$, where ${\left\{A_k\right\}}_{k=1}^K$ is a family of sets of indices associated either with marginal or conditional distributions involving some $y_j$, $j\in \{1,\dots ,m\}$ and $w_k$, $k=1,\dots ,K$ are non-negative and known weights. If the weights are all equal, then they can be ignored. In this case, all the statistical procedures give equivalent results. The composite maximum likelihood estimator (CMLE), ${\widehat{\mathit{\theta }}}_c$, is obtained by maximizing, in respect to $\mathit{\theta }\in \mathsf{\Theta }$, the expression (1).

The CMLE is consistent and asymptotically normal and, based on it, we can establish hypothesis testing procedures in a similar way to the classical likelihood ratio test, Wald test or Rao’s score test. A development of the asymptotic theory of the CMLE including its application to obtain the composite ratio statistics, Wald-type tests and Rao score tests in the context of composite likelihood can be seen in [10]. However, in [11,12,13] is shown that the CMLE and the derived testing procedures present an important lack of robustness. In this sense, [11,12,13] derived some new distance-based estimators and tests with good robustness behaviour without an important loss of efficiency. In this paper, we are going to consider the composite minimum density power divergence estimator (CMDPDE), introduced in [12], in order to present a model selection criterion in a composite likelihood framework.

Model selection criteria, for summarizing data evidence in favor of a model, is a very well studied subject in statistical literature, overall in the context of full likelihood. The construction of such criteria requires a measure of similarity between two models, which are typically described in terms of their distributions. This can be achieved if an unbiased estimator of the expected overall discrepancy is found, which measures the statistical distance between the true, but unknown model, and the entertained model. Therefore, the model with the smallest value of the criterion is the most preferable model. The use of divergence measures, in particular Kullback–Leibler divergence ([14]), to measure this discrepancy, is the main idea of some of the most known criteria: Akaike Information Criterion (AIC, [15,16]), the criterion proposed by Takeuchi (TIC, [17]) and other modifications of AIC [18]. DIC criterion, based on the density power divergence (DPD), was presented in [19] and, recently, [20] presented a local BHHJ power divergence information criterion following [21]. In the context of the composite likelihood there are some criteria based on Kullback–Leibler divergence, see for instance [22,23,24] and references therein. To the best of our knowledge only Kullback–Leibler divergence was used to develop model selection criteria in a composite likelihood framework. To fill this gap, our interest is now focused on DPD.

In this paper, we present a new information criterion for model selection in the framework of composite likelihood based on DPD measure. This divergence measure, introduced and studied in the case of complete likelihood by [25], has been considered previously in [12,13] in the context of composite likelihood. In these papers, a new estimator, the CMDPDE, was introduced and its robustness in relation to the CMLE as well as the robustness of some families of test statistics were studied, but the problem of model selection was not considered. This problem is considered in this paper. The criterion introduced in this paper will be called composite likelihood DIC criterion (CLDIC). The motivation of considering a criterion based on DPD instead of Kullback–Leibler divergence is due to the robustness of the procedures based on DPD in statistical inference, not only in the context of full likelihood [25,26], but also in the context of composite likelihood [12,13]. In Section 2, the CMDPDE is presented and some properties of this estimator are discussed. The new model selection criterion, CLDIC, based on CMDPDE is introduced in Section 3 and some of its asymptotic properties are studied. A simulation study is carried out in Section 4 and some numerical examples are presented in Section 5. Finally, some concluding remarks are presented in Section 6.

2. Composite Minimum Density Power Divergence Estimator

Given two probability density functions g and f, associated with two m-dimensional random variables respectively, the DPD ([25]) measures a statistical distance between g and f by

$$\begin{array}{c}\hfill d_{\alpha }(g,f)={\int }_{{\mathbb{R}}^m}\left\{f{\left(\mathit{y}\right)}^{1+\alpha }-\left(1+\frac{1}{\alpha }\right)f{\left(\mathit{y}\right)}^{\alpha }g\left(\mathit{y}\right)+\frac{1}{\alpha }g{\left(\mathit{y}\right)}^{1+\alpha }\right\}d\mathit{y},\end{array}$$

(2)

for $\alpha >0,$ while for $\alpha =0$ it is defined by

$$d_0(g,f)=\underset{\alpha \to 0^{+}}{lim}{d}_{\alpha }(g,f)=d_{KL}(g,f),$$

where $d_{KL}(g,f)$ is the Kullback–Leibler divergence (see, for example, [26]). For $\alpha =1$, the expression (2) leads to the $L_2$ distance $L_2(g,f)={\int }_{{\mathbb{R}}^m}{\left(f\left(\mathit{y}\right)-g\left(\mathit{y}\right)\right)}^{2}d\mathit{y}.$ It is also interesting to note that (2) is a special case of the so-called Bregman divergence

$${\int }_{{\mathbb{R}}^m}\left[T\left(g\left(\mathit{y}\right)\right)-T\left(f\left(\mathit{y}\right)\right)-\{g\left(\mathit{y}\right)-f\left(\mathit{y}\right\}{T}^{\prime }\left(f\left(\mathit{y}\right)\right)\right]d\mathit{y}.$$

(3)

If we consider $T\left(l\right)=\frac{1}{\alpha }{l}^{1+\alpha }$ in (3), we get $d_{\alpha }(g,f)$. The parameter $\alpha$ controls the trade-off between robustness and asymptotic efficiency of the parameter estimates which are the minimizers of this family of divergences. For more details about this family of divergence measures we refer to [27].

Let now ${\mathit{Y}}_1,\dots ,{\mathit{Y}}_n$ be independent and identically distributed replications of $\mathit{Y}$ which are characterized by the true but unknown distribution g. Taking into account that the true model g is unknown, suppose that $\mathsf{\Xi }=\{f(·;\mathit{\theta }),\mathit{\theta }\in \mathsf{\Theta }\subseteq {\mathbb{R}}^p,p\ge 1\}$ is a parametric identifiable family of candidate distributions to describe the observations ${\mathit{y}}_1,\dots ,{\mathit{y}}_n$. Then, the DPD between the true model g and the composite likelihood function, $\mathcal{CL}(\mathit{\theta },·),$ associated to the parametric model $f(·;\mathit{\theta })$ is defined as

$$\begin{array}{cc}\hfill d_{\alpha }(g\left(·\right),\mathcal{CL}(\mathit{\theta },·))=& {\int }_{{\mathbb{R}}^m}\left\{\mathcal{CL}{(\mathit{\theta },\mathit{y})}^{1+\alpha }-\left(1+\frac{1}{\alpha }\right)\mathcal{CL}{(\mathit{\theta },\mathit{y})}^{\alpha }g\left(\mathit{y}\right)+\frac{1}{\alpha }g{\left(\mathit{y}\right)}^{1+\alpha }\right\}d\mathit{y},\hfill \end{array}$$

(4)

for $\alpha >0$, while for $\alpha =0$ we have $d_{KL}(g\left(·\right),\mathcal{CL}(\mathit{\theta },·))$, which is defined by

$$d_{KL}(g\left(·\right),\mathcal{CL}(\mathit{\theta },·))={\int }_{{\mathbb{R}}^m}g\left(\mathit{y}\right)log\frac{g\left(\mathit{y}\right)}{\mathcal{CL}(\mathit{\theta },\mathit{y})}d\mathit{y}.$$

(5)

In Section 3, we are going to introduce and study the CLDIC criterion based on (4).

Let

$${\left\{M_k\right\}}_{k\in \left\{1,\dots ,\ell \right\}}$$

(6)

be a family of candidate models to govern the observations ${\mathit{Y}}_1,\dots ,{\mathit{Y}}_n$. We shall assume that the true model is included in ${\left\{M_k\right\}}_{k\in \left\{1,\dots ,\ell \right\}}.$ For a specific $k=1,\dots ,\ell$, the parametric model $M_k$ is described by the composite likelihood function

$$\mathcal{CL}(\mathit{\theta },·),\mathit{\theta }\in {\mathsf{\Theta }}_k\subset {\mathbb{R}}^k.$$

In this setting, it is quite clear that the most suitable candidate model to describe the observations is the model that minimizes the DPD in (4). However, the unknown parameter $\mathit{\theta }$ is included in it, so it is not possible to use directly this measure for the choice of the most suitable model. A way to overcome this problem is to plug-in, in (4), the unknown parameter $\mathit{\theta }$ by an estimator which is desirable to obey some nice properties, like consistency and asymptotic normality. Based on this point, the CMDPDE, introduced in [12], can be used. This estimator is described in the sequel for the sake of completeness.

If we denote the kernel of (4) as

$$W_{\alpha }\left(\mathit{\theta }\right)={\int }_{{\mathbb{R}}^m}\mathcal{CL}{(\mathit{\theta },\mathit{y})}^{1+\alpha }d\mathit{y}-\left(1+\frac{1}{\alpha }\right){\int }_{{\mathbb{R}}^m}\mathcal{CL}{(\mathit{\theta },\mathit{y})}^{\alpha }g\left(\mathit{y}\right)d\mathit{y},$$

(7)

we can write

$$d_{\alpha }(g\left(·\right),\mathcal{CL}(\mathit{\theta },·))=W_{\alpha }\left(\mathit{\theta }\right)+\frac{1}{\alpha }{\int }_{{\mathbb{R}}^m}g{\left(\mathit{y}\right)}^{1+\alpha }d\mathit{y}$$

and the term $\frac{1}{\alpha }{\int }_{{\mathbb{R}}^m}g{\left(\mathit{y}\right)}^{1+\alpha }d\mathit{y}$ does not depend on $\mathit{\theta }$ and could be ignored in (9). A natural estimator of $W_{\alpha }\left(\mathit{\theta }\right)$, given in (7), can be obtained by observing that the last integral in (7), can be expressed in the form ${\int }_{{\mathbb{R}}^m}\mathcal{CL}{(\mathit{\theta },\mathit{y})}^{\alpha }dG\left(\mathit{y}\right)$, for G the distribution function corresponding to g. Hence, if the empirical distribution function of ${\mathit{Y}}_1,\dots ,{\mathit{Y}}_n$ will be exploited, this last integral is approximated by $\frac{1}{n}{\displaystyle \sum _{i=1}^n}\mathcal{CL}{(\mathit{\theta },{\mathit{Y}}_i)}^{\alpha }$, i.e.,

$$W_{n,\alpha }\left(\mathit{\theta }\right)={\int }_{{\mathbb{R}}^m}\mathcal{CL}{(\mathit{\theta },\mathit{y})}^{\alpha +1}d\mathit{y}-\left(1+\frac{1}{\alpha }\right)\frac{1}{n}\sum _{i=1}^n\mathcal{CL}{(\mathit{\theta },{\mathit{Y}}_i)}^{\alpha }.$$

(8)

Definition 1.

The CMDPDE of θ, ${\widehat{\mathit{\theta }}}_c^{\alpha },$ is defined, for $\alpha >0$, by

$${\widehat{\mathit{\theta }}}_c^{\alpha }=arg\underset{\mathit{\theta }\in \mathsf{\Theta }}{min}{W}_{n,\alpha }\left(\mathit{\theta }\right).$$

(9)

We shall denote the score of the composite likelihood by

$$\mathit{u}(\mathit{\theta },\mathit{y})=\frac{\partial log\mathcal{CL}(\mathit{\theta },\mathit{y})}{\partial \mathit{\theta }}.$$

(10)

Let ${\mathit{\theta }}_0$ be the true value of the parameter $\mathit{\theta }$. In [12], it was shown that the asymptotic distribution of ${\widehat{\mathit{\theta }}}_c^{\alpha }$ is given by

$$\sqrt{n}({\widehat{\mathit{\theta }}}_c^{\alpha }-{\mathit{\theta }}_0)\underset{n\to \infty }{\stackrel{\mathcal{L}}{⟶}}\mathcal{N}\left({\mathbf{0}}_p,{\mathit{H}}_{\alpha }{\left({\mathit{\theta }}_0\right)}^{-1}{\mathit{J}}_{\alpha }\left({\mathit{\theta }}_0\right){\mathit{H}}_{\alpha }{\left({\mathit{\theta }}_0\right)}^{-1}\right),$$

being

$${\mathit{H}}_{\alpha }\left(\mathit{\theta }\right)={\int }_{{\mathbb{R}}^m}\mathcal{CL}{(\mathit{\theta },\mathit{y})}^{\alpha +1}\mathit{u}(\mathit{\theta },\mathit{y})\mathit{u}{(\mathit{\theta },\mathit{y})}^{T}d\mathit{y}$$

(11)

and

$$\begin{array}{cc}\hfill {\mathit{J}}_{\alpha }\left(\mathit{\theta }\right)& ={\int }_{{\mathbb{R}}^m}\mathcal{CL}{(\mathit{\theta },\mathit{y})}^{2\alpha +1}\mathit{u}(\mathit{\theta },\mathit{y})\mathit{u}{(\mathit{\theta },\mathit{y})}^{T}d\mathit{y}\hfill \\ \hfill & -{\int }_{{\mathbb{R}}^m}\mathcal{CL}{(\mathit{\theta },\mathit{y})}^{\alpha +1}\mathit{u}(\mathit{\theta },\mathit{y})d\mathit{y}{\int }_{{\mathbb{R}}^m}\mathit{u}{(\mathit{\theta },\mathit{y})}^T\mathcal{CL}{(\mathit{\theta },\mathit{y})}^{1+\alpha }d\mathit{y}.\hfill \end{array}$$

(12)

Remark 1.

For $\alpha =0$ we get the CMLE of θ

$${\widehat{\mathit{\theta }}}_c=arg\underset{\mathit{\theta }\in \mathsf{\Theta }}{min}\left(-\frac{1}{n}\sum _{i=1}^n\log \mathcal{CL}(\mathit{\theta },{\mathit{y}}_i)\right).$$

(13)

At the same time it is well-known that

$$\sqrt{n}({\widehat{\mathit{\theta }}}_c-\mathit{\theta })\underset{n\to \infty }{\stackrel{\mathcal{L}}{⟶}}\mathcal{N}\left({\mathbf{0}}_p,{\mathit{G}}_{\ast }{\left(\mathit{\theta }\right)}^{-1}\right),$$

where ${\mathit{G}}_{\ast }\left(\mathit{\theta }\right)$ denotes the Godambe information matrix defined by ${\mathit{G}}_{\ast }\left(\mathit{\theta }\right)=\mathit{H}\left(\mathit{\theta }\right)\mathit{J}{\left(\mathit{\theta }\right)}^{-1}\mathit{H}\left(\mathit{\theta }\right)$, with $\mathit{H}\left(\mathit{\theta }\right)$ being the sensitivity or Hessian matrix and $\mathit{J}\left(\mathit{\theta }\right)$ being the variability matrix, defined, respectively, by

$$\begin{array}{c}\hfill \mathit{H}\left(\mathit{\theta }\right)=E_{\mathit{\theta }}\left[-\frac{\partial }{\partial \mathit{\theta }}\mathit{u}{(\mathit{\theta },\mathit{y})}^T\right],\mathit{J}\left(\mathit{\theta }\right)=E_{\mathit{\theta }}\left[\mathit{u}(\mathit{\theta },\mathit{y})\mathit{u}{(\mathit{\theta },\mathit{y})}^T\right].\end{array}$$

3. A New Model Selection Criterion

In order to describe the CLDIC criterion we consider the model $M_k$ given in (6). Following standard methodology (cf. [28], pp. 240), the most suitable candidate model to describe the data ${\mathit{Y}}_1,\dots ,{\mathit{Y}}_n$ is the model that minimizes the expected estimated DPD

$$E_{{\mathit{Y}}_1,\dots ,{\mathit{Y}}_n}\left[d_{\alpha }(g\left(·\right),\mathcal{CL}({\widehat{\mathit{\theta }}}_c^{\alpha },·))\right],$$

(14)

subject to the assumption that the unknown model g is belonging to $\mathsf{\Xi }$, i.e., the true model is included in ${\left\{M_s\right\}}_{s\in \left\{1,\dots ,\ell \right\}}$ and taking into account that ${\widehat{\mathit{\theta }}}_c^{\alpha }$, defined in (9), is a consistent and asymptotic normally distributed estimator of $\mathit{\theta }$. However, this expected value is still depending on the unknown parameter $\mathit{\theta }$. So, as a criterion, it should be used an asymptotically unbiased estimator of (14), for $g\in \mathsf{\Xi }$.

The most appropriate model to select is the model which minimizes the expected value

$$E_{{\mathit{Y}}_1,\dots ,{\mathit{Y}}_n}\left[W_{\alpha }\left({\widehat{\mathit{\theta }}}_c^{\alpha }\right)\right].$$

This expected value is still depending on the unknown parameter $\mathit{\theta }$. So, an asymptotically unbiased estimator of the above expected value could be the basis of a selection criterion, for $g\in \mathsf{\Xi }$. In order to proceed with the derivation of such an asymptotically unbiased estimator of $E_{{\mathit{Y}}_1,\dots ,{\mathit{Y}}_n}\left[W_{\alpha }\left({\widehat{\mathit{\theta }}}_c^{\alpha }\right)\right]$. The empirical version of $W_{\alpha }\left(\mathit{\theta }\right)$, in (7), is $W_{n,\alpha }\left(\mathit{\theta }\right)$, given in (8), and plays a central role in the development of the model selection criterion on the basis of the next theorem which expresses the expected value $E_{{\mathit{Y}}_1,\dots ,{\mathit{Y}}_n}\left[W_{\alpha }\left({\widehat{\mathit{\theta }}}_c^{\alpha }\right)\right]$ by means of the respective expected value of $W_{n,\alpha }\left({\widehat{\mathit{\theta }}}_c^{\alpha }\right)$, in an asymptotically equivalent way.

Theorem 1.

If the true distribution g belongs to the parametric family $\mathsf{\Xi }$ and ${\mathit{\theta }}_0$ denotes the true value of the parameter θ, then we have

$$E_{{\mathit{Y}}_1,\dots ,{\mathit{Y}}_n}\left[W_{\alpha }\left({\widehat{\mathit{\theta }}}_c^{\alpha }\right)\right]=E_{{\mathit{Y}}_1,\dots ,{\mathit{Y}}_n}\left[W_{n,\alpha }\left({\widehat{\mathit{\theta }}}_{\alpha }\right)+\frac{\alpha +1}{n}trace\left({\mathit{J}}_{\alpha }\left({\mathit{\theta }}_0\right){\mathit{H}}_{\alpha }{\left({\mathit{\theta }}_0\right)}^{-1}\right)\right]+o_p\left(1\right)$$

with ${\mathit{H}}_{\alpha }\left(\mathit{\theta }\right)$ and ${\mathit{J}}_{\alpha }\left(\mathit{\theta }\right)$ given in (11) and (12), respectively.

Based on the above theorem, the proof of which is presented in a full detail in the Appendix A, an asymptotic unbiased estimator of $E_{{\mathit{Y}}_1,\dots .,{\mathit{Y}}_n}\left[W_{\alpha }\left({\widehat{\mathit{\theta }}}_c^{\alpha }\right)\right]$ is given by

$$W_{n,\alpha }\left({\widehat{\mathit{\theta }}}_c^{\alpha }\right)+\frac{\alpha +1}{n}trace\left({\mathit{J}}_{\alpha }\left({\widehat{\mathit{\theta }}}_c^{\alpha }\right){\mathit{H}}_{\alpha }{\left({\widehat{\mathit{\theta }}}_c^{\alpha }\right)}^{-1}\right).$$

This ascertainment is the basis and a strong motivation for the next definition which introduces the model selection criterion.

Definition 2.

Let ${\left\{M_k\right\}}_{k\in \left\{1,\dots ,\ell \right\}}$ be candidate models for the observations ${\mathit{Y}}_1,\dots ,{\mathit{Y}}_n$. The selected model $M^{\ast }$ verifies

$$M^{\ast }=\underset{k\in \{1,\dots ,,\ell \}}{min}CLDI{C}_{\alpha }\left(M_k\right),$$

where

$$CLDI{C}_{\alpha }\left(M_k\right)=W_{n,\alpha }\left({\widehat{\mathit{\theta }}}_c^{\alpha }\right)+\frac{\alpha +1}{n}trace\left({\mathit{J}}_{\alpha }\left({\widehat{\mathit{\theta }}}_c^{\alpha }\right){\mathit{H}}_{\alpha }{\left({\widehat{\mathit{\theta }}}_c^{\alpha }\right)}^{-1}\right),$$

$W_{n,\alpha }\left(\mathit{\theta }\right)$ was given in (8) and ${\mathit{J}}_{\alpha }\left(\mathit{\theta }\right)$ and ${\mathit{H}}_{\alpha }\left(\mathit{\theta }\right)$ were defined in (11) and (12), respectively.

The next remark summarizes the model selection criterion in the case $\alpha =0$ and it therefore extends, in a sense, the pioneer and classic AIC.

Remark 2.

For $\alpha =0$ we have,

$$d_{KL}(g(·),\mathcal{CL}(\mathit{\theta },·))=W_0\left(\mathit{\theta }\right)+{\int }_{{\mathbb{R}}^n}g\left(\mathit{y}\right)logg\left(\mathit{y}\right)d\mathit{y}$$

with $W_0\left(\mathit{\theta }\right)=-{\int }_{{\mathbb{R}}^n}log\mathcal{CL}(\mathit{\theta },\mathit{y})g\left(\mathit{y}\right)d\mathit{y}$. Therefore, the most appropriate model which should be selected, is the model which minimizes the expected value

$$E_{{\mathit{Y}}_1,\dots ,{\mathit{Y}}_n}\left[W_0\left({\widehat{\mathit{\theta }}}_c\right)\right],$$

(15)

where ${\widehat{\mathit{\theta }}}_c$ is the CMLE of ${\mathit{\theta }}_0$ defined in (9).

The expected value (15) is still depending on the unknown parameter θ. A natural estimator of $W_0\left({\widehat{\mathit{\theta }}}_c\right)$ can be obtained by replacing the distribution function G, of g, by the empirical distribution function based on ${\mathit{Y}}_1,\dots ,{\mathit{Y}}_n$,

$$W_{n,0}\left(\mathit{\theta }\right)=-\frac{1}{n}\sum _{i=1}^{n}log\mathcal{CL}(\mathit{\theta },{\mathit{y}}_i).$$

Based on it, we select the model $M^{\ast }$ that verifies

$$M^{\ast }=\underset{k\in \{1,\dots ,\ell \}}{min}CLDI{C}_0\left(M_k\right),$$

with

$$CLDI{C}_0\left(M_k\right)=H_{n,0}\left({\widehat{\mathit{\theta }}}_c\right)+\frac{1}{n}trace\left(\mathit{J}\left({\widehat{\mathit{\theta }}}_c\right)\mathit{H}{\left({\widehat{\mathit{\theta }}}_c\right)}^{-1}\right),$$

where $\mathit{J}\left({\widehat{\mathit{\theta }}}_c\right)$ and $\mathit{H}\left({\widehat{\mathit{\theta }}}_c\right)$ are defined in Remark 1. In a manner, quite similar to that of the previous theorem, it can be established that $CLDI{C}_0\left(M_k\right)$ is an asymptotic unbiased estimator of $E_{{\mathit{Y}}_1,\dots ,{\mathit{Y}}_n}\left[W_0\left({\widehat{\mathit{\theta }}}_c\right)\right].$

This would be the model selection criterion in a composite likelihood framework based on Kullback–Leibler divergence. We can observe that this criterion coincides with the criterion given in [22] as a generalization of the classical criterion of Akaike, which will be referred from now as Composite Akaike Information Criterion (CAIC).

4. Numerical Simulations

4.1. Scenario 1: Two-Component Mixed Model

We are starting with a simulation example, which is motivated and follows ideas from the paper [29] and the Example 4.1 in [20] which will compare the behaviour of the proposed criteria with the CAIC criterion, for $\alpha =0$ (see Remark 2).

Consider the random vector $\mathit{Y}={(Y_1,Y_2,Y_3,Y_4)}^T$ from an unknown density g and let now ${\mathit{Y}}_1,\dots ,{\mathit{Y}}_n$ be independent and identically distributed replications of $\mathit{Y}$ which are described by the true but unknown distribution g. Taking into account that the true model g is unknown, suppose that $\{f(·;\mathit{\theta }),\mathit{\theta }\in \mathsf{\Theta }\subseteq {\mathbb{R}}^p,p\ge 1\}$ is a parametric identifiable family of candidate distributions to describe the observations ${\mathit{y}}_1,\dots ,{\mathit{y}}_n$. Let also $\mathcal{CL}(\mathit{\theta },\mathit{y})$ denotes the composite likelihood function associated to the parametric model $f(·;\mathit{\theta })$.

We consider the problem of choosing (on the basis of n independent and identically distributed replications ${\mathit{y}}_1,\dots ,{\mathit{y}}_n$ of $\mathit{Y}={(Y_1,Y_2,Y_3,Y_4)}^T$) between a 4-variate normal distribution, $\mathcal{N}\left({\mathit{\mu }}^N,\mathsf{\Sigma }\right)$, with ${\mathit{\mu }}^N={({\mu }_1^N,{\mu }_2^N,{\mu }_3^N,{\mu }_4^N)}^T$ and

$$\mathsf{\Sigma }=\left(\begin{array}{cccc}1& \rho & 2\rho & 2\rho \\ \rho & 1& 2\rho & 2\rho \\ 2\rho & 2\rho & 1& \rho \\ 2\rho & 2\rho & \rho & 1\end{array}\right),$$

and a 4-variate t-distribution with $\nu$ degrees of freedom, $t_{\nu }\left({\mathit{\mu }}^{t_{\nu }},{\mathsf{\Sigma }}^{\ast }\right)$, with different location parameters ${\mathit{\mu }}^{t_{\nu }}={({\mu }_1^{t_{\nu }},{\mu }_2^{t_{\nu }},{\mu }_3^{t_{\nu }},{\mu }_4^{t_{\nu }})}^T$ and same variance-covariance matrix $\mathsf{\Sigma }$, and density,

$$C_m{\left|{\mathsf{\Sigma }}^{\ast }\right|}^{-1/2}{\left[1+\frac{1}{\nu }{(\mathit{y}-{\mu }^{t_{\nu }})}^T{\left({\mathsf{\Sigma }}^{\ast }\right)}^{-1}(\mathit{y}-{\mathit{\mu }}^{t_{\nu }})\right]}^{-(\nu +m)/2},$$

with ${\mathsf{\Sigma }}^{\ast }=\frac{\nu -2}{\nu }\mathsf{\Sigma }$, $C_m={\left(\pi \nu \right)}^{-m/2}\frac{\mathsf{\Gamma }\left[\right(\nu +m)/2]}{\mathsf{\Gamma }(\nu /2)}$ and $m=4$.

Consider the composite likelihood function,

$$\mathcal{CL}N(\rho ,\mathit{y})=f_{A_1}^N(\mathit{y};\rho )f_{A_2}^N(\mathit{y};\rho ),$$

with $f_{A_1}^N(\mathit{y};\rho )=f_{12}^N(y_1,y_2;{\mu }_1^N,{\mu }_2^N;\rho )$ and $f_{A_2}^N(\mathit{y};\rho )=f_{34}^N(y_3,y_4;{\mu }_3^N,{\mu }_4^N;\rho )$, where $f_{12}^N$ and $f_{34}^N$ are the densities of the marginals of $\mathit{Y}$, i.e., bivariate normal distributions with mean vectors ${({\mu }_1^N,{\mu }_2^N)}^T$ and ${({\mu }_3^N,{\mu }_4^N)}^T$, respectively, and common variance-covariance matrix

$${\mathsf{\Sigma }}_0=\left(\begin{array}{cc}1& \rho \\ \rho & 1\end{array}\right).$$

In a similar manner consider the composite likelihood

$$\mathcal{CL}{t}_{\nu }(\rho ,\mathit{y})=f_{A_1}^{t_{\nu }}(\mathit{y};\rho )f_{A_2}^{t_{\nu }}(\mathit{y};\rho ),$$

with $f_{A_1}^{t_{\nu }}(\mathit{y};\rho )=f_{12}^{t_{\nu }}(y_1,y_2;{\mu }_1^{t_{\nu }},{\mu }_2^{t_{\nu }};\rho )$ and $f_{A_2}^{t_{\nu }}(\mathit{y};\rho )=f_{34}^{t_{\nu }}(y_3,y_4;{\mu }_3^{t_{\nu }},{\mu }_4^{t_{\nu }};\rho )$, where $f_{12}^{t_{\nu }}$ and $f_{34}^{t_{\nu }}$ are the densities of the marginals of $\mathit{Y}$, i.e., bivariate t-distributions with mean vectors ${({\mu }_1^{t_{\nu }},{\mu }_2^{t_{\nu }})}^T$ and ${({\mu }_3^{t_{\nu }},{\mu }_4^{t_{\nu }})}^T$, respectively, and common variance-covariance matrix

$${\mathsf{\Sigma }}_0=\left(\begin{array}{cc}1& \rho \\ \rho & 1\end{array}\right).$$

Under this formulation, the simulation study follows in the next two scenarios.

4.1.1. Scenario 1a

Following Example 4.1 in [20], the steps of the simulation study are the following:

• Generate 1000 samples of size $n=5,7,10,20,40,50,70,100$ from a two component mixture of two 4-variate distributions, namely, a 4-variate normal and a 4-variate t-distribution,

  $$h_{\omega }\left(\mathit{y}\right)=\omega N\left({\mathit{\mu }}^N,\mathsf{\Sigma }\right)+(1-\omega )t_{\nu }\left({\mathit{\mu }}^{t_{\nu }},{\mathsf{\Sigma }}^{\ast }\right),0\le \omega \le 1,$$

  with ${\mathit{\mu }}^N=(0,0,0.5,0)$ and ${\mathit{\mu }}^{t_{\nu }}=(3.2,1.5,0.5,2)$, for $\omega =0,0.25,0.45,0.5,0.55,0.75,1$, $\nu =5,10,30$ degrees of freedom and with specific values of $\rho =-0.15,-0.10,0.10$. As pointed out in [29], taking into account that $\mathsf{\Sigma }$ should be semi-positive definite, the following condition is imposed: $-\frac{1}{5}\le \rho \le \frac{1}{3}$.
• Estimate the common parameter $\rho$, separately in each model, by using the CMDPDE estimator for different values of the tuning parameter $\alpha =0,0.3$. The composite density which corresponds to the mixture $h_{\omega }\left(\mathit{y}\right)$ is defined by

  $$\mathcal{CL}(\rho ,\mathit{y})=\omega \mathcal{CL}N(\rho ,\mathit{y})+(1-\omega )\mathcal{CL}{t}_{\nu }(\rho ,\mathit{y}),0\le \omega \le 1,$$

  and it is used to obtain the CMDPDE estimator, $\widehat{\rho }$, of $\rho$.
• Define the mixture composite likelihood function

  $$\mathcal{CL}(\widehat{\rho },\mathit{y})=\omega \mathcal{CL}N(\widehat{\rho },\mathit{y})+(1-\omega )\mathcal{CL}{t}_{\nu }(\widehat{\rho },\mathit{y}),0\le \omega \le 1.$$
• Calculate $CLDI{C}_{\alpha }\left(M_k\right)$, the value of the model selection criterion considered in this paper, for the two candidate models, with

  $$CLDI{C}_{\alpha }\left(M_k\right)=W_{n,\alpha }\left(\widehat{\rho }\right)+\frac{\alpha +1}{n}trace\left({\mathit{J}}_{\alpha }\left(\widehat{\rho }\right){\mathit{H}}_{\alpha }{\left(\widehat{\rho }\right)}^{-1}\right).$$
  An explanation of how to obtain this value for the both candidate models is given in Appendix B.
• Compute the times that the 4-variate normal model was selected.

Results are summarized in

Table 1. Main results, Scenario 1a.

	$\mathit{\alpha }=0$ (CAIC)							$\mathit{\alpha }=0.3$
$\mathit{\omega }$	0	0.25	0.45	0.5	0.55	0.75	1	0	0.25	0.45	0.5	0.55	0.75	1
$\nu =5,\rho =-0.15$
n = 5	0	1	269	499	713	996	1000	0	0	273	498	712	1000	1000
7	0	1	246	504	758	998	1000	0	1	220	511	738	999	1000
10	0	0	202	482	775	1000	1000	0	0	185	467	771	1000	1000
20	0	0	114	486	871	1000	1000	0	0	112	473	866	1000	1000
40	0	0	41	459	947	1000	1000	0	0	54	496	954	1000	1000
50	0	0	21	475	964	1000	1000	0	0	41	556	986	1000	1000
70	0	0	9	461	985	1000	1000	0	0	48	656	995	1000	1000
100	0	0	5	472	992	1000	1000	0	0	142	885	1000	1000	1000
$\nu =10,\rho =-0.15$
5	0	3	222	445	688	996	1000	0	3	218	433	688	997	1000
7	0	1	191	439	720	1000	1000	0	0	179	431	690	999	1000
10	0	0	163	432	747	1000	1000	0	0	152	402	725	1000	1000
20	0	0	59	399	819	1000	1000	0	0	49	361	773	1000	1000
40	0	0	19	336	912	1000	1000	0	0	12	326	899	1000	1000
50	0	0	6	362	936	1000	1000	0	0	10	334	925	1000	1000
70	0	0	1	292	960	999	1000	0	0	2	356	973	1000	1000
100	0	0	0	301	983	1000	1000	0	0	1	531	992	1000	1000
$\nu =30,\rho =-0.15$
5	0	4	237	423	677	997	1000	0	2	235	413	656	996	1000
7	0	0	155	394	689	1000	1000	0	0	141	379	677	999	1000
10	0	0	144	413	719	1000	1000	0	0	134	393	701	1000	1000
20	0	0	57	351	801	1000	1000	0	0	40	311	764	1000	1000
40	0	0	11	296	904	1000	1000	0	0	8	263	882	1000	1000
50	0	0	6	271	918	1000	1000	0	0	3	253	903	1000	1000
70	0	0	1	225	942	1000	1000	0	0	0	229	941	1000	1000
100	0	0	0	208	978	1000	1000	0	0	0	303	989	1000	1000
$\nu =10,\rho =-0.10$
5	0	4	242	464	680	996	1000	0	3	238	459	682	999	1000
7	0	0	187	461	733	997	1000	0	0	199	457	731	998	1000
10	0	0	162	445	738	1000	1000	0	0	165	407	713	1000	1000
20	0	0	62	378	807	1000	1000	0	0	59	354	789	1000	1000
40	0	0	19	357	902	999	1000	0	0	14	333	895	1000	1000
50	0	0	6	325	932	1000	1000	0	0	8	325	931	1000	1000
70	0	0	2	305	954	1000	1000	0	0	6	367	967	1000	1000
100	0	0	0	307	979	1000	1000	0	0	2	507	993	1000	1000
$\nu =10,\rho =0.10$
5	0	11	268	459	669	991	1000	1	11	268	478	680	993	1000
7	0	1	211	456	720	999	1000	0	3	207	464	716	998	1000
10	0	0	168	423	704	1000	1000	0	0	162	403	702	1000	1000
20	0	0	86	360	789	1000	999	0	0	89	357	786	1000	1000
40	0	0	35	367	893	1000	1000	0	0	38	398	896	1000	1000
50	0	0	19	331	886	1000	1000	0	0	19	360	913	1000	1000
70	0	0	11	311	933	1000	1000	0	0	16	379	963	1000	1000
100	0	0	2	276	969	1000	1000	0	0	7	490	985	1000	1000

. Extreme values of $\omega =0,1$ represent the times that the 4-variate normal model was selected under the 4-variate t-distribution and 4-variate normal distribution, respectively. This means that, for $\omega =1$, the perfect discrimination will be achieved when 1000 of the 1000 simulated samples are correctly assigned, while for $\omega =0$, the more near to 0, the better discrimination of the criterion. $\omega =0.5$ means that each sample was generated both from the normal and t-distribution in the same proportion.

4.1.2. Scenario 1b

Same Scenario is evaluated under the more-closed means ${\mathit{\mu }}^N=(0,1.5,0.5,-0.75)$ and ${\mathit{\mu }}^{t_{\nu }}=(0,1.5,0.5,2)$ for moderate-large sample sizes and $\alpha \in \{0,0.2,0.4\}$. Here $\nu =5$ and $\rho =-0.15$. Results are shown in

Table 2. Main results, Scenario 1b.

	$\mathit{\alpha }=0$ (CAIC)				$\mathit{\alpha }=0.2$				$\mathit{\alpha }=0.4$
	0	0.25	0.75	1	0	0.25	0.75	1	0	0.25	0.75	1
n = 40	0	0	39	731	0	0	537	961	0	0	580	949
50	0	0	24	732	0	0	859	990	0	0	944	994
60	0	0	14	772	0	0	999	1000	0	1	999	1000
70	0	0	9	734	0	0	999	1000	0	27	999	1000
80	0	0	5	770	0	1	1000	1000	0	326	1000	1000
90	0	0	4	782	0	23	1000	1000	2	794	1000	1000
100	0	0	4	802	0	173	1000	1000	26	978	1000	1000

. In this case, the models under consideration are more similar, so it would be understandable that the CLDIC criterion did not discriminate in such as good way.

4.2. Scenario 2: Three-Component Mixed Model

Now, we consider a mixed model composed on two 4-variate normal distributions and a 4-variate t-distribution with $\nu =10$ degrees of freedom. The three distributions have common variance-covariance matrix, as in the previous scenario, with unknown $\rho =-0.15$ and different but known means ${\mathit{\mu }}_1^N=(0,0,0.5,0)$, ${\mathit{\mu }}_2^N=(0,1.5,0.5,0)$ and ${\mathit{\mu }}^t=(0,1.5,0.5,2)$. The model is defined by

$$\omega \mathcal{N}({\mathit{\mu }}_1^N,\mathsf{\Sigma })+\lambda \mathcal{N}({\mathit{\mu }}_2^N,\mathsf{\Sigma })+(1-\omega -\lambda )t_{\nu =10}({\mathit{\mu }}^t,{\mathsf{\Sigma }}^{\ast }),0\le \omega ,\lambda ,\omega +\lambda \le 1,$$

with $\mathsf{\Sigma }$ being again a common variance-covariance matrix with unknown parameter $\rho$ of the form

$$\mathsf{\Sigma }=\left(\begin{array}{cccc}1& \rho & 2\rho & 2\rho \\ \rho & 1& 2\rho & 2\rho \\ 2\rho & 2\rho & 1& \rho \\ 2\rho & 2\rho & \rho & 1\end{array}\right).$$

Following the same steps that in the first scenario, we generate 1000 samples of the three-component mixture for different sample sizes $n=5,7,10,20,40,50,70,100$ and different values of $\omega$ and $\lambda$. Then, we consider the problem of choosing among one of the two 4-variate normal distributions and the 4-variate t-distribution through the CLDIC criterion, for different values of the tuning parameter $\alpha =0,0.3,0.5,0.7$. See

Table 3. Main results, Scenario 2.

	$\mathit{\alpha }=0$ (CAIC)			$\mathit{\alpha }=0.3$			$\mathit{\alpha }=0.5$			$\mathit{\alpha }=0.7$
Model ${}^{\ast }$	N1	N2	MT	N1	N2	MT	N1	N2	MT	N1	N2	MT
True model: $\mathcal{N}({\mathit{\mu }}_1^N,\mathsf{\Sigma })$
n = 5	957	24	19	950	16	34	939	23	38	936	28	36
7	970	19	11	966	13	24	961	13	26	950	22	28
10	993	3	4	986	4	10	979	6	15	971	6	23
20	1000	0	0	1000	0	0	998	0	2	997	0	3
40	1000	0	0	1000	0	0	1000	0	0	1000	0	0
50	1000	0	0	1000	0	0	1000	0	0	1000	0	0
70	1000	0	0	1000	0	0	1000	0	0	1000	0	0
100	1000	0	0	1000	0	0	1000	0	0	999	0	0
True model: $\mathcal{N}({\mathit{\mu }}_2^N,\mathsf{\Sigma })$
5	29	638	333	34	610	356	38	639	323	50	646	304
7	15	622	363	13	589	398	17	599	384	28	627	345
10	6	610	384	5	540	455	5	540	455	11	586	403
20	1	612	387	1	518	481	1	472	527	1	527	472
40	0	566	434	0	650	350	0	590	410	0	614	386
50	0	561	439	0	804	196	0	797	203	0	835	165
70	0	584	416	0	987	13	0	994	6	0	998	2
100	0	520	480	0	1000	0	0	1000	0	0	1000	0
True model: $t_{\nu =10}({\mathit{\mu }}^t,\mathsf{\Sigma })$
5	2	15	983	1	6	993	1	8	991	3	15	982
7	0	3	997	0	1	999	2	2	996	0	4	996
10	0	1	999	0	2	998	0	2	998	0	3	997
20	0	0	1000	0	0	1000	0	0	1000	0	0	1000
40	0	0	1000	0	0	1000	0	0	1000	0	0	1000
50	0	0	1000	0	0	1000	0	0	1000	0	0	1000
70	0	0	1000	0	0	1000	0	0	1000	0	0	1000
100	0	0	1000	0	0	1000	0	4	996	0	296	704
True model: $0.7\mathcal{N}({\mathit{\mu }}_2^N,\mathsf{\Sigma })+0.3t_{\nu =10}({\mathit{\mu }}^t,\mathsf{\Sigma })$
5	6	384	610	6	375	619	4	401	595	11	452	537
7	1	331	668	1	294	705	1	317	682	1	373	626
10	1	261	738	1	218	781	1	253	746	1	306	693
20	0	109	891	0	101	899	0	107	893	0	141	859
40	0	26	974	0	126	874	0	122	878	0	166	834
50	0	13	987	0	311	689	0	345	655	0	445	555
70	0	6	994	0	948	52	0	982	18	0	994	6
100	0	2	998	0	1000	0	0	1000	0	0	999	1
True model: $\frac{1}{3}\mathcal{N}({\mathit{\mu }}_1^N,\mathsf{\Sigma })+\frac{1}{3}\mathcal{N}({\mathit{\mu }}_2^N,\mathsf{\Sigma })+\frac{1}{3}{t}_{\nu =10}({\mathit{\mu }}^t,\mathsf{\Sigma })$
5	127	377	496	121	363	516	107	392	501	107	424	469
7	87	357	556	70	339	591	66	356	578	63	396	541
10	69	326	605	61	314	625	56	330	614	45	381	574
20	37	259	704	25	298	677	17	337	646	15	349	636
40	7	145	848	9	452	539	4	508	488	1	469	530
50	2	122	876	5	744	251	3	814	183	3	853	144
70	0	99	901	4	996	0	4	996	0	4	996	0
100	0	36	964	355	645	0	645	355	0	856	144	0

${}^{\ast }$ Here the model candidates are expressed as N1, N2, MT to denote $\mathcal{N}({\mathit{\mu }}_1^N,\mathsf{\Sigma })$, $\mathcal{N}({\mathit{\mu }}_2^N,\mathsf{\Sigma })$ and $t_{10}({\mathit{\mu }}^t,\mathsf{\Sigma })$, respectively.

for results. Here, the normal models are denoted by N1 and N2, respectively, while the 4-variate t-distribution is denoted by MT. The first three cases evaluate the selected model under these multivariate distributions. In the last two scenarios, a mixed model is considered as the true distribution.

4.3. Discussion of Results

In Scenario 1a, two well-differentiated multivariate models are considered. In this case CLDIC criterion works in a very efficient way, with an almost-perfect discrimination for extreme values of $\omega$. The good behaviour is also observed for not so extreme values of $\omega$, such as $\omega =0.55$ or $0.45$. We can not observe a significant difference in the choice of $\alpha$.

In Scenario 1b we consider closer models, which affect the discrimination power of the CLDIC. However, in this case, we do observe great differences when considering different $\alpha$. While the discrimination power of CLDIC for $\alpha =0$ (CAIC) and $\omega =1$ is around $75\%$, for $\alpha =0.2$ or $\alpha =0.4$ the behaviour is excellent. This happens also for large but not extreme values of $\omega$, such as $\omega =0.75$. However, a medium value of $\alpha$ turns into a worse discrimination for low values of $\omega$.

Scenario 2 deals with three different models, two multivariate normal and one multivariate t (N1, N2 and MT, respectively). The second normal distribution is closer to MT in terms of means. While CLDIC criterion discriminate well between N1 and N2 and between N1 and MT, it has difficulties in distinguishing N2 an MT distributions, overall for small samples sizes and $\alpha =0$.

It seems, therefore, that when we have well-discriminated models, CLDIC criterion works very well, independently of the sample size and the tuning parameter $\alpha$ considered. Dealing with closer models leads, as expected, to worst results, overall for $\alpha =0$ (CAIC).

Note that the behaviour of Wald-type and Rao tests based on CMDPDEs was studied in [12,13] through extensive simulation studies.

5. Numerical Examples

5.1. Choice of the Tuning Parameter

In the previous sections, we have seen that CLDIC criterion works generally very well, independently of $\alpha$, but that some values present a better behaviour, overall when distinguishing similar models. In these situations, it appears that values close to $0.2$ or $0.3$ work well, while CAIC criterion presents a worse behaviour. A data-driven approach for the choice of the tuning parameter which would be helpful in practice. The approach of [30] was adapted In [13], for the choice of the optimum $\alpha$ in CMDPDEs. This approach consisted on minimizing the estimated mean squared error by means of a pilot estimator, ${\mathit{\theta }}^P$. This approximation is given by

$${\widehat{MSE}}_{\alpha }={({\widehat{\mathit{\theta }}}_c^{\alpha }-{\mathit{\theta }}^P)}^T({\widehat{\mathit{\theta }}}_c^{\alpha }-{\mathit{\theta }}^P)+\frac{1}{n}\mathrm{Trace}\left({\mathit{H}}_{\alpha }^{-1}\left({\widehat{\mathit{\theta }}}_c^{\alpha }\right){\mathit{J}}_{\alpha }\left({\widehat{\mathit{\theta }}}_c^{\alpha }\right){\mathit{H}}_{\alpha }^{-1}\left({\widehat{\mathit{\theta }}}_c^{\alpha }\right)\right),$$

(16)

where ${\mathit{H}}_{\alpha }\left(\mathit{\theta }\right)$ and ${\mathit{J}}_{\alpha }\left(\mathit{\theta }\right)$ are given in (11) and (12). The optimum $\alpha$ will be the one that minimizes expression (16). The choice of the pilot estimator is probably one of the major drawbacks of this approach, as it may lead to a choice of $\alpha$ too close to that used for the pilot estimator. A pilot estimator with $\alpha \approx 0.4$, was proposed in [13] after some simulations, in concordance with [30], where the initial choice of a pilot is suggested to be a robust one in order to obtain the best results in terms of robustness.

5.2. Iris Data

The Iris data (Fisher, [31]) includes 3 categories of 50 sample values each, where each category refers to a type of iris plant: setosa, versicolor and virginica. Each plant is categorized in its class and described by other 4 variables: (1) sepal length, (2) sepal width, (3) petal length and (4) petal width. This is one of the most known data sets for discriminant analysis. [32] proposed the use of a Gaussian finite mixture for modeling Iris data, in which each known class is modeled by a single Gaussian term with the same variance-covariance matrix. The resulting model is as follows

$$\begin{array}{c}\hfill f\left(\mathit{x}\right)=\frac{1}{3}\mathcal{N}({\mathit{\mu }}_1,\mathsf{\Sigma })+\frac{1}{3}\mathcal{N}({\mathit{\mu }}_2,\mathsf{\Sigma })+\frac{1}{3}\mathcal{N}({\mathit{\mu }}_3,\mathsf{\Sigma }),\end{array}$$

(17)

with

$${\mathit{\mu }}_1={({\mu }_{11},{\mu }_{12},{\mu }_{13},{\mu }_{14})}^T,{\mathit{\mu }}_2={({\mu }_{21},{\mu }_{22},{\mu }_{23},{\mu }_{24})}^T,{\mathit{\mu }}_3={({\mu }_{31},{\mu }_{32},{\mu }_{33},{\mu }_{34})}^T$$

and

$$\begin{array}{c}\hfill \mathsf{\Sigma }=\left(\begin{array}{cccc}{\sigma }_1^2& {\sigma }_{12}& {\sigma }_{13}& {\sigma }_{14}\\ {\sigma }_{21}& {\sigma }_2^2& {\sigma }_{23}& {\sigma }_{24}\\ {\sigma }_{31}& {\sigma }_{32}& {\sigma }_3^2& {\sigma }_{34}\\ {\sigma }_{41}& {\sigma }_{42}& {\sigma }_{43}& {\sigma }_4^2\end{array}\right).\end{array}$$

Exact values can be obtained by MclustDA() function of mclust package in R Software ([32]).

We propose a composite likelihood approach to modeling (17) where we suppose independence between the two first and two last variables. This is

$$f_{CL}\left(\mathit{y}\right)=\frac{1}{3}CL{N}_1+\frac{1}{3}CL{N}_2+\frac{1}{3}CL{N}_3,$$

(18)

with

$$\begin{array}{c}\hfill CL{N}_i=f_{A_{i1}}^N({\rho }_{12},\mathit{y})f_{A_{i2}}^N({\rho }_{34},\mathit{y}),\end{array}$$

where $f_{A_{i1}}^N({\rho }_{12},\mathit{y})=f_{A_{i1}}^N({\rho }_{12},{\mu }_{i1},{\mu }_{i2},{\mathsf{\Sigma }}_{A_1},\mathit{y})$ and $f_{A_{i2}}^N({\rho }_{34},\mathit{y})=f_{A_{i2}}^N({\rho }_{34},{\mu }_{i3},{\mu }_{i4},{\mathsf{\Sigma }}_{A_2},\mathit{y})$, $i=1,2,3$ are bivariate normals with variance-covariance matrices

$$\begin{array}{c}\hfill {\mathsf{\Sigma }}_{A_1}=\left(\begin{array}{cc}{\sigma }_1^2& {\rho }_{12}{\sigma }_1{\sigma }_2\\ {\rho }_{12}{\sigma }_1{\sigma }_2& {\sigma }_2^2\end{array}\right),{\mathsf{\Sigma }}_{A_2}=\left(\begin{array}{cc}{\sigma }_3^2& {\rho }_{34}{\sigma }_3{\sigma }_4\\ {\rho }_{34}{\sigma }_3{\sigma }_4& {\sigma }_4^2\end{array}\right).\end{array}$$

We are going to evaluate the behavior of the CLDIC criterion proposed in previous sections. After estimating parameters ${\rho }_{12}$ and ${\rho }_{34}$ in (18), we consider 10 different subsets of the IRIS data:

• SE subset: 50 first observations, corresponding to Setosa plants ($n=50$).
• VE subset: 50 second observations, corresponding to Versicolor plants ($n=50$).
• VI subset: 50 last observations, corresponding to Virginica plants ($n=50$).
• SE(VE) subset: SE subset with 2 first observations of VE subset ($n=52$).
  Equivalently: SE(VI), VE(SE), VE(VI), VI(SE) and VI(VE).
• VI(SE+VE) subset: VI subset with 2 first observations of SE and VE subsets ($n=54$).

In

Table 4. Selected model in each of the subsets. Iris data.

$\mathit{\alpha }$	SE	VE	VI	SE(VE)	SE(VI)	VE(SE)	VE(VI)	VI(SE)	VI(VE)	VI(SE+VE)
0 (CAIC)	CN1	CN2	CN3	CN1	CN1	CN1${}^{\ast }$	CN2	CN1${}^{\ast }$	CN3	CN3
0.2	CN1	CN2	CN3	CN1	CN1	CN2	CN2	CN3	CN3	CN3
0.3	CN1	CN2	CN3	CN1	CN1	CN2	CN2	CN3	CN3	CN3
0.4	CN1	CN2	CN3	CN1	CN1	CN2	CN2	CN1${}^{\ast }$	CN3	CN3
0.5	CN1	CN2	CN3	CN1	CN1	CN2	CN2	CN1${}^{\ast }$	CN3	CN3
0.8	CN1	CN2	CN3	CN1	CN1	CN2	CN2	CN1${}^{\ast }$	CN3	CN3
0.22	CN1	CN2	CN3	CN1	CN1	CN2	CN2	CN3	CN3	CN3

, chosen models for each one of the subsets are obtained by the proposed CLDIC criterion. When a “pure” subset is considered, all the tuning parameters lead to optimal decisions, but when a “contaminated” subset is under consideration, only $\alpha =0.2,0.3$ have an optimal response in all the cases.

We now apply the ad hoc approach presented in Section 5.1 for selecting the tuning parameter $\alpha$ in a composite likelihood framework. Applying this procedure to our data set though a grid search of length 100 and by means of a pilot estimator with $\alpha =0.4$ leads to the optimal tuning parameter $\alpha =0.22$, what is in concordance with the obtained results (see

Table 5. Selected $\alpha$ for different pilot estimators, ad-hoc tuning parameter selection procedure. Iris and Wine data.

	${\mathit{\alpha }}_{\mathbf{pilot}}$	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
Iris	${\alpha }_{opt}$	0.31	0.17	0.20	0.21	0.22	0.23	0.24	0.24	0.25	0.25	0.25
Wine	${\alpha }_{opt}$	0.45	0.46	0.47	0.49	0.51	0.53	0.55	0.56	0.56	0.56	0.57

). We can see that the use of other pilot estimators would not affect very much to the final decission.

5.3. Wine Data

We now work with Wine data ([33]), which contain a chemical analysis of 178 Italian wines from three different cultivars (Barolo, Grignolino, Barbera) yielded 13 measurements. In order to illustrate our criterion, we will work with only first four explanatory variables: Alcohol, Malic, Ash and Alkalinity. As in the previous section, we adjust a Gaussian mixture model with weights, in this case: $59/178$, $72/178$ and $47/178$ corresponding to Barolo, Grignolino and Barbera classes, respectively. We now consider these 10 different subsets of the Wine data:

• BO subset: 20 first observations of Barolo wines ($n=20$).
• GR subset: 20 first observations of Grignolino wines ($n=20$).
• BA subset: 20 first observations of Barbera wines ($n=20$).
• BO(GR) subset: BO subset with 5 first observations of GR subset ($n=25$).
  Equivalently: BO(BA), GR(BO), GR(BA), BA(BO) and BA(GR).
• BA(BO+GR) subset: BA subset with 3 first observations of BO and GR subsets ($n=26$).

We can observe how, for medium values of $\alpha$, the discrimination is perfect (see

Table 6. Selected model in each of the subsets. Wine data.

$\mathit{\alpha }$	BO	GR	BA	BO(GR)	BO(BA)	GR(BO)	GR(BA)	BA(BO)	BA(GR)	BA(BO+GR)
0 (CAIC)	CN1	CN2	CN3	CN1	CN1	CN2	CN2	CN3	CN3	CN2${}^{\ast }$
0.2	CN1	CN2	CN3	CN1	CN1	CN2	CN2	CN3	CN3	CN3
0.3	CN1	CN2	CN3	CN1	CN1	CN2	CN2	CN3	CN3	CN3
0.4	CN1	CN2	CN3	CN1	CN1	CN2	CN2	CN3	CN3	CN3
0.5	CN1	CN2	CN3	CN1	CN1	CN2	CN2	CN3	CN3	CN3
0.8	CN1	CN2	CN3	CN1	CN1	CN2	CN2	CN2${}^{\ast }$	CN2${}^{\ast }$	CN3
0.51	CN1	CN2	CN3	CN1	CN1	CN2	CN2	CN3	CN3	CN3

). Applying ad-hoc tuning parameter choice procedure we obtain ${\alpha }_{opt}\approx 0.51$, with a perfect discrimination again (Table 5).

6. Conclusions and Future Research

In this paper, we have addressed the problem of model selection in the framework of composite likelihood methodology, on the basis of the DPD as a measure of the closeness of the composite density and the true model that drives the data. In this context, an information criterion is introduced and studied which is defined by means of composite minimum distance type estimators of the unknown parameters, well-known for having nice robustness properties. Thanks to a simulation study, we have shown that the proposed here model selection criterion works well in practice and mainly that the use of CMDPDE makes the criterion more robust than the criteria based on the classic CMLE and the Kullback–Leibler divergence, given in [22]. The analysis of two real data examples of the literature illustrate on how the model selection criterion, presented here, can be applied in practical cases. This paper is a part of a series of papers by the authors where composite likelihood ideas and methods are harmonically weaved with divergence theoretic methods in order to develop statistical inference (estimation and testing of hypotheses) and model selection criteria, as well. We envision future work in some directions. The development of change point methodology on the basis of composite density with CMDPDE and divergence measures would be maybe an appealing problem for a future research on the topic. However, all the information theoretic methods developed on the basis of the composite likelihood depend on the choice of the family of sets ${\left\{A_k\right\}}_{k=1}^K$, appeared in Formula (1). A question is raised at this point: how the information theoretic procedures developed on the basis of the composite likelihood are affected by this family of sets? It is an appealing problem which deserves also investigation in a future work.