Penalized EM algorithm and copula skeptic graphical models for inferring networks for mixed variables

Graphical models are efficient tools for studying multivariate distributions through a compact, graphical representation of the joint probability distribution of the underlying random variables. The treatment of graphical models simplifies significantly, when one focuses on normally distributed variables. Let the random vector $Y=(Y_{1},\ldots,Y_{p})^{T}$ be assumed to be Gaussian with a positive-definite covariance matrix $\Sigma$ of dimension $p\times p$. A graphical model $G=(V,E)$, where $V$ corresponds to the set of nodes with $p$ elements and $E\subset V\times V$ of ordered pairs of distinct nodes called the edges of $G$, for $\mathcal{N}_{p}(0,\Sigma)$ is called a Gaussian graphical model, if on the graph $G$, the edges $E$ represent conditional dependence among the random variables. Absence of an edge between any pair of random variables or zero value of a precision matrix, $\Theta=\Sigma^{-1}$, corresponds to conditional independence of the two random variables given the remaining ones.

In practice, we encounter both discrete and continuous variables that may not be Gaussian. Thus, the assumption of multivariate normal distribution would be too restrictive. To relax the normality requirement, we use the copula framework to construct multivariate distributions for arbitrary marginals as discussed in Section 1 above. In order to use the properties of the Gaussian graphical model, we consider the Gaussian copula. The Gaussian copula with correlation matrix $\Gamma$ of dimension $p\times p$ having $p(p-1)/2$ free parameters is given by:

$$C(u_{1},\ldots,u_{p}\mid\Gamma)=\Phi_{p}(\Phi^{-1}(u_{1}),\ldots,\Phi^{-1}(u_{p})\mid\Gamma),$$

and the corresponding Gaussisn copula-based distribution function is

$\displaystyle H(y_{1},\ldots,y_{p}|\Gamma,F_{1},\ldots,F_{p})=\Phi_{p}(\Phi^{-1}(F_{1}(y_{1})),\ldots,\Phi^{-1}(F_{p}(y_{p}))\mid\Gamma).$

(2.1)

Here $\Phi(\cdot)$ represents the CDF of the standard normal distribution and $\Phi_{p}(\cdot\mid\Gamma)$ is the CDF of multivariate normal distribution $\mathcal{N}_{P}(0,\Gamma)$.

We note that under the Gaussian copula the correlation matrix $\Gamma$ is the matrix of correlation coefficients among the transformed variables $Z_{j}=\Phi^{-1}(F_{j}(y_{j}))$, $j=1,\ldots,p$, which represent the maximum pairwise correlations among the $Y_{j}$s, $j=1,\ldots,p$. If the univariate marginal distributions are normal, then entries of the correlation matrix represent exactly the pairwise correlation coefficients of the variables.

We now focus on graphical modeling for observed variables $Y$ of mixed type, i.e. in the case that they represent a collection of continuous, binary, ordinal or count variables. Suppose the $j$-th variable $Y_{j}$ has marginal distribution $F_{j}$ with its pseudo-inverse $F_{j}^{-1}$. In copula modeling, the marginals are treated as nuisance parameters and estimated nonparametrically mainly using the rescaled empirical distribution: $\hat{F}_{j}(y)=\frac{1}{n+1}\sum_{i=1}^{n}\mathcal{I}\{Y_{ji}\leq y\}$, $j=1,\ldots,p$. A copula graphical model, discussed above, can be constructed by introducing a vector of latent variables $Z\sim\mathcal{N_{p}}(0,\Gamma)$ that are related to the observed variables $Y$ as $Y_{j}=\hat{F}_{j}^{-1}(\Phi(Z_{j}))$, $j=1,\ldots,p$. In the case of mixed variables, the graphical structure, i.e. the conditional independence implied by the graph structure, is assumed to hold exclusively on the latent variable $Z$.

The aim of the inference procedure is to infer the graphical structure $G$, defined by the latent variable $Z$. Though the $Z$s are not observable, the observed $Y_{j}$s do provide a limited amount of information about them. Since the $\hat{F}_{j}$s are nondecreasing, observing $Y_{k_{1}}<Y_{k_{2}}$ implies that $Z_{k_{1}}<Z_{k_{2}}$. More generally, observing the n-dimensional vector $Y_{i}$ tells us that $Z_{i}$ lies in

$$D(Y_{i})=\{z\in\Re^{n}|a_{i}(Y_{ij})<Z_{j}\leq b_{i}(Y_{ij}))\},$$

(2.2)

where $a_{i}(y)=\inf\{z|\hat{F}_{j}^{-1}(\Phi(z))=y\}$ and $b_{i}(y)=\sup\{z|\hat{F}_{j}^{-1}(\Phi(z))=y\}$. In fact, for every ordinal $Y_{j}$, we can identify a set of thresholds $\tau_{j}=\left(\tau_{j0},\tau_{j1},\ldots,\tau_{j{n_{j}}}\right)$ with

$$-\infty=\tau_{j0}<\tau_{j1}<\cdots<\tau_{j{n_{j}}}=\infty,$$

such that for some ordered set of values $\left\{c_{j1}<\cdots<c_{j{w_{j}}}\right\}$,

$\displaystyle Y_{j}=\sum^{n_{j}}_{r=1}c_{jr}\times I\left\{{\tau_{j,{r-1}}<z_{j}\leq\tau_{jr}}\right\}.$

It follows that the mapping of the ordered values of $Y_{j}$ into some defined intervals $(\tau_{jr},\tau_{j,r+1}]$ of the latent variable $Z_{j}$ relies on the following relationship

$$\tau_{jr}=\Phi^{-1}(\hat{F}_{j}(c_{jr}))),~~r=1,\ldots,n_{j}-1.$$

The collection of these intervals is the set $D(Y)=D(Y_{1},\ldots,Y_{n})$ in (2.2). In case of missing observations, we consider data missing in this study as missing completely at random so that the missing values are easily determined from the latent variable distribution defined on the interval $(-\infty,\infty)$. Then the occurrence of event $Z\in D(Y)$ is taken as the observed data to infer about the copula parameter or the graph structure separately from the marginal distributions. Such inference approach is similar to the extended rank likelihood in Hoff (2007) and copula Gaussian graphical modeling in Dobra and Lenkoski (2011).

The Expectation-Maximization (EM) algorithm (Dempster, Laird, and Rubin, 1977) is a popular method for maximum likelihood estimation in the case of incomplete data, which naturally occur in our setting as a result of the latent nature of $Z$ as discussed in the previous section. In this section we consider EM algorithm with penalized likelihood. Green (1990) studied convergence properties of the EM algorithm for penalized likelihood.

The marginal likelihood of $Y$ where we consider $F_{1},\ldots,F_{p}$ as nuisance parameters; see also Hoff (2007), is

$\displaystyle L_{Y}(\Theta)=\int_{D(Y)}p(z\mid\Theta)dz$

(3.1)

For large sample sizes the precision matrix $\Theta$ can be estimated by maximizing the log-likelihood $l(\Theta)$ as a function of $\Theta$. Whereas for high-dimensional data $(n<<p)$, we add an $l_{1}$-norm penalty to encourage sparsity in the precision matrix and the identification of the underlying graph. The $l_{1}$ penalized log-likelihood is given by

$\displaystyle\ell_{\lambda}(\Theta)=\log L_{Y}(\Theta)-\lambda\left\|\Theta\right\|_{1},$

(3.2)

where the scalar parameter $\lambda\geq 0$ controls the size of the penalty.

Due to the complexity of maximizing the marginal log-likelihood $l_{Y}(\Theta)$ in (3.1) and (3.2), we employ the EM algorithm. We have discussed in the previous section that the observed data $Y$ provide some information on the latent variables $Z$ such that the occurence of the event $Z\in D(Y)$ is taken as the observed data to infer on $\Theta$. We now recall from standard EM algorithm setting that the loglikelihood of the observed data can be expressed as

$$\ell(\Theta)=Q(\Theta\mid\Theta^{(m)})-H(\Theta\mid\Theta^{(m)}),$$

(3.3)

where $\Theta^{(m)}$ an estimate of $\Theta$ from the previous step of the algorithm,

$$Q(\Theta\mid\Theta^{(m)})=E\left[\log L_{Z,Z\in D(Y)}(\Theta)\mid Z\in D(Y),\Theta^{(m)}\right],$$

(3.4)

and

$$H(\Theta\mid\Theta^{(m)})=E\left[\log L_{Z\mid Z\in D(Y)}(\Theta)\mid Z\in D(Y),\Theta^{(m)}\right].$$

(3.5)

The penalized log-likelihood takes the form

$\displaystyle\ell_{\lambda}(\Theta)=Q(\Theta\mid\Theta^{(m)})-H(\Theta\mid\Theta^{(m)})-\lambda\left\|\Theta\right\|_{1},$

(3.6)

such that

$\displaystyle\Theta_{\lambda}^{({m+1})}=\textrm{argmax}_{\Theta}\left\{Q(\Theta\mid\Theta^{(m)})-H(\Theta\mid\Theta^{(m)})-\lambda\left\|\Theta\right\|_{1}\right\}.$

(3.7)

Further, from standard EM approach, $H(\Theta\mid\Theta^{(m)})\leq H(\Theta^{(m)}\mid\Theta^{(m)})$ for any $\Theta$ in the parameter space. Thus, obtaining an updated estimate of the parameter by maximizing the $l_{1}$ penalized log-likelihood in (3.7) reduces to

$\displaystyle\Theta_{\lambda}^{({m+1})}=\textrm{argmax}_{\Theta}\left\{Q(\Theta\mid\Theta^{(m)})-\lambda\left\|\Theta\right\|_{1}\right\}.$

(3.8)

The EM optimization strategy alternates iteratively between the E-step, computing conditional expectation of the complete log-likelihood $Q(\Theta\mid\Theta^{(m)})$ and the M-step, maximizing $Q(\Theta\mid\Theta^{(m)})$, with a sparsity penalty $\lambda\left\|\Theta\right\|_{1}$, over $\Theta$.

E-step: The complete data likelihood depends on the joint distribution of $(Z,Z\in D(Y))$ given by

$\displaystyle p(Z,Z\in D(Y)\mid\Theta)=\left\{\begin{array}[pos]{ll}p(Z\mid\Theta)~~~~~~~~~~~~~~~~~~~~Z\in D(Y)\\ 0~~~~~~~~~~~~~~~~~~~~~~~~~~~~~Z\notin D(Y)\end{array}\right.$

where $p(Z\mid\Theta)$ is the multivariate normal density with mean zero and variance $\Sigma=\Theta^{-1}$. Then the complete log likelihood of $(Z,Z\in D(y))$ for a random sample of size $n$ after ignoring constants with respect to $\Theta$ is given by

	$\displaystyle l_{Z}(\Theta)$	$\displaystyle=$	$\displaystyle\sum_{i=1}^{n}\log\left(p(Z_{i}\mid\Theta)\right)I_{\{Z_{i}\in D(Y_{i})\}}$		(3.10)
		$\displaystyle=$	$\displaystyle-\frac{np}{2}\log(2\pi)+\frac{n}{2}\log\det(\Theta)-\frac{1}{2}\sum^{n}_{i=1}Z_{i}^{T}\Theta Z_{i}I_{\{Z_{i}\in D(Y_{i})\}}$

Using the complete log likelihood in (3.10), it follows that

	$\displaystyle Q(\Theta\mid\Theta^{(m)})$	$\displaystyle=$	$\displaystyle E\left[l_{Z}(\Theta)\mid Z\in D(Y),\Theta^{(m)}\right]$		(3.11)
		$\displaystyle=$	$\displaystyle\frac{n}{2}\left\{\log\det(\Theta)-\frac{1}{n}\sum_{i=1}^{n}\left(E\left[Z_{i}^{T}\Theta Z_{i}\mid Z_{i}\in D(Y_{i}),\Theta^{(m)}\right]\right)\right\}$
		$\displaystyle=$	$\displaystyle\frac{n}{2}\left\{\log\det(\Theta)-\mbox{Tr}\left(\Theta\frac{1}{n}\sum_{i=1}^{n}E\left[Z_{i}Z_{i}^{T}\mid Z_{i}\in D(Y_{i}),\Theta^{(m)}\right]\right)\right\}$
		$\displaystyle=$	$\displaystyle\frac{n}{2}\left\{\log\det(\Theta)-\mbox{Tr}\left(\Theta\bar{R}\right)\right\},$

where Tr stands for the trace of a matrix and $\bar{R}$ is the expected empirical covariance function of the latent variables given $Z\in D(Y)$:

$$\bar{R}=\frac{1}{n}\sum_{i=1}^{n}E\left[Z_{i}Z_{i}^{T}\mid z_{i}\in D(Y_{i}),\Theta^{(m)}\right],$$

where

	$\displaystyle E\left[Z_{i}Z_{i}^{T}\mid Z_{i}\in D(Y_{i}),\Theta^{(m)}\right]$	$\displaystyle=$	$\displaystyle E\left[Z_{i}\mid Z_{i}\in D(y_{i}),\Theta^{(m)}\right]E\left[Z_{i}\mid Z_{i}\in D(y_{i}),\Theta^{(m)}\right]^{T}$		(3.12)
			$\displaystyle+~~cov\left[Z_{i}\mid Z_{i}\in D(y_{i}),\Theta^{(m)}\right]$

Note that the conditional latent random variable $\left\{Z\mid Z\in D(Y)\right\}$ follows a truncated multivariate normal distribution. Analytical expressions for the computations of moments for truncated multivariate normal distribution are given in Wilhelm and Manjunath (2010) and references therein. However, due to the computational complexity, obtaining analytical solutions is only feasible for very few variables. Another approach is to simulate a large sample from the truncated multivariate normal distribution and calculate the sample conditional covariance matrix and sample conditional mean to estimate $E\left[Z_{i}Z_{i}^{T}\mid Z_{i}\in D(y_{i}),\Theta^{(m)}\right]$ using (3.12).

Alternatively, towards a computational efficient approach instead of mapping all mixed variables to a latent space as discussed above we may partition the mixed variables into two as continuous denoted by $Y_{c}$, and ordered variables that includes ordinal, binary and counts denoted by $Y_{o}$. We also partition the correlation matrix along with the variables grouping as

$\displaystyle\Sigma=\begin{bmatrix}\Sigma_{cc}&\Sigma_{co}\\ \Sigma_{oc}&\Sigma_{oo}\end{bmatrix}.$

We then take $Z=(Z_{c},Z_{o})\sim N\left(0,\Theta\right)$, where $\Theta=\Sigma^{-1}$ with $Z_{c}=\Phi^{-1}(\hat{F}(Y_{c}))$ is transformed normal scores of observed continuous variables using the rescaled empirical distribution based on the natural Gaussian copula semiparametric approach and $Z_{o}\in D(Y_{o})$ is the latent normal score corresponding to the ordered observed variables $Y_{o}$ obtained in a similar way as discussed in Section 2.1 .

With a similar argument as above the complete data likelihood depends on the joint distribution: $p(Z_{c},Z_{o},Z_{o}\in D(Y_{o})\mid\Theta)=p(Z_{c},Z_{o})$, for $Z_{o}\in D(Y_{o})$. Such that the complete data loglikelihood after ignoring constants is given by

	$\displaystyle l_{Z_{c},Z_{o}}(\Theta)$	$\displaystyle=$	$\displaystyle\sum_{i=1}^{n}\log\left(p(Z_{c_{i}},Z_{0_{i}}\mid\Theta)\right)I_{\{Z_{o}\in D(Y_{o})\}}$		(3.13)
		$\displaystyle=$	$\displaystyle\frac{n}{2}\log\det(\Theta)-\frac{1}{2}\sum^{n}_{i=1}\left[Z_{c_{i}},Z_{o_{i}}\right]^{T}\Theta\left[Z_{c_{i}},Z_{o_{i}}\right]I_{\{Z_{o_{i}}\in D(Y_{o_{i}})\}}$

Using this complete data log-likelihood and after ignoring constants with respect to $\Theta$, it follows that

	$\displaystyle Q(\Theta\mid\Theta^{(m)})$	$\displaystyle=$	$\displaystyle E_{Z_{o}}\left[l_{Z_{c},Z_{o}}(\Theta)\mid Y_{c},Z_{o}\in D(Y_{o}),\Theta^{(m)}\right]$
		$\displaystyle=$	$\displaystyle\frac{n}{2}\left\{\log\det(\Theta)\right.$
			$\displaystyle~~\left.-\mbox{Tr}\left(\Theta\frac{1}{n}\sum_{i=1}^{n}E_{Z_{o}}\left[\left[Z_{c_{i}},Z_{o_{i}}\right]\left[Z_{c_{i}},Z_{o_{i}}\right]^{T}\mid Y_{c_{i}},Z_{o_{i}}\in D(Y_{o_{i}}),\Theta^{(m)}\right]\right)\right\}$
		$\displaystyle=$	$\displaystyle\frac{n}{2}\left\{\log\det(\Theta)-\mbox{Tr}\left(\Theta\tilde{R}\right)\right\},$		(3.15)

where

$\displaystyle\tilde{R}=\frac{1}{n}\sum_{i=1}^{n}E_{Z_{o}}\left[\left[Z_{c_{i}},Z_{o_{i}}\right]\left[Z_{c_{i}},Z_{o_{i}}\right]^{T}\mid Y_{c_{i}},Z_{o_{i}}\in D(Y_{o_{i}}),\Theta^{(m)}\right].$

(3.16)

An estimate of $\tilde{R}$ can be obtained after evaluating the expectations as follows.

	$\displaystyle E_{Z_{o}}\left[Z_{c_{i}}Z_{c_{i}}^{T}\mid Y_{c_{i}},Z_{o_{i}}\in D(Y_{o_{i}}),\Theta^{(m)}\right]=Z_{c_{i}}Z_{c_{i}}^{T}$
	$\displaystyle E_{Z_{o}}\left[Z_{c_{i}}Z_{o_{i}}^{T}\mid Y_{c_{i}},Z_{o_{i}}\in D(Y_{o_{i}}),\Theta^{(m)}\right]=Z_{c_{i}}\hat{Z}_{o_{i}}^{T}$
	$\displaystyle E_{Z_{o}}\left[Z_{o_{i}}Z_{o_{i}}^{T}\mid Y_{c_{i}},Z_{o_{i}}\in D(Y_{o_{i}}),\Theta^{(m)}\right]=\hat{Z}_{o_{i}}\hat{Z}_{o_{i}}^{T}+\Theta^{{(m)}^{-1}}_{oo}.$

where $\hat{Z}_{o_{i}}^{T}=E_{Z_{o}}\left[Z_{o_{i}}\mid Z_{o_{i}}\in D(Y_{o_{i}}),\Theta^{(m)}\right]-\Theta^{{(m)}^{-1}}_{oo}\Theta^{(m)}_{oc}Z_{c_{i}}$, is a conditional expectation defined on the distribution $p(Z_{o}\mid Y_{c})$ for $Z_{o}\in D(Y_{o})$.

M-step: This involves updating the parameter $\Theta$ using the $l_{1}$ penalized log-likelihood given by

$\displaystyle\Theta_{\lambda}^{({m+1})}=\textrm{argmax}_{\Theta}\left\{Q(\Theta\mid\Theta^{(m)})-\lambda\left\|\Theta\right\|_{1}\right\}.$

(3.17)

We next substitute the $Q$ function by (3.11) or (3.15) from the E-step to obtain

$\displaystyle\Theta_{\lambda}^{({m+1})}=\textrm{argmax}_{\Theta}\left\{\log\det(\Theta)-\mbox{Tr}(\Theta\bar{R})-\lambda\left\|\Theta\right\|_{1}\right\}.$

(3.18)

The maximization problem in (3.18) takes the form of $l_{1}$ penalized likelihood for Gaussian graphical models and computation is done efficiently using the graphic lasso algorithm (Friedman et al., 2008). This algorithm is fast and allows the re-use of the estimate under one value of the tuning parameter as a “warm”’ start for the next value. The determination of a value for $\lambda$ in case of penalized inference with EM algorithm is discussed in Section 3.3.

The copula EM glasso approach discussed in the previous section, though it is a natural approach, it is computationally expensive, since it calls MCMC in the E-step and glasso in the M-step. In particular, in a very high dimensional setting the computational issue requires further attention. One approach is to seek a one-to-one mapping of the observed and latent variables that avoids the Monte Carlo EM algorithm resulting from a one-to-many mapping. The semiparametric copula approach of estimating marginals through the rescaled empirical distribution is a one-to-one mapping or transformation of the observed data. Instead of directly using the transformed data to estimate $\Theta$, a sample correlation matrix can be compute from pairwise rank correlations. In this regard, Liu et al. (2012) considered a nonparanormal skeptic approach to obtain sparse estimates of $\Theta$ for binary and continuous variables using one step glasso based on the estimated correlation matrix.

We note that the use a one step glasso approach is computationally efficient. Further, we note that rank correlations like Kendall’s tau and Spearman’s rho can be better approximated by a carefully chosen parametric bivariate copula model that takes in to account the underlying bivariate dependence structure. The vast literature on copulas deals with bivariate copula models and has demonstrated their potential to capture various types of dependence structure.

It is known, for example, that the population version of Kendall’s tau is related to parametric copula models parametrized by $\gamma_{ij}$ via

$\displaystyle\tau_{ij}=4\int^{1}_{0}\int^{1}_{0}C(u,v\mid\theta_{ij})dC(u,v\mid\gamma_{ij})-1.$

For commonly used copula models, there is closed form representation of the Kendall’s tau using the bivariate copula parameter, see for example Nelsen (2006). An estimate of Kendall’s tau is obtained using

$$\hat{\tau}_{ij}=\begin{cases}4\int^{1}_{0}\int^{1}_{0}C(u,v\mid\hat{\theta}_{ij})dC(u,v\mid\hat{\gamma}_{ij})-1&\text{for}\quad i\neq j\\ 1&\text{for}\quad i=j\end{cases}$$

or its closed form representation. Further Kendall’s tau is related to the correlation coefficient, $\Gamma$, by

$$\hat{\Gamma}_{ij}=\begin{cases}\sin\left(\frac{\pi}{2}\hat{\tau}_{ij}\right)&\text{for}\quad i\neq j\\ 1&\text{for}\quad i=j\end{cases}$$

Then to obtain sparse estimates, glasso can be implemented that uses the estimated correlation matrix $\hat{\Gamma}$ in the direct optimization of the objective function:

$\displaystyle\hat{\Theta}_{\lambda}=\textrm{argmax}_{\Theta}\left\{\log\det(\Theta)-\mbox{Tr}(\Theta\hat{\Gamma})-\lambda\left\|\Theta\right\|_{1}\right\}.$

For high dimensional data, the empirical covariance matrix is singular and poses computational problems. However, our $l_{1}$ penalized approach guarantees with probability one a positive definite precision matrix with the additional property of being sparse. Note that sparseness refers to the property that all parameters that are zero are actually estimated as zero with probability tending to one. This helps to assess conditional independence based on entries of the precision matrix (Dempster, 1972).

Under the $l_{1}$ penalized maximum likelihood setting the sparsity of the estimated precision matrix is controlled by the penalty parameter $\lambda$ in (3.18). We follow information based criteria in order to obtain reasonably sparse precision matrix. One could also use cross-validation for the choice of $\lambda$ which we have not consider in this article.

We now consider (3.3) that suggests the log likelihood of the observed data can be computed at EM convergence, see for example Ibrahim et al. (2008). Let the estimate $\widehat{\Theta}_{\lambda}$ is obtained at EM convergence for a given value of $\lambda$. The log likelihood of the observed data is

$$\log L_{Y}(\widehat{\Theta}_{\lambda})=Q(\widehat{\Theta}_{\lambda}\mid\widehat{\Theta}_{\lambda})-H(\widehat{\Theta}_{\lambda}\mid\widehat{\Theta}_{\lambda}).$$

Thus a model selection criterion is defined by

	$\displaystyle IC(\lambda)$	$\displaystyle=$	$\displaystyle-2\log L_{Y}(\widehat{\Theta}_{\lambda})+pen(\widehat{\Theta}_{\lambda})$
		$\displaystyle=$	$\displaystyle-2Q(\widehat{\Theta}_{\lambda}\mid\widehat{\Theta}_{\lambda})+2H(\widehat{\Theta}_{\lambda}\mid\widehat{\Theta}_{\lambda})+pen(\widehat{\Theta}_{\lambda}),$

where $pen(\widehat{\Theta}_{\lambda})$ refers to a penalty term. Different forms of $pen(\widehat{\Theta}_{\lambda})$ lead to different model selection criteria. Let $d$ denotes the number of non-zero upper or lower off-diagonal elements of $\widehat{\Theta}_{\lambda}$. Thus we define AIC and BIC as follows:

	$\displaystyle AIC(\lambda)$	$\displaystyle=$	$\displaystyle-2Q(\widehat{\Theta}_{\lambda}\mid\widehat{\Theta}_{\lambda})+2H(\widehat{\Theta}_{\lambda}\mid\widehat{\Theta}_{\lambda})+2d,$
	$\displaystyle BIC(\lambda)$	$\displaystyle=$	$\displaystyle-2Q(\widehat{\Theta}_{\lambda}\mid\widehat{\Theta}_{\lambda})+2H(\widehat{\Theta}_{\lambda}\mid\widehat{\Theta}_{\lambda})+\log(n)d.$

Then we choose the optimal value of the penalty parameter as the one that minimizes these criteria on a grid of candidate values for $\lambda$.

We carried out simulation studies with a variety of data structures to compare how well competing methods recover the true graph structure. Though our EM approach is computationally expensive, we noticed that in our simulations the EM algorithm converges very quickly with a maximum of ten iterations and 100 MCMC samples for hundreds of variables.

For the purpose of comparison we considered the following approaches:

• 1.

  Proposed copula EM glasso (CopulaEM).

• 2.

  Proposed copula skeptic glasso (CopulaTau)

• 3.

  Nonparanormal normal-score based estimation with truncation presented in Liu et al. (2012) (NPNscore)

• 4.

  Nonparanormal skeptic using Kendall’s tau presented in Liu et al. (2012) (NPNtau)

In our simulation we consider sample sizes (n=200) and number of variables(p=100) which are of mixed types that include binary(10), ordinal(10), count (10), nonnormal (eg. Chisquare (10)), and the remaining 60 are normal variables with outliers (none , 1% , 20%). In case of outliers, observations are replaced by a value either 5 or -5 with probability 0.6, see also Liu et al. (2012). ROC curves are used to compare performance of the different approaches in recovering the true graph.

Figures 1 and 2 displays ROC curves based on averages of true positive rates and false positive rates computed from 100 times repeated simulations at each of 10 grid points of the tuning parameter. For mixed data with no and low level outliers, we see that the difference in the performance of recovering the true graph based on the various methods is negligible though our copula EM glasso shows slightly better performance in case of no outliers.

In case of high level of outliers with mixed variables, the performance of the proposed copula skeptic and nonparanormal skeptic are comparable but the proposed copula skeptic out performs the nonparanormal skeptic. This suggests that a careful choice of parametric bivariate copulas results in better performance over the nonparametric approaches.