---
title: "heterocop: an R package for Gaussian copula semi-parametric inference for heterogeneous data"
author: "Ekaterina Tomilina"
date: "`r Sys.Date()`"
bibliography: biblio.bib
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{heterocop: an R package for Gaussian copula semi-parametric inference for heterogeneous data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(heterocop)
require(dplyr)
require(kableExtra)
require(knitr)
```

This package enables the user to quantify dependencies between mixed (continuous, discrete, binary) variables in the framework of a Gaussian copula model by estimating the correlation matrix of the copula (@tmj) as well as its conditional correlation structure.

## 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). 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 (@ICGC) 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".
```{r,warning=FALSE,message=FALSE}
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.
```{r,echo=FALSE,warning=FALSE,message=FALSE}
knitr::kable(head(R[,1:6]), digits = 2) 
```

The precision matrix, $\Omega=\Sigma^{-1}$ that encodes the latent conditional correlation network, can also be estimated via gLASSO penalized inversion (@friedman). 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.
```{r, warning=FALSE,message=FALSE}
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.
```{r,echo=FALSE,warning=FALSE,message=FALSE}
knitr::kable(head(O[[2]][,1:6]), digits = 2) 
```

The minimum $\lambda$ and $HBIC$ are returned by the package.
```{r,echo=FALSE,warning=FALSE,message=FALSE}
O[[3]]
O[[4]]
```

Finally, it is possible to graphically represent the evolution of the HBIC depending on the $\lambda$ by running the following code:
```{r,warning=FALSE,message=FALSE}
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.
```{r}
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.

```{r,warning=FALSE,message=FALSE}
cor_network_graph(O[[2]],TS=0,legend=c(rep("RNAseq",5),rep("Proteins",5),rep("Mutations",5)))
```

## Simulation

Our package is also able to simulate data distributed according to the Gaussian copula model. Two functions enable us to generate two types of correlation matrices: block-wise and sparse. The diag_block_matrix function enables the user to get a block-wise correlation matrix. It takes as arguments a vector containing the block sizes and a vector containing the coefficients of each block. An example is shown below.

```{r}
R <- diag_block_matrix(c(3,2),c(0.4,0.8))
```
```{r,echo=FALSE,warning=FALSE,message=FALSE}
knitr::kable(R, digits = 2, col.names = NULL)
```

The matrix_gen function enables the user to generate a sparse correlation matrix of initial sparsity parameter $\gamma$, which has to be specified. It is based on the Cholesky decomposition where a lower triangular matrix of sparsity $\gamma$ is generated before being multiplied by its transpose in order to obtain the final matrix. Note that the initial parameter is not equal to the final parameter, which is also returned by the function. In the example below, the first element of the list is the resulting matrix, and the second element of the list is the final sparsity parameter.
```{r}
R <- matrix_gen(5,0.81)
```

```{r,echo=FALSE,warning=FALSE,message=FALSE}
knitr::kable(R, digits = 2, col.names = NULL)
```

The CopulaSim function enables the user to generate a data set which CDF can be expressed as a Gaussian copula of correlation matrix R (to be specified). In the example below, we first generate a block diagonal correlation matrix R and then generate the data set. Then, CopulaSim takes as arguments the number of observations, the correlation matrix of the copula, a vector containing the probability distributions and their parameters, the number of repetitions of each distribution, and enables the user to randomize their order. It returns a list of three elements: the data frame containing the generated data, the correlation matrix, and the permutation realized on the rows and columns of R order after randomization.

```{r}
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)
```

Additionally, the gauss_gen function, which is used in CopulaSim, generates the latent Gaussian variables linked by the correlation matrix R. Its only arguments are the correlation matrix R and the number of observations.

```{r}
latent_data <- gauss_gen(R,10)
```

```{r,echo=FALSE,warning=FALSE,message=FALSE}
knitr::kable(latent_data, digits = 2)%>%
kableExtra::kable_styling(font_size = 8, full_width = FALSE)
```

## References


