Compute (conditional) mutual information

Description

For discrete or categorical variables, the (conditional) mutual information is computed using the empirical frequencies minus a complexity cost (computed as BIC or with the Normalized Maximum Likelihood). When continuous variables are present, each continuous variable is discretized for each mutual information estimate so as to maximize the mutual information minus the complexity cost (see Cabeli 2020).

Usage

computeMutualInfo(
  x,
  y,
  df_conditioning = NULL,
  maxbins = NULL,
  cplx = c("nml", "bic"),
  n_eff = -1,
  sample_weights = NULL,
  is_continuous = NULL,
  plot = FALSE
)

Arguments

Details

For a pair of continuous variables X and Y, the mutual information I(X;Y) will be computed iteratively. In each iteration, the algorithm optimizes the partitioning of X and then of Y, in order to maximize

Ik(X_{d};Y_{d}) = I(X_{d};Y_{d}) - cplx(X_{d};Y_{d})

where cplx(X_{d}; Y_{d}) is the complexity cost of the corresponding partitioning (see Cabeli 2020). Upon convergence, the information terms I(X_{d};Y_{d}) and Ik(X_{d};Y_{d}), as well as the partitioning of X_{d} and Y_{d} in terms of cutpoints, are returned.

For conditional mutual information with a conditioning set U, the computation is done based on

Ik(X;Y|U) = 0.5*(Ik(X_{d};Y_{d},U_{d}) - Ik(X_{d};U_{d}) + Ik(Y_{d};X_{d},U_{d}) - Ik(Y_{d};U_{d})),

where each of the four summands is estimated separately.

Value

A list that contains :

• cutpoints1: Only when X is continuous, a vector containing the cutpoints for the partitioning of X.

• cutpoints2: Only when Y is continuous, a vector containing the cutpoints for the partitioning of Y.

• n_iterations: Only when at least one of the input variables is continuous, the number of iterations it takes to reach the convergence of the estimated information.

• iteration1, iteration2, ... Only when at least one of the input variables is continuous, the list of vectors of cutpoints of each iteration.

• info: The estimation of (conditional) mutual information without the complexity cost.

• infok: The estimation of (conditional) mutual information with the complexity cost (Ik = I - cplx).

• plot: Only when 'plot == TRUE', the plot object.

References

• Cabeli et al., PLoS Comput. Biol. 2020, Learning clinical networks from medical records based on information estimates in mixed-type data

• Affeldt et al., UAI 2015, Robust Reconstruction of Causal Graphical Models based on Conditional 2-point and 3-point Information

Examples

library(miic)
N <- 1000
# Dependence, conditional independence : X <- Z -> Y
Z <- runif(N)
X <- Z * 2 + rnorm(N, sd = 0.2)
Y <- Z * 2 + rnorm(N, sd = 0.2)
res <- computeMutualInfo(X, Y, plot = FALSE)
message("I(X;Y) = ", res$info)
res <- computeMutualInfo(X, Y, df_conditioning = matrix(Z, ncol = 1), plot = FALSE)
message("I(X;Y|Z) = ", res$info)


# Conditional independence with categorical conditioning variable : X <- Z -> Y
Z <- sample(1:3, N, replace = TRUE)
X <- -as.numeric(Z == 1) + as.numeric(Z == 2) + 0.2 * rnorm(N)
Y <- as.numeric(Z == 1) + as.numeric(Z == 2) + 0.2 * rnorm(N)
res <- miic::computeMutualInfo(X, Y, cplx = "nml")
message("I(X;Y) = ", res$info)
res <- miic::computeMutualInfo(X, Y, matrix(Z, ncol = 1), is_continuous = c(TRUE, TRUE, FALSE))
message("I(X;Y|Z) = ", res$info)


# Independence, conditional dependence : X -> Z <- Y
X <- runif(N)
Y <- runif(N)
Z <- X + Y + rnorm(N, sd = 0.1)
res <- computeMutualInfo(X, Y, plot = TRUE)
message("I(X;Y) = ", res$info)
res <- computeMutualInfo(X, Y, df_conditioning = matrix(Z, ncol = 1), plot = TRUE)
message("I(X;Y|Z) = ", res$info)

Compute (conditional) three-point information

Description

Three point information is defined and computed as the difference of mutual information and conditional mutual information, e.g.

I(X;Y;Z|U) = I(X;Y|U) - Ik(X;Y|U,Z)

For discrete or categorical variables, the three-point information is computed with the empirical frequencies minus a complexity cost (computed as BIC or with the Normalized Maximum Likelihood).

Usage

computeThreePointInfo(
  x,
  y,
  z,
  df_conditioning = NULL,
  maxbins = NULL,
  cplx = c("nml", "bic"),
  n_eff = -1,
  sample_weights = NULL,
  is_continuous = NULL
)

Arguments

Details

For variables X, Y, Z and a set of conditioning variables U, the conditional three point information is defined as

Ik(X;Y;Z|U) = Ik(X;Y|U) - Ik(X;Y|U,Z)

where Ik is the shifted or regularized conditional mutual information. See computeMutualInfo for the definition of Ik.

Value

A list that contains :

• i3: The estimation of (conditional) three-point information without the complexity cost.

• i3k: The estimation of (conditional) three-point information with the complexity cost (i3k = i3 - cplx).

• i2: For reference, the estimation of (conditional) mutual information I(X;Y|U) used in the estimation of i3.

• i2k: For reference, the estimation of regularized (conditional) mutual information Ik(X;Y|U) used in the estimation of i3k.

References

• Cabeli et al., PLoS Comput. Biol. 2020, Learning clinical networks from medical records based on information estimates in mixed-type data

• Affeldt et al., UAI 2015, Robust Reconstruction of Causal Graphical Models based on Conditional 2-point and 3-point Information

Examples

library(miic)
N <- 1000
# Dependence, conditional independence : X <- Z -> Y
Z <- runif(N)
X <- Z * 2 + rnorm(N, sd = 0.2)
Y <- Z * 2 + rnorm(N, sd = 0.2)
res <- computeThreePointInfo(X, Y, Z)
message("I(X;Y;Z) = ", res$i3)
message("Ik(X;Y;Z) = ", res$i3k)


# Independence, conditional dependence : X -> Z <- Y
X <- runif(N)
Y <- runif(N)
Z <- X + Y + rnorm(N, sd = 0.1)
res <- computeThreePointInfo(X, Y, Z)
message("I(X;Y;Z) = ", res$i3)
message("Ik(X;Y;Z) = ", res$i3k)

Genomic and ploidy alterations in breast tumors

Description

The dataset contains 807 samples without predisposing Brca1/2 germline mutations and includes 204 somatic mutations (from whole exome sequencing) and expression level information for 91 genes.

Usage

data(cosmicCancer)

Format

A data.frame object.

References

Forbes SA, Beare D, Gunasekaran P, Leung K, Bindal N, et al. (2015) Nucleic Acids Res 43:D805–D811. (PubMed link)

Genomic and ploidy alterations in breast tumors

Description

The dataset contains 807 samples without predisposing Brca1/2 germline mutations and includes 204 somatic mutations (from whole exome sequencing) and expression level information for 91 genes, category order file.

Usage

data(cosmicCancer_stateOrder)

Format

A data.frame object.

References

Forbes SA, Beare D, Gunasekaran P, Leung K, Bindal N, et al. (2015) Nucleic Acids Res 43:D805–D811. (PubMed link)

Covid cases

Description

Demo dataset of chronological series to be used in temporal mode of miic. Evolution of Covid cases on a subset of EU countries from 12/31/2019 to 06/18/2020. Source of the data : European Centre for Disease Prevention and Control.

Usage

data(covidCases)

Format

A data.frame object.

References

ECDC (ECDC link)

Discretize a real valued distribution

Description

This function performs minimum description length (MDL)-optimal histogram density estimation as described in Kontkanen and Myllymäki (2007) and returns the cutpoints found to give the best model according to the MDL principle.

Usage

discretizeMDL(x = NULL, max_bins = 20)

Arguments

Value

A list containing the cutpoints of the best discretization.

References

• Kontkanen P, Myllymäki P. MDL histogram density estimation. Artificial Intelligence and Statistics 2007 Mar 11 (pp. 219-226).

Examples

library(miic)
# Bimodal normal distribution
N <- 300
modes <- sample(1:2, size = N, replace = TRUE)
x <- as.numeric(modes == 1) * rnorm(N, mean = 0, sd = 1) +
     as.numeric(modes == 2) * rnorm(N, mean = 5, sd = 2)
MDL_disc <- discretizeMDL(x)
hist(x, breaks = MDL_disc$cutpoints)

N <- 2000
modes <- sample(1:2, size = N, replace = TRUE)
x <- as.numeric(modes == 1) * rnorm(N, mean = 0, sd = 1) +
     as.numeric(modes == 2) * rnorm(N, mean = 5, sd = 2)
MDL_disc <- discretizeMDL(x)
hist(x, breaks = MDL_disc$cutpoints)

Iterative dynamic programming for (conditional) mutual information through optimized discretization.

Description

This function chooses cutpoints in the input distributions by maximizing the mutual information minus a complexity cost (computed as BIC or with the Normalized Maximum Likelihood). The (conditional) mutual information computed on the optimized discretized distributions effectively estimates the mutual information of the original continuous variables.

Usage

discretizeMutual(
  x,
  y,
  matrix_u = NULL,
  maxbins = NULL,
  cplx = "nml",
  n_eff = NULL,
  sample_weights = NULL,
  is_continuous = NULL,
  plot = TRUE
)

Arguments

Details

For a pair of continuous variables X and Y, the algorithm will iteratively choose cutpoints on X then on Y, maximizing I(X_{d};Y_{d}) - cplx(X_{d};Y_{d}) where cplx(X_{d};Y_{d}) is the complexity cost of the considered discretizations of X and Y (see Cabeli 2020). Upon convergence, the discretization scheme of X_{d} and Y_{d} is returned as well as I(X_{d};Y_{d}) and I(X_{d};Y_{d})-cplx(X_{d};Y_{d}).

With a set of conditioning variables U, the discretization scheme maximizes each term of the sum I(X;Y|U) \sim 0.5*(I(X_{d};Y_{d}, U_{d}) - I(X_{d};U_{d}) + I(Y_{d};X_{d}, U_{d}) - I(Y_{d};U_{d})).

Discrete variables can be passed as factors and will be used "as is" to maximize each term.

Value

A list that contains :

• two vectors containing the cutpoints for each variable : cutpoints1 corresponds to x, cutpoints2 corresponds to y.

• n_iterations is the number of iterations performed before convergence of the (C)MI estimation.

• iteration1, iteration2, ..., lists containing the cutpoint vectors for each iteration.

• info and infok, the estimated (C)MI value and (C)MI minus the complexity cost.

• if plot == TRUE, a plot object (requires ggplot2 and gridExtra).

References

• Cabeli et al., PLoS Comput. Biol. 2020, Learning clinical networks from medical records based on information estimates in mixed-type data

Examples

library(miic)
N <- 1000
# Dependence, conditional independence : X <- Z -> Y
Z <- runif(N)
X <- Z * 2 + rnorm(N, sd = 0.2)
Y <- Z * 2 + rnorm(N, sd = 0.2)
res <- discretizeMutual(X, Y, plot = FALSE)
message("I(X;Y) = ", res$info)
res <- discretizeMutual(X, Y, matrix_u = matrix(Z, ncol = 1), plot = FALSE)
message("I(X;Y|Z) = ", res$info)


# Conditional independence with categorical conditioning variable : X <- Z -> Y
Z <- sample(1:3, N, replace = TRUE)
X <- -as.numeric(Z == 1) + as.numeric(Z == 2) + 0.2 * rnorm(N)
Y <- as.numeric(Z == 1) + as.numeric(Z == 2) + 0.2 * rnorm(N)
res <- miic::discretizeMutual(X, Y, cplx = "nml")
message("I(X;Y) = ", res$info)
res <- miic::discretizeMutual(X, Y, matrix(Z, ncol = 1), is_continuous = c(TRUE, TRUE, FALSE))
message("I(X;Y|Z) = ", res$info)


# Independence, conditional dependence : X -> Z <- Y
X <- runif(N)
Y <- runif(N)
Z <- X + Y + rnorm(N, sd = 0.1)
res <- discretizeMutual(X, Y, plot = TRUE)
message("I(X;Y) = ", res$info)
res <- discretizeMutual(X, Y, matrix_u = matrix(Z, ncol = 1), plot = TRUE)
message("I(X;Y|Z) = ", res$info)

Estimation of the temporal causal discovery parameters

Description

This function estimates the number of layers and number of time steps between each layer that are needed to cover the dynamic of a temporal dataset when reconstructing a temporal causal graph. Using autocorrelation decay, the function computes the average relaxation time of the variables and, based on a maximum number of nodes, deduces the number of layers and number of time steps between each layer to be used.

Usage

estimateTemporalDynamic(
  input_data,
  state_order = NULL,
  mov_avg = NULL,
  max_nodes = 50,
  verbose_level = 1
)

Arguments

Value

A named list with two items:

• n_layers: the number of layers

• delta_t: the number of time steps between the layers

Export miic result for plotting (with igraph)

Description

This function creates an object built from the result returned by miic that is ready to be fed to the plotting method.

Usage

export(
  miic_obj,
  method = "igraph",
  pcor_palette = NULL,
  display = "compact",
  show_self_loops = TRUE
)

Arguments

Details

The behavior depends on the method used for the export.

For igraph, edge attributes are passed to the igraph graph and can be accessed with e.g. E(g)$partial_correlation. See miic for more details on edge parameters. By default, edges are colored according to the partial correlation between two nodes conditioned on the conditioning set (negative is blue, null is gray and positive is red) and their width is based on the conditional mutual information minus the complexity cost.

Value

A graph object adapted to the method.

Examples

library(miic)
data(hematoData)

# execute MIIC (reconstruct graph)
miic_obj <- miic(
  input_data = hematoData, latent = "yes",
  n_shuffles = 10, conf_threshold = 0.001
)

# Using igraph
if(require(igraph)) {
g = export(miic_obj, "igraph")
plot(g) # Default visualisation, calls igraph::plot.igraph()

# Specifying layout (see ?igraph::layout_)
l <-layout_with_kk(g)
plot(g, layout=l)

# Override some graphical parameters
plot(g, edge.curved = .2)
plot(g, vertex.shape="none", edge.color="gray85", vertex.label.color="gray10")
}

# In temporal mode, execute MIIC
data(covidCases)
tmiic_obj <- miic(input_data = covidCases, mode = "TS", n_layers = 3, delta_t = 1, mov_avg = 14)

# Plot by default the compact display of the temporal network using igraph
if(require(igraph)) {
g = export (tmiic_obj)
plot(g)

# Plot the raw temporal network using igraph
g = export(tmiic_obj, display="raw")
plot(g)

# Plot the complete temporal network using igraph (completed by stationarity)
g = export(tmiic_obj, display="lagged")
plot(g)

# Specifying layout (see ?igraph::layout_)
l <- layout_on_grid(g, width = 5, height = 3, dim = 2)
plot(g, layout=l)

# For compact temporal display, please be aware that the rendering of
# igraph::plot.igraph() is not optimal when the graph contains
# multiple edges between the same nodes.
# So, the recommend way to plot a compact graph is to use tmiic plotting:
plot(tmiic_obj)
}

Early blood development: single cell binary gene expression data

Description

Binarized expression data of 33 transcription factors involved in early differentiation of primitive erythroid and endothelial cells (3934 cells).

Usage

data(hematoData)

Format

A data.frame object.

References

Moignard et al. (2015) Nat Biotechnol 33(3):269-76 (PubMed link)

MIIC, causal network learning algorithm including latent variables

Description

MIIC (Multivariate Information-based Inductive Causation) combines constraint-based and information-theoretic approaches to disentangle direct from indirect effects amongst correlated variables, including cause-effect relationships and the effect of unobserved latent causes.

Usage

miic(
  input_data,
  state_order = NULL,
  true_edges = NULL,
  black_box = NULL,
  n_threads = 1,
  cplx = "nml",
  orientation = TRUE,
  ort_proba_ratio = 1,
  ort_consensus_ratio = NULL,
  propagation = FALSE,
  latent = "orientation",
  n_eff = -1,
  n_shuffles = 0,
  conf_threshold = 0,
  sample_weights = NULL,
  test_mar = TRUE,
  consistent = "no",
  max_iteration = 100,
  consensus_threshold = 0.8,
  negative_info = FALSE,
  mode = "S",
  n_layers = NULL,
  delta_t = NULL,
  mov_avg = NULL,
  keep_max_data = FALSE,
  max_nodes = 50,
  verbose = FALSE
)

Arguments

Details

Starting from a complete graph, the method iteratively removes dispensable edges, by uncovering significant information contributions from indirect paths, and assesses edge-specific confidences from randomization of available data. The remaining edges are then oriented based on the signature of causality in observational data. Miic distinguishes genuine causal edges (with both reliable arrow heads and tails) from putative causal edges (with one reliable arrow head only) and latent causal edges (with both reliable arrow heads). (see Ribeiro-Dantas 2024)

In temporal mode, miic reorganizes the dataset using the n_layers and delta_t parameters to transform the time steps into lagged samples. As starting point, a lagged graph is created with only edges having at least one node laying on the last time step. Then, miic standard algorithm is applied to remove dispensable edges. The remaining edges are then duplicated to ensure time invariance (stationary dynamic) and oriented using the temporality and the signature of causality in observational data. The use of temporal mode is presented in Simon 2024.

The method relies on information theoretic principles which replace (conditional) independence tests as described in Affeldt 2015, Cabeli 2020, Cabeli 2021 and Ribeiro-Dantas 2024. It deals with both categorical and continuous variables by performing optimal context-dependent discretization. As such, the input data frame may contain both numerical columns which will be treated as continuous, or character / factor columns which will be treated as categorical. For further details on the optimal discretization method and the conditional independence test, see the function discretizeMutual. The user may also choose to run miic with scheme presented in Li 2019 and Ribeiro-Dantas 2024 to improve the end result's interpretability by ensuring consistent separating sets.

Value

A miic-like object that contains:

• summary: a data frame with information about the relationship between relevant pair of variables.

  As returning the information on all possible pairs of variables could lead to an huge data frame, by convention, the summary does not include pair of variables not sharing information at all (I'(x,y) <= 0). However, as exception to this convention, when a ground truth is supplied (using the true_edges parameter), the edges that are not retained by MIIC because the variables does not share information at all but are present in the true edges will be included in the summary to report correctly all the false negative edges.

  So, the summary contains these categories of edges:

  • edges retained

  • edges not retained after conditioning on some contributor(s)

  • edges not retained without conditioning but present in true edges

  while these edges are not considered as relevant and are not included:

  • edges not retained without conditioning and not in true edges

  Information available in the summary are:

  • x: X node name

  • y: Y node name

  • type: contains 'N' if the edge has been removed or 'P' for retained edges. If the true graph is supplied in the true_edges parameter, 'P' becomes 'TP' (True Positive) or 'FP' (False Positive), while 'N' becomes 'TN' (True Negative) or 'FN' (False Negative). Note that, as the summary does not contain all the removed edges, edges not present have to be considered as 'N' and, if the true graph is supplied, as 'TN'.

  • ai: the contributing nodes found by the method which contribute to the mutual information between x and y, and possibly separate them.

  • raw_contributions: describes the share of total mutual information between x and y explained by each contributor, measured by I'(x;y;ai|{aj}) / I'(x;y), where {aj} is the separating set before adding ai.

  • contributions: describes the share of remaining mutual information between x and y explained by each successive contributors, measured by I'(x;y;ai|{aj}) / I'(x;y|{aj}), where {aj} is the separating set before adding ai.

  • info: the mutual information I(x;y) times n_xy, the number of samples without missing or NA values for both x and y.

  • n_xy: gives the number of samples on which the information without conditioning has been computed. If the input dataset has no missing value, the number of samples is the same for all pairs and corresponds to the total number of samples.

  • info_cond: the conditional mutual information I(x;y|ai) times the number of samples without NA n_xy_ai used in the computation. info_cond is equal to info when ai is an empty set.

  • cplx: the complexity term for the pair (x, y) taking into account the contributing nodes ai.

  • n_xy_ai: the number of samples without NA in x, y and all nodes in ai on which the information and the complexity terms are computed. If the input dataset has no missing value, the number of samples is the same for all pairs and corresponds to the total number of samples.

  • info_shifted: value equal to info_cond - cplx. Used to decide whether the edge is retained (when positive), or removed (when zero or possibly negative when the parameter negative_info is set to TRUE).

  • ort_inferred: the orientation of the edge (x, y). 0: edge removed, 1: un-directed, 2: directed from X to Y, -2: directed from Y to X, 6: bi-directed.
    When the consistent option is turned on and there is more than one graph in the consistent cycle, this is the inferred orientation of the edge in the last graph in the cycle.

  • ort_ground_truth: the orientation of the edge (x, y) in the ground truth graph when true edges are provided.

  • is_inference_correct: indicates if the inferred orientation agrees with the provided ground truth. TRUE: agrees, FALSE: disagrees and set to NA when no ground truth is supplied.

  • is_causal: boolean value indicating the causal nature of the arrow tips of an edge, based on the probabilities given in the columns p_y2x and p_x2y. TRUE: when the edges is directed and both the head and the tail are set with high confidence (adjustable with the ort_consensus_ratio parameter), FALSE otherwise or NA if the edge is not retained. More formally, an oriented edge is marked as genuine causal when (1 - p_{head}) / p_{head} < ort_consensus_ratio and p_{tail} / (1 - p_{tail}) < ort_consensus_ratio.
    A directed edge not marked as genuine causal indicates that only the head is set with high confidence, while the tail remains uncertain. This corresponds to a putative causal edge, which could either be a genuine causal edge or a bi-directed edge from a latent confounder.
    Note that the genuine causality is deducible only when latent variables are allowed and propagation is not allowed.

  • ort_consensus: Not computed (NAs) when consistency is not activated or, when consistency is on, if there is only one graph returned (no cycle). When computed, indicates the consensus orientation of the edge determined from the consensus skeleton and the ort_consensus_ratio threshold on averaged orientation probabilities over the cycle of graphs. Possible values are 0: not connected, 1: un-directed, -2 or 2: directed and 6: bi-directed (latent variable).

  • is_causal_consensus: Not computed (NAs) when consistency is not activated or, when consistency is on, if there is only one graph returned (no cycle). When computed, work in the same way as is_causal but on the consensus graph.

  • edge_stats: Not computed (NAs) when consistency is not activated or, when consistency is on, if there is only one graph returned (no cycle). When computed, contains the frequencies of all ort_inferred values present in the cycle of graphs for the edge (x, y), in the format [percentage(orientation)], separated by ";". e.g. In a cycle of 4 graphs, if an edge is three times marked as 2 (directed) and one time marked as 1 (un-directed), edge_stats will contain "75%(2);25%(1)".

  • sign: the sign of the partial correlation between variables x and y, conditioned on the contributing nodes ai.

  • partial_correlation: value of the partial correlation for the edge (x, y) conditioned on the contributing nodes ai.

  • p_y2x: probability of the arrowhead from y to x, of the inferred orientation, derived from the three-point mutual information (see Verny 2017 and Ribeiro-Dantas 2024). NA if the edge is removed.

  • p_x2y: probability of the arrowhead from x to y, of the inferred orientation, derived from the three-point mutual information (see Verny 2017 and Ribeiro-Dantas 2024). NA if the edge is removed.

  • confidence: computed only when the confidence cut is activated, NA otherwise. When computed, it corresponds to a measure of the strength of the retained edges: it is the ratio between the probability to reject the edge exp(-info_shifted(x;y|ai)) in the original dataset and the mean probability to do the same in n_shuffles number of randomized datasets. Edges with confidence > conf_threshold will be filtered out from the graph. (see parameters n_shuffles and conf_threshold)

• edges: a data frame with the raw edges output coming from the C++ core function. This data frame is used internally by MIIC to produce the summary and contains all pairs of variables (x, y).

• triples: this data frame lists the orientation probabilities of the two edges of all unshielded triples of the reconstructed network with the structure: node1 – mid-node – node2:

  • node1: node at the end of the unshielded triplet

  • p1: probability of the arrowhead node1 <- mid-node

  • p2: probability of the arrowhead node1 -> mid-node

  • mid-node: node at the center of the unshielded triplet

  • p3: probability of the arrowhead mid-node <- node2

  • p4: probability of the arrowhead mid-node -> node2

  • node2: node at the end of the unshielded triplet

  • ni3: 3 point (conditional) mutual information * N

  • conflict: indicates if there is a conflict between the computed probabilities and the ni3 value

• adj_matrix: the adjacency matrix is a square matrix used to represent the inferred graph. The entries of the matrix indicate whether pairs of vertices are adjacent or not in the graph. The matrix can be read as a (row, column) set of couples where the row represents the source node and the column the target node. Since miic can reconstruct mixed networks (including directed, un-directed and bi-directed edges), we will have a different digit for each case:

  • 1: (x, y) edge is un-directed

  • 2: (x, y) edge is directed as x -> y

  • -2: (x, y) edge is directed as x <- y

  • 6: (x, y) edge is bi-directed

• proba_adj_matrix: the probability adjacency matrix is a square matrix used to represent the orientation probabilities associated to the edges tips. The value at ("row", "column") is the probability, for the edge between "row" and "column" nodes, of the edge tip on the "row" side. A probability less than 0.5 is an indication of a possible tail (cause) and a probability greater than 0.5 a possible head (effect).

• adj_matrices: present only when consistency is activated. The list of the adjacency matrices, one for each graph which is part of the resulting cycle of graphs. Each item is a square matrix with the same layout as adj_matrix.

• proba_adj_matrices: present only when consistency is activated. The list of the probability adjacency matrices, one for each graph which is part of the resulting cycle of graphs. Each item is a square matrix with the same layout as proba_adj_matrix.

• proba_adj_average: present only when consistency is activated. The average probability adjacency matrix is a square matrix used to represent the orientation probabilities associated to the edges tips of the consensus graph. Its layout is the same as proba_adj_matrix and it contains the averaged probability of edges tips over the resulting cycle of graphs.

• is_consistent: present only when consistency is activated. TRUE if the returned graph is consistent, FALSE otherwise.

• time: execution time of the different steps and total run-time of the causal graph reconstruction by MIIC.

• interrupted: TRUE if causal graph reconstruction has been interrupted, FALSE otherwise.

• scores: present only when true edges have been supplied. Contains the scores of the returned graph in regard of the ground truth:

  • tp: number of edges marked as True Positive

  • fp: number of edges marked as False Positive

  • fn: number of edges marked as False Negative

  • precision: Precision

  • recall: Recall

  • fscore: F1-Score

• params: the list of parameters used for the network reconstruction. The parameters not supplied are initialized to their default values. Otherwise, the parameters are checked and corrected if necessary.

• state_order: the state order used for the network reconstruction. If no state order is supplied, it is generated by using default values. Otherwise, it is the state order checked and corrected if necessary.

• black_box: present only if a black box has been supplied, the black box, checked and corrected if necessary, used for the network reconstruction.

• true_edges: present only if the true edges have been supplied, the true edges, checked and corrected if necessary, used for the network evaluation.

• tmiic: present only in temporal mode. Named list containing the full list of edges completed by stationarity, the lagged state order and, if a black box or true edges have been supplied, the lagged versions of these inputs.

References

• Simon et al., eLife 2024, CausalXtract: a flexible pipeline to extract causal effects from live-cell time-lapse imaging data

• Ribeiro-Dantas et al., iScience 2024, Learning interpretable causal networks from very large datasets, application to 400,000 medical records of breast cancer patients

• Cabeli et al., NeurIPS 2021, Reliable causal discovery based on mutual information supremum principle for finite dataset

• Cabeli et al., PLoS Comput. Biol. 2020, Learning clinical networks from medical records based on information estimates in mixed-type data

• Li et al., NeurIPS 2019, Constraint-based causal structure learning with consistent separating sets

• Verny et al., PLoS Comput. Biol. 2017, Learning causal networks with latent variables from multivariate information in genomic data

• Affeldt et al., UAI 2015, Robust Reconstruction of Causal Graphical Models based on Conditional 2-point and 3-point Information

See Also

discretizeMutual for optimal discretization and (conditional) independence test.

Examples

library(miic)

# EXAMPLE HEMATOPOIESIS
data(hematoData)

# execute MIIC (reconstruct graph)
miic_obj <- miic(
  input_data = hematoData[1:1000,], latent = "yes",
  n_shuffles = 10, conf_threshold = 0.001
)

# plot graph
if(require(igraph)) {
 plot(miic_obj, method="igraph")
}


# write graph to graphml format. Note that to correctly visualize
# the network we created the miic style for Cytoscape (http://www.cytoscape.org/).

writeCytoscapeNetwork(miic_obj, file = file.path(tempdir(), "temp"))

# EXAMPLE CANCER
data(cosmicCancer)
data(cosmicCancer_stateOrder)
# execute MIIC (reconstruct graph)
miic_obj <- miic(
  input_data = cosmicCancer, state_order = cosmicCancer_stateOrder, latent = "yes",
  n_shuffles = 100, conf_threshold = 0.001
)

# plot graph
if(require(igraph)) {
 plot(miic_obj)
}

# write graph to graphml format. Note that to correctly visualize
# the network we created the miic style for Cytoscape (http://www.cytoscape.org/).
writeCytoscapeNetwork(miic_obj, file = file.path(tempdir(), "temp"))

# EXAMPLE COVID CASES (time series demo)
data(covidCases)
# execute MIIC (reconstruct graph in temporal mode)
tmiic_obj <- miic(input_data = covidCases, mode = "TS", n_layers = 3, delta_t = 1, mov_avg = 14)

# to plot the default graph (compact)
if(require(igraph)) {
 plot(tmiic_obj)
}

# to plot the raw temporal network
if(require(igraph)) {
  plot(tmiic_obj, display="raw")
}

# to plot the full temporal network
if(require(igraph)) {
  plot(tmiic_obj, display="lagged")
}

Basic plot function of a miic network inference result

Description

This function calls export to build a plottable object from the result returned by miic and plot it.

Usage

## S3 method for class 'miic'
plot(x, method = "igraph", pcor_palette = NULL, ...)

Arguments

Details

See the documentation of export for further details.

See Also

export for graphical exports, igraph.plotting

Basic plot function of a temporal miic (tmiic) network inference result

Description

This function calls export to build a plottable object from the result returned by miic in temporal mode and plot it.

Usage

## S3 method for class 'tmiic'
plot(
  x,
  method = "igraph",
  pcor_palette = NULL,
  display = "compact",
  show_self_loops = TRUE,
  positioning_for_grid = "greedy",
  orientation_for_grid = "L",
  ...
)

Arguments

Details

See the documentation of export for further details.

See Also

export for graphical exports, igraph.plotting

Examples

library(miic)

#' # EXAMPLE COVID CASES (time series demo)
data(covidCases)
# execute MIIC (reconstruct graph in temporal mode)
tmiic_obj <- miic(input_data = covidCases, mode = "TS", n_layers = 3, delta_t = 1, mov_avg = 14)

# to plot the default compact graph
if(require(igraph)) {
  plot(tmiic_obj)
}

# to plot the raw temporal network
if(require(igraph)) {
  plot(tmiic_obj, display="raw")
}

# to plot the full temporal network
if(require(igraph)) {
  plot(tmiic_obj, display="lagged")
}

GraphML converting function for miic graph

Description

Convert miic graph to GraphML format.

Usage

writeCytoscapeNetwork(miic_obj, file, layout = NULL)

Arguments

Value

None

Style writing function for the miic network

Description

This function writes the miic style for a correct visualization using the cytoscape tool (http://www.cytoscape.org/).

Usage

writeCytoscapeStyle(file)

Arguments

Details

The style is written in the xml file format.

Value

None