Estimating Gaussian Copulas with Missing Data with and without Expert Knowledge

Abstract

In this work, we present a rigorous application of the Expectation Maximization algorithm to determine the marginal distributions and the dependence structure in a Gaussian copula model with missing data. We further show how to circumvent a priori assumptions on the marginals with semiparametric modeling. Further, we outline how expert knowledge on the marginals and the dependency structure can be included. A simulation study shows that the distribution learned through this algorithm is closer to the true distribution than that obtained with existing methods and that the incorporation of domain knowledge provides benefits.

1. Introduction

Even though the amount of data is increasing due to new technologies, big data are by no means good data. For example, missing values are ubiquitous in various fields, from the social sciences [1] to manufacturing [2]. For explanatory analysis or decision making, one is often interested in the joint distribution of a multivariate dataset, and its estimation is a central topic in statistics [3]. At the same time, there exists background knowledge in many domains that can help to compensate for the potential shortcomings of datasets. For instance, domain experts have an understanding of the causal relationships in the data generation process [4]. It is the scope of this paper to unify expert knowledge and datasets with missing data to derive approximations of the underlying joint distribution.

To estimate the multivariate distribution, we use copulas, where the dependence structure is assumed to belong to a parametric family, while the marginals are estimated nonparametrically. Genest et al. [5] showed that for complete datasets, a two-step approach consisting of the estimation of the marginals with an empirical cumulative distribution function (ecdf) and subsequent derivation of the dependence structure is consistent. This idea is even transferable to high dimensions [6].

In the case of missing values, the situation becomes more complex. Here, nonparametric methods do not scale well with the number of dimensions [7]. On the other hand, assuming that the distribution belongs to a parametric family, it can often be derived by using the EM algorithm [8]. However, this assumption is, in general, restrictive. Due to the encouraging results for complete datasets, there have been several works that have investigated the estimation of the joint distribution under a copula model. The authors of [9,10] even discussed the estimation in a missing-not-at-random (MNAR) setting. While MNAR is less restrictive than missing at random (MAR), it demands the explicit modeling of the missing mechanism [11]. On the contrary, the authors of [12,13] provided results in cases in which data were missing completely at random (MCAR). This strong assumption is rarely fulfilled in practice. Therefore, we assume an MAR mechanism in what follows [11].

Another interesting contribution [14] assumed external covariates, such that the probability of a missing value depended exclusively on them and not on the variables under investigation. They applied inverse probability weighting (IPW) and the two-step approach of [5]. While they proved a consistent result, it is unclear how this approach can be adapted to a setting without those covariates. IPW for general missing patterns is computationally demanding, and no software exists [15,16]. Thus, IPW is mostly applied with monotone missing patterns that appear, for example, in longitudinal studies [17]. The popular work of [18] proposed an EM algorithm in order to derive the joint distribution in a Gaussian copula model with data MAR [11]. However, their approach had weaknesses:

• The presented algorithm was inexact. Among other things, the algorithm simplified by assuming that the marginals and the copula could be estimated separately (compare Equation (6) in [18] and Equation (11) in this paper).

• If there was no a priori knowledge of the parametric family of all marginals, Ref. [18] proposed using the ecdf of the observed data points. Afterwards, they exclusively derived the parameters of the copula. This estimator of the marginals was biased [19,20], which is often overlooked in the copula literature, e.g., [21] (Section 4.3), [22] (Section 3), [23] (Section 3), or [24] (Section 3).

• The description of the simulation study was incomplete and the results were not reproducible.

The aim of this paper is to close these gaps, and our contributions are the following:

• We give a rigorous derivation of the EM algorithm under a Gaussian copula model. Similarly to [5], it consists of two separate steps, which estimate the marginals and the copula, respectively. However, these two steps alternate.

• We show how prior knowledge about the marginals and the dependency structure can be utilized in order to achieve better results.

• We propose a flexible parametrization of the marginals when a priori knowledge is absent. This allows us to learn the underlying marginal distributions; see Figure 1.

• We provide a Python library that implements the proposed algorithm.

The structure of this paper is as follows. In Section 2, we review some background information about the Gaussian copula. We proceed by presenting the method (Section 3). In Section 4, we investigate its performance and the effect of domain knowledge in simulation studies. We conclude in Section 5. All technical aspects and proofs in this paper are given in Appendix A and Appendix B.

2. The Gaussian Copula Model

2.1. Notation and Assumptions

In the following, we consider a p-dimensional dataset $\{{\mathbf{x}}^1,\dots ,{\mathbf{x}}^N\}\subset {\mathbb{R}}^p$ of size N, where ${\mathbf{x}}^1,\dots ,{\mathbf{x}}^N$ are i.i.d. samples from a p-dimensional random vector $\mathbf{X}=\left(X_1,\dots ,X_p\right)$ with a joint distribution function F and marginal distribution functions $F_1,\dots ,F_p$. We denote the entries of ${\mathbf{x}}^{\ell }$ by ${\mathbf{x}}^{\ell }=\left(x_1^{\ell },\dots ,x_p^{\ell }\right)\forall \ell =1,\dots ,N$. The parameters of the marginals are represented by $\theta =\left({\theta }_1,\dots ,{\theta }_p\right)$, where ${\theta }_j$ is the parameter of $F_j$, so we write $F_j^{{\theta }_j}$, where ${\theta }_j$ can be a vector itself.

For $\ell \in \{1,\dots ,p\}$, we define $\mathbf{obs}(\ell )\subset \{1,\dots ,p\}$ as the index set of the observed and $\mathbf{mis}(\ell )\subset \{1,\dots ,p\}$ as the index set of the missing columns of ${\mathbf{x}}^{\ell }$. Hence, $\mathbf{mis}(\ell )\cup \mathbf{obs}(\ell )=\{1,\dots ,p\}$ and $\mathbf{mis}(\ell )\cap \mathbf{obs}(\ell )=\varnothing$. $\mathbf{R}=\left(R_1,\dots ,R_p\right)\in {\{0,1\}}^p$ is a random vector for which $R_i=0$ if $X_i$ is missing and $R_i=1$ if $X_i$ can be observed. Further, we define $\varphi$ to be the density function and $\Phi$ to be the distribution function of the one-dimensional standard normal distribution. ${\Phi }_{\mathit{\mu },\Sigma }$ stands for the distribution function of a p-variate normal distribution with covariance $\Sigma \in {\mathbb{R}}^{p\times p}$ and mean $\mathit{\mu }\in {\mathbb{R}}^p$. To simplify the notation, we define ${\Phi }_{\Sigma }:={\Phi }_{0,\Sigma }$. For a matrix $A\in {\mathbb{R}}^{p\times p}$, the entry of the i-th row and the j-th column is denoted by $A_{ij}$, while for index sets $\mathbf{S},\mathbf{T}\subset \{1,\dots ,p\}$, $A_{\mathbf{S},\mathbf{T}}$ is the submatrix of A with the row number in $\mathbf{S}$ and column number in $\mathbf{T}$. For a (random) vector $\mathbf{x}$ ($\mathbf{X}$), ${\mathbf{x}}_{\mathbf{S}}$ (${\mathbf{X}}_{\mathbf{S}}$) is the subvector containing entries with the index in $\mathbf{S}$.

Throughout, we assume F to be strictly increasing and continuous in every component. Therefore, $F_j$ is strictly increasing and continuous for all $j\in \{1,\dots ,p\}$, and so is the existing inverse function $F_j^{-1}$. For $\mathbf{S}=\{s_1,\dots ,s_k\}\subset \{1,\dots ,p\}$, we define $F_{\mathbf{S}}:{\mathbf{R}}^{\left|S\right|}\to {\mathbf{R}}^{\left|S\right|}$ by

$F_{\mathbf{S}}(x_{s_1},\dots ,x_{s_k})=\left(F_{s_1}\left(x_{s_1}\right),\dots ,F_{s_k}\left(x_{s_k}\right)\right).$

This work assumes that data are Missing at Random (MAR), as defined by [11], i.e.,

(1)${\mathbb{P}}_{\mathbf{X},\mathbf{R}}\left(\mathbf{R}=\mathbf{r}|{\mathbf{X}}_{\mathbf{r}}={\mathbf{x}}_{-\mathbf{r}},{\mathbf{X}}_{\mathbf{r}}={\mathbf{x}}_{\mathbf{r}}\right)={\mathbb{P}}_{\mathbf{X},\mathbf{R}}\left(\mathbf{R}=\mathbf{r}|{\mathbf{X}}_{\mathbf{r}}={\mathbf{x}}_{\mathbf{r}}\right),$

where ${\mathbf{X}}_{\mathbf{r}}:={\mathbf{X}}_{\{i:{\mathbf{r}}_i=1\}}$ are the observed and ${\mathbf{X}}_{-\mathbf{r}}:={\mathbf{X}}_{\{i:{\mathbf{r}}_i=0\}}$ are the missing entries of $\mathbf{X}$.

2.2. Properties

Sklar’s theorem [25] decomposes F into its marginals $F_1,\dots ,F_p$ and its dependency structure C with

(2)$F(x_1,\dots ,x_p)=C\left(F_1\left(x_1\right),\dots ,F_p\left(x_p\right)\right).$

Here, C is a copula, which means it is a p-dimensional distribution function with support ${[0,1]}^p$ whose marginal distributions are uniform. In this paper, we focus on Gaussian copulas, where

(3)$C_{\Sigma }(u_1,\dots ,u_p)={\Phi }_{\Sigma }\left({\Phi }^{-1}\left(u_1\right),\dots ,{\Phi }^{-1}\left(u_p\right)\right)$

and $\Sigma$ is a covariance matrix with ${\Sigma }_{jj}=1\forall j\in \{1,\dots ,p\}$. Beyond all multivariate normal distributions, there are distributions with non-normal marginals whose copula is Gaussian. Hence, the Gaussian copula model provides an extension of the normality assumption. Consider a random vector $\mathbf{X}$ whose copula is $C_{\Sigma }$. Under the transformation

$\mathbf{Z}:={\Phi }^{-1}\circ F\left(\mathbf{X}\right):=\left({\Phi }^{-1}\circ F_1\left(X_1\right),\dots ,{\Phi }^{-1}\circ F_p\left(X_p\right)\right),$

it holds that

(4)$\begin{array}{cc}\hfill F_{\mathbf{Z}}(z_1\dots ,z_p)& =\mathbb{P}\left(Z_1\le z_1,\dots ,Z_p\le z_p\right)\hfill \\ \hfill & =\mathbb{P}\left(X_1\le F_1^{-1}\left(\Phi \left(z_1\right)\right),\dots ,X_p\le F_p^{-1}\left(\Phi \left(z_p\right)\right)\right)\hfill \\ \hfill & ={\Phi }_{\Sigma }\left({\Phi }^{-1}\left(F_1\left(F_1^{-1}\left(\Phi \left(z_1\right)\right)\right)\right),\dots ,{\Phi }^{-1}\left(F_p\left(F_p^{-1}\left(\Phi \left(z_p\right)\right)\right)\right)\right)\hfill \\ \hfill & ={\Phi }_{\Sigma }(z_1,\dots ,z_p)\hfill \end{array}$

and hence, $\mathbf{Z}$ is normally distributed with mean 0 and covariance $\Sigma$. The two-step approaches given in [5,6] use this property and apply the following scheme:

• Find consistent estimates $\widehat{F_1},\dots ,\widehat{F_p}$ for the marginal distributions $F_1,\dots ,F_p$.

• Find $\Sigma$ by estimating the covariance of the random vector

  $\mathbf{Z}=\left({\Phi }^{-1}\left(\widehat{F_1}\left(X_1\right)\right),\dots ,{\Phi }^{-1}\left(\widehat{F_p}\left(X_p\right)\right)\right).$

From now on, we assume that the marginals of $\mathbf{X}$ have existing density functions $f_1,\dots ,f_p$. Then, by using Equation (4) and a change of variables, we can derive the joint density function

(5)$f_{F_1,\dots ,F_p,\Sigma }(x_1,\dots ,x_p)=f(x_1,\dots ,x_p)={|\Sigma |}^{-\frac{1}{2}}exp\left(-\frac{1}{2}{\mathbf{z}}^T\left({\Sigma }^{-1}-I\right)\mathbf{z}\right)\prod _{j=1}^p{f}_j\left(x_j\right),$

where $\mathbf{z}:=\left({\Phi }^{-1}\left(F_1\left(x_1\right)\right),\dots ,{\Phi }^{-1}\left(F_p\left(x_p\right)\right)\right)$. As for the multivariate normal distribution, we can identify the conditional independencies ([6]) from the inverse of the covariance matrix $K:={\Sigma }^{-1}$ by using the property

(6)$K_{jk}=K_{kj}=0\iff X_j\perp X_k|\left\{X_i:i\in \{1,\dots ,p\}\setminus \{j,k\}\right\}.$

K is called the precision matrix. In order to slim down the notation, we define

${\Phi }^{-1}\left(F_{\mathbf{S}}\left({\mathbf{x}}_{\mathbf{S}}\right)\right):=\left({\Phi }^{-1}\left(F_{s_1}\left(x_{s_1}\right)\right),\dots ,{\Phi }^{-1}\left(F_{s_k}\left(x_{s_k}\right)\right)\right)$

and similarly

$F_{\mathbf{S}}^{-1}\left(\Phi \left({\mathbf{z}}_{\mathbf{S}}\right)\right):=\left(F_{s_1}^{-1}\left(\Phi \left(z_{s_1}\right)\right),\dots ,F_{s_1}^{-1}\left(\Phi \left(z_{s_k}\right)\right)\right).$

The former function transforms the data of a Gaussian copula distribution to be normally distributed. The latter mapping takes multivariate normally distributed data and returns data following a Gaussian copula distribution with marginals $F_{s_1},\dots ,F_{s_k}$. The conditional density functions have a closed form.

Proposition 1

(Conditional Distribution of Gaussian Copula). Let $\mathbf{S}=\{s_1,\dots ,s_k\}$ and $\mathbf{T}=\{t_1,\dots ,t_{k^{\prime }}\}$ be such that $\mathbf{T}\dot{\cup }\mathbf{S}=\{1,\dots ,p\}$.

• The conditional density of ${\mathbf{X}}_{\mathbf{T}}|{\mathbf{X}}_{\mathbf{S}}={\mathbf{x}}_{\mathbf{S}}$ is given by

  $f\left({\mathbf{x}}_{\mathbf{T}}\right|{\mathbf{X}}_{\mathbf{S}}={\mathbf{x}}_{\mathbf{S}})={\left|{\Sigma }^{\prime }\right|}^{-\frac{1}{2}}exp\left(-\frac{1}{2}{({\mathbf{z}}_{\mathbf{T}}-\mathit{\mu })}^T{\Sigma }^{\prime -1}({\mathbf{z}}_{\mathbf{T}}-\mathit{\mu })\right)exp\left(\frac{1}{2}{\mathbf{z}}_{\mathbf{T}}^T{\mathbf{z}}_{\mathbf{T}}\right)\prod _{j\in \mathbf{T}}{f}_j\left(x_j\right),$

  where $\mathit{\mu }={\Sigma }_{\mathbf{T},\mathbf{S}}{\Sigma }_{\mathbf{S},\mathbf{S}}^{-1}{\mathbf{z}}_{\mathbf{S}}$, ${\Sigma }^{\prime }={\Sigma }_{\mathbf{T},\mathbf{T}}-{\Sigma }_{\mathbf{T},\mathbf{S}}{\Sigma }_{\mathbf{S},\mathbf{S}}^{-1}{\Sigma }_{\mathbf{S},\mathbf{T}}$, ${\mathbf{z}}_{\mathbf{T}}={\Phi }^{-1}\left(F_{\mathbf{T}}\left({\mathbf{x}}_{\mathbf{T}}\right)\right)$  and ${\mathbf{z}}_{\mathbf{S}}={\Phi }^{-1}\left(F_{\mathbf{S}}\left({\mathbf{x}}_{\mathbf{S}}\right)\right).$

• ${\Phi }^{-1}\left(F_{\mathbf{T}}\left({\mathbf{X}}_{\mathbf{T}}\right)\right)|{\mathbf{X}}_{\mathbf{S}}={\mathbf{x}}_{\mathbf{s}}$ is normally distributed with mean μ and covariance ${\Sigma }^{\prime }$.

• The expectation of $h\left({\mathbf{X}}_{\mathbf{T}}\right)$ with respect to the density $f\left({\mathbf{x}}_{\mathbf{T}}\right|{\mathbf{X}}_{\mathbf{S}}={\mathbf{x}}_{\mathbf{S}})$ can be expressed by

  $\int h\left({\mathbf{x}}_{\mathbf{T}}\right)f\left({\mathbf{x}}_{\mathbf{T}}\right|{\mathbf{X}}_{\mathbf{S}}={\mathbf{x}}_{\mathbf{S}})d{\mathbf{x}}_{\mathbf{T}}=\int h\left(F_{\mathbf{T}}^{-1}\left(\Phi \left({\mathbf{z}}_{\mathbf{T}}\right)\right)\right){\varphi }_{\mathit{\mu },{\Sigma }^{\prime }}\left({\mathbf{z}}_{\mathbf{T}}\right)d{\mathbf{z}}_{\mathbf{T}}.$

Proposition 1 shows that the conditional distribution’s copula is Gaussian as well. More importantly, we can derive an algorithm for sampling from the conditional distribution.

Algorithm 1:Sampling from the conditional distribution of a Gaussian copula

Input: ${\mathbf{x}}_{\mathbf{S}},\Sigma ,F_1,\dots ,F_p$ Result: m samples of ${\mathbf{X}}_{\mathbf{T}}|{\mathbf{X}}_{\mathbf{S}}={\mathbf{x}}_{\mathbf{S}}$ Calculate ${\mathbf{z}}_{\mathbf{S}}:={\Phi }^{-1}\left(F_{\mathbf{S}}\left({\mathbf{x}}_{\mathbf{S}}\right)\right)$ Calculate $\mathit{\mu }$ and ${\Sigma }^{\prime }$ as in Proposition 1 using ${\mathbf{z}}_{\mathbf{S}}$ and $\Sigma$ Draw samples $\{{\mathbf{z}}^1,\dots ,{\mathbf{z}}^m\}$ from $\mathcal{N}(\mathit{\mu },{\Sigma }^{\prime })$ return {$F_{\mathbf{T}}^{-1}\left(\Phi \left({\mathbf{z}}^1\right)\right),\dots ,F_{\mathbf{T}}^{-1}\left(\Phi \left({\mathbf{z}}^m\right)\right)\}$

The very last step follows with Proposition 1, as it holds for any measurable $A\subset {\mathbb{R}}^{k^{\prime }}$:

$\begin{array}{cc}\hfill \mathbb{P}\left({\mathbf{X}}_{\mathbf{T}}\in A|{\mathbf{X}}_{\mathbf{S}}={\mathbf{x}}_{\mathbf{S}}\right)& =\int 1_A\left({\mathbf{x}}_{\mathbf{T}}\right)f\left({\mathbf{x}}_{\mathbf{T}}\right|{\mathbf{X}}_{\mathbf{S}}={\mathbf{x}}_{\mathbf{S}})d{\mathbf{x}}_{\mathbf{T}}=\int 1_A\left(F_{\mathbf{T}}^{-1}\left(\Phi \left({\mathbf{z}}_{\mathbf{T}}\right)\right)\right){\varphi }_{\mathit{\mu },{\Sigma }^{\prime }}\left({\mathbf{z}}_{\mathbf{T}}\right)d{\mathbf{z}}_{\mathbf{T}}.\hfill \end{array}$

3. The EM Algorithm in the Gaussian Copula Model

3.1. The EM Algorithm

Let $\{{\mathbf{y}}^1,\dots ,{\mathbf{y}}^N\}\subset {\mathbb{R}}^p$ be a dataset following a distribution with parameter $\psi$ and corresponding density function $g_{\psi }\left(·\right)$, where observations are MAR. The EM algorithm [8] finds a local optimum of the log-likelihood function

$\begin{array}{cc}\hfill \sum _{\ell =1}^{N}ln\left(g_{\psi }\left({\mathbf{y}}_{\mathbf{obs}(\ell )}^{\ell }\right)\right)& =\sum _{\ell =1}^N\int ln\left(g_{\psi }\left(\left({\mathbf{y}}_{\mathbf{obs}(\ell )}^{\ell },{\mathbf{y}}_{\mathbf{mis}(\ell )}\right)\right)\right)\hfill \\ \hfill & g_{\psi }\left({\mathbf{y}}_{\mathbf{mis}(\ell )}|{\mathbf{Y}}_{\mathbf{obs}(\ell )}^{\ell }={\mathbf{y}}_{\mathbf{obs}(\ell )}^{\ell }\right)d{\mathbf{y}}_{\mathbf{mis}(\ell )}\hfill \\ \hfill & =\sum _{\ell =1}^N{\mathbb{E}}_{\psi }\left(ln\left(g_{\psi }\left(\left({\mathbf{y}}_{\mathbf{obs}(\ell )},{\mathbf{y}}_{\mathbf{mis}(\ell )}\right)\right)\right)|{\mathbf{Y}}_{\mathbf{obs}}^{\ell }={\mathbf{y}}_{\mathbf{obs}(\ell )}^{\ell }\right).\hfill \end{array}$

After choosing a start value ${\psi }^0$, it does so by iterating the following two steps.

• E-Step: Calculate

  (7)$\begin{array}{cc}\hfill \lambda \left(\psi \right|{\mathbf{y}}^1,\dots ,{\mathbf{y}}^N,{\psi }^t)& :=\sum _{\ell =1}^N{\mathbb{E}}_{{\psi }^t}\left(ln\left(g_{\psi }\left(\left({\mathbf{y}}_{\mathbf{obs}(\ell )},{\mathbf{y}}_{\mathbf{mis}(\ell )}\right)\right)\right)|{\mathbf{Y}}_{\mathbf{obs}}^{\ell }={\mathbf{y}}_{\mathbf{obs}(\ell )}^{\ell }\right)\hfill \\ \hfill & =\sum _{\ell =1}^N\lambda \left(\psi \right|{\mathbf{y}}^{\ell },{\psi }^t).\hfill \end{array}$

• M-Step: Set

  (8)${\psi }^{t+1}=\underset{\psi }{argmax}\lambda \left(\psi \right|{\mathbf{y}}^1,\dots ,{\mathbf{y}}^N,{\psi }^t)$

  and $t=t+1$.

For our purposes, there are two extensions of interest:

• If there is no closed formula for the right-hand side of Equation (7), one can apply Monte Carlo integration [26] as an approximation. This is called the Monte Carlo EM algorithm.

• If $\psi =\left({\psi }_1,\dots ,{\psi }_v\right)$ and the joint maximization of (8) with respect to $\psi$ is not feasible, Ref. [27] proposed a sequential maximization. Thus, we optimize (8) with respect to ${\psi }_i$ while holding ${\psi }_1={\psi }_1^{t+1},\dots ,{\psi }_{i-1}={\psi }_{i-1}^{t+1},{\psi }_{i+1}={\psi }_{i+1}^t,\dots ,{\psi }_v={\psi }_v^t$ fixed before we continue with ${\psi }_{i+1}$. This is called the Expectation Conditional Maximization (ECM) algorithm.

3.2. Applying the ECM Algorithm on the Gaussian Copula Model

As we need a full parametrization of the Gaussian copula model for the EM algorithm, we assume parametric marginal distributions $F_1^{{\theta }_1},\dots ,F_p^{{\theta }_p}$ with densities $f_1^{{\theta }_1},\dots ,f_p^{{\theta }_p}$. According to Equation (5), the joint density with respect to the parameters $\theta =\left({\theta }_1,\dots ,{\theta }_p\right)$ and $\Sigma$ has the form

(9)$f_{\theta ,\Sigma }(x_1,\dots ,x_p)={|\Sigma |}^{-\frac{1}{2}}exp\left(-\frac{1}{2}{\mathbf{z}}_{\theta }^T\left({\Sigma }^{-1}-I\right){\mathbf{z}}_{\theta }\right)\prod _{j=1}^p{f}_j^{{\theta }_j}\left(x_j\right),$

where ${\mathbf{z}}_{\theta }:=\left({\Phi }^{-1}\left(F_1^{{\theta }_1}\left(x_1\right)\right),\dots ,{\Phi }^{-1}\left(F_p^{{\theta }_p}\left(x_p\right)\right)\right)$. Section 3.3 will describe how we can keep the flexibility for the marginals despite the parametrization. However, first, we outline the EM algorithm for general parametric marginal distributions.

3.2.1. E-Step

Set $K:={\Sigma }^{-1}$ and $K^t:={\Sigma }^{t^{-1}}$. For simplicity, we pick one summand in Equation (7). By Equation (7) and (9), it holds with $\psi =\left(\theta ,\Sigma \right)$ and ${\mathbf{x}}^{\ell }$ taking the role of ${\mathbf{y}}^{\ell }$:

(10)$\begin{array}{cc}\hfill \lambda (\theta ,\Sigma |{\mathbf{x}}^{\ell },{\theta }^t,{\Sigma }^t)& ={\mathbb{E}}_{{\theta }^t,{\Sigma }^t}\left(ln\left(f_{\theta ,\Sigma }\left(\left({\mathbf{x}}_{\mathbf{obs}(\ell )},{\mathbf{x}}_{\mathbf{mis}(\ell )}\right)\right)\right)|{\mathbf{X}}_{\mathbf{obs}(\ell )}={\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell }\right)\hfill \\ \hfill & =-\frac{1}{2}ln\left(|\Sigma |\right)\hfill \\ \hfill & -\frac{1}{2}{\mathbb{E}}_{{\Sigma }^t,{\theta }^t}\left({{\mathbf{z}}_{\theta }}^T\left(K-I\right){\mathbf{z}}_{\theta }|{\mathbf{X}}_{\mathbf{obs}(\ell )}={\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell }\right)\hfill \\ \hfill & +\sum _{j=1}^p{\mathbb{E}}_{{\Sigma }^t,{\theta }^t}\left(ln\left(f_j^{{\theta }_j}\left(x_j\right)\right)|{\mathbf{X}}_{\mathbf{obs}(\ell )}={\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell }\right).\hfill \end{array}$

The first and last summand depend only on $\Sigma$ and $\theta$, respectively. Thus, of special interest is the second summand, for which we obtain the following with Proposition 1:

(11)${\mathbb{E}}_{{\Sigma }^t,{\theta }^t}\left({\mathbf{z}}_{\theta }^T\left(K-I\right){\mathbf{z}}_{\theta }|{\mathbf{X}}_{\mathbf{obs}(\ell )}={\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell }\right)=\int \left({\mathbf{z}}_{\theta ,{\theta }^t}^T\left(K-I\right){\mathbf{z}}_{\theta ,{\theta }^t}\right){\varphi }_{\mathit{\mu },{{\Sigma }^t}^{\prime }}\left({\mathbf{q}}_{\mathbf{mis}(\ell )}\right)d{\mathbf{q}}_{\mathbf{mis}(\ell )},$

where

${\mathbf{z}}_{\theta ,{\theta }^t}:=\left({\Phi }^{-1}\left(F_1^{{\theta }_1}\left(F_1^{{{\theta }_1^t}^{-1}}\left(\Phi \left(q_1\right)\right)\right)\right),\dots ,{\Phi }^{-1}\left(F_p^{{\theta }_p}\left(F_p^{{{\theta }_p^t}^{-1}}\left(\Phi \left(q_p\right)\right)\right)\right)\right).$

Here,

$\mathit{\mu }={\Sigma }_{\mathbf{mis}(\ell ),\mathbf{obs}(\ell )}{\Sigma }_{\mathbf{obs}(\ell ),\mathbf{obs}(\ell )}^{-1}{\Phi }^{-1}\left(F_{\mathbf{obs}(\ell )}^{{\theta }^t}\left(x_{\mathbf{obs}(\ell )}^{\ell }\right)\right)$

and

${{\Sigma }^t}^{\prime }={\Sigma }_{\mathbf{mis}(\ell ),\mathbf{mis}(\ell )}^t-{\Sigma }_{\mathbf{mis}(\ell ),\mathbf{obs}(\ell )}^t{\left({\Sigma }_{\mathbf{obs}(\ell ),\mathbf{obs}(\ell )}^t\right)}^{-1}{\Sigma }_{\mathbf{obs}(\ell ),\mathbf{mis}(\ell )}^t.$

At this point, the authors of [18] neglected that, in general,

$F_j^{{\theta }_j^t}\ne F_j^{{\theta }_j},j=1,\dots ,p$

holds, and hence, (11) depends not only on $\Sigma$, but also on $\theta$. This let us reconsider their approach, as we describe below.

3.2.2. M-Step

The joint optimization with respect to $\theta$ and $\Sigma$ is difficult, as there is no closed form for Equation (10). We circumvent this problem by sequentially optimizing with respect to $\Sigma$ and $\theta$ by applying the ECM algorithm. The maximization routine is the following.

• Set ${\Sigma }^{t+1}={argmax}_{\Sigma }{\sum }_{l=1}^N\lambda ({\theta }^t,\Sigma |{\mathbf{x}}^{\ell },{\theta }^t,{\Sigma }^t)$.

• Set ${\theta }^{t+1}={argmax}_{\theta }{\sum }_{l=1}^N\lambda (\theta ,{\Sigma }^{t+1}|{\mathbf{x}}^{\ell },{\theta }^t,{\Sigma }^t)$.

This is a two-step approach consisting of estimating the copula first and the marginals second. However, both steps are executed iteratively, which is typical for the EM algorithm.

Estimating $\Sigma$

As we are maximizing Equation (10) with respect to $\Sigma$ with a fixed $\theta ={\theta }^t$, the last summand can be neglected. By a change-of-variables argument, we show the following in Theorem A1:

${\mathbb{E}}_{{\Sigma }^t,{\theta }^t}\left({{\mathbf{z}}_{{\theta }^t}}^T\left(K-I\right){\mathbf{z}}_{{\theta }^t}|{\mathbf{X}}_{\mathbf{obs}(\ell )}={\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell }\right)=tr\left({\Sigma }^{-1}{V}_{\ell }\right),$

where $V_{\ell }$ depends on ${\Sigma }^t$ and ${\mathbf{z}}_{{\theta }^t,\mathbf{obs}(\ell )}={\Phi }^{-1}\left(F_{\mathbf{obs}(\ell )}^{{\theta }^t}\left({\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell }\right)\right)$. Thus, considering all observations, we search for

(12)$\begin{array}{cc}\hfill {\Sigma }^{t+1}& =\underset{\Sigma ,{\Sigma }_{\ell \ell }=1\forall \ell =1,\dots ,p}{argmax}\frac{1}{N}\sum _{l=1}^N\lambda ({\theta }^t,\Sigma |{\mathbf{x}}^{\ell },{\theta }^t,{\Sigma }^t)\hfill \\ \hfill & =\underset{\Sigma ,{\Sigma }_{\ell \ell }=1\forall \ell =1,\dots ,p}{argmax}\frac{1}{N}\sum _{\ell =1}^N-\frac{1}{2}ln\left(\right|\Sigma \left|\right)-\frac{1}{2}tr\left({\Sigma }^{-1}{V}_{\ell }\right)\hfill \\ \hfill & =\underset{\Sigma ,{\Sigma }_{\ell \ell }=1\forall \ell =1,\dots ,p}{argmax}-\frac{1}{2}ln\left(\right|\Sigma \left|\right)-\frac{1}{2}tr\left({\Sigma }^{-1}\frac{1}{N}\sum _{\ell =1}^N{V}_{\ell }\right),\hfill \end{array}$

which only depends on the statistic $S:=\frac{1}{N}{\sum }_{\ell =1}^N{V}_{\ell }$. Generally, this maximization can be formalized as a convex optimization problem that can be solved by a gradient descent. However, the properties of this estimator are not understood (for example, a scaling of S by $a\in {\mathbb{R}}_{>0}$ leads to a different solution; see Appendix A.3). To overcome this issue, we instead approximate the solution with the correlation matrix

$\underset{\Sigma ,{\Sigma }_{\ell \ell }=1\forall \ell =1,\dots ,p}{argmax}-\frac{1}{2}ln\left(\right|\Sigma \left|\right)-\frac{1}{2}tr\left({\Sigma }^{-1}S\right)\approx PSP,$

where $P\in {\mathbb{R}}^p$ is the diagonal matrix with entries $P_{jj}=\frac{1}{\sqrt{S_{jj}}},\forall j=1,\dots ,p$. This was also proposed in [28] (Section 2.2).

In cases in which there is expert knowledge on the dependency structure of the underlying distribution, one can adapt Equation (12) accordingly. We discuss this in more detail in Section 4.4.

Estimating $\theta$

We now focus on finding ${\theta }^{t+1}$, which is the maximizer of

$\begin{array}{cc}\hfill \sum _{\ell =1}^N\lambda (\theta ,{\Sigma }^{t+1}|{\mathbf{x}}^{\ell },{\theta }^t,{\Sigma }^t)& =\sum _{\ell =1}^N{\mathbb{E}}_{{\theta }^t,{\Sigma }^t}\left(ln\left(f_{\theta ,{\Sigma }^{t+1}}\left({\mathbf{x}}_{\mathbf{obs}(\ell )},{\mathbf{x}}_{\mathbf{mis}(\ell )}\right)\right)|{\mathbf{X}}_{\mathbf{obs}(\ell )}={\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell }\right)\hfill \\ \hfill & =\sum _{\ell =1}^N\int ln\left(f_{\theta ,{\Sigma }^{t+1}}\left({\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell },{\mathbf{x}}_{\mathbf{mis}(\ell )}\right)\right)\hfill \\ & f_{{\theta }^t,{\Sigma }^t}\left({\mathbf{x}}_{\mathbf{mis}(\ell )}|{\mathbf{X}}_{\mathbf{obs}(\ell )}={\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell }\right)d{\mathbf{x}}_{\mathbf{mis}(\ell )}\hfill \end{array}$

with respect to $\theta$. As there is, in general, no closed formula for the right-hand side, we use Monte Carlo integration. Again, we start by considering a single observation ${\mathbf{x}}^{\ell }$ to simplify terms. Employing Algorithm 1, we receive M samples ${\mathbf{x}}_{\mathbf{mis}(\ell ),1}^{\ell },\dots ,{\mathbf{x}}_{\mathbf{mis}(\ell ),M}^{\ell }$ from the distribution of $X_{\mathbf{mis}(\ell )}|X_{\mathbf{obs}(\ell )}=x_{\mathbf{obs}(\ell )}^{\ell }$ given the parameters ${\theta }^t$ and ${\Sigma }^t$. We set ${\mathbf{x}}_{\mathbf{obs}(\ell ),m}^{\ell }={\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell }\forall m=1,\dots ,M$. Then, by Equation (9),

(13)$\begin{array}{cc}\hfill \lambda (\theta ,{\Sigma }^{t+1}|{\mathbf{x}}^{\ell },{\theta }^t,{\Sigma }^t)\approx C+\frac{1}{M}\sum _{m=1}^M& -\frac{1}{2}{\left({\Phi }^{-1}\left(F_1^{{\theta }_1}\left(x_{1,m}^{\ell }\right)\right),\dots ,{\Phi }^{-1}\left(F_p^{{\theta }_p}\left(x_{p,m}^{\ell }\right)\right)\right)}^T\hfill \\ \hfill & \left(K^{t+1}-I\right)\hfill \\ \hfill & \left({\Phi }^{-1}\left(F_1^{{\theta }_1}\left(x_{1,m}^{\ell }\right)\right),\dots ,{\Phi }^{-1}\left(F_p^{{\theta }_p}(x_{p,m}^{\ell }\right)\right)\hfill \\ \hfill & +\sum _{j=1}^{p}ln\left(f_j^{{\theta }_j}\left(x_{j,m}^{\ell }\right)\right).\hfill \end{array}$

Hence, considering all observations, we set

(14)$\begin{array}{cc}\hfill {\theta }^{t+1}=\underset{\theta }{argmax}\frac{1}{M}\sum _{\ell =1}^N\sum _{m=1}^M& -\frac{1}{2}{\left({\Phi }^{-1}\left(F_1^{{\theta }_1}\left(x_{1,m}^{\ell }\right)\right),\dots ,{\Phi }^{-1}\left(F_p^{{\theta }_p}\left(x_{p,m}^{\ell }\right)\right)\right)}^T\hfill \\ \hfill & \left(K^{t+1}-I\right)\hfill \\ \hfill & \left({\Phi }^{-1}\left(F_1^{{\theta }_1}\left(x_{1,m}^{\ell }\right)\right),\dots ,{\Phi }^{-1}\left(F_p^{{\theta }_p}(x_{p,m}^{\ell }\right)\right)\hfill \\ \hfill & +\sum _{j=1}^{p}ln\left(f_j^{{\theta }_j}\left(x_{j,m}^{\ell }\right)\right).\hfill \end{array}$

Note that we only use the Monte Carlo samples to update the parameters of the marginal distributions $\theta$. We would also like to point out some interesting aspects about Equations (13) and (14):

• The summand ${\sum }_{\ell =1}^N{\sum }_{m=1}^{M}ln\left(f_j^{{\theta }_j}\left(x_{j,m}^{\ell }\right)\right)$ describes how well the marginal distributions fit the (one-dimensional) data.

• The estimations of the marginals are interdependent. Hence, in order to maximize with respect to ${\theta }_j$, we have to take into account all other components of $\theta$.

• The first summand adjusts for the dependence structure in the data. If all observations at step $t+1$ are assumed to be independent, then $K^{t+1}=I$, and this term is 0.

• More generally, the derivative $\frac{\partial \lambda (\theta ,{\Sigma }^{t+1}|{\mathbf{x}}^{\ell },{\theta }^t,{\Sigma }^t)}{\partial {\theta }_j}$ depends on ${\theta }_k$ if and only if $K_{jk}^{t+1}\ne 0$. This means that if ${\Sigma }^{t+1}$ implies the conditional independence of column j and k given all other columns (Equation (6)), the optimal ${\theta }_j$ can be found without considering ${\theta }_k$. This, e.g., is the case if we set entries of the precision matrix to 0. Thus, the incorporation of prior knowledge reduces the complexity of the identification of the marginal distributions.

The intuition behind the derived EM algorithm is simple. Given a dataset with missing values, we estimate the dependency structure. With the identified dependency structure, we can derive likely locations of the missing values. Again, these locations help us to find a better dependency structure. This leads to the proposed cyclic approach. The framework of the EM algorithm guarantees the convergence of this procedure to a local maximum for $M\to \infty$ in Equation (14).

3.3. Modelling with Semiparametric Marginals

In the case in which the missing mechanism is MAR, the estimation of the marginal distribution using only complete observations is biased. Even worse, any moment of the distribution can be distorted. Thus, one needs a priori knowledge in order to identify the parametric family of the marginals [19,20]. If their family is known, one can directly apply the algorithm of Section 3.2. If this is not the case, we propose the use of a mixture model parametrization of the form

(15)$F_j^{{\theta }_j}\left(x_j\right)=\frac{1}{g}\sum _{k=1}^g\Phi \left(\frac{x_j-{\theta }_{jk}}{{\sigma }_j}\right),{\theta }_{j1}\le \dots \le {\theta }_{jg},\forall j=1,\dots ,p,$

where ${\sigma }_j$ is a hyperparameter and the ordering of the ${\theta }_{jk}$ ensures the identifiability.

Using mixture models for density estimation is a well-known idea (e.g., [29,30,31]). As the authors of [31] noted, mixture models vary between being parametric and being non-parametric, where flexibility increases with g. It is reasonable to choose Gaussian mixture models, as their density functions are dense in the set of all density functions with respect to the $L^1$-norm [29] (Section 3.2). This flexibility and the provided parametrization make the mixture models a natural choice.

3.4. A Blueprint of the Algorithm

The complete algorithm is summarized in Algorithm 2. For the Monte Carlo EM algorithm, Ref. [26] proposed the stabilization of the parameters with a rather small number of samples M and to increase this number substantially in the latter steps of the algorithm. This seems to be reasonable for line 2 of Algorithm 2 as well.

If there is no a priori knowledge about the marginals, we propose that we follow Section 3.3. We choose the initial ${\theta }^0$ such that the cumulative distribution function of the mixture model fits the ecdf of the observed data points. For an empirical analysis of the role of g, see Section 4.3.3. For ${\sigma }_1,\dots ,{\sigma }_p$, we use a rule of thumb inspired by [3] and set

${\sigma }_j=1.06\frac{\widehat{{\sigma }_j}}{g^{1/5}},$

where $\widehat{{\sigma }_j}$ is the standard deviation of the observed data points in the j-th component.

Algorithm 2:Blueprint for the EM algorithm for the Gaussian copula model

4. Simulation Study

We analyze the performance of the proposed estimator in two studies. First, we consider scenarios for two-dimensional datasets and check the potential of the algorithm. In the second part, we explore how expert knowledge can be incorporated and how this affects the behavior and performance. The proposed procedure, which is indexed with EM in the figures below, is compared with:

• Standard COPula Estimator (SCOPE): The marginal distributions are estimated by the ecdf of the observed data points. This was proposed by [18] if the parametric family is unknown, and it is the state-of-the art approach. Thus, we apply an EM algorithm to determine the correlation structure on the mapped data points

  $z_j^{\ell }={\Phi }^{-1}\left(\widehat{F_j}\left(x_j^{\ell }\right)\right),\ell =1,\dots ,N,j=1,\dots ,p,$

  where $\widehat{F_j}$ is the ecdf of the observed data points in column j. Its corresponding results are indexed with SCOPE in the figures and tables.

• Known marginals: The distribution of the marginals is completely known. The idea is to eliminate the difficulty of finding them. Here, we apply the EM algorithm for the correlation structure on

  $z_j^{\ell }={\Phi }^{-1}\left(F_j\left(x_j^{\ell }\right)\right),\ell =1,\dots ,N,j=1,\dots ,p,$

  where $F_j$ is the real marginal distribution function. Its corresponding results are indexed with a 0 in the figures and tables.

• Markov chain–Monte Carlo (MCMC) approach [21]: The author proposed an MCMC scheme to estimate the copula in a Bayesian fashion. Therefore, Ref. [21] derived the distribution of the multivariate ranks. The marginals are treated as nuisance parameters. We employed the R package sbgcop, which is available on CRAN, as it provides not only a posterior distribution of the correlation matrix $\Sigma$, but also imputations for missing values. In order to compare the approach with the likelihood-based methods, we set

  ${\widehat{\Sigma }}_{MCMC}=\frac{1}{M}\sum _{m=1}^M{\Sigma }^m,$

  where $\{{\Sigma }^m:m=1,\dots ,M\}$ are samples of the posterior distribution of the correlation matrix. For the marginals, we defined

  ${\widehat{F}}_{j,MCMC}\left(x\right)=\frac{1}{MN}\sum _{\ell =1}^N\sum _{m=1}^M{1}_{\{x_{j,m}^{\ell }\le x\}},$

  where $x_{j,m}^{\ell }$ is the m-th of the total of M imputations for $x_j^{\ell }$ and $x_{j,m}^{\ell }=x_j^{\ell }\forall m=1,\dots ,M$ if $x_j^{\ell }$ can be observed. The samples were drawn from the posterior distribution. The corresponding results were indexed with the MCMC approach in the figures and tables.

Sklar’s theorem shows that the joint distribution can be decomposed into the marginals and the copula. Thus, we analyze them separately.

4.1. Adapting the EM Algorithm

In Section 4.3 and Section 4.4, we chose $g=15$, for which we saw a sufficient flexibility. A sensitivity analysis of the procedure with respect to g can be found in Section 4.3.3. The initial ${\theta }^0$ was chosen by fitting the marginals to the existing observations, and ${\Sigma }^0$ was the identity matrix. For the number of Monte Carlo samples M, we observed that with $M=20$, $\theta$ stabilized after around 10 steps. Cautiously, we ran 20 steps before we increased M to 1000, for which we run another five steps. We stopped the algorithm when the condition $\parallel {\Sigma }^{t+1}-{\Sigma }^t{\parallel }_1<{10}^{-5}$ was fulfilled.

4.2. Data Generation

We considered a two-dimensional dataset (we would have liked to include the setup of the simulation study of [18]; however, neither could the missing mechanism be extracted from the paper nor did the authors provide it on request) with a priori unknown marginals $F_1$ and $F_2$, whose copula was Gaussian with the correlation parameter $\rho \in [-1,1]$. The marginals were chosen to be ${\chi }^2$ with six and seven degrees of freedom. The data matrix $D\in {\mathbb{R}}^{N\times 2}$ kept N (complete) observations of the random vector. We enforced the following MAR mechanism:

• Remove every entry in D with probability $0\le p_{MCAR}<1$. We denote the resulting data matrix (with missing entries) as $D^{MCAR}={\left(D_{\ell j}^{MCAR}\right)}_{\ell =1,\dots ,N,j=1,2}$.

• If $D_{\ell 1}^{MCAR}$ and $D_{\ell 2}^{MCAR}$ are observed, remove $D_{\ell 2}^{MCAR}$ with probability

  $\begin{array}{cc}\hfill \mathbb{P}\left(R_2=0|X_1=D_{\ell 1},X_2=D_{\ell 2}\right)& =\mathbb{P}\left(R_2=0|X_1=D_{\ell 1}\right)\hfill \\ \hfill & =\frac{1}{1+exp\left(-{\beta }_0-{\beta }_1{\Phi }^{-1}\left(F_1\left(D_{\ell 1}\right)\right)\right)}\hfill \end{array}$

  We call the resulting data matrix $D^{MAR}$.

The missing patterns were non-monotone. Aside from $p_{MCAR}$, the parameters ${\beta }_0$ and ${\beta }_1$ controlled how many entries were absent in the final dataset. Assuming that $\rho >0$, ${\beta }_1>0$, and $|{\beta }_0|$ was not too large, the ecdf of the observed values of $X_2$ was shifted to the left compared to the true distribution function (changing the signs of ${\beta }_1$ and/or $\rho$ may change the direction of the shift, but the situation is analogous). This can be seen in Figure 1, where we chose $N=200$, $\rho =0.5$, $\beta =({\beta }_0,{\beta }_1)=(0,2)$. The marginal distribution of $X_1$ could be estimated well by the ecdf of the observed data.

4.3. Results

This subsection explores how different specifications of the data-generating process presented in Section 4.2 influenced the estimation of the joint distribution. First, we investigate the influence of the share of missing values (controlled via $\beta$) and the dependency (controlled via $\rho$) by fixing the number of observations (denoted by N) to 100. Then, we vary N to study the behavior of the algorithms for larger sample sizes. Afterwards, we carry out a sensitivity analysis of the EM algorithm with respect to g, the number of mixtures. Finally, we study the computational demands of the algorithms.

4.3.1. The Effects of Dependency and Share of Missing Values

We investigate two different choices for the setup in Section 4.2 by setting the parameters to $\rho =0.1$, $\beta =(-1,1)$ and $\rho =0.5$, $\beta =(0,2)$. For both, we draw 1000 datasets with $N=100$ each and apply the estimators. To evaluate the methods, we look at two different aspects.

First, we compare the estimators for $\rho$ with respect to bias and standard deviations. The results are depicted in the corresponding third columns of Table 1 and are summarized as boxplots in Figure A1 in Appendix B.3. We see that no method is clearly superior. While the EM algorithm has a stronger bias for $\rho =0.5$ than that of SCOPE, it also has a smaller standard deviation. The MCMC approach shows the largest bias. As even known marginals (${\rho }_0$) do not lead to substantially better estimators compared to SCOPE (${\rho }_{SCOPE}$) or the proposed (${\rho }_{EM}$) approach, we deduce that (at least in this setting) the estimators for the marginals are almost negligible. MCMC performs notably worse.

Second, we investigate the Cramer–von Mises statistics $\omega$ between the estimated and the true marginal distribution (${\omega }^1$ statistic for the first marginal, ${\omega }^2$ statistic for the second marginal). The results are shown in Table 1 (corresponding first two columns) and are summarized as boxplots in Figure A2 in Appendix B.3. While for $\rho =0.1$, the proposed estimator behaves only slightly better than SCOPE, we see that the benefit becomes larger in the case of high correlation and more missing values, especially when estimating the second marginal. This is in line with the intuition that if the correlation is vanishing, the two random variables $X_1$ and $X_2$ become independent. Thus, $R_2$, the missing value indicator, and $X_2$ become independent. (Note that there is a difference from the case in which $\rho \ne 0$, and hence, the missingness probability $R_2$ isconditionally independent from $X_2$ given $X_1$.) In that case, we can estimate the marginal of $X_2$ using the ecdf of the observed data points. Hence, SCOPE’s estimates of the marginals should be good for small values of $\rho$. An illustration can be found in Figure 2. Again, the MCMC approach performs the worst.

4.3.2. Varying the Sample Size N

To investigate the behavior of the methods for larger sample sizes, we repeat the experiment from Section 4.2 with $\rho =0.5,\beta =(0,2)$ for the sample sizes $N=100,200,500,1000$. The results are depicted in Table 2 and Figure A3, Figure A4 and Figure A5 in Appendix B.3. The bias of SCOPE and EM algorithm for $\rho$ seem to vanish for large N, while the MCMC approach remains biased. Studying the estimation of the true marginals, the approximation of the second marginal via MCMC and SCOPE improves only slowly and is still poor for the largest sample sizes $N=1000$. In contrast, the EM algorithm performs best in small sample sizes, and the mean (of ${\omega }^1$ and ${\omega }^2$) and standard deviations (of all three values) move towards 0 for increasing N.

4.3.3. The Impacts of Varying the Number of Mixtures g

The proposed EM algorithm relies on the hyperparameter g, the number of mixtures in Equation (15). To analyze the behavior of the EM algorithm with respect to g, we additionally run the EM algorithm with $g=5$ and $g=30$ on the 1000 datasets of Section 4.2 for $\rho =0.5$, $\beta =(0,2)$, and $N=100$. We did not adjust the number of steps in the EM algorithm to keep the results comparable. The results can be found in Table 3. We see that the choice of g does not have a large effect on the estimation of $\rho$. However, an increased g leads to better estimates for $X_1$. This is in line with the intuition that the ecdf of the first components is an unbiased estimate for the distribution function of $X_1$, and setting g to the number of samples corresponds to the kernel density estimator. On the other hand, the estimator for $X_2$ benefits slightly from $g=5$, as ${\omega }_{EM}^2$ has a lower mean and standard deviation compared to the choice $g=30$. However, this effect is small and almost non-existent when we compare $g=5$ with $g=15$. As the choice $g=15$ leads to better estimates of the first marginal compared to $g=5$, we see this choice as a good compromise for our setting. For applications without prior knowledge, we recommend considering g as additional tuning parameter (via cross-validation).

4.3.4. Run Time

We analyze the computational demands of the different algorithms by comparing their run times in the study of Section 4.3.1 with $\rho =0.5$ and $\beta =(0,2)$ (the settings $\rho =0.1$ and $\beta =(-1,1)$ lead to similar results and are omitted). The run times of all presented algorithms depend not only on the dataset, but also on the parameters (e.g., convergence criterion and ${\Sigma }^0$ for SCOPE). Thus, we do not aim for an extensive study, but focus on the magnitudes. We compare the proposed EM algorithm with a varying number of mixtures ($g=5,15,30$) with MCMC and SCOPE. The results are shown in Table 4. We see that the EM algorithm has the longest run time, which depends on the number of mixtures g. The MCMC approach and the proposed EM algorithm have a higher computational demand than SCOPE, as they are trying to model the interaction between the copula and the marginals. As mentioned in the onset, we could reduce the run time of the EM algorithm by going down to only 10 steps instead of 20.

4.4. Inclusion of Expert Knowledge

In the presence of prior knowledge on the dependency structure, the presented EM algorithm is highly flexible. While information on the marginals can be used to parametrize the copula model, expert knowledge on the dependency structure can be incorporated by adapting Equation (12). In the case of soft constraints on the covariance or precision matrix, one can replace Equation (12) with a penalized covariance estimation, where the penalty reflects the expert assessment [32,33]. Similarly, one can define a prior distribution on the covariance matrices and set ${\Sigma }^{t+1}$ as the mode of the posterior distribution (the MAP estimate) of $\Sigma$ given the statistic S of Equation (12).

Another possibility could be that we are aware of conditional independencies in the data-generating process. This is, for example, the case when causal relationships are known [4]. To exemplify the latter, we consider a three-dimensional dataset $\mathbf{X}$ with the Gaussian copula $C_{\Sigma }$ and marginals $X_1,X_2,X_3$, which are ${\chi }^2$ distributed with six, seven, and five degrees of freedom. The precision is set to

$K={\Sigma }^{-1}={\Delta }^{1/2}\left(\begin{array}{ccc}1& 0.5& 0.5\\ 0.5& 1& 0\\ 0.5& 0& 1\end{array}\right){\Delta }^{1/2},$

where ${\Delta }^{1/2}$ is a diagonal matrix, which ensures that the diagonal elements of $\Sigma$ are 1. We see that $X_2$ and $X_3$ are conditionally independent given $X_1$. The missing mechanism is similar to the one in Section 4.2. The missingness of $X_3$ depends on $X_1$ and $X_2$, while the probability of a missing $X_1$ or $X_2$ is independent of the others. The mechanism is, again, MAR. Details can be found in Appendix B.2. We compare the proposed method with prior knowledge on the zeros in the precision matrix (indexed by KP, EM in the figures) with the EM, SCOPE, and MCMC algorithms without background knowledge. We again sample 1000 datasets with 50 observations each from the real distribution. The background knowledge on the precision is used by restricting the non-zero elements in Equation (12). Therefore, we apply the procedure presented in [34] (Chapter 17.3.1) to find ${\Sigma }^{t+1}$. The means and standard deviations of the estimates are presented in Table 5.

First, we evaluate the estimated dependency structures by calculating the Frobenius norm of the estimation error $\Sigma -\widehat{\Sigma }$. The EM algorithm with background knowledge (KP, EM) performs best and is more stable than its competitors. Apart from MCMC, the other procedures behave similarly, which indicates again that the exact knowledge of the marginal distributions is not too relevant for identifying the dependency structure. MCMC performs the worst.

Second, we see that the proposed EM estimators return marginal distributions that are closer to the truth, while the estimate with background knowledge (KP, EM) performs the best. Thus, the background knowledge on the copula also transfers into better estimates for the marginal distribution—in particular, for $X_3$. This is due to Equation (14) and the comments thereafter. The zeros in the precision structure indicate which other marginals are relevant in order to identify the parameter of a marginal. In our case, $X_2$ provides no additional information for $X_3$. This information is provided to the EM algorithm through the restriction of the precision matrix.

Finally, we compare the EM estimates of the joint distribution. The relative entropy or Kullback–Leibler divergence is a popular tool for estimating the difference between two distributions [35,36], where one of them is absolutely continuous with respect to the other. A lower number indicates a higher similarity. Due to the discrete structure of the marginals of SCOPE and MCMC, we cannot calculate their relative entropy with respect to the truth. However, we would like to analyze how the estimate of the proposed procedure improves if we include expert knowledge. The results are depicted in Table 6. Again, we observe that the incorporation of extra knowledge improves the estimates. This is in line with Table 5, as the estimation of all components in the joint distribution of Equation (3) is improved by the domain knowledge.

5. Discussion

In this paper, we investigated the estimation of the Gaussian copula and the marginals with an incomplete dataset, for which we derived a rigorous EM algorithm. The procedure iteratively searches for the marginal distributions and the copula. It is, hence, similar to known methods for complete datasets. We saw that if the data are missing at random, a consistent estimate of a marginal distribution depends on the copula and other marginals.

The EM algorithm relies on a complete parametrization of the marginals. The parametric family of the marginals is, in general, a priori unknown and cannot be identified through the observed data points. For this case, we presented a novel idea of employing mixture models. Although this is practically always a misspecification, our simulation study revealed that the combination of our EM algorithm and marginal mixture models delivers better estimates for the joint distribution than currently used procedures do. In principle, uncertainty quantification of the parameters derived by the proposed EM algorithm can be achieved by bootstrapping [37].

There are different possibilities for incorporating expert knowledge. Information on the parametric family of the marginals can be used for their parametrization. However, causal and structural understandings of the data-generating process can also be utilized [4,38,39]. For example, this can be achieved by restricting the correlation matrix or its inverse, the precision matrix. We presented how one can restrict the non-zero elements of the precision, which enforces conditional independencies. Our simulation study showed that this leads not only to an improved estimate for the dependency structure, but also to better estimates for the marginals. This translates into a lower relative entropy between the real distribution and the estimate. We also discussed how soft constraints on the dependency structure can be included.

We note that the focus of this paper is on estimating the joint distribution without precise specification of its subsequent use. Therefore, we did not discuss imputation methods (see, e.g., [40,41,42,43]). However, Gaussian copula models were employed as a device for multiple imputation (MI) with some success [22,24,44]. The resulting complete datasets can be used for inference. All approaches that we are aware of estimate the marginals by using the ecdf of the observed data points. The findings in Section 4 translate into better draws for the missing values.

Additionally, the joint distribution can be utilized for regressing a potentially multivariate $\mathbf{Y}$ on $\mathbf{Z}$ even if data are missing. By applying the EM algorithm on $\mathbf{X}:=\left(\mathbf{Y},\mathbf{Z}\right)$ and by Proposition 1, one even obtains the whole conditional distribution of $\mathbf{Y}$ given $\mathbf{Z}=\mathbf{z}$.

We have shown how to incorporate a causal understanding of the data-generating process. However, in the potential outcome framework of [45], the derivation of a causal relationship can also be interpreted as a missing data problem in which the missing patterns are “misaligned” [46]. Our algorithm is applicable for this.

figuresection tablesection

Author Contributions

Conceptualization, M.K.; methodology, M.K.; software, M.K.; validation, M.P. and M.K.; formal analysis, M.K.; investigation, M.K.; resources, M.P. and M.K.; data curation, M.K.; writing—original draft preparation, M.K.; writing—review and editing, M.P. and M.K; visualization, M.K.; supervision, M.P.; project administration, M.P. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data generation procedures of the simulation studies and the proposed algorithm are availabe at https://github.com/mkrtl/misscop, accessed on 26 October 2022.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Estimates of the proposed EM algorithm ([Forumla omitted. See PDF.], orange line), the Standard Copula Estimator ([Forumla omitted. See PDF.], blue line, corresponds to ecdf), the Markov chain–Monte Carlo approach ([Forumla omitted. See PDF.], purple line) for the marginals [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and the truth ([Forumla omitted. See PDF.], green line) of a two-dimensional example dataset generated as described in Section 4.2 with [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.].

Figure 2. Dependency graph for [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.]. [Forumla omitted. See PDF.] is independent of [Forumla omitted. See PDF.] if either [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] are independent ([Forumla omitted. See PDF.]) or if [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] are independent ([Forumla omitted. See PDF.]).

Table 1

Comparison of the algorithms with respect to the Cramer–von Mises distance between the estimated and the true first (${\omega }^1$) and true second marginal distributions (${\omega }^2$), as well as the correlation ($\rho$). Shown are the mean and standard deviation of the proposed EM algorithm (EM), the method based on known marginals (0), the Standard Copula Estimator (SCOPE), and the Markov chain–Monte Carlo approach (MCMC) for 1000 datasets generated as described in Section 4.3.1.

		Mean			Standard Deviation
Setting	Method	${\mathbf{\omega }}^{\mathbf{1}}$	${\mathbf{\omega }}^{\mathbf{2}}$	$\mathbf{\rho }$	${\mathbf{\omega }}^{\mathbf{1}}$	${\mathbf{\omega }}^{\mathbf{2}}$	$\mathbf{\rho }$
$\rho =0.1,\beta =(-1,1)$	EM	8.55	10.41	0.107	9.30	11.67	0.139
	0	-	-	0.109	-	-	0.144
	SCOPE	9.13	12.25	0.105	8.47	11.00	0.144
	MCMC	18.21	24.99	0.094	16.62	21.89	0.127
$\rho =0.5,\beta =(0,2)$	EM	8.03	16.48	0.455	8.68	19.47	0.139
	0	-	-	0.498	-	-	0.138
	SCOPE	9.06	45.25	0.486	8.25	36.11	0.143
	MCMC	17.90	59.34	0.393	16.13	57.15	0.131

Table 2

Comparison of the algorithms with respect to the Cramer–von Mises distance between the estimated and the true first (${\omega }^1$) and true second marginal distributions (${\omega }^2$), as well as the correlation ($\rho$). Shown are the mean and standard deviation of the proposed EM algorithm (EM), the method based on known marginals (0), the Standard Copula Estimator (SCOPE), and the Markov chain–Monte Carlo approach (MCMC) for 1000 datasets generated as described in Section 4.2 with $\rho =0.5$ and $\beta =(0,2)$ and varying sample sizes $N=100,200,500,1000$.

		Mean			Standard Deviation
N	Method	${\mathbf{\omega }}^{\mathbf{1}}$	${\mathbf{\omega }}^{\mathbf{2}}$	$\mathbf{\rho }$	${\mathbf{\omega }}^{\mathbf{1}}$	${\mathbf{\omega }}^{\mathbf{2}}$	$\mathbf{\rho }$
$N=100$	EM	8.03	16.48	0.455	8.68	19.47	0.139
	0	-	-	0.498	-	-	0.138
	SCOPE	9.06	45.25	0.486	8.25	36.11	0.143
	MCMC	17.90	59.34	0.393	16.13	57.15	0.131
$N=200$	EM	4.91	8.53	0.469	5.46	8.88	0.098
	0	-	-	0.500	-	-	0.094
	SCOPE	4.76	37.38	0.493	4.18	25.35	0.096
	MCMC	9.27	42.91	0.370	8.01	36.23	0.089
$N=500$	EM	3.01	3.83	0.480	2.92	3.59	0.063
	0	-	-	0.499	-	-	0.060
	SCOPE	2.05	31.92	0.497	1.85	14.95	0.060
	MCMC	4.01	31.41	0.0360	3.49	20.51	0.051
$N=1000$	EM	2.25	2.74	0.486	1.92	2.40	0.047
	0	-	-	0.500	-	-	0.042
	SCOPE	1.08	30.60	0.499	0.93	11.13	0.043
	MCMC	1.99	28.13	0.365	1.84	14.49	0.037

Table 3

Comparison of the proposed EM algorithm with respect to the Cramer–von Mises distance between the estimated and the true first (${\omega }^1$) and true second marginal distributions (${\omega }^2$), as well as the correlation ($\rho$), for different numbers of mixtures g in Equation (15). Shown are the mean and standard deviation for $g=5,15,30$ and for 1000 datasets generated as described in Section 4.2 with $\rho =0.5$ and $\beta =(0,2)$.

	Mean			Standard Deviation
# Mixtures	${\mathbf{\omega }}^{\mathbf{1}}$	${\mathbf{\omega }}^{\mathbf{2}}$	$\mathbf{\rho }$	${\mathbf{\omega }}^{\mathbf{1}}$	${\mathbf{\omega }}^{\mathbf{2}}$	$\mathbf{\rho }$
$g=5$	13.82	16.38	0.469	14.17	19.69	0.145
$g=15$	8.03	16.48	0.455	8.68	19.47	0.139
$g=30$	7.17	18.73	0.454	7.48	20.98	0.140

Table 4

Comparison of the algorithms with respect to the run time in seconds. Shown are the mean and standard deviation of the proposed EM algorithm (EM) with the number of mixtures g set to $5,15,30$, the Standard Copula Estimator (SCOPE), and the Markov chain–Monte Carlo approach (MCMC) for 1000 datasets generated as described in Section 4.2 with $\rho =0.5$ and $\beta =(0,2)$.

	Run Time in Seconds
Method	Mean	Standard Deviation
EM ($g=5$)	21.78	3.27
EM ($g=15$)	55.94	11.39
EM ($g=30$)	161.57	38.00
SCOPE	0.45	0.11
MCMC	12.98	0.87

Table 5

Comparison of the algorithms with respect to the Cramer–von Mises distance between the estimated and the true first marginal distribution (${\omega }^1$), true second marginal distribution (${\omega }^2$), and true third marginal distribution (${\omega }^3$), as well as the correlation ($\rho$). Shown are the mean and standard deviation of the proposed EM algorithm (EM), the proposed EM algorithm with prior knowledge on the conditional independencies (KP, EM), the method based on known marginals (0), the Standard Copula Estimator (SCOPE), and the Markov chain–Monte Carlo approach (MCMC) for 1000 datasets generated as described in Section 4.4.

	Mean				Standard Deviation
Method	${\mathbf{\omega }}^{\mathbf{1}}$	${\mathbf{\omega }}^{\mathbf{2}}$	${\mathbf{\omega }}^{\mathbf{3}}$	$\left|\right|\widehat{\Sigma }-\Sigma {\left|\right|}_{\mathbf{2}}$	${\mathbf{\omega }}^{\mathbf{1}}$	${\mathbf{\omega }}^{\mathbf{2}}$	${\mathbf{\omega }}^{\mathbf{3}}$	$\left|\right|\widehat{\Sigma }-\Sigma {\left|\right|}_{\mathbf{2}}$
EM	12.12	13.38	21.15	0.229	13.89	14.25	22.44	0.113
KP, EM	12.04	13.28	19.66	0.182	13.93	14.37	20.88	0.111
0	-	-	-	0.227	-	-	-	0.108
SCOPE	17.57	17.55	26.69	0.232	16.75	15.55	24.84	0.113
MCMC	36.85	35.70	80.22	0.263	32.82	33.24	78.57	0.140

Table 6

Comparison of the algorithms with respect to the Kullback–Leibler divergence ($D_{KL}$) between the true joint distribution (F) and the estimates. Shown are the mean and standard deviation of the proposed EM algorithm (EM) and the proposed EM algorithm with prior knowledge on the conditional independencies (KP, EM) for 1000 datasets generated as described in Section 4.4.

	Mean(${\mathit{D}}_{\mathbf{KL}}(\mathit{F},\mathit{·})$)	Standard Deviation(${\mathit{D}}_{\mathbf{KL}}(\mathit{F},\mathit{·})$)
EM	1.37	0.53
KP, EM	1.26	0.32

Appendix A. Technical Results

Appendix A.1. Proof of Conditional Distribution

Proof of Proposition 1.

We prove in the order of the proposition, which is a multivariate generalization of [47].

• 1.

  We inspect the conditional density function: $$\begin{array}{cc}\hfill f\left({\mathbf{x}}_{\mathbf{T}}\right|{\mathbf{X}}_{\mathbf{S}}={\mathbf{x}}_{\mathbf{S}})& =\frac{{|\Sigma |}^{-\frac{1}{2}}exp\left(-\frac{1}{2}{\mathbf{z}}^T\left({\Sigma }^{-1}-I\right)\mathbf{z}\right){\prod }_{j=1}^p{f}_j\left(x_j\right)}{|{\Sigma }_{\mathbf{S},\mathbf{S}}{|}^{-\frac{1}{2}}exp\left(-\frac{1}{2}{\mathbf{z}}_{\mathbf{S}}^T\left({\Sigma }_{\mathbf{S},\mathbf{S}}^{-1}-I\right){\mathbf{z}}_{\mathbf{S}}\right){\prod }_{j\in \mathbf{S}}{f}_j\left(x_j\right)}\hfill \\ \hfill & =\frac{{|\Sigma |}^{-\frac{1}{2}}exp\left(-\frac{1}{2}{\mathbf{z}}^T{\Sigma }^{-1}\mathbf{z}\right)exp\left(\frac{1}{2}{\mathbf{z}}^T\mathbf{z}\right){\prod }_{j=1}^p{f}_j\left(x_j\right)}{|{\Sigma }_{\mathbf{S},\mathbf{S}}{|}^{-\frac{1}{2}}exp\left(-\frac{1}{2}{\mathbf{z}}_{\mathbf{S}}^T{\Sigma }_{\mathbf{S},\mathbf{S}}^{-1}{\mathbf{z}}_{\mathbf{S}}\right)exp\left(\frac{1}{2}{\mathbf{z}}_{\mathbf{S}}^T{\mathbf{z}}_{\mathbf{S}}\right){\prod }_{j\in \mathbf{S}}{f}_j\left(x_j\right)}\hfill \\ \hfill & =\frac{{|\Sigma |}^{-\frac{1}{2}}exp\left(-\frac{1}{2}{\mathbf{z}}^T{\Sigma }^{-1}\mathbf{z}\right)exp\left(\frac{1}{2}{\mathbf{z}}_{\mathbf{T}}^T{\mathbf{z}}_{\mathbf{T}}\right){\prod }_{j\in \mathbf{T}}{f}_j\left(x_j\right)}{|{\Sigma }_{\mathbf{S},\mathbf{S}}{|}^{-\frac{1}{2}}exp\left(-\frac{1}{2}{\mathbf{z}}_{\mathbf{S}}^T{\Sigma }_{\mathbf{S},\mathbf{S}}^{-1}{\mathbf{z}}_{\mathbf{S}}\right)}\hfill \end{array}$$

  Using well-known factorization lemmas and using the Schur complement (see, for example, [48] (Section 4.3.4)) applied on ${\Sigma }^{-1}$, we encounter (A1)$$f\left({\mathbf{x}}_{\mathbf{T}}\right|{\mathbf{X}}_{\mathbf{S}}={\mathbf{x}}_{\mathbf{S}})={\left|{\Sigma }^{\prime }\right|}^{-\frac{1}{2}}exp\left(-\frac{1}{2}{({\mathbf{z}}_{\mathbf{T}}-\mathit{\mu })}^T{\Sigma }^{\prime -1}({\mathbf{z}}_{\mathbf{T}}-\mathit{\mu })\right)exp\left(\frac{1}{2}{\mathbf{z}}_{\mathbf{T}}^T{\mathbf{z}}_{\mathbf{T}}\right)\prod _{j\in \mathbf{T}}{f}_j\left(x_j\right).$$

• 2.

  The distribution of $${\Phi }^{-1}\left(F_{\mathbf{T}}\left({\mathbf{X}}_{\mathbf{T}}\right)\right)|{\mathbf{X}}_{\mathbf{S}}={\mathbf{x}}_{\mathbf{S}}$$ follows with a change-of-variable argument. Using Equation (A1), we observe for any measurable set A that $$\begin{array}{cc}\hfill & \mathbb{P}\left(\left({\Phi }^{-1}\left(F_{\mathbf{T}}\left({\mathbf{X}}_{\mathbf{T}}\right)\right)|{\mathbf{X}}_{\mathbf{S}}={\mathbf{x}}_{\mathbf{s}}\right)\in A\right)\hfill \\ \hfill & ={\int }_{F^{-1}\left(\Phi \left(A\right)\right)}{\left|{\Sigma }^{\prime }\right|}^{-\frac{1}{2}}exp\left(-\frac{1}{2}{({\mathbf{z}}_{\mathbf{T}}-\mathit{\mu })}^T{\Sigma }^{\prime -1}({\mathbf{z}}_{\mathbf{T}}-\mathit{\mu })\right)exp\left(\frac{1}{2}{\mathbf{z}}_{\mathbf{T}}^T{\mathbf{z}}_{\mathbf{T}}\right)\prod _{j\in \mathbf{T}}{f}_j\left(x_j\right)d{\mathbf{x}}_{\mathbf{T}}\hfill \\ \hfill & ={\int }_A{\varphi }_{\mathit{\mu },{\Sigma }^{\prime }}\left({\mathbf{q}}_{\mathbf{T}}\right)d{\mathbf{q}}_{\mathbf{T}},\hfill \end{array}$$ where, in the second equation, we used the transformation ${\mathbf{q}}_{\mathbf{T}}={\Phi }^{-1}\left(F_{\mathbf{T}}\left({\mathbf{x}}_{\mathbf{T}}\right)\right)$ and the fact that $$\left|D\left({\varphi }^{-1}\left(F_{\mathbf{T}}\left({\mathbf{x}}_{\mathbf{T}}\right)\right)\right)\right|=2{\pi }^{\frac{\left|\mathbf{T}\right|}{2}}exp\left(\frac{1}{2}{\left({\Phi }^{-1}\left(F_{\mathbf{T}}\left({\mathbf{x}}_{\mathbf{T}}\right)\right)\right)}^T\left({\Phi }^{-1}\left(F_{\mathbf{T}}\left({\mathbf{x}}_{\mathbf{T}}\right)\right)\right)\right)\prod _{j\in \mathbf{T}}{f}_j\left(x_j\right).$$

• 3.

  This proof is analogous to the one above, and we finally obtain $$\begin{array}{c}\hfill \int h\left({\mathbf{x}}_{\mathbf{T}}\right)f\left({\mathbf{x}}_{\mathbf{T}}\right|{\mathbf{X}}_{\mathbf{S}}={\mathbf{x}}_{\mathbf{S}})d{\mathbf{x}}_{\mathbf{T}}=\int h\left(F^{-1}\left(\Phi \left({\mathbf{z}}_{\mathbf{T}}\right)\right)\right){\varphi }_{\mathit{\mu },{\Sigma }^{\prime }}\left({\mathbf{z}}_{\mathbf{T}}\right)d{\mathbf{z}}_{\mathbf{T}}.\end{array}$$

The result can be generalized to the case in which $\mathbf{S}\cup \mathbf{T}\ne \{1,\dots ,p\}$. □

Appendix A.2. Closed-Form Solution of the E-Step for θ = θ t

Theorem A1.

We assume w.l.o.g. that ${\mathbf{x}}^{\ell }=({\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell },{\mathbf{x}}_{\mathbf{mis}(\ell )}^{\ell })$ and set $${\mathbf{z}}_{{\theta }^t}:=\left({\mathbf{z}}_{\mathbf{obs}(\ell ),{\theta }^t},{\mathbf{z}}_{\mathbf{mis}(\ell )}\right):=\left({\Phi }^{-1}\left(F_{\mathbf{obs}(\ell )}^{{\theta }^t}\left({\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell }\right)\right),{\mathbf{z}}_{\mathbf{mis}(\ell )}\right).$$

Then, it holds that $${\mathbb{E}}_{{\Sigma }^t,{\theta }^t}\left({{\mathbf{z}}_{{\theta }^t}}^T{\Sigma }^{-1}{\mathbf{z}}_{{\theta }^t}|{\mathbf{X}}_{\mathbf{obs}(\ell )}={\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell }\right)=tr\left({\Sigma }^{-1}V\right),$$

where $V=\left(\begin{array}{cc}{\mathbf{z}}_{\mathbf{obs}(\ell ),{\theta }^t}{{\mathbf{z}}_{\mathbf{obs}(\ell ),{\theta }^t}}^T& {\mathbf{z}}_{\mathbf{obs}(\ell ),{\theta }^t}{\mathit{\mu }}^T\\ \mathit{\mu }{{\mathbf{z}}_{\mathbf{obs}(\ell ),{\theta }^t}}^T& {\Sigma }^{\prime }+\mathit{\mu }{\mathit{\mu }}^T\end{array}\right),$ $\mathit{\mu }={\Sigma }_{\mathbf{mis}(\ell ),\mathbf{obs}(\ell )}^t{{\Sigma }_{\mathbf{obs}(\ell ),\mathbf{obs}(\ell )}^t}^{-1}{\mathbf{z}}_{\mathbf{obs}(\ell ),{\theta }^t}$ and ${\Sigma }^{\prime }={\Sigma }_{\mathbf{mis}(\ell ),\mathbf{mis}(\ell )}^t-{\Sigma }_{\mathbf{mis}(\ell ),\mathbf{obs}(\ell )}^t{{\Sigma }_{\mathbf{obs}(\ell ),\mathbf{obs}(\ell )}^t}^{-1}{\Sigma }_{\mathbf{obs}(\ell ),\mathbf{mis}(\ell )}^t$.

Proof. 

We define $F_{{\theta }^t}\left({\mathbf{x}}_{\mathbf{mis}(\ell )}\right):=F_{{\theta }^t}({\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell },{\mathbf{x}}_{\mathbf{mis}(\ell )})$. Then, $$\begin{array}{cc}\hfill & {\mathbb{E}}_{{\Sigma }^t,{\theta }^t}\left({\mathbf{z}}_{{\theta }^t}^T{\Sigma }^{-1}{\mathbf{z}}_{{\theta }^t}|{\mathbf{X}}_{\mathbf{obs}(\ell )}={\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell }\right)\hfill \\ \hfill & ={\mathbb{E}}_{{\Sigma }^t,{\theta }^t}\left({\left({\Phi }^{-1}\left(F_{{\theta }^t}\left({\mathbf{x}}_{\mathbf{mis}(\ell )}\right)\right)\right)}^T{\Sigma }^{-1}\left({\Phi }^{-1}\left(F_{{\theta }^t}\left({\mathbf{x}}_{\mathbf{mis}(\ell )}\right)\right)\right)|{\mathbf{X}}_{\mathbf{obs}(\ell )}={\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell }\right)\hfill \\ \hfill & =\int {\left({\Phi }^{-1}\left(F_{{\theta }^t}\left({\mathbf{x}}_{\mathbf{mis}(\ell )}\right)\right)\right)}^T{\Sigma }^{-1}\left({\Phi }^{-1}\left(F_{{\theta }^t}\left({\mathbf{x}}_{\mathbf{mis}(\ell )}\right)\right)\right)\hfill \\ \hfill & f_{{\theta }^t,{\Sigma }^t}\left({\mathbf{x}}_{\mathbf{mis}(\ell )}|{\mathbf{X}}_{\mathbf{obs}(\ell )}={\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell }\right)d{\mathbf{x}}_{\mathbf{mis}(\ell )}.\hfill \end{array}$$

We now apply Proposition 1 and encounter $$\begin{array}{cc}\hfill & \int {\left({\Phi }^{-1}\left(F_{{\theta }^t}\left({\mathbf{x}}_{\mathbf{mis}(\ell )}\right)\right)\right)}^T{\Sigma }^{-1}{\Phi }^{-1}\left(F_{{\theta }^t}\left({\mathbf{x}}_{\mathbf{mis}(\ell )}\right)\right)\hfill \\ \hfill & f_{{\theta }^t,{\Sigma }^t}\left({\mathbf{x}}_{\mathbf{mis}(\ell )}|{\mathbf{X}}_{\mathbf{obs}(\ell )}={\mathbf{x}}_{\mathbf{obs}(\ell )}^{\ell }\right)d{\mathbf{x}}_{\mathbf{mis}(\ell )}\hfill \\ \hfill & =\int {\mathbf{z}}_{{\theta }^t}^T{\Sigma }^{-1}{\mathbf{z}}_{{\theta }^t}{\varphi }_{{\Sigma }^{\prime },\mathit{\mu }}\left({\mathbf{z}}_{\mathbf{mis}(\ell )}\right)d{\mathbf{z}}_{\mathbf{mis}(\ell )}\hfill \\ \hfill & =\int tr\left({\mathbf{z}}_{{\theta }^t}{\mathbf{z}}_{{\theta }^t}^T{\Sigma }^{-1}\right){\varphi }_{{\Sigma }^{\prime },\mathit{\mu }}\left({\mathbf{z}}_{\mathbf{mis}(\ell )}\right)d{\mathbf{z}}_{\mathbf{mis}(\ell )}\hfill \\ \hfill & =tr\left({\Sigma }^{-1}\int {\mathbf{z}}_{{\theta }^t}{\mathbf{z}}_{{\theta }^t}^T{\varphi }_{{\Sigma }^{\prime },\mathit{\mu }}\left({\mathbf{z}}_{\mathbf{mis}(\ell )}\right)d{\mathbf{z}}_{\mathbf{mis}(\ell )}\right).\hfill \end{array}$$

The last integral is understood element-wise. By the first and second moment of ${\Phi }_{{\Sigma }^{\prime },\mathit{\mu }}$, it follows that $$\begin{array}{cc}\hfill \int {\mathbf{z}}_{{\theta }^t}{\mathbf{z}}_{{\theta }^t}^T{\varphi }_{{\Sigma }^{\prime },\mathit{\mu }}\left({\mathbf{z}}_{\mathbf{mis}(\ell ),{\theta }^t}\right)d{\mathbf{z}}_{\mathbf{mis}(\ell ),{\theta }^t}& =\int \left({\mathbf{z}}_{\mathbf{obs}(\ell ),{\theta }^t},{\mathbf{z}}_{\mathbf{mis}(\ell ),{\theta }^t}\right){\left({\mathbf{z}}_{\mathbf{obs}(\ell ),{\theta }^t},{\mathbf{z}}_{\mathbf{mis}(\ell ),{\theta }^t}\right)}^T\hfill \\ \hfill & {\varphi }_{{\Sigma }^{\prime },\mathit{\mu }}\left({\mathbf{z}}_{\mathbf{mis}(\ell ),{\theta }^t}\right)d{\mathbf{z}}_{\mathbf{mis}(\ell ),{\theta }^t}\hfill \\ \hfill & =\left(\begin{array}{cc}{\mathbf{z}}_{\mathbf{obs}(\ell ),{\theta }^t}{\mathbf{z}}_{\mathbf{obs}(\ell ),{\theta }^t}^T& {\mathbf{z}}_{\mathbf{obs}(\ell ),{\theta }^t}{\mathit{\mu }}^T\\ \mathit{\mu }{\mathbf{z}}_{\mathbf{obs}(\ell ),{\theta }^t}^T& {\Sigma }^{\prime }+\mathit{\mu }{\mathit{\mu }}^T\end{array}\right).\hfill \end{array}$$

□

Appendix A.3. Maximizer of argmax Σ,Σ jj =1∀j=1,…,p λ(θ t, Σ|θ t, Σ t )

We are interested in $$\underset{{\Sigma }_{jj}=1\forall j=1,\dots ,p}{argmax}l(\Sigma ):=\underset{{\Sigma }_{jj}=1\forall j=1,\dots ,p}{argmax}-log\left(|\Sigma |\right)-tr\left({\Sigma }^{-1}S\right),$$ where $\Sigma ,S\in {\mathbb{R}}^{p\times p}$ are positive definite matrices. Clearly, $${\Sigma }_{jj}=1\iff 1=e_j^T\Sigma e_j=tr\left(e_j^T\Sigma e_j\right)=tr\left(e_j{e}_j^T\Sigma \right).$$

Using the Lagrangian, we obtain the following function to optimize $$L(\Sigma ,\lambda )=-log\left(|\Sigma |\right)-tr\left({\Sigma }^{-1}S\right)+\sum _{j=1}^p{\lambda }_j\left(tr\left(e_j{e}_j^T\Sigma \right)-1\right).$$

Applying the identities $\frac{\partial tr\left(AX\right)}{\partial X}=A$, $\frac{\partial tr\left(A{X}^{-1}\right)}{\partial X}=-X^{-1}A{X}^{-1}$, and $\frac{\partial log\left(\right|X\left|\right)}{\partial X}=X^{-1}$, we obtain the derivative with respect to $\Sigma$:$$\frac{\partial L}{\partial \Sigma }=-{\Sigma }^{-1}+{\Sigma }^{-1}S{\Sigma }^{-1}-\left(\sum _{j=1}^p{\lambda }_j\left(e_j{e}_j^T\right)\right)\stackrel{!}{=}0.$$

This is equivalent to $$-K+KSK=D_{\lambda },$$ where $D_{\lambda }$ is the diagonal matrix with entries $\lambda =\left({\lambda }_1,\dots ,{\lambda }_p\right)$ and $K:={\Sigma }^{-1}$. We see that the scaling of S by $a\in {\mathbb{R}}_{>0}$ leads, in general, to a different solution K, and hence, the estimator is not invariant under strictly monotone linear transformations of S.

We can also formulate the task as a convex optimization problem:$$\underset{{\left(K^{-1}\right)}_{ii}=1\forall i=1,\dots ,p}{argmin}-log\left(\left|K\right|\right)+tr\left(KS\right).$$

Appendix B. Details of the Simulation Studies

Appendix B.1. Drawing Samples from the Joint Distributions

Appendix B.1.1. Estimators of the Percentile Function

• In the case of SCOPE, consider the marginal observed data points, which we assume to be ordered $y_1\le \dots \le y_N$. We use the following linearly interpolated estimator for the percentile function: $$\widehat{F^{-1}}\left(u\right)=\left\{\begin{array}{cc}{y}_1\hfill & \hfill \mathrm{for}u\le \frac{1}{N+1}\\ y_N,\hfill & \hfill \mathrm{for}u>\frac{N}{N+1}\\ \frac{u-\frac{i}{N+1}}{\frac{i+1}{N+1}-\frac{i}{N+1}}(y_{i+1}-y_i)+y_i,\hfill & \hfill \mathrm{for}u\in \left(\frac{i}{N+1},\frac{i+1}{N+1}\right]\end{array}\right\}$$

• To estimate the percentile function for the mixture models, we choose with equal probability (all Gaussians have equal weight) one component of the mixture and then draw a random number with its mean ${\theta }_{jk}$ and standard deviation ${\sigma }_j$, $j=1,\dots ,p,k=1,\dots ,g$. In this manner, we generate $N^{\prime }$ samples $y_1^{\prime },\dots ,y_{N^{\prime }}^{\prime }$. The estimator for the percentile function is then chosen to be analogous to the one above. A higher $N^{\prime }$ leads to a more exact result. We choose $N^{\prime }$ to be 10,000.

Appendix B.1.2. Sampling

Given an estimator $\widehat{\rho }$ and estimators for the percentile functions $\widehat{F_1^{-1}},\widehat{F_2^{-1}}$, we obtain M samples from the learned joint distribution with $$y_{\ell }=\left(y_{\ell 1},y_{\ell 2}\right)=\left(\widehat{F_1^{-1}}\left(u_{\ell 1}\right),\widehat{F_2^{-1}}\left(u_{\ell 2}\right)\right)=\left(\widehat{F_1^{-1}}\left(\Phi \left(z_{\ell 1}\right)\right),\widehat{F_2^{-1}}\left(\Phi \left(z_{\ell 2}\right)\right)\right),\ell =1,\dots ,M,$$ where $z_{\ell }=\left(z_{\ell 1},z_{\ell 2}\right),\ell =1,\dots ,M$ are draws from a bivariate normal distribution with mean 0 and covariance $\left(\begin{array}{cc}1& \widehat{\rho }\\ \widehat{\rho }& 1\end{array}\right)$. In the case of the gold standard, we set $\widehat{F_j^{-1}}=F_j^{-1},j=1,2$. We obtain samples of the real underlying distribution by using the correct percentile functions, as in the gold standard, and, additionally, $\widehat{\rho }=\rho$. The procedure for three dimensions is analogous.

Appendix B.2. Missing Mechanism for Section 4.4

The missing mechanism is similar to the two-dimensional case. The marginals are chosen to be ${\chi }^2$ with six, seven, and five degrees of freedom. The data matrix $D\in {\mathbb{R}}^{N\times 3}$ keeps N (complete) observations of the random vector. We enforce the following missing data mechanism:

• Again, we remove every entry in the data matrix D with probability $0\le p_{MCAR}<1$. The resulting data matrix (with missing entries) is denoted as $$D^{MCAR}={\left(D_{\ell j}^{MCAR}\right)}_{\ell =1,\dots ,N,j=1,2,3}.$$

• If $D_{\ell 1}^{MCAR}$, $D_{\ell 2}^{MCAR}$, and $D_{\ell 3}^{MCAR}$ are observed, we remove $D_{\ell 3}^{MCAR}$ with probability $$\mathbb{P}\left(R_3=0|X_1=D_{\ell 1},X_2=D_{\ell 2}\right)=h(D_{\ell 1},D_{\ell 2};\beta ),$$ where $$h(D_{\ell 1},D_{\ell 2};\beta )=\frac{1}{1+exp\left(-\left({\beta }_0+{\beta }_1{\Phi }^{-1}\left(F_1\left(D_{\ell 1}\right)\right)+{\beta }_2{\Phi }^{-1}\left(F_2\left(D_{\ell 2}\right)\right)\right)\right)}$$ and $\beta =\left({\beta }_0,{\beta }_1,{\beta }_2\right)$.

We call the resulting data matrix $D^{MAR}$. Its missing patterns are, again, non-monotone, and the data are MAR, but not MCAR. In Section 4.4, we set $\beta =(0,2,2)$.

Appendix B.3. Complementary Figures

Figure A1. Comparison of the algorithms with respect to the correlation [Forumla omitted. See PDF.]. Shown are the boxplots for the Standard Copula Estimator (SCOPE), the proposed EM algorithm (EM), the method based on known marginals (0), and the Markov chain–Monte Carlo approach (MCMC) for 1000 datasets generated as described in Section 4.2, where [Forumla omitted. See PDF.] are depicted in the left canvas and [Forumla omitted. See PDF.] are depicted in the right canvas. The true correlations are indicated by the dashed line.

Figure A2. Comparison of the algorithms with respect to the Cramer–von Mises distance between the estimated and the true first ([Forumla omitted. See PDF.]) and true second marginal distributions ([Forumla omitted. See PDF.]). Shown are the boxplots on a logarithmic scale for the proposed EM algorithm (EM), the Standard Copula Estimator (SCOPE), and the Markov chain–Monte Carlo approach (MCMC) for 1000 datasets generated as described in Section 4.2, where [Forumla omitted. See PDF.] are depicted in the left canvas and [Forumla omitted. See PDF.] are depicted in the right canvas.

Figure A3. Comparison of the algorithms with respect to the correlation ([Forumla omitted. See PDF.]). Shown are the mean (upper canvas) and standard deviation (lower canvas) of the Standard Copula Estimator (SCOPE), the proposed EM algorithm (EM), and the Markov chain–Monte Carlo approach (MCMC) for 1000 datasets generated as described in Section 4.2 with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] and for varying sample sizes [Forumla omitted. See PDF.], where the true [Forumla omitted. See PDF.] is [Forumla omitted. See PDF.] (dashed line in the upper canvas).

Figure A4. Comparison of the algorithms with respect to the Cramer–von Mises statistic [Forumla omitted. See PDF.] between the estimated and the true first marginal distribution. Shown are the mean (upper canvas) and standard deviation (lower canvas) of the Standard Copula Estimator (SCOPE), the proposed EM algorithm (EM), and the Markov chain–Monte Carlo approach (MCMC) for 1000 datasets generated as described in Section 4.2 with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] and for varying sample sizes of [Forumla omitted. See PDF.].

Figure A5. Comparison of the algorithms with respect to the Cramer–von Mises statistic [Forumla omitted. See PDF.] between the estimated and the true second marginal distribution. Shown are the mean (upper canvas) and standard deviation (lower canvas) of the Standard Copula Estimator (SCOPE), the proposed EM algorithm (EM), and the Markov chain–Monte Carlo approach (MCMC) for 1000 datasets generated as described in Section 4.2 with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] and for varying sample sizes of N = 100, 200, 500, 1000.

References

1. Thurow, M.; Dumpert, F.; Ramosaj, B.; Pauly, M. Imputing missings in official statistics for general tasks–our vote for distributional accuracy. Stat. J. IAOS; 2021; 37, pp. 1379-1390. [DOI: https://dx.doi.org/10.3233/SJI-210798]

2. Liu, Y.; Dillon, T.; Yu, W.; Rahayu, W.; Mostafa, F. Missing value imputation for industrial IoT sensor data with large gaps. IEEE Internet Things J.; 2020; 7, pp. 6855-6867. [DOI: https://dx.doi.org/10.1109/JIOT.2020.2970467]

3. Silverman, B. Density Estimation for Statistics and Data Analysis; Routledge: London, UK, 2018.

4. Kertel, M.; Harmeling, S.; Pauly, M. Learning causal graphs in manufacturing domains using structural equation models. arXiv; 2022; [DOI: https://dx.doi.org/10.48550/ARXIV.2210.14573] arXiv: 2210.14573

5. Genest, C.; Ghoudi, K.; Rivest, L.P. A semiparametric estimation procedure of dependence parameters in multivariate families of distributions. Biometrika; 1995; 82, pp. 543-552. [DOI: https://dx.doi.org/10.1093/biomet/82.3.543]

6. Liu, H.; Han, F.; Yuan, M.; Lafferty, J.; Wasserman, L. High-dimensional semiparametric gaussian copula graphical models. Ann. Stat.; 2012; 40, pp. 2293-2326. [DOI: https://dx.doi.org/10.1214/12-AOS1037]

7. Titterington, D.; Mill, G. Kernel-based density estimates from incomplete data. J. R. Stat. Soc. Ser. B Methodol.; 1983; 45, pp. 258-266. [DOI: https://dx.doi.org/10.1111/j.2517-6161.1983.tb01249.x]

8. Dempster, A.; Laird, N.; Rubin, D. Maximum likelihood from incomplete data via the EM algorithm. J. R. Stat. Soc. Ser. B Methodol.; 1977; 39, pp. 1-22.

9. Shen, C.; Weissfeld, L. A copula model for repeated measurements with non-ignorable non-monotone missing outcome. Stat. Med.; 2006; 25, pp. 2427-2440. [DOI: https://dx.doi.org/10.1002/sim.2355]

10. Gomes, M.; Radice, R.; Camarena Brenes, J.; Marra, G. Copula selection models for non-Gaussian outcomes that are missing not at random. Stat. Med.; 2019; 38, pp. 480-496. [DOI: https://dx.doi.org/10.1002/sim.7988]

11. Rubin, D.B. Inference and missing data. Biometrika; 1976; 63, pp. 581-592. [DOI: https://dx.doi.org/10.1093/biomet/63.3.581]

12. Cui, R.; Groot, P.; Heskes, T. Learning causal structure from mixed data with missing values using Gaussian copula models. Stat. Comput.; 2019; 29, pp. 311-333. [DOI: https://dx.doi.org/10.1007/s11222-018-9810-x]

13. Wang, H.; Fazayeli, F.; Chatterjee, S.; Banerjee, A. Gaussian copula precision estimation with missing values. Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics; Reykjavik, Iceland, 22–25 April 2014; PMLR: Reykjavik, Iceland, 2014; Volume 33, pp. 978-986.

14. Hamori, S.; Motegi, K.; Zhang, Z. Calibration estimation of semiparametric copula models with data missing at random. J. Multivar. Anal.; 2019; 173, pp. 85-109. [DOI: https://dx.doi.org/10.1016/j.jmva.2019.02.003]

15. Robins, J.M.; Gill, R.D. Non-response models for the analysis of non-monotone ignorable missing data. Stat. Med.; 1997; 16, pp. 39-56. [DOI: https://dx.doi.org/10.1002/(SICI)1097-0258(19970115)16:1<39::AID-SIM535>3.0.CO;2-D]

16. Sun, B.; Tchetgen, E.J.T. On inverse probability weighting for nonmonotone missing at random data. J. Am. Stat. Assoc.; 2018; 113, pp. 369-379. [DOI: https://dx.doi.org/10.1080/01621459.2016.1256814] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30034062]

17. Seaman, S.R.; White, I.R. Review of inverse probability weighting for dealing with missing data. Stat. Methods Med. Res.; 2013; 22, pp. 278-295. [DOI: https://dx.doi.org/10.1177/0962280210395740]

18. Ding, W.; Song, P. EM algorithm in gaussian copula with missing data. Comput. Stat. Data Anal.; 2016; 101, pp. 1-11. [DOI: https://dx.doi.org/10.1016/j.csda.2016.01.008]

19. Efromovich, S. Adaptive nonparametric density estimation with missing observations. J. Stat. Plan. Inference; 2013; 143, pp. 637-650. [DOI: https://dx.doi.org/10.1016/j.jspi.2012.10.008]

20. Dubnicka, S.R. Kernel density estimation with missing data and auxiliary variables. Aust. N. Z. J. Stat.; 2009; 51, pp. 247-270. [DOI: https://dx.doi.org/10.1111/j.1467-842X.2009.00541.x]

21. Hoff, P. Extending the rank likelihood for semiparametric copula estimation. Ann. Appl. Stat.; 2007; 1, pp. 265-283. [DOI: https://dx.doi.org/10.1214/07-AOAS107]

22. Hollenbach, F.; Bojinov, I.; Minhas, S.; Metternich, N.; Ward, M.; Volfovsky, A. Multiple imputation using gaussian copulas. Sociol. Methods Res.; 2021; 50, pp. 1259-1283. [DOI: https://dx.doi.org/10.1177/0049124118799381]

23. Di Lascio, F.; Giannerini, S.; Reale, A. Exploring copulas for the imputation of complex dependent data. Stat. Methods Appl.; 2015; 24, pp. 159-175. [DOI: https://dx.doi.org/10.1007/s10260-014-0287-2]

24. Houari, R.; Bounceur, A.; Kechadi, T.; Tari, A.; Euler, R. A new method for estimation of missing data based on sampling methods for data mining. Adv. Intell. Syst. Comput.; 2013; 225, pp. 89-100. [DOI: https://dx.doi.org/10.1007/978-3-319-00951-3_9]

25. Sklar, A. Fonctions de repartition an dimensions et leurs marges. Publ. Inst. Statist. Univ. Paris; 1959; 8, pp. 229-231.

26. Wei, G.; Tanner, M. A monte carlo implementation of the EM algorithm and the poor man’s data augmentation algorithms. J. Am. Stat. Assoc.; 1990; 85, pp. 699-704. [DOI: https://dx.doi.org/10.1080/01621459.1990.10474930]

27. Meng, X.L.; Rubin, D. Maximum likelihood estimation via the ECM algorithm: A general framework. Biometrika; 1993; 80, pp. 267-278. [DOI: https://dx.doi.org/10.1093/biomet/80.2.267]

28. Guo, J.; Levina, E.; Michailidis, G.; Zhu, J. Graphical models for ordinal data. J. Comput. Graph. Stat.; 2015; 24, pp. 183-204. [DOI: https://dx.doi.org/10.1080/10618600.2014.889023]

29. McLachlan, G.; Lee, S.; Rathnayake, S. Finite mixture models. Annu. Rev. Stat. Its Appl.; 2019; 6, pp. 355-378. [DOI: https://dx.doi.org/10.1146/annurev-statistics-031017-100325]

30. Hwang, J.; Lay, S.; Lippman, A. Nonparametric multivariate density estimation: A comparative study. IEEE Trans. Signal Process.; 1994; 42, pp. 2795-2810. [DOI: https://dx.doi.org/10.1109/78.324744]

31. Scott, D.; Sain, S. Multidimensional density estimation. Handb. Stat.; 2005; 24, pp. 229-261.

32. Zuo, Y.; Cui, Y.; Yu, G.; Li, R.; Ressom, H. Incorporating prior biological knowledge for network-based differential gene expression analysis using differentially weighted graphical LASSO. BMC Bioinform.; 2017; 18, 99. [DOI: https://dx.doi.org/10.1186/s12859-017-1515-1]

33. Li, Y.; Jackson, S.A. Gene network reconstruction by integration of prior biological knowledge. G3 Genes Genomes Genet.; 2015; 5, pp. 1075-1079. [DOI: https://dx.doi.org/10.1534/g3.115.018127] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25823587]

34. Hastie, T.; Tibshirani, R.; Friedman, J. The Elements of Statistical Learning: Data Mining, Inference, and Prediction; Springer Series in Statistics; Springer: Berlin/Heidelberg, Germany, 2009.

35. Joyce, J.M. Kullback-Leibler divergence. International Encyclopedia of Statistical Science; Springer: Berlin/Heidelberg, Germany, 2011; pp. 720-722.

36. Contreras-Reyes, J.E.; Arellano-Valle, R.B. Kullback–Leibler divergence measure for multivariate skew-normal distributions. Entropy; 2012; 14, pp. 1606-1626. [DOI: https://dx.doi.org/10.3390/e14091606]

37. Honaker, J.; King, G.; Blackwell, M. Amelia II: A program for missing data. J. Stat. Softw.; 2011; 45, pp. 1-47. [DOI: https://dx.doi.org/10.18637/jss.v045.i07]

38. Holzinger, A.; Langs, G.; Denk, H.; Zatloukal, K.; Müller, H. Causability and explainability of artificial intelligence in medicine. WIREs Data Min. Knowl. Discov.; 2019; 9, e1312. [DOI: https://dx.doi.org/10.1002/widm.1312] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32089788]

39. Dinu, V.; Zhao, H.; Miller, P.L. Integrating domain knowledge with statistical and data mining methods for high-density genomic SNP disease association analysis. J. Biomed. Inform.; 2007; 40, pp. 750-760. [DOI: https://dx.doi.org/10.1016/j.jbi.2007.06.002]

40. Rubin, D. Multiple imputation after 18+ years. J. Am. Stat. Assoc.; 1996; 91, pp. 473-489. [DOI: https://dx.doi.org/10.1080/01621459.1996.10476908]

41. Van Buuren, S. Flexible Imputation of Missing Data; CRC Press: Boca Raton, FL, USA, 2018.

42. Ramosaj, B.; Pauly, M. Predicting missing values: A comparative study on non-parametric approaches for imputation. Comput. Stat.; 2019; 34, pp. 1741-1764. [DOI: https://dx.doi.org/10.1007/s00180-019-00900-3]

43. Ramosaj, B.; Amro, L.; Pauly, M. A cautionary tale on using imputation methods for inference in matched-pairs design. Bioinformatics; 2020; 36, pp. 3099-3106. [DOI: https://dx.doi.org/10.1093/bioinformatics/btaa082]

44. Zhao, Y.; Udell, M. Missing value imputation for mixed data via gaussian copula. Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining; Virtual Event, CA, USA, 6–10 July 2020; Association for Computing Machinery: New York, NY, USA, 2020; pp. 636-646. [DOI: https://dx.doi.org/10.1145/3394486.3403106]

45. Rubin, D.B. Causal Inference Using Potential Outcomes: Design, Modeling, Decisions. J. Am. Stat. Assoc.; 2005; 100, pp. 322-331. [DOI: https://dx.doi.org/10.1198/016214504000001880]

46. Ding, P.; Li, F. Causal inference: A missing data perspective. Stat. Sci.; 2017; 33, [DOI: https://dx.doi.org/10.1214/18-STS645]

47. Käärik, E.; Käärik, M. Modeling dropouts by conditional distribution, a copula-based approach. J. Stat. Plan. Inference; 2009; 139, pp. 3830-3835. [DOI: https://dx.doi.org/10.1016/j.jspi.2009.05.020]

48. Murphy, K. Machine Learning: A Probabilistic Perspective; The MIT Press: Cambridge, MA, USA, 2012.