Context

When working with \(d\) mixed variables \(X_1, ..., X_d\), it can be complicated to infer a correlation network as well-known statistical methods do not work in case of mixed data. Indeed, most classical network inference methods such as gLasso or Gaussian graphical models rely on the assumption that the data follow a Gaussian distribution. Recently, the Non-paranormal distribution was introduced for network inference for continuous, non-Gaussian data (Liu (2009)). It consists in a transformation of the cumulative distribution functions via a Gaussian copula and provides results for non-Gaussian but continuous variables. We propose an extension of this model to the case of mixed variables.

The Model

Let \(X_1, ..., X_d\) be mixed variables (continuous or discrete). Let \(F_1, ..., F_d\) denote their marginal CDFs, \(\Phi^{-1}\) the inverse of the standard normal CDF and \(\Phi_\Sigma\) the Gaussian CDF of correlation matrix \(\Sigma\). We assume that the multivariate CDF of the vector \((X_1,\dots,X_d)\) is given by: \[F(X_1, ..., X_d)=C_\Sigma(F_1(X_1), ..., F_d(X_d)):=\Phi_\Sigma(\Phi^{-1}(F_1(X_1)), ..., \Phi^{-1}(F_d(X_d))).\]

Estimation

Our package enables the estimation of \(\Sigma\), the correlation matrix of the copula, and \(\Omega=\Sigma^{-1}\), the conditional covariance matrix of the copula..

In order to estimate \(\Sigma\), the rho_estim function uses a semiparametric pairwise maximum likelihood estimator. It returns the estimated correlation matrix of the copula and takes as arguments the data set and the variable types in a vector. In the example below, we have used a subset of the ICGC data set (Zhang, Bajari, and D. (2019)) which contains 5 RNA-seq, 5 protein and 5 mutation variables. We have specified the variable types, where a “C” stands for “continuous” and a “D” for “discrete”.

data(icgc_data)
R <- rho_estim(icgc_data,c(rep("D",5),rep("C",5),rep("D",5)))

A \(6\times6\) subset of the obtained copula correlation matrix is represented below.

The precision matrix, \(\Omega=\Sigma^{-1}\) that encodes the latent conditional correlation network, can also be estimated via gLASSO penalized inversion (Friedman, Hastie, and Tibshirani (2008)). Because it automatically sets some coefficients to zero, no posterior thresholding is needed. Our function omega_estim takes as arguments the data set and the type of variables (or the correlation matrix if it has already been estimated), a grid of penalization parameters \(\lambda\), and the number of observations in the data if a correlation matrix has been entered as the first parameter.

O <- omega_estim(R, lambda=seq(0.4,0.5,0.01), n=250)

It returns a list containing: the correlation matrix, the optimal precision matrix according to the HBIC criterion, the optimal corresponding \(\lambda\), all tested values of \(\lambda\) and all corresponding values of the \(HBIC\). A \(6\times6\) subset of the obtained copula penalized conditional correlation matrix is represented below.

The minimum \(\lambda\) and \(HBIC\) are returned by the package.

#> [1] "The minimal lambda is 0.4"
#> [1] "The minimal HBIC is 11.0917177044302"

Finally, it is possible to graphically represent the evolution of the HBIC depending on the \(\lambda\) by running the following code:

plot(O[[5]],O[[6]],xlab="Lambda",ylab="HBIC",type="l")

Graphical representation

The cor_network_graph function enables to visualize the obtained network. For a correlation network, it takes as arguments the dataset, the correlation matrix, the threshold and a legend.

cor_network_graph(R,TS=0.3,legend=c(rep("RNAseq",5),rep("Proteins",5),rep("Mutations",5)))

To visualize the conditional correlation network that is encoded in \(\Omega\), it is sufficient to specify the threshold parameter to exactly zero as the estimated null coefficients have been shrunk to zero by penalization.

cor_network_graph(O[[2]],TS=0,legend=c(rep("RNAseq",5),rep("Proteins",5),rep("Mutations",5)))

R <- diag_block_matrix(c(3,2),c(0.4,0.8))

R <- matrix_gen(5,0.81)

R <- diag_block_matrix(c(3,5,2),c(0.7,0.3,0.5))
CopulaSim(5,R,c(rep("qnorm(0,1)",5),rep("qexp(0.5)",3),rep("qbinom(4,0.8)",2)),random=TRUE)
#> [[1]]
#>            X1           X2         X3         X4         X5         X6
#> 1 -0.31432935 -0.003397378 -0.3042208 -0.1121517 -0.3293876 1.85493071
#> 2 -0.04569702 -0.157287732 -1.8411259  0.3596431  1.3008398 8.30431146
#> 3 -0.63241425 -0.571101700  1.1884076 -0.9861470 -0.7240118 0.02478069
#> 4 -0.34398010 -0.284438037  0.4719280 -1.1694105 -1.3806585 0.30800447
#> 5 -0.98796694 -0.692130436 -2.1990930 -1.2800052  1.2636802 4.57877686
#>          X7        X8 X9 X10
#> 1 3.1266226 2.8760999  3   3
#> 2 1.1930258 0.7823285  3   3
#> 3 0.6112743 0.1576969  1   3
#> 4 1.5508764 0.7029956  3   4
#> 5 0.3770914 1.3759504  4   3
#> 
#> [[2]]
#>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#>  [1,]  1.0  0.0  0.5  0.0  0.0  0.0  0.0  0.0  0.0   0.0
#>  [2,]  0.0  1.0  0.0  0.0  0.0  0.0  0.7  0.0  0.7   0.0
#>  [3,]  0.5  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
#>  [4,]  0.0  0.0  0.0  1.0  0.3  0.3  0.0  0.3  0.0   0.3
#>  [5,]  0.0  0.0  0.0  0.3  1.0  0.3  0.0  0.3  0.0   0.3
#>  [6,]  0.0  0.0  0.0  0.3  0.3  1.0  0.0  0.3  0.0   0.3
#>  [7,]  0.0  0.7  0.0  0.0  0.0  0.0  1.0  0.0  0.7   0.0
#>  [8,]  0.0  0.0  0.0  0.3  0.3  0.3  0.0  1.0  0.0   0.3
#>  [9,]  0.0  0.7  0.0  0.0  0.0  0.0  0.7  0.0  1.0   0.0
#> [10,]  0.0  0.0  0.0  0.3  0.3  0.3  0.0  0.3  0.0   1.0
#> 
#> [[3]]
#>  [1]  3  2 10  4  8  5  1  6  7  9

latent_data <- gauss_gen(R,10)