From 346caeae3cda7122477a477f9a44b6e803b8b956 Mon Sep 17 00:00:00 2001 From: santikka Date: Fri, 18 Dec 2020 16:37:18 +0200 Subject: [PATCH] removed ggm dependency --- .Rbuildignore | 2 + DESCRIPTION | 23 ++++--- NAMESPACE | 32 +++++++++- NEWS | 4 ++ R/ancestors.R | 12 +++- R/c.components.R | 12 ++-- R/causal.effect.R | 15 ++--- R/causal.parents.R | 2 +- R/children.R | 10 +++- R/compare.graphs.R | 10 ++-- R/connected.R | 4 +- R/dSep.R | 110 ++++++++++++++++++++++++++++++++++ R/descendants.R | 4 +- R/descendent.sets.R | 3 +- R/exclusion.restrictions.R | 6 +- R/gather.R | 14 ----- R/generalize.R | 14 ++--- R/generalize_ex.R | 14 ++--- R/id.R | 8 +-- R/idc.R | 9 ++- R/identify.R | 8 +-- R/independence.restrictions.R | 4 +- R/insert.R | 18 ++---- R/irrelevant.R | 4 +- R/join.R | 20 ++----- R/latent.projection.R | 24 ++++---- R/meta.transport.R | 6 +- R/observed.graph.R | 4 +- R/parents.R | 10 +++- R/parse.expression.R | 12 ++-- R/parse.graphml.R | 13 ++-- R/parse.graphml.internal.R | 26 ++++---- R/parse.graphml.standard.R | 49 ++++++++------- R/pid.R | 18 +++--- R/powerset.R | 7 +++ R/q.constraints.R | 8 +-- R/rc.R | 12 ++-- R/recover.R | 20 +++---- R/sc.components.R | 14 ++--- R/simplify.R | 6 +- R/soid.R | 10 ++-- R/surrogate.outcome.R | 18 +++--- R/tr.target.R | 2 +- R/transport.R | 10 ++-- R/trmz.R | 23 ++++--- R/trmz_ex.R | 32 +++++----- R/trso.R | 19 +++--- R/unobserved.graph.R | 14 ++--- R/verma.constraints.R | 28 +++------ R/wrap.dSep.R | 8 +-- R/zzaux.effect.R | 4 +- 51 files changed, 440 insertions(+), 319 deletions(-) create mode 100644 .Rbuildignore create mode 100644 R/dSep.R delete mode 100644 R/gather.R create mode 100644 R/powerset.R diff --git a/.Rbuildignore b/.Rbuildignore new file mode 100644 index 0000000..af21c6d --- /dev/null +++ b/.Rbuildignore @@ -0,0 +1,2 @@ +^.*\.git$ +^README.md$ \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index 6c45606..8f8acdc 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,11 +1,18 @@ Package: causaleffect -Version: 1.3.11 -Date: 2020-05-14 +Version: 1.3.12 +Date: 2020-12-18 Title: Deriving Expressions of Joint Interventional Distributions and Transport Formulas in Causal Models -Author: Santtu Tikka -Maintainer: Santtu Tikka -Imports: ggm, igraph, XML -Suggests: R.rsp +Authors@R: person(given = "Santtu", + family = "Tikka", + role = c("aut", "cre"), + email = "santtuth@gmail.com", + comment = c(ORCID = "0000-0003-4039-4342")) +BugReports: https://github.com/santikka/causaleffect/issues +Description: Functions for identification and transportation of causal effects. Provides a conditional causal effect identification algorithm (IDC) by Shpitser, I. and Pearl, J. (2006) , an algorithm for transportability from multiple domains with limited experiments by Bareinboim, E. and Pearl, J. (2014) and a selection bias recovery algorithm by Bareinboim, E. and Tian, J. (2015) . All of the previously mentioned algorithms are based on a causal effect identification algorithm by Tian , J. (2002) . +License: GPL (>= 2) +Imports: igraph +Suggests: R.rsp, XML VignetteBuilder: R.rsp -Description: Functions for identification and transportation of causal effects. Provides a conditional causal effect identification algorithm (IDC) by Shpitser, I. and Pearl, J. (2006) , an algorithm for transportability from multiple domains with limited experiments by Bareinboim, E. and Pearl, J. (2014) and a selection bias recovery algorithm by Bareinboim, E. and Tian, J. (2015) . All of the previously mentioned algorithms are based on a causal effect identification algorithm by Tian , J. (2002) . -License: GPL-2 \ No newline at end of file +NeedsCompilation: no +Author: Santtu Tikka [aut, cre] () +Maintainer: Santtu Tikka \ No newline at end of file diff --git a/NAMESPACE b/NAMESPACE index a17d79e..543ee85 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,31 @@ -export(aux.effect, causal.effect, generalize, get.expression, meta.transport, parse.graphml, recover, surrogate.outcome, transport, verma.constraints) -import(igraph, XML) -importFrom(ggm, dSep, powerset) +export(aux.effect) +export(generalize) +export(causal.effect) +export(get.expression) +export(meta.transport) +export(parse.graphml) +export(transport) +export(recover) +export(surrogate.outcome) +export(verma.constraints) +importFrom(igraph, get.vertex.attribute) +importFrom(igraph, set.edge.attribute) +importFrom(igraph, vertex.attributes) +importFrom(igraph, topological.sort) +importFrom(igraph, induced.subgraph) +importFrom(igraph, decompose.graph) +importFrom(igraph, edge.attributes) +importFrom(igraph, subgraph.edges) +importFrom(igraph, get.adjacency) +importFrom(igraph, neighborhood) +importFrom(igraph, read.graph) +importFrom(igraph, get.edges) +importFrom(igraph, vertices) +importFrom(igraph, vcount) +importFrom(igraph, is.dag) +importFrom(igraph, edges) +importFrom(igraph, `%->%`) +importFrom(igraph, V) +importFrom(igraph, E) importFrom(stats, setNames) importFrom(utils, combn) \ No newline at end of file diff --git a/NEWS b/NEWS index e85f56d..b752a8c 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,7 @@ +Changes from version 1.3.11 to 1.3.12 + * The package no longer depends on the 'ggm' package. + * The package no longer requires the 'XML' package, now suggests instead. + Changes from version 1.3.10 to 1.3.11 * Fixed inconsistency with function arguments when computing causal effects with surrogate experiments using 'aux.effect'. * Fixed a rare issue with simplification. diff --git a/R/ancestors.R b/R/ancestors.R index 7f46ee0..59b8207 100644 --- a/R/ancestors.R +++ b/R/ancestors.R @@ -1,6 +1,12 @@ ancestors <- function(node, G, topo) { - an.ind <- unique(unlist(neighborhood(G, order = vcount(G), nodes = node, mode = "in"))) - an <- V(G)[an.ind]$name + an.ind <- unique(unlist(igraph::neighborhood(G, order = igraph::vcount(G), nodes = node, mode = "in"))) + an <- igraph::V(G)[an.ind]$name an <- an %ts% topo - return (an) + return(an) +} + +ancestors_unsrt <- function(node, G) { + an.ind <- unique(unlist(igraph::neighborhood(G, order = igraph::vcount(G), nodes = node, mode = "in"))) + an <- igraph::V(G)[an.ind]$name + return(an) } diff --git a/R/c.components.R b/R/c.components.R index 737531e..8c00dd5 100644 --- a/R/c.components.R +++ b/R/c.components.R @@ -1,18 +1,18 @@ c.components <- function(G, topo) { - A <- as.matrix(get.adjacency(G)) - v <- get.vertex.attribute(G, "name") + A <- as.matrix(igraph::get.adjacency(G)) + v <- igraph::get.vertex.attribute(G, "name") indices <- which(A >= 1 & t(A) >= 1, arr.ind = TRUE) bidirected <- NULL - e <- E(G) + e <- igraph::E(G) if (nrow(indices) > 0) { bidirected <- unlist(apply(indices, 1, function(x) { e[v[x[1]] %->% v[x[2]]] })) } - G.bidirected <- subgraph.edges(G, bidirected, delete.vertices = FALSE) - subgraphs <- decompose.graph(G.bidirected) + G.bidirected <- igraph::subgraph.edges(G, bidirected, delete.vertices = FALSE) + subgraphs <- igraph::decompose.graph(G.bidirected) cc <- lapply(subgraphs, function(x) { - v.sub <- get.vertex.attribute(x, "name") + v.sub <- igraph::get.vertex.attribute(x, "name") return(v.sub %ts% topo) }) cc.rank <- order(sapply(cc, function(x) { diff --git a/R/causal.effect.R b/R/causal.effect.R index 301716a..f58d8e9 100644 --- a/R/causal.effect.R +++ b/R/causal.effect.R @@ -1,11 +1,11 @@ causal.effect <- function(y, x, z = NULL, G, expr = TRUE, simp = FALSE, steps = FALSE, primes = FALSE, prune = FALSE, stop_on_nonid = TRUE) { - if (length(edge.attributes(G)) == 0) { - G <- set.edge.attribute(G, "description", 1:length(E(G)), NA) + if (length(igraph::edge.attributes(G)) == 0) { + G <- igraph::set.edge.attribute(G, "description", 1:length(igraph::E(G)), NA) } G.obs <- observed.graph(G) - if (!is.dag(G.obs)) stop("Graph 'G' is not a DAG") - topo <- topological.sort(G.obs) - topo <- get.vertex.attribute(G, "name")[topo] + if (!igraph::is.dag(G.obs)) stop("Graph 'G' is not a DAG") + topo <- igraph::topological.sort(G.obs) + topo <- igraph::get.vertex.attribute(G, "name")[topo] if (length(setdiff(y, topo)) > 0) stop("Set 'y' contains variables not present in the graph.") if (length(setdiff(x, topo)) > 0) stop("Set 'x' contains variables not present in the graph.") if (length(z) > 0 && !identical(z, "")) { @@ -40,11 +40,8 @@ causal.effect <- function(y, x, z = NULL, G, expr = TRUE, simp = FALSE, steps = if (res$tree$call$id) { if (simp) { G.unobs <- unobserved.graph(G) - G.adj <- as.matrix(get.adjacency(G.unobs)) - topo.u <- topological.sort(G.unobs) - topo.u <- get.vertex.attribute(G.unobs, "name")[topo.u] res.prob <- deconstruct(res.prob, probability(), topo) - res.prob <- parse.expression(res.prob, topo, G.adj, G, G.obs) + res.prob <- parse.expression(res.prob, topo, G.unobs, G, G.obs) } attr(res.prob, "algorithm") <- algo attr(res.prob, "query") <- list(y = y, x = x, z = z) diff --git a/R/causal.parents.R b/R/causal.parents.R index 16953ef..7ba3c13 100644 --- a/R/causal.parents.R +++ b/R/causal.parents.R @@ -1,5 +1,5 @@ causal.parents <- function(node, vi, G, G.obs, topo) { - G.vi <- induced.subgraph(G, vi) + G.vi <- igraph::induced.subgraph(G, vi) cc <- c.components(G.vi, topo) t <- Find(function(x) node %in% x, cc) pa.t <- parents(t, G.obs) diff --git a/R/children.R b/R/children.R index 854a214..fa29831 100644 --- a/R/children.R +++ b/R/children.R @@ -1,6 +1,12 @@ children <- function(node, G, topo) { - ch.ind <- unique(unlist(neighborhood(G, order = 1, nodes = node, mode = "out"))) - ch <- V(G)[ch.ind]$name + ch.ind <- unique(unlist(igraph::neighborhood(G, order = 1, nodes = node, mode = "out"))) + ch <- igraph::V(G)[ch.ind]$name ch <- ch %ts% topo return(ch) } + +children_unsrt <- function(node, G) { + ch.ind <- unique(unlist(igraph::neighborhood(G, order = 1, nodes = node, mode = "out"))) + ch <- igraph::V(G)[ch.ind]$name + return(ch) +} \ No newline at end of file diff --git a/R/compare.graphs.R b/R/compare.graphs.R index 48f6ac7..b17cc3b 100644 --- a/R/compare.graphs.R +++ b/R/compare.graphs.R @@ -1,8 +1,8 @@ compare.graphs <- function(G1, G2) { - e1 <- as.data.frame(get.edges(G1, E(G1))) - e1[ ,3] <- edge.attributes(G1) - e2 <- as.data.frame(get.edges(G2, E(G2))) - e2[ ,3] <- edge.attributes(G2) + e1 <- as.data.frame(igraph::get.edges(G1, igraph::E(G1))) + e1[ ,3] <- igraph::edge.attributes(G1) + e2 <- as.data.frame(igraph::get.edges(G2, igraph::E(G2))) + e2[ ,3] <- igraph::edge.attributes(G2) n1 <- nrow(e1) n2 <- nrow(e2) if (n1 != n2) return(FALSE) @@ -10,6 +10,6 @@ compare.graphs <- function(G1, G2) { if (ncol(e2) == 2) e2$description <- "O" e1[which(is.na(e1[,3])), 3] <- "O" e2[which(is.na(e2[,3])), 3] <- "O" - if (all(duplicated(rbind(e1,e2))[(n1+1):(2*n1)])) return(TRUE) + if (all(duplicated(rbind(e1, e2))[(n1+1):(2*n1)])) return(TRUE) return(FALSE) } \ No newline at end of file diff --git a/R/connected.R b/R/connected.R index 495442c..3a3621f 100644 --- a/R/connected.R +++ b/R/connected.R @@ -1,6 +1,6 @@ connected <- function(y, G, topo) { - connected <- unique(unlist(neighborhood(G, order = vcount(G), nodes = y, mode = "all"))) - co <- V(G)[connected]$name + connected <- unique(unlist(igraph::neighborhood(G, order = igraph::vcount(G), nodes = y, mode = "all"))) + co <- igraph::V(G)[connected]$name co <- co %ts% topo return(co) } diff --git a/R/dSep.R b/R/dSep.R new file mode 100644 index 0000000..78c204e --- /dev/null +++ b/R/dSep.R @@ -0,0 +1,110 @@ +# Implements relevant path separation (rp-separation) for testing d-separation. For details, see: +# +# Relevant Path Separation: A Faster Method for Testing Independencies in Bayesian Networks +# Cory J. Butz, Andre E. dos Santos, Jhonatan S. Oliveira; +# Proceedings of the Eighth International Conference on Probabilistic Graphical Models, +# PMLR 52:74-85, 2016. +# +# Note that the roles of Y and Z have been reversed from the paper, meaning that +# we are testing whether X is separated from Y given Z in G. + +dSep <- function(G, x, y, z) { + an_z <- ancestors_unsrt(z, G) + an_xyz <- ancestors_unsrt(union(union(x, y), z), G) + stack_top <- length(x) + stack_size <- max(stack_top, 64) + stack <- rep(FALSE, stack_size) + stack[1:stack_top] <- TRUE + names(stack)[1:stack_top] <- x + visited_top <- 0 + visited_size <- 64 + visited <- rep(FALSE, visited_size) + names(visited) <- rep(NA, visited_size) + is_visited <- FALSE + while (stack_top > 0) { + is_visited <- FALSE + el <- stack[stack_top] + el_name <- names(el) + stack_top <- stack_top - 1 + if (visited_top > 0) { + for (i in 1:visited_top) { + if (el == visited[i] && identical(el_name, names(visited[i]))) { + is_visited <- TRUE + break + } + } + } + if (!is_visited) { + if (el_name %in% y) return(FALSE) + visited_top <- visited_top + 1 + if (visited_top > visited_size) { + visited_old <- visited + visited_size_old <- visited_size + visited_size <- visited_size * 2 + visited <- rep(FALSE, visited_size) + visited[1:visited_size_old] <- visited_old + names(visited[1:visited_size_old]) <- names(visited_old) + } + visited[visited_top] <- el + names(visited)[visited_top] <- el_name + if (el && !(el_name %in% z)) { + visitable_parents <- intersect(setdiff(parents_unsrt(el_name, G), el_name), an_xyz) + visitable_children <- intersect(setdiff(children_unsrt(el_name, G), el_name), an_xyz) + n_vis_pa <- length(visitable_parents) + n_vis_ch <- length(visitable_children) + if (n_vis_pa + n_vis_ch > 0) { + while (n_vis_pa + n_vis_ch + stack_top > stack_size) { + stack_old <- stack + stack_size_old <- stack_size + stack_size <- stack_size * 2 + stack <- rep(FALSE, stack_size) + stack[1:stack_size_old] <- stack_old + names(stack[1:stack_size_old]) <- names(stack_old) + } + stack_add <- stack_top + n_vis_pa + n_vis_ch + stack[(stack_top + 1):(stack_add)] <- c(rep(TRUE, n_vis_pa), rep(FALSE, n_vis_ch)) + names(stack)[(stack_top + 1):(stack_add)] <- c(visitable_parents, visitable_children) + stack_top <- stack_add + } + } else if (!el) { + if (!(el_name %in% z)) { + visitable_children <- intersect(setdiff(children_unsrt(el_name, G), el_name), an_xyz) + n_vis_ch <- length(visitable_children) + if (n_vis_ch > 0) { + while (n_vis_ch + stack_top > stack_size) { + stack_old <- stack + stack_size_old <- stack_size + stack_size <- stack_size * 2 + stack <- rep(FALSE, stack_size) + stack[1:stack_size_old] <- stack_old + names(stack[1:stack_size_old]) <- names(stack_old) + } + stack_add <- stack_top + n_vis_ch + stack[(stack_top + 1):(stack_add)] <- rep(FALSE, n_vis_ch) + names(stack)[(stack_top + 1):(stack_add)] <- visitable_children + stack_top <- stack_add + } + } + if (el_name %in% an_z) { + visitable_parents <- intersect(setdiff(parents_unsrt(el_name, G), el_name), an_xyz) + n_vis_pa <- length(visitable_parents) + if (n_vis_pa > 0) { + while (n_vis_pa + stack_top > stack_size) { + stack_old <- stack + stack_size_old <- stack_size + stack_size <- stack_size * 2 + stack <- rep(FALSE, stack_size) + stack[1:stack_size_old] <- stack_old + names(stack[1:stack_size_old] <- stack_old) + } + stack_add <- stack_top + n_vis_pa + stack[(stack_top + 1):(stack_add)] <- rep(TRUE, n_vis_pa) + names(stack)[(stack_top + 1):(stack_add)] <- visitable_parents + stack_top <- stack_add + } + } + } + } + } + return(TRUE) +} \ No newline at end of file diff --git a/R/descendants.R b/R/descendants.R index 832be7c..8641c78 100644 --- a/R/descendants.R +++ b/R/descendants.R @@ -1,6 +1,6 @@ descendants <- function(node, G, topo) { - de.ind <- unique(unlist(neighborhood(G, order = vcount(G), nodes = node, mode = "out"))) - de <- V(G)[de.ind]$name + de.ind <- unique(unlist(igraph::neighborhood(G, order = igraph::vcount(G), nodes = node, mode = "out"))) + de <- igraph::V(G)[de.ind]$name de <- de %ts% topo return(de) } diff --git a/R/descendent.sets.R b/R/descendent.sets.R index 0b5bbf1..304fabf 100644 --- a/R/descendent.sets.R +++ b/R/descendent.sets.R @@ -15,13 +15,12 @@ descendent.sets <- function(node, s, G.s.obs, topo) { desc <- desc[!vapply(desc, is.null, logical(1))] n.desc <- length(desc) if (n.desc > 0) { - desc.pow <- powerset(1:n.desc, nonempty = TRUE) + desc.pow <- powerset(1:n.desc)[-1] n.sets <- 2^n.desc - 1 D <- vector(mode = "list", length = n.sets) for (i in 1:n.sets) { D[[i]] <- Reduce(union, desc[desc.pow[[i]]]) } - # cat("Descendant sets of ", s, " not containing ", node, " are" , as.character(unique(D)), "\n") return(unique(D)) } return(list()) diff --git a/R/exclusion.restrictions.R b/R/exclusion.restrictions.R index 48bb324..66192fa 100644 --- a/R/exclusion.restrictions.R +++ b/R/exclusion.restrictions.R @@ -1,12 +1,12 @@ exclusion.restrictions <- function(G) { G.obs <- observed.graph(G) - topo <- topological.sort(G.obs) - v <- get.vertex.attribute(G, "name")[topo] + topo <- igraph::topological.sort(G.obs) + v <- igraph::get.vertex.attribute(G, "name")[topo] ex <- lapply(v, function(y) { pa <- setdiff(parents(y, G.obs, topo), y) Z <- setdiff(v, union(y, pa)) if (length(Z) > 0) { - Z.pow <- powerset(setdiff(v, union(y, pa)), nonempty = TRUE) + Z.pow <- powerset(setdiff(v, union(y, pa)))[-1] return(list(pa = pa, Z = Z.pow)) } else return(NULL) }) diff --git a/R/gather.R b/R/gather.R deleted file mode 100644 index d9ce8a4..0000000 --- a/R/gather.R +++ /dev/null @@ -1,14 +0,0 @@ -gather <- function(P, i, G.Adj) { - ind <- c() - vari <- P$children[[i]]$var - for (k in (1:length(P$children))[-i]) { - if (vari %in% P$children[[k]]$cond) { - if (!dSep(G.Adj, P$children[[k]]$var, vari, setdiff(P$children[[k]]$cond, vari))) { - ind <- c(ind, k) - } else { - P$children[[k]]$cond <- setdiff(P$children[[k]]$cond, vari) - } - } - } - return(list(ind, P)) -} \ No newline at end of file diff --git a/R/generalize.R b/R/generalize.R index ca0d7a6..a840f3f 100644 --- a/R/generalize.R +++ b/R/generalize.R @@ -1,21 +1,21 @@ generalize <- function(y, x, Z, D, expr = TRUE, simp = FALSE, steps = FALSE, primes = FALSE, stop_on_nonid = TRUE) { d <- length(D) z <- length(Z) - v <- get.vertex.attribute(D[[1]], "name") - s <- v[which(vertex.attributes(D[[1]])$description == "S")] + v <- igraph::get.vertex.attribute(D[[1]], "name") + s <- v[which(igraph::vertex.attributes(D[[1]])$description == "S")] if (length(s) > 0) stop("The causal diagram cannot contain selection variables.") if (d != z) stop("Number of available experiments does not match number of domains.") if (length(intersect(x, y)) > 0) stop("Sets 'x' and 'y' are not disjoint.") - topo <- lapply(D, function(k) topological.sort(observed.graph(k))) - topo <- lapply(1:d, function(k) get.vertex.attribute(D[[k]], "name")[topo[[k]]]) + topo <- lapply(D, function(k) igraph::topological.sort(observed.graph(k))) + topo <- lapply(1:d, function(k) igraph::get.vertex.attribute(D[[k]], "name")[topo[[k]]]) D <- lapply(D, function(k) { - if (length(edge.attributes(k)) == 0) { - k <- set.edge.attribute(k, "description", 1:length(E(k)), NA) + if (length(igraph::edge.attributes(k)) == 0) { + k <- igraph::set.edge.attribute(k, "description", 1:length(igraph::E(k)), NA) } return(k) }) for (i in 1:d) { - if (!is.dag(observed.graph(D[[i]]))) { + if (!igraph::is.dag(observed.graph(D[[i]]))) { if (i > 1) stop("Selection diagram 'D[", i, "]' is not a DAG.") else stop("Causal diagram 'D[", i, "]' is not a DAG.") } diff --git a/R/generalize_ex.R b/R/generalize_ex.R index 8637097..0afc583 100644 --- a/R/generalize_ex.R +++ b/R/generalize_ex.R @@ -1,21 +1,21 @@ generalize_ex <- function(y, x, Z, D, expr = TRUE, simp = FALSE, steps = FALSE, primes = FALSE, prioritize = FALSE, stop_on_nonid = FALSE) { d <- length(D) z <- length(Z) - v <- get.vertex.attribute(D[[1]], "name") - s <- v[which(vertex.attributes(D[[1]])$description == "S")] + v <- igraph::get.vertex.attribute(D[[1]], "name") + s <- v[which(igraph::vertex.attributes(D[[1]])$description == "S")] if (length(s) > 0) stop("The causal diagram cannot contain selection variables.") if (d != z) stop("Number of available experiments does not match number of domains.") if (length(intersect(x, y)) > 0) stop("Sets 'x' and 'y' are not disjoint.") - topo <- lapply(D, function(k) topological.sort(observed.graph(k))) - topo <- lapply(1:d, function(k) get.vertex.attribute(D[[k]], "name")[topo[[k]]]) + topo <- lapply(D, function(k) igraph::topological.sort(observed.graph(k))) + topo <- lapply(1:d, function(k) igraph::get.vertex.attribute(D[[k]], "name")[topo[[k]]]) D <- lapply(D, function(k) { - if (length(edge.attributes(k)) == 0) { - k <- set.edge.attribute(k, "description", 1:length(E(k)), NA) + if (length(igraph::edge.attributes(k)) == 0) { + k <- igraph::set.edge.attribute(k, "description", 1:length(igraph::E(k)), NA) } return(k) }) for (i in 1:d) { - if (!is.dag(observed.graph(D[[i]]))) { + if (!igraph::is.dag(observed.graph(D[[i]]))) { if (i > 1) stop("Selection diagram 'D[", i, "]' is not a DAG.") else stop("Causal diagram 'D[", i, "]' is not a DAG.") } diff --git a/R/id.R b/R/id.R index 1110195..d56ad54 100644 --- a/R/id.R +++ b/R/id.R @@ -23,7 +23,7 @@ id <- function(y, x, P, G, G.obs, v, topo, tree) { # line 2 if (length(setdiff(v, an)) != 0) { - G.an <- induced.subgraph(G, an) + G.an <- igraph::induced.subgraph(G, an) G.an.obs <- observed.graph(G.an) if (P$product | P$fraction) { P$sumset <- union(setdiff(v, an), P$sumset) %ts% topo @@ -40,7 +40,7 @@ id <- function(y, x, P, G, G.obs, v, topo, tree) { } # line 3 - G.xbar <- subgraph.edges(G, E(G)[!(to(x) | (from(x) & (description == "U" & !is.na(description))))], delete.vertices = FALSE) + G.xbar <- igraph::subgraph.edges(G, igraph::E(G)[!(to(x) | (from(x) & (description == "U" & !is.na(description))))], delete.vertices = FALSE) an.xbar <- ancestors(y, observed.graph(G.xbar), topo) w <- setdiff(setdiff(v, x), an.xbar) w.len <- length(w) @@ -55,7 +55,7 @@ id <- function(y, x, P, G, G.obs, v, topo, tree) { } # line 4 - G.remove.x <- induced.subgraph(G, v[!(v %in% x)]) + G.remove.x <- igraph::induced.subgraph(G, v[!(v %in% x)]) s <- c.components(G.remove.x, topo) if (length(s) > 1) { tree$call$line <- 4 @@ -129,7 +129,7 @@ id <- function(y, x, P, G, G.obs, v, topo, tree) { tree$call$s.prime <- s s.len <- length(s) ind <- which(v %in% s) - G.s <- induced.subgraph(G, s) + G.s <- igraph::induced.subgraph(G, s) G.s.obs <- observed.graph(G.s) product.list <- vector(mode = "list", length = s.len) P.prod <- probability() diff --git a/R/idc.R b/R/idc.R index 45e48e9..be10b00 100644 --- a/R/idc.R +++ b/R/idc.R @@ -6,14 +6,13 @@ idc <- function(y, x, z, P, G, G.obs, v, topo, tree, prune) { offset <- (prune) * 2 from <- NULL G.xz <- unobserved.graph(G) - edges.to.x <- E(G.xz)[to(x)] - edges.from.z <- E(G.xz)[from(z)] - G.xz <- subgraph.edges(G.xz, E(G.xz)[setdiff(E(G.xz), union(edges.to.x, edges.from.z))], delete.vertices = FALSE) - A <- as.matrix(get.adjacency(G.xz)) + edges.to.x <- igraph::E(G.xz)[to(x)] + edges.from.z <- igraph::E(G.xz)[from(z)] + G.xz <- igraph::subgraph.edges(G.xz, igraph::E(G.xz)[setdiff(igraph::E(G.xz), union(edges.to.x, edges.from.z))], delete.vertices = FALSE) nxt <- list() for (node in z) { cond <- setdiff(z, node) - if (dSep(A, y, node, union(x, cond))) { + if (wrap.dSep(G.xz, y, node, union(x, cond))) { tree$call$line <- 9 + offset tree$call$z.prime <- node nxt <- idc(y, union(x, node) %ts% topo, cond, P, G, G.obs, v, topo, list(), prune) diff --git a/R/identify.R b/R/identify.R index fd5a9f2..e3e5d97 100644 --- a/R/identify.R +++ b/R/identify.R @@ -1,9 +1,9 @@ identify <- function(C, T, Q, G, topo, tree) { - v <- get.vertex.attribute(G, "name") - s <- v[which(vertex.attributes(G)$description == "S")] + v <- igraph::get.vertex.attribute(G, "name") + s <- v[which(igraph::vertex.attributes(G)$description == "S")] v <- v %ts% topo G.obs <- observed.graph(G) - G.T <- induced.subgraph(G, T) + G.T <- igraph::induced.subgraph(G, T) G.T.obs <- observed.graph(G.T) tree$call <- list(y = C, x = setdiff(v, C), C = C, T = T, P = activate.selection.variable(Q, s), G = G.T, line = "", v = v, alg = "Identify", id = FALSE) anc.c <- ancestors(C, G.T.obs, topo) @@ -33,7 +33,7 @@ identify <- function(C, T, Q, G, topo, tree) { # iii) if (all(C %in% A) && all(A %in% T)) { - G.A <- induced.subgraph(G, A) + G.A <- igraph::induced.subgraph(G, A) cc <- c.components(G.A, topo) T.prime <- Find(function(x) all(C %in% x), cc) T.one <- intersect(T.prime, A) diff --git a/R/independence.restrictions.R b/R/independence.restrictions.R index a052814..5d58783 100644 --- a/R/independence.restrictions.R +++ b/R/independence.restrictions.R @@ -1,7 +1,7 @@ independence.restrictions <- function(G) { G.obs <- observed.graph(G) - topo <- topological.sort(G.obs) - v <- get.vertex.attribute(G, "name")[topo] + topo <- igraph::topological.sort(G.obs) + v <- igraph::get.vertex.attribute(G, "name")[topo] cc <- c.components(G, v) cc.len <- length(cc) if (cc.len > 1) { diff --git a/R/insert.R b/R/insert.R index 932feaf..b509801 100644 --- a/R/insert.R +++ b/R/insert.R @@ -1,30 +1,20 @@ -insert <- function(J, D, M, cond, S, O, G.adj, G, G.obs, topo) { +insert <- function(J, D, M, cond, S, O, G.unobs, G, G.obs, topo) { mis.ind <- which(M %in% D) if (length(mis.ind) == 0) return(list(J, D)) mis <- M[mis.ind] M <- mis[length(mis)] if (M %in% cond) return(list(J, D)) - # non.mis <- min(which(J %in% O)) - # V.prev <- J[non.mis] - # ind <- which(topo == V.prev) - # V.pi <- topo[0:(ind-1)] - # anc <- ancestors(M, G.obs, topo) - # cond.diff <- setdiff(V.pi, anc) - # ds <- powerset(cond.diff, nonempty = FALSE) J.min <- min(which(J %in% topo)) V.prev <- J[J.min] ind <- which(topo == V.prev) V.pi <- topo[0:(ind-1)] - ds <- powerset(setdiff(V.pi, M), nonempty = FALSE) + ds <- powerset(setdiff(V.pi, M)) n <- length(ds) for (i in 1:n) { - # add <- union(anc, ds[[i]]) - # add.M <- setdiff(add, M) - # a.set <- union(setdiff(add, D), setdiff(D, add)) add <- union(ds[[i]], M) a.set <- union(setdiff(add, D), setdiff(D, add)) - if (wrap.dSep(G.adj, J, a.set, setdiff(D, a.set)) && - wrap.dSep(G.adj, M, S, setdiff(ds[[i]], S))) { + if (wrap.dSep(G.unobs, J, a.set, setdiff(D, a.set)) && + wrap.dSep(G.unobs, M, S, setdiff(ds[[i]], S))) { J.new <- union(J, M) D.new <- ds[[i]] return(list(J.new, D.new, M, ds[[i]])) diff --git a/R/irrelevant.R b/R/irrelevant.R index 9c57b58..7258a6e 100644 --- a/R/irrelevant.R +++ b/R/irrelevant.R @@ -1,10 +1,10 @@ -irrelevant <- function(children, vari, sumset, G.adj) { +irrelevant <- function(children, vari, sumset, G.unobs) { n <- length(children) irrel <- NULL for (i in 1:n) { ch <- children[[i]] if (!(ch$var %in% sumset)) { - if (wrap.dSep(G.adj, ch$var, vari, setdiff(ch$cond, vari))) irrel <- c(irrel, i) + if (wrap.dSep(G.unobs, ch$var, vari, setdiff(ch$cond, vari))) irrel <- c(irrel, i) } } return(irrel) diff --git a/R/join.R b/R/join.R index a4e461b..8ba5729 100644 --- a/R/join.R +++ b/R/join.R @@ -1,4 +1,4 @@ -join <- function(J, D, vari, cond, S, M, O, G.adj, G, G.obs, topo) { +join <- function(J, D, vari, cond, S, M, O, G.unobs, G, G.obs, topo) { J.new <- character() D.new <- character() if (length(J) == 0) { @@ -6,35 +6,25 @@ join <- function(J, D, vari, cond, S, M, O, G.adj, G, G.obs, topo) { D.new <- cond return(list(J.new, D.new)) } - # non.mis <- min(which(J %in% O)) - # V.prev <- J[non.mis] - # ind <- which(topo == V.prev) - # anc <- ancestors(vari, G.obs, topo) - # cond.diff <- setdiff(V.pi, anc) - # ds <- powerset(cond.diff, nonempty = FALSE) J.min <- min(which(J %in% topo)) V.prev <- J[J.min] ind <- which(topo == V.prev) V.pi <- topo[0:(ind-1)] - ds <- powerset(setdiff(V.pi, vari), nonempty = FALSE) + ds <- powerset(setdiff(V.pi, vari)) n <- length(ds) for (i in 1:n) { - # add <- union(anc, ds[[i]]) - # add.v <- setdiff(add, vari) - # a.set <- union(setdiff(add, D), setdiff(D, add)) - # b.set <- union(setdiff(add.v, cond), setdiff(cond, add.v)) add <- union(ds[[i]], vari) a.set <- union(setdiff(add, D), setdiff(D, add)) b.set <- union(setdiff(ds[[i]], cond), setdiff(cond, ds[[i]])) - if (wrap.dSep(G.adj, J, a.set, setdiff(D, a.set)) && - wrap.dSep(G.adj, vari, b.set, setdiff(cond, b.set))) { + if (wrap.dSep(G.unobs, J, a.set, setdiff(D, a.set)) && + wrap.dSep(G.unobs, vari, b.set, setdiff(cond, b.set))) { J.new <- union(J, vari) D.new <- ds[[i]] return(list(J.new, D.new)) } } if (any(M %in% D)) { - joint <- insert(J, D, M, cond, S, O, G.adj, G, G.obs, topo) + joint <- insert(J, D, M, cond, S, O, G.unobs, G, G.obs, topo) if (length(joint[[1]]) > length(J)) { return(joint) } diff --git a/R/latent.projection.R b/R/latent.projection.R index ac56fe7..80bee45 100644 --- a/R/latent.projection.R +++ b/R/latent.projection.R @@ -3,33 +3,33 @@ latent.projection <- function(G, l) { from <- NULL description <- NULL for (i in 1:length(l)) { - e <- E(G) - v <- get.vertex.attribute(G, "name") + e <- igraph::E(G) + v <- igraph::get.vertex.attribute(G, "name") inc.edges <- e[to(l[i]) & (is.na(description) | description != "U")] out.edges <- e[from(l[i]) & (is.na(description) | description != "U")] unobs.edges <- e[to(l[i]) & description == "U" & !is.na(description)] - inc.ind <- get.edges(G, inc.edges)[ ,1] - out.ind <- get.edges(G, out.edges)[ ,2] - unobs.ind <- setdiff(get.edges(G, unobs.edges)[ ,1], out.ind) + inc.ind <- igraph::get.edges(G, inc.edges)[ ,1] + out.ind <- igraph::get.edges(G, out.edges)[ ,2] + unobs.ind <- setdiff(igraph::get.edges(G, unobs.edges)[ ,1], out.ind) inc.len <- length(inc.ind) out.len <- length(out.ind) unobs.len <- length(unobs.ind) if (inc.len > 0 & out.len > 0) { obs.new <- t(as.matrix(expand.grid(inc.ind, out.ind))) - G <- G + edges(v[c(obs.new)], description = rep(NA, ncol(obs.new))) # replace path v_1 -> L -> v_2 with v_1 -> v_2 + G <- G + igraph::edges(v[c(obs.new)], description = rep(NA, ncol(obs.new))) # replace path v_1 -> L -> v_2 with v_1 -> v_2 } if (out.len > 1) { unobs.new <- combn(out.ind, 2) - G <- G + edges(v[c(unobs.new, unobs.new[2:1, ])], description = rep("U", 2 * ncol(unobs.new))) # replace path v_1 <- L -> v_2 with v_1 <-> v_2 + G <- G + igraph::edges(v[c(unobs.new, unobs.new[2:1, ])], description = rep("U", 2 * ncol(unobs.new))) # replace path v_1 <- L -> v_2 with v_1 <-> v_2 } if (unobs.len > 0 & out.len > 0) { unobs.old <- t(as.matrix(expand.grid(unobs.ind, out.ind))) - G <- G + edges(v[c(unobs.old, unobs.old[2:1, ])], description = rep("U", 2 * ncol(unobs.old))) # replace path v_1 <-> L -> v_2 with v_1 <-> v_2 + G <- G + igraph::edges(v[c(unobs.old, unobs.old[2:1, ])], description = rep("U", 2 * ncol(unobs.old))) # replace path v_1 <-> L -> v_2 with v_1 <-> v_2 } - G <- induced.subgraph(G, setdiff(v, l[i])) - e.dat <- as.data.frame(get.edges(G, E(G))) - e.dat[ ,3] <- edge.attributes(G) - G <- subgraph.edges(G, which(!duplicated(e.dat)), delete.vertices = FALSE) + G <- igraph::induced.subgraph(G, setdiff(v, l[i])) + e.dat <- as.data.frame(igraph::get.edges(G, E(G))) + e.dat[ ,3] <- igraph::edge.attributes(G) + G <- igraph::subgraph.edges(G, which(!duplicated(e.dat)), delete.vertices = FALSE) } return(G) } \ No newline at end of file diff --git a/R/meta.transport.R b/R/meta.transport.R index 41dc50f..1bffa15 100644 --- a/R/meta.transport.R +++ b/R/meta.transport.R @@ -1,8 +1,8 @@ meta.transport <- function(y, x, D, expr = TRUE, simp = TRUE, steps = FALSE, primes = FALSE, stop_on_nonid = TRUE) { - v <- get.vertex.attribute(D[[1]], "name") - s <- v[which(vertex.attributes(D[[1]])$description == "S")] + v <- igraph::get.vertex.attribute(D[[1]], "name") + s <- v[which(igraph::vertex.attributes(D[[1]])$description == "S")] interventions <- setdiff(v, union(y, s)) - D.causal <- induced.subgraph(D[[1]], v[!(v %in% s)]) + D.causal <- igraph::induced.subgraph(D[[1]], v[!(v %in% s)]) D.all <- list() D.all[[1]] <- D.causal D.all[2:(length(D)+1)] <- D diff --git a/R/observed.graph.R b/R/observed.graph.R index 0e0c27c..545207d 100644 --- a/R/observed.graph.R +++ b/R/observed.graph.R @@ -1,5 +1,5 @@ observed.graph <- function(G) { - obs.edges <- setdiff(E(G), E(G)[which(edge.attributes(G)$description == "U")]) - G.obs <- subgraph.edges(G, E(G)[obs.edges], delete.vertices = FALSE) + obs.edges <- setdiff(igraph::E(G), igraph::E(G)[which(igraph::edge.attributes(G)$description == "U")]) + G.obs <- igraph::subgraph.edges(G, igraph::E(G)[obs.edges], delete.vertices = FALSE) return(G.obs) } diff --git a/R/parents.R b/R/parents.R index e353d45..98b61f2 100644 --- a/R/parents.R +++ b/R/parents.R @@ -1,6 +1,12 @@ parents <- function(node, G.obs, topo) { - pa.ind <- unique(unlist(neighborhood(G.obs, order = 1, nodes = node, mode = "in"))) - pa <- V(G.obs)[pa.ind]$name + pa.ind <- unique(unlist(igraph::neighborhood(G.obs, order = 1, nodes = node, mode = "in"))) + pa <- igraph::V(G.obs)[pa.ind]$name pa <- pa %ts% topo return(pa) +} + +parents_unsrt <- function(node, G.obs) { + pa.ind <- unique(unlist(igraph::neighborhood(G.obs, order = 1, nodes = node, mode = "in"))) + pa <- igraph::V(G.obs)[pa.ind]$name + return(pa) } \ No newline at end of file diff --git a/R/parse.expression.R b/R/parse.expression.R index dfb1383..0f26c7a 100644 --- a/R/parse.expression.R +++ b/R/parse.expression.R @@ -1,8 +1,8 @@ -parse.expression <- function(P, topo, G.adj, G, G.obs) { +parse.expression <- function(P, topo, G.unobs, G, G.obs) { if (P$fraction) { P <- cancel.out(P) if (P$fraction) { - P$den <- parse.expression(P$den, topo, G.adj, G, G.obs) + P$den <- parse.expression(P$den, topo, G.unobs, G, G.obs) if (length(P$den) == 0) { sum_p <- P$sumset P <- P$num @@ -23,7 +23,7 @@ parse.expression <- function(P, topo, G.adj, G, G.obs) { P$sumset <- setdiff(P$sumset, nodep) %ts% topo } } - P$num <- parse.expression(P$num, topo, G.adj, G, G.obs) + P$num <- parse.expression(P$num, topo, G.unobs, G, G.obs) P <- cancel.out(P) } return(P) @@ -35,7 +35,7 @@ parse.expression <- function(P, topo, G.adj, G, G.obs) { parse_children <- P$children[non_atomic] P$children <- P$children[!non_atomic] for (i in 1:length(parse_children)) { - P.parse <- parse.expression(parse_children[[i]], topo, G.adj, G, G.obs) + P.parse <- parse.expression(parse_children[[i]], topo, G.unobs, G, G.obs) if (!is.null(P.parse$collapse)) { P$children <- c(P$children, P.parse$children) } else { @@ -58,7 +58,7 @@ parse.expression <- function(P, topo, G.adj, G, G.obs) { ord.sum <- order(sapply(P$sumset, FUN = function(x) which(topo == x)), decreasing = TRUE) P$children <- P$children[ord.children] P$sumset <- P$sumset[ord.sum] - P <- simplify(P, topo, G.adj, G, G.obs) + P <- simplify(P, topo, G.unobs, G, G.obs) if (length(P$children) == 0) return(NULL) } P.parse <- probability(product = TRUE, children = list()) @@ -82,7 +82,7 @@ parse.expression <- function(P, topo, G.adj, G, G.obs) { sum_p <- P$sumset P <- P$children[[1]] P$sumset <- union(sum_p, P$sumset) %ts% topo - P <- parse.expression(P, topo, G.adj, G, G.obs) + P <- parse.expression(P, topo, G.unobs, G, G.obs) } } if (length(P$children) > 0) P.parse$children[[j + 1]] <- P diff --git a/R/parse.graphml.R b/R/parse.graphml.R index 2023240..0ec3dec 100644 --- a/R/parse.graphml.R +++ b/R/parse.graphml.R @@ -1,7 +1,10 @@ parse.graphml <- function(file, format = c("standard", "internal"), nodes = c(), use.names = TRUE) { - format <- match.arg(format) - res <- switch(format, standard = parse.graphml.standard(file, nodes, use.names), - internal = parse.graphml.internal(file, nodes, use.names), - stop(paste("Unknown file format:", format))) - return(res) + if (!requireNamespace("XML", quietly = TRUE)) message("Package 'XML' is required for importing GraphML files.") + else { + format <- match.arg(format) + res <- switch(format, standard = parse.graphml.standard(file, nodes, use.names), + internal = parse.graphml.internal(file, nodes, use.names), + stop(paste("Unknown file format:", format))) + return(res) + } } diff --git a/R/parse.graphml.internal.R b/R/parse.graphml.internal.R index a6d8c53..cbb46ab 100644 --- a/R/parse.graphml.internal.R +++ b/R/parse.graphml.internal.R @@ -1,24 +1,24 @@ parse.graphml.internal <- function(file, nodes, use.names) { - doc <- xmlParse(file, useInternalNodes = TRUE) - top <- xmlRoot(doc) + doc <- XML::xmlParse(file, useInternalNodes = TRUE) + top <- XML::xmlRoot(doc) graph <- top[["graph"]] if (!use.names) { - if (xmlSize(which(xmlSApply(graph, xmlName) == "node")) != length(nodes)) stop("Incorrect number of node names") + if (XML::xmlSize(which(XML::xmlSApply(graph, XML::xmlName) == "node")) != length(nodes)) stop("Incorrect number of node names") } ns <- c(ns = "http://graphml.graphdrawing.org/xmlns") if (use.names) { - node.data <- getNodeSet(doc, "//ns:data[contains(@key,'d6')]/*", ns) - for (i in 1:length(node.data)) nodes[i] <- xmlValue(node.data[[i]]["NodeLabel"]$NodeLabel[1]$text) + node.data <- XML::getNodeSet(doc, "//ns:data[contains(@key,'d6')]/*", ns) + for (i in 1:length(node.data)) nodes[i] <- XML::xmlValue(node.data[[i]]["NodeLabel"]$NodeLabel[1]$text) } - all <- getNodeSet(doc, "//*[position() > 1]", ns) - keep <- getNodeSet(doc, "//ns:edge | //ns:node | //ns:graph | //ns:key[@attr.name = 'description'] | //ns:data[contains(@key,'d9')]", ns) - remove <- getNodeSet(doc, "//ns:key[@for='port'] | //ns:key[@for='graphml'] | //ns:data[@key='d4'] | //ns:data[@key='d6'] + all <- XML::getNodeSet(doc, "//*[position() > 1]", ns) + keep <- XML::getNodeSet(doc, "//ns:edge | //ns:node | //ns:graph | //ns:key[@attr.name = 'description'] | //ns:data[contains(@key,'d9')]", ns) + remove <- XML::getNodeSet(doc, "//ns:key[@for='port'] | //ns:key[@for='graphml'] | //ns:data[@key='d4'] | //ns:data[@key='d6'] | //ns:data[@key='d7'] | //ns:data[@key='d8'] | //ns:data[@key='d10']" , ns) - removeNodes(union(setdiff(all, keep), remove)) + XML::removeNodes(union(setdiff(all, keep), remove)) temp <- tempfile(fileext = ".graphml") - temp.xml <- saveXML(doc, file = temp) - free(doc) - igrph <- read.graph(temp.xml, format = "graphml") - igrph <- set.vertex.attribute(igrph, "name", value = nodes) + temp.xml <- XML::saveXML(doc, file = temp) + XML::free(doc) + igrph <- igraph::read.graph(temp.xml, format = "graphml") + igrph <- igraph::set.vertex.attribute(igrph, "name", value = nodes) return(igrph) } diff --git a/R/parse.graphml.standard.R b/R/parse.graphml.standard.R index ba04c86..791239a 100644 --- a/R/parse.graphml.standard.R +++ b/R/parse.graphml.standard.R @@ -1,23 +1,23 @@ parse.graphml.standard <- function(file, nodes, use.names) { - doc <- xmlParse(file, useInternalNodes = TRUE) - top <- xmlRoot(doc) + doc <- XML::xmlParse(file, useInternalNodes = TRUE) + top <- XML::xmlRoot(doc) graph <- top[["graph"]] if (!use.names) { - if (xmlSize(which(xmlSApply(graph, xmlName) == "node")) != length(nodes)) stop("Incorrect number of node names") + if (XML::xmlSize(which(XML::xmlSApply(graph, XML::xmlName) == "node")) != length(nodes)) stop("Incorrect number of node names") } ns <- c(ns = "http://graphml.graphdrawing.org/xmlns") - edges <- which(xmlSApply(graph, xmlName) == "edge") + edges <- which(XML::xmlSApply(graph, XML::xmlName) == "edge") remove.id <- 0 removals <- list() for (i in 1:length(edges)) { current.edge <- graph[[edges[i]]] - edge.attr <- xmlAttrs(current.edge) - datas <- which(xmlSApply(current.edge, xmlName) == "data") + edge.attr <- XML::xmlAttrs(current.edge) + datas <- which(XML::xmlSApply(current.edge, XML::xmlName) == "data") for (j in 1:length(datas)) { current.data <- current.edge[[datas[j]]] - if (xmlAttrs(current.data) == "d10") { + if (XML::xmlAttrs(current.data) == "d10") { arc <- current.data[[1]] - arw.attr <- xmlAttrs(arc[["Arrows"]]) + arw.attr <- XML::xmlAttrs(arc[["Arrows"]]) src <- edge.attr[["source"]] trgt <- edge.attr[["target"]] two.arrows <- !identical(arw.attr[["source"]], "none") & !identical(arw.attr[["target"]], "none") @@ -25,29 +25,28 @@ parse.graphml.standard <- function(file, nodes, use.names) { if (two.arrows | no.arrows) { remove.id <- remove.id + 1 removals[[remove.id]] <- current.edge - e1 <- newXMLNode("edge", parent = doc, attrs = c(id = "e", source = src, target = trgt), - newXMLNode("data", attrs = c(key = "d9"), cdata = TRUE, "U" )) - e2 <- newXMLNode("edge", parent = doc, attrs = c(id = "e", source = trgt, target = src), - newXMLNode("data", attrs = c(key = "d9"), cdata = TRUE, "U" )) - addChildren(graph, kids = list(e1, e2)) + e1 <- XML::newXMLNode("edge", parent = doc, attrs = c(id = "e", source = src, target = trgt), + XML::newXMLNode("data", attrs = c(key = "d9"), cdata = TRUE, "U" )) + e2 <- XML::newXMLNode("edge", parent = doc, attrs = c(id = "e", source = trgt, target = src), + XML::newXMLNode("data", attrs = c(key = "d9"), cdata = TRUE, "U" )) + XML::addChildren(graph, kids = list(e1, e2)) } } } } - removeNodes(removals) + XML::removeNodes(removals) if (use.names) { - node.data <- getNodeSet(doc, "//ns:data[contains(@key,'d6')]/*", ns) - for (i in 1:length(node.data)) nodes[i] <- xmlValue(node.data[[i]]["NodeLabel"]$NodeLabel[1]$text) + node.data <- XML::getNodeSet(doc, "//ns:data[contains(@key,'d6')]/*", ns) + for (i in 1:length(node.data)) nodes[i] <- XML::xmlValue(node.data[[i]]["NodeLabel"]$NodeLabel[1]$text) } - all <- getNodeSet(doc, "//*[position() > 1]", ns) - keep <- getNodeSet(doc, "//ns:edge | //ns:node | //ns:graph | //ns:key[@attr.name = 'description'] | //ns:data[contains(@key,'d9')]", ns) - remove <- getNodeSet(doc, "//ns:key[@for='port'] | //ns:key[@for='graphml'] | //ns:data[@key='d4'] | //ns:data[@key='d6'] - | //ns:data[@key='d7'] | //ns:data[@key='d8'] | //ns:data[@key='d10']" , ns) - removeNodes(union(setdiff(all, keep), remove)) + all <- XML::getNodeSet(doc, "//*[position() > 1]", ns) + keep <- XML::getNodeSet(doc, "//ns:edge | //ns:node | //ns:graph | //ns:key[@attr.name = 'description'] | //ns:data[contains(@key,'d9')]", ns) + remove <- XML::getNodeSet(doc, "//ns:key[@for='port'] | //ns:key[@for='graphml'] | //ns:data[@key='d4'] | //ns:data[@key='d6'] | //ns:data[@key='d7'] | //ns:data[@key='d8'] | //ns:data[@key='d10']" , ns) + XML::removeNodes(union(setdiff(all, keep), remove)) temp <- tempfile(fileext = ".graphml") - temp.xml <- saveXML(doc, file = temp) - free(doc) - igrph <- read.graph(temp.xml, format = "graphml") - igrph <- set.vertex.attribute(igrph, "name", value = nodes) + temp.xml <- XML::saveXML(doc, file = temp) + XML::free(doc) + igrph <- igraph::read.graph(temp.xml, format = "graphml") + igrph <- igraph::set.vertex.attribute(igrph, "name", value = nodes) return(igrph) } diff --git a/R/pid.R b/R/pid.R index 54e8b09..7a48e97 100644 --- a/R/pid.R +++ b/R/pid.R @@ -22,7 +22,7 @@ pid <- function(y, x, P, G, G.obs, v, topo, tree) { # line 2 if (length(setdiff(v, an)) != 0) { - G.an <- induced.subgraph(G, an) + G.an <- igraph::induced.subgraph(G, an) G.an.obs <- observed.graph(G.an) if (P$product | P$fraction) { P$sumset <- union(setdiff(v, an), P$sumset) %ts% topo @@ -38,12 +38,12 @@ pid <- function(y, x, P, G, G.obs, v, topo, tree) { } # line 3 - G.xbar <- subgraph.edges(G, E(G)[!(to(x) | (from(x) & (description == "U" & !is.na(description))))], delete.vertices = FALSE) + G.xbar <- igraph::subgraph.edges(G, igraph::E(G)[!(to(x) | (from(x) & (description == "U" & !is.na(description))))], delete.vertices = FALSE) co <- connected(y, G.xbar, topo) z <- setdiff(an, co) if (length(z) > 0) { v.new <- setdiff(v, z) %ts% topo - G.z <- induced.subgraph(G, v.new) + G.z <- igraph::induced.subgraph(G, v.new) G.z.obs <- observed.graph(G.z) G.z.l <- latent.projection(G, z) if (compare.graphs(G.z, G.z.l)) { @@ -68,13 +68,13 @@ pid <- function(y, x, P, G, G.obs, v, topo, tree) { t <- c() for (i in 1:length(v.x)) { r <- setdiff(ancestors(v.x[i], G.xbar.obs, topo), de) - G.vbar <- subgraph.edges(G, E(G)[!(to(v.x[i]) | (from(v.x[i]) & (description == "U" & !is.na(description))))], delete.vertices = FALSE) + G.vbar <- igraph::subgraph.edges(G, igraph::E(G)[!(to(v.x[i]) | (from(v.x[i]) & (description == "U" & !is.na(description))))], delete.vertices = FALSE) co.vi <- connected(setdiff(v, r), G.vbar, topo) t <- union(t, setdiff(r, co.vi)) } if (length(t) > 0) { v.new <- setdiff(v, t) %ts% topo - G.t <- induced.subgraph(G, v.new) + G.t <- igraph::induced.subgraph(G, v.new) G.t.obs <- observed.graph(G.t) if (P$product | P$fraction) { P$sumset <- union(t, P$sumset) %ts% topo @@ -93,7 +93,7 @@ pid <- function(y, x, P, G, G.obs, v, topo, tree) { if (length(x) == 1) { cc <- c.components(G, topo) s <- Find(function(z) x %in% z, cc) - G.s <- induced.subgraph(G, s) + G.s <- igraph::induced.subgraph(G, s) G.s.obs <- observed.graph(G.s) ch <- setdiff(children(x, G.s.obs, topo), x) if (length(ch) == 0) { @@ -105,7 +105,7 @@ pid <- function(y, x, P, G, G.obs, v, topo, tree) { G.prime <- latent.projection(G, v.xy[i]) cc.prime <- c.components(G.prime, topo) s.prime <- Find(function(z) x %in% z, cc.prime) - G.s.prime <- induced.subgraph(G.prime, s.prime) + G.s.prime <- igraph::induced.subgraph(G.prime, s.prime) G.s.prime.obs <- observed.graph(G.s.prime) ch.prime <- setdiff(children(x, G.s.prime.obs, topo), x) if (length(ch.prime) == 0) { @@ -149,7 +149,7 @@ pid <- function(y, x, P, G, G.obs, v, topo, tree) { } # line 7 - G.remove.x <- induced.subgraph(G, v[!(v %in% x)]) + G.remove.x <- igraph::induced.subgraph(G, v[!(v %in% x)]) s <- c.components(G.remove.x, topo) if (length(s) > 1) { tree$call$line <- 7 @@ -223,7 +223,7 @@ pid <- function(y, x, P, G, G.obs, v, topo, tree) { tree$call$s.prime <- s s.len <- length(s) ind <- which(v %in% s) - G.s <- induced.subgraph(G, s) + G.s <- igraph::induced.subgraph(G, s) G.s.obs <- observed.graph(G.s) product.list <- vector(mode = "list", length = s.len) P.prod <- probability() diff --git a/R/powerset.R b/R/powerset.R new file mode 100644 index 0000000..f0279a1 --- /dev/null +++ b/R/powerset.R @@ -0,0 +1,7 @@ +powerset <- function(set) { + n <- length(set) + if (n == 0) return(list(c())) + indices <- sapply(0:(2^n-1), function(p) as.logical(intToBits(p)[1:n]), simplify = FALSE) + subsets <- lapply(indices, function(i) set[i]) + return(subsets) +} \ No newline at end of file diff --git a/R/q.constraints.R b/R/q.constraints.R index f096a25..8e45402 100644 --- a/R/q.constraints.R +++ b/R/q.constraints.R @@ -1,6 +1,6 @@ q.constraints <- function(s, node, G, G.obs, G.unobs, topo, topo.u, constraints) { - G.s <- induced.subgraph(G, s) - v <- get.vertex.attribute(G, "name") + G.s <- igraph::induced.subgraph(G, s) + v <- igraph::get.vertex.attribute(G, "name") v <- v %ts% topo G.s.obs <- observed.graph(G.s) desc.sets <- descendent.sets(node, s, G.s.obs, topo) @@ -65,8 +65,8 @@ q.constraints <- function(s, node, G, G.obs, G.unobs, topo, topo.u, constraints) ))) } d.prime <- s_d - G.d <- induced.subgraph(G.s, d.prime) - v <- get.vertex.attribute(G.d, "name") + G.d <- igraph::induced.subgraph(G.s, d.prime) + v <- igraph::get.vertex.attribute(G.d, "name") v <- v %ts% topo cc.d <- c.components(G.d, topo) if (length(cc.d) > 1) e <- Find(function(x) node %in% x, cc.d) diff --git a/R/rc.R b/R/rc.R index 0e7ed33..af4a181 100644 --- a/R/rc.R +++ b/R/rc.R @@ -1,8 +1,8 @@ rc <- function(D, P, G, topo, tree) { - v.s <- get.vertex.attribute(G, "name") - s <- v.s[which(vertex.attributes(G)$description == "S")] - G.causal <- induced.subgraph(G, v.s[!(v.s %in% s)]) - v <- get.vertex.attribute(G.causal, "name") + v.s <- igraph::get.vertex.attribute(G, "name") + s <- v.s[which(igraph::vertex.attributes(G)$description == "S")] + G.causal <- igraph::induced.subgraph(G, v.s[!(v.s %in% s)]) + v <- igraph::get.vertex.attribute(G.causal, "name") v <- v %ts% topo G.obs <- observed.graph(G.causal) G.s.obs <- observed.graph(G) @@ -19,7 +19,7 @@ rc <- function(D, P, G, topo, tree) { } else { P$var <- anc.union } - nxt <- rc(D, P, induced.subgraph(G, anc.union), topo, list()) + nxt <- rc(D, P, igraph::induced.subgraph(G, anc.union), topo, list()) tree$call$line <- 2 tree$call$id <- nxt$tree$call$id tree$call$anc.d <- anc.d @@ -72,7 +72,7 @@ rc <- function(D, P, G, topo, tree) { } else { P.new$den <- product.list[[1]] } - nxt <- rc(D, P.new, induced.subgraph(G, setdiff(v.s, c.set)), topo, list()) + nxt <- rc(D, P.new, igraph::induced.subgraph(G, setdiff(v.s, c.set)), topo, list()) tree$call$line <- 6 tree$call$id <- nxt$tree$call$id tree$call$c.set <- c.set diff --git a/R/recover.R b/R/recover.R index af58cbe..54867ab 100644 --- a/R/recover.R +++ b/R/recover.R @@ -1,28 +1,28 @@ recover <- function(y, x, G, expr = TRUE, simp = TRUE, steps = FALSE, primes = FALSE, stop_on_nonid = TRUE) { - if (length(edge.attributes(G)) == 0) { - G <- set.edge.attribute(G, "description", 1:length(E(G)), NA) + if (length(igraph::edge.attributes(G)) == 0) { + G <- igraph::set.edge.attribute(G, "description", 1:length(igraph::E(G)), NA) } - if (!is.dag(observed.graph(G))) stop("Graph 'G' is not a DAG") - sel <- which(vertex.attributes(G)$description == "S") + if (!igraph::is.dag(observed.graph(G))) stop("Graph 'G' is not a DAG") + sel <- which(igraph::vertex.attributes(G)$description == "S") if (length(sel) == 0) stop("No selection variables present in the diagram.") if (length(sel) > 1) stop("Multiple selection variables are not supported.") G.obs <- observed.graph(G) - topo <- topological.sort(G.obs) - v.s <- get.vertex.attribute(G, "name") + topo <- igraph::topological.sort(G.obs) + v.s <- igraph::get.vertex.attribute(G, "name") topo <- v.s[topo] if (length(setdiff(y, topo)) > 0) stop("Set 'y' contains variables not present in the graph.") if (length(setdiff(x, topo)) > 0) stop("Set 'x' contains variables not present in the graph.") if (length(intersect(x, y)) > 0) stop("Sets 'x' and 'y' are not disjoint. ") s <- v.s[sel] s <- topo[which(topo %in% s)] - G.causal <- induced.subgraph(G, v.s[!(v.s %in% s)]) - v <- get.vertex.attribute(G.causal, "name") + G.causal <- igraph::induced.subgraph(G, v.s[!(v.s %in% s)]) + v <- igraph::get.vertex.attribute(G.causal, "name") v <- topo[which(topo %in% v)] T.prime <- setdiff(v, x) - G.T.prime <- induced.subgraph(G.causal, T.prime) + G.T.prime <- igraph::induced.subgraph(G.causal, T.prime) G.T.prime.obs <- observed.graph(G.T.prime) D <- ancestors(y, G.T.prime.obs, topo) - G.D <- induced.subgraph(G.causal, D) + G.D <- igraph::induced.subgraph(G.causal, D) cc <- c.components(G.D, topo) cg <- length(cc) res <- probability() diff --git a/R/sc.components.R b/R/sc.components.R index 010720e..de9be5f 100644 --- a/R/sc.components.R +++ b/R/sc.components.R @@ -1,9 +1,9 @@ sc.components <- function(D, topo) { from <- NULL - A <- as.matrix(get.adjacency(D)) - v <- get.vertex.attribute(D, "name") - s <- v[which(vertex.attributes(D)$description == "S")] - e <- E(D) + A <- as.matrix(igraph::get.adjacency(D)) + v <- igraph::get.vertex.attribute(D, "name") + s <- v[which(igraph::vertex.attributes(D)$description == "S")] + e <- igraph::E(D) bidirected <- NULL selection <- e[from(s)] indices <- which(A >= 1 & t(A) >= 1, arr.ind = TRUE) @@ -12,10 +12,10 @@ sc.components <- function(D, topo) { e[v[x[1]] %->% v[x[2]]] })) } - D.bidirected <- subgraph.edges(D, union(bidirected, selection), delete.vertices = FALSE) - subgraphs <- decompose.graph(D.bidirected) + D.bidirected <- igraph::subgraph.edges(D, union(bidirected, selection), delete.vertices = FALSE) + subgraphs <- igraph::decompose.graph(D.bidirected) cc.s <- lapply(subgraphs, function(x) { - v.sub <- get.vertex.attribute(x, "name") + v.sub <- igraph::get.vertex.attribute(x, "name") return(topo[which(topo %in% v.sub)]) }) cc.s.rank <- order(sapply(cc.s, function(x) { diff --git a/R/simplify.R b/R/simplify.R index 39198b0..a967d64 100644 --- a/R/simplify.R +++ b/R/simplify.R @@ -1,4 +1,4 @@ -simplify <- function(P, topo, G.adj, G, G.obs) { +simplify <- function(P, topo, G.unobs, G, G.obs) { j <- 0 while (j < length(P$sumset)) { P.orig <- P @@ -14,7 +14,7 @@ simplify <- function(P, topo, G.adj, G, G.obs) { J <- character() D <- character() if (i > 1) { - irrel <- irrelevant(P$children[1:(i-1)], P$sumset[j], P$sumset, G.adj) + irrel <- irrelevant(P$children[1:(i-1)], P$sumset[j], P$sumset, G.unobs) irl.len <- length(irrel) if (irl.len > 0) { i <- i - irl.len @@ -26,7 +26,7 @@ simplify <- function(P, topo, G.adj, G, G.obs) { M <- topo[!(topo %in% vars)] O <- topo[(topo %in% vars)] while (k <= i) { - joint <- join(J, D, P$children[[k]]$var, P$children[[k]]$cond, P$sumset[j], M, O, G.adj, G, G.obs, topo) + joint <- join(J, D, P$children[[k]]$var, P$children[[k]]$cond, P$sumset[j], M, O, G.unobs, G, G.obs, topo) if (length(joint[[1]]) <= length(J)) { J <- character() D <- character() diff --git a/R/soid.R b/R/soid.R index bb4f812..c897c48 100644 --- a/R/soid.R +++ b/R/soid.R @@ -21,7 +21,7 @@ soid <- function(y, x, P, G, G.obs, I, J, v, topo, tree) { # line 2 if (length(setdiff(v, an)) != 0) { - G.an <- induced.subgraph(G, an) + G.an <- igraph::induced.subgraph(G, an) G.an.obs <- observed.graph(G.an) if (P$product | P$fraction) { P$sumset <- union(setdiff(v, an), P$sumset) %ts% topo @@ -51,7 +51,7 @@ soid <- function(y, x, P, G, G.obs, I, J, v, topo, tree) { P.new <- P P.new$do <- i$x v.l <- v.l[!(v.l %in% i$x)] - G.lz <- induced.subgraph(G.l, v.l) + G.lz <- igraph::induced.subgraph(G.l, v.l) G.lz.obs <- observed.graph(G.lz) nxt <- soid(y, setdiff(x, i$x) %ts% topo, P, G.lz, G.lz.obs, I, i$x, v.l, topo, list()) nxt$P <- activate.interventions(nxt$P, 0, i$x) @@ -72,7 +72,7 @@ soid <- function(y, x, P, G, G.obs, I, J, v, topo, tree) { } # line 3 - G.xbar <- subgraph.edges(G, E(G)[!(to(x) | (from(x) & (description == "U" & !is.na(description))))], delete.vertices = FALSE) + G.xbar <- igraph::subgraph.edges(G, igraph::E(G)[!(to(x) | (from(x) & (description == "U" & !is.na(description))))], delete.vertices = FALSE) an.xbar <- ancestors(y, observed.graph(G.xbar), topo) w <- setdiff(setdiff(v, x), an.xbar) w.len <- length(w) @@ -86,7 +86,7 @@ soid <- function(y, x, P, G, G.obs, I, J, v, topo, tree) { } # line 4 - G.remove.x <- induced.subgraph(G, v[!(v %in% x)]) + G.remove.x <- igraph::induced.subgraph(G, v[!(v %in% x)]) s <- c.components(G.remove.x, topo) if (length(s) > 1) { tree$call$line <- 4 @@ -156,7 +156,7 @@ soid <- function(y, x, P, G, G.obs, I, J, v, topo, tree) { tree$call$s.prime <- s s.len <- length(s) ind <- which(v %in% s) - G.s <- induced.subgraph(G, s) + G.s <- igraph::induced.subgraph(G, s) G.s.obs <- observed.graph(G.s) product.list <- vector(mode = "list", length = s.len) P.prod <- probability() diff --git a/R/surrogate.outcome.R b/R/surrogate.outcome.R index 71f345c..5e5a2cb 100644 --- a/R/surrogate.outcome.R +++ b/R/surrogate.outcome.R @@ -1,12 +1,12 @@ surrogate.outcome <- function(y, x, S, G, expr = TRUE, steps = FALSE, primes = FALSE, stop_on_nonid = TRUE) { d <- length(S) + 1 - v <- get.vertex.attribute(G, "name") + v <- igraph::get.vertex.attribute(G, "name") if (length(intersect(x, y)) > 0) stop("Sets 'x' and 'y' are not disjoint.") D <- vector(mode = "list", length = d) Z <- vector(mode = "list", length = d) D[[1]] <- G G.obs <- observed.graph(G) - topo.obs <- topological.sort(G.obs) + topo.obs <- igraph::topological.sort(G.obs) Z[[1]] <- character(0) for (i in 1:(d-1)) { Di <- G @@ -16,23 +16,23 @@ surrogate.outcome <- function(y, x, S, G, expr = TRUE, steps = FALSE, primes = F if ((tr.len <- length(target)) > 0) { for (j in 1:tr.len) { tr.node <- paste0(c("S_{",i,",",j,"}"), collapse = "") - Di <- Di + vertex(tr.node, description = "S") - Di <- Di + edge(tr.node, target[j]) + Di <- Di + igraph::vertex(tr.node, description = "S") + Di <- Di + igraph::edge(tr.node, target[j]) } } D[[i+1]] <- Di Z[[i+1]] <- Zi } - topo <- lapply(D, function(k) topological.sort(observed.graph(k))) - topo <- lapply(1:d, function(k) get.vertex.attribute(D[[k]], "name")[topo[[k]]]) + topo <- lapply(D, function(k) igraph::topological.sort(observed.graph(k))) + topo <- lapply(1:d, function(k) igraph::get.vertex.attribute(D[[k]], "name")[topo[[k]]]) D <- lapply(D, function(k) { - if (length(edge.attributes(k)) == 0) { - k <- set.edge.attribute(k, "description", 1:length(E(k)), NA) + if (length(igraph::edge.attributes(k)) == 0) { + k <- igraph::set.edge.attribute(k, "description", 1:length(igraph::E(k)), NA) } return(k) }) for (i in 1:d) { - if (!is.dag(observed.graph(D[[i]]))) { + if (!igraph::is.dag(observed.graph(D[[i]]))) { if (i > 1) stop("Selection diagram 'D[", i, "]' is not a DAG.") else stop("Causal diagram 'D[", i, "]' is not a DAG.") } diff --git a/R/tr.target.R b/R/tr.target.R index ffa12f2..f90f458 100644 --- a/R/tr.target.R +++ b/R/tr.target.R @@ -3,7 +3,7 @@ tr.target <- function(Z, W, G.obs, G, topo.obs) { from <- NULL description <- NULL target <- setdiff(descendants(Z, G.obs, topo.obs), W) - G.xbar <- subgraph.edges(G, E(G)[!(to(Z) | (from(Z) & (description == "U" & !is.na(description))))], delete.vertices = FALSE) # remove id nonid + G.xbar <- igraph::subgraph.edges(G, igraph::E(G)[!(to(Z) | (from(Z) & (description == "U" & !is.na(description))))], delete.vertices = FALSE) # remove id nonid G.xbar.obs <- observed.graph(G.xbar) nontarget <- setdiff(ancestors(W, G.xbar.obs, topo.obs), target) cc <- c.components(G, topo.obs) diff --git a/R/transport.R b/R/transport.R index 2850055..43debf3 100644 --- a/R/transport.R +++ b/R/transport.R @@ -1,11 +1,11 @@ transport <- function(y, x, z = NULL, D, expr = TRUE, simp = TRUE, steps = FALSE, primes = FALSE, stop_on_nonid = TRUE) { - if (length(edge.attributes(D)) == 0) { - D <- set.edge.attribute(D, "description", 1:length(E(D)), NA) + if (length(igraph::edge.attributes(D)) == 0) { + D <- igraph::set.edge.attribute(D, "description", 1:length(igraph::E(D)), NA) } - v <- get.vertex.attribute(D, "name") - s <- v[which(vertex.attributes(D)$description == "S")] + v <- igraph::get.vertex.attribute(D, "name") + s <- v[which(igraph::vertex.attributes(D)$description == "S")] if (is.null(z)) z <- setdiff(v, union(y, s)) - D.causal <- induced.subgraph(D, v[!(v %in% s)]) + D.causal <- igraph::induced.subgraph(D, v[!(v %in% s)]) res <- generalize(y = y, x = x, Z = list(character(0), z), D = list(D.causal, D), expr = FALSE, simp = simp, steps = TRUE, primes = primes, stop_on_nonid = stop_on_nonid) res.prob <- res$P attr(res.prob, "algorithm") <- "trz" diff --git a/R/trmz.R b/R/trmz.R index 0603d9d..0ea1d8c 100644 --- a/R/trmz.R +++ b/R/trmz.R @@ -3,11 +3,11 @@ trmz <- function(y, x, P, J, domain, w.index, D, Z, topo, tree) { from <- NULL description <- NULL d <- length(D) - v.s <- lapply(D, function(x) get.vertex.attribute(x, "name")) - s <- lapply(1:d, function(x) v.s[[x]][which(vertex.attributes(D[[x]])$description == "S")]) + v.s <- lapply(D, function(x) igraph::get.vertex.attribute(x, "name")) + s <- lapply(1:d, function(x) v.s[[x]][which(igraph::vertex.attributes(D[[x]])$description == "S")]) s <- lapply(1:d, function(x) topo[[x]][which(topo[[x]] %in% s[[x]])]) D.causal <- D[[1]] - v <- get.vertex.attribute(D.causal, "name") + v <- igraph::get.vertex.attribute(D.causal, "name") v <- topo[[1]][which(topo[[1]] %in% v)] D.obs <- observed.graph(D.causal) D.s.obs <- lapply(D, function(x) observed.graph(x)) @@ -34,7 +34,7 @@ trmz <- function(y, x, P, J, domain, w.index, D, Z, topo, tree) { # line 2 if (length(setdiff(v, anc)) != 0) { - anc.graph <- lapply(1:d, function(x) induced.subgraph(D[[x]], anc.s[[x]])) + anc.graph <- lapply(1:d, function(x) igraph::induced.subgraph(D[[x]], anc.s[[x]])) if (P$product | P$fraction | P$sum) { P$sumset <- union(setdiff(v, anc), P$sumset) P <- simplify.expression(P, NULL) @@ -50,7 +50,7 @@ trmz <- function(y, x, P, J, domain, w.index, D, Z, topo, tree) { } # line 3 - D.x.overbar <- subgraph.edges(D.causal ,E(D.causal)[!(to(x) | (from(x) & (description == "U" & !is.na(description))))], delete.vertices = FALSE) + D.x.overbar <- igraph::subgraph.edges(D.causal, igraph::E(D.causal)[!(to(x) | (from(x) & (description == "U" & !is.na(description))))], delete.vertices = FALSE) anc.xbar <- ancestors(y, observed.graph(D.x.overbar), topo[[1]]) w <- setdiff(setdiff(v, x), anc.xbar) if (length(w) != 0) { @@ -64,7 +64,7 @@ trmz <- function(y, x, P, J, domain, w.index, D, Z, topo, tree) { } # line 4 - D.remove.x <- induced.subgraph(D.causal, v[!(v %in% x)]) + D.remove.x <- igraph::induced.subgraph(D.causal, v[!(v %in% x)]) cc <- c.components(D.remove.x, topo[[1]]) cc.len <- length(cc) if (cc.len > 1) { @@ -144,7 +144,7 @@ trmz <- function(y, x, P, J, domain, w.index, D, Z, topo, tree) { tree$call$line <- 8 tree$call$cprime <- cc ind <- which(v %in% cc) - cc.graph <- lapply(1:d, function(x) induced.subgraph(D[[x]], cc.s[[x]])) + cc.graph <- lapply(1:d, function(x) igraph::induced.subgraph(D[[x]], cc.s[[x]])) kappa <- c() if (cc.len > 1) { for (i in 1:cc.len) { @@ -193,14 +193,13 @@ trmz <- function(y, x, P, J, domain, w.index, D, Z, topo, tree) { ind <- 0 W.new <- w.index + 1 for (i in 1:length(D)) { - A <- unobserved.graph(D[[i]]) - A <- subgraph.edges(A, E(A)[!to(x)], delete.vertices = FALSE) - A <- as.matrix(get.adjacency(A)) - if (wrap.dSep(A, s[[i]], y, x) & (length(intersect(Z[[i]], x)) != 0)) { + D.unobs <- unobserved.graph(D[[i]]) + D.unobs <- igraph::subgraph.edges(D.unobs, igraph::E(D.unobs)[!to(x)], delete.vertices = FALSE) + if (wrap.dSep(D.unobs, s[[i]], y, x) & (length(intersect(Z[[i]], x)) != 0)) { P.new <- P P.new$domain <- i xcapz <- intersect(Z[[i]], x) - D.remove.xcapz <- lapply(1:d, function(x) induced.subgraph(D[[x]], v.s[[x]][!(v.s[[x]] %in% xcapz)])) + D.remove.xcapz <- lapply(1:d, function(x) igraph::induced.subgraph(D[[x]], v.s[[x]][!(v.s[[x]] %in% xcapz)])) nxt <- trmz(y, setdiff(x, Z[[i]]), P, xcapz, i, W.new, D.remove.xcapz, Z, topo, list()) if (nxt$tree$call$id) { ind <- ind + 1 diff --git a/R/trmz_ex.R b/R/trmz_ex.R index 891000a..0a352fd 100644 --- a/R/trmz_ex.R +++ b/R/trmz_ex.R @@ -3,11 +3,11 @@ trmz_ex <- function(y, x, P, J, domain, w.index, D, Z, topo, tree, prioritize) { from <- NULL description <- NULL d <- length(D) - v.s <- lapply(D, function(x) get.vertex.attribute(x, "name")) - s <- lapply(1:d, function(x) v.s[[x]][which(vertex.attributes(D[[x]])$description == "S")]) + v.s <- lapply(D, function(x) igraph::get.vertex.attribute(x, "name")) + s <- lapply(1:d, function(x) v.s[[x]][which(igraph::vertex.attributes(D[[x]])$description == "S")]) s <- lapply(1:d, function(x) topo[[x]][which(topo[[x]] %in% s[[x]])]) D.causal <- D[[1]] - v <- get.vertex.attribute(D.causal, "name") + v <- igraph::get.vertex.attribute(D.causal, "name") v <- topo[[1]][which(topo[[1]] %in% v)] D.obs <- observed.graph(D.causal) D.s.obs <- lapply(D, function(x) observed.graph(x)) @@ -34,7 +34,7 @@ trmz_ex <- function(y, x, P, J, domain, w.index, D, Z, topo, tree, prioritize) { # line 2 if (length(setdiff(v, anc)) != 0) { - anc.graph <- lapply(1:d, function(x) induced.subgraph(D[[x]], anc.s[[x]])) + anc.graph <- lapply(1:d, function(x) igraph::induced.subgraph(D[[x]], anc.s[[x]])) if (P$product | P$fraction | P$sum) { P$sumset <- union(setdiff(v, anc), P$sumset) P <- simplify.expression(P, NULL) @@ -50,7 +50,7 @@ trmz_ex <- function(y, x, P, J, domain, w.index, D, Z, topo, tree, prioritize) { } # line 3 - D.x.overbar <- subgraph.edges(D.causal ,E(D.causal)[!(to(x) | (from(x) & (description == "U" & !is.na(description))))], delete.vertices = FALSE) + D.x.overbar <- igraph::subgraph.edges(D.causal, igraph::E(D.causal)[!(to(x) | (from(x) & (description == "U" & !is.na(description))))], delete.vertices = FALSE) anc.xbar <- ancestors(y, observed.graph(D.x.overbar), topo[[1]]) w <- setdiff(setdiff(v, x), anc.xbar) if (length(w) != 0) { @@ -64,7 +64,7 @@ trmz_ex <- function(y, x, P, J, domain, w.index, D, Z, topo, tree, prioritize) { } # line 4 - D.remove.x <- induced.subgraph(D.causal, v[!(v %in% x)]) + D.remove.x <- igraph::induced.subgraph(D.causal, v[!(v %in% x)]) cc <- c.components(D.remove.x, topo[[1]]) cc.len <- length(cc) if (cc.len > 1) { @@ -99,14 +99,13 @@ trmz_ex <- function(y, x, P, J, domain, w.index, D, Z, topo, tree, prioritize) { ind <- 0 W.new <- w.index + 1 for (i in 1:length(D)) { - A <- unobserved.graph(D[[i]]) - A <- subgraph.edges(A, E(A)[!to(x)], delete.vertices = FALSE) - A <- as.matrix(get.adjacency(A)) - if (wrap.dSep(A, s[[i]], y, x) & (length(intersect(Z[[i]], x)) != 0)) { + D.unobs <- unobserved.graph(D[[i]]) + D.unobs <- igraph::subgraph.edges(D.unobs, igraph::E(D.unobs)[!to(x)], delete.vertices = FALSE) + if (wrap.dSep(D.unobs, s[[i]], y, x) & (length(intersect(Z[[i]], x)) != 0)) { P.new <- P P.new$domain <- i xcapz <- intersect(Z[[i]], x) - D.remove.xcapz <- lapply(1:d, function(x) induced.subgraph(D[[x]], v.s[[x]][!(v.s[[x]] %in% xcapz)])) + D.remove.xcapz <- lapply(1:d, function(x) igraph::induced.subgraph(D[[x]], v.s[[x]][!(v.s[[x]] %in% xcapz)])) nxt <- trmz_ex(y, setdiff(x, Z[[i]]), P, xcapz, i, W.new, D.remove.xcapz, Z, topo, list(), prioritize) if (nxt$tree$call$id) { ind <- ind + 1 @@ -188,7 +187,7 @@ trmz_ex <- function(y, x, P, J, domain, w.index, D, Z, topo, tree, prioritize) { tree$call$line <- 8 tree$call$cprime <- cc ind <- which(v %in% cc) - cc.graph <- lapply(1:d, function(x) induced.subgraph(D[[x]], cc.s[[x]])) + cc.graph <- lapply(1:d, function(x) igraph::induced.subgraph(D[[x]], cc.s[[x]])) kappa <- c() if (cc.len > 1) { for (i in 1:cc.len) { @@ -245,14 +244,13 @@ trmz_ex <- function(y, x, P, J, domain, w.index, D, Z, topo, tree, prioritize) { ind <- 0 W.new <- w.index + 1 for (i in 1:length(D)) { - A <- unobserved.graph(D[[i]]) - A <- subgraph.edges(A, E(A)[!to(x)], delete.vertices = FALSE) - A <- as.matrix(get.adjacency(A)) - if (wrap.dSep(A, s[[i]], y, x) & (length(intersect(Z[[i]], x)) != 0)) { + D.unobs <- unobserved.graph(D[[i]]) + D.unobs <- igraph::subgraph.edges(D.unobs, igraph::E(D.unobs)[!to(x)], delete.vertices = FALSE) + if (wrap.dSep(D.unobs, s[[i]], y, x) & (length(intersect(Z[[i]], x)) != 0)) { P.new <- P P.new$domain <- i xcapz <- intersect(Z[[i]], x) - D.remove.xcapz <- lapply(1:d, function(x) induced.subgraph(D[[x]], v.s[[x]][!(v.s[[x]] %in% xcapz)])) + D.remove.xcapz <- lapply(1:d, function(x) igraph::induced.subgraph(D[[x]], v.s[[x]][!(v.s[[x]] %in% xcapz)])) nxt <- trmz_ex(y, setdiff(x, Z[[i]]), P, xcapz, i, W.new, D.remove.xcapz, Z, topo, list(), prioritize) if (nxt$tree$call$id) { ind <- ind + 1 diff --git a/R/trso.R b/R/trso.R index cbcc9e6..920b3bd 100644 --- a/R/trso.R +++ b/R/trso.R @@ -3,11 +3,11 @@ trso <- function(y, x, P, J, domain, D, Z, topo, tree) { from <- NULL description <- NULL d <- length(D) - v.s <- lapply(D, function(x) get.vertex.attribute(x, "name")) - s <- lapply(1:d, function(x) v.s[[x]][which(vertex.attributes(D[[x]])$description == "S")]) + v.s <- lapply(D, function(x) igraph::get.vertex.attribute(x, "name")) + s <- lapply(1:d, function(x) v.s[[x]][which(igraph::vertex.attributes(D[[x]])$description == "S")]) s <- lapply(1:d, function(x) topo[[x]][which(topo[[x]] %in% s[[x]])]) D.causal <- D[[1]] - v <- get.vertex.attribute(D.causal, "name") + v <- igraph::get.vertex.attribute(D.causal, "name") v <- topo[[1]][which(topo[[1]] %in% v)] D.obs <- observed.graph(D.causal) D.s.obs <- lapply(D, function(x) observed.graph(x)) @@ -34,7 +34,7 @@ trso <- function(y, x, P, J, domain, D, Z, topo, tree) { # line 2 if (length(setdiff(v, anc)) != 0) { - anc.graph <- lapply(1:d, function(x) induced.subgraph(D[[x]], anc.s[[x]])) + anc.graph <- lapply(1:d, function(x) igraph::induced.subgraph(D[[x]], anc.s[[x]])) if (P$product | P$fraction | P$sum) { P$sumset <- union(setdiff(v, anc), P$sumset) P <- simplify.expression(P, NULL) @@ -50,7 +50,7 @@ trso <- function(y, x, P, J, domain, D, Z, topo, tree) { } # line 3 - D.x.overbar <- subgraph.edges(D.causal ,E(D.causal)[!(to(x) | (from(x) & (description == "U" & !is.na(description))))], delete.vertices = FALSE) + D.x.overbar <- igraph::subgraph.edges(D.causal, igraph::E(D.causal)[!(to(x) | (from(x) & (description == "U" & !is.na(description))))], delete.vertices = FALSE) anc.xbar <- ancestors(y, observed.graph(D.x.overbar), topo[[1]]) w <- setdiff(setdiff(v, x), anc.xbar) if (length(w) != 0) { @@ -64,7 +64,7 @@ trso <- function(y, x, P, J, domain, D, Z, topo, tree) { } # line 4 - D.remove.x <- induced.subgraph(D.causal, v[!(v %in% x)]) + D.remove.x <- igraph::induced.subgraph(D.causal, v[!(v %in% x)]) cc <- c.components(D.remove.x, topo[[1]]) cc.len <- length(cc) if (cc.len > 1) { @@ -92,10 +92,9 @@ trso <- function(y, x, P, J, domain, D, Z, topo, tree) { tree$call$line <- 6 ind <- 0 for (i in 2:length(D)) { - A <- unobserved.graph(D[[i]]) - A <- subgraph.edges(A, E(A)[!to(x)], delete.vertices = FALSE) - A <- as.matrix(get.adjacency(A)) - if (wrap.dSep(A, s[[i]], y, x) & (length(intersect(Z[[i]], x)) > 0)) { + D.unobs <- unobserved.graph(D[[i]]) + D.unobs <- igraph::subgraph.edges(D.unobs, igraph::E(D.unobs)[!to(x)], delete.vertices = FALSE) + if (wrap.dSep(D.unobs, s[[i]], y, x) & (length(intersect(Z[[i]], x)) > 0)) { P.new <- P P.new$domain <- i xcapz <- intersect(Z[[i]], x) diff --git a/R/unobserved.graph.R b/R/unobserved.graph.R index 00a1e95..b9f73ac 100644 --- a/R/unobserved.graph.R +++ b/R/unobserved.graph.R @@ -1,16 +1,16 @@ unobserved.graph <- function(G) { - unobs.edges <- which(edge.attributes(G)$description == "U") + unobs.edges <- which(igraph::edge.attributes(G)$description == "U") u <- length(unobs.edges) if (u > 0) { - e <- get.edges(G, unobs.edges) + e <- igraph::get.edges(G, unobs.edges) e <- e[e[ ,1] > e[ ,2], , drop = FALSE] e.len <- nrow(e) new.nodes <- paste0("u_{", 1:e.len, "}") - G <- G + vertices(new.nodes, description = rep("U", e.len)) - v <- get.vertex.attribute(G, "name") - G <- G + edges(c(rbind(new.nodes, v[e[ ,1]]), rbind(new.nodes, v[e[ ,2]]))) - obs.edges <- setdiff(E(G), E(G)[unobs.edges]) - G.unobs <- subgraph.edges(G, E(G)[obs.edges], delete.vertices = FALSE) + G <- G + igraph::vertices(new.nodes, description = rep("U", e.len)) + v <- igraph::get.vertex.attribute(G, "name") + G <- G + igraph::edges(c(rbind(new.nodes, v[e[ ,1]]), rbind(new.nodes, v[e[ ,2]]))) + obs.edges <- setdiff(igraph::E(G), igraph::E(G)[unobs.edges]) + G.unobs <- igraph::subgraph.edges(G, igraph::E(G)[obs.edges], delete.vertices = FALSE) return(G.unobs) } return(G) diff --git a/R/verma.constraints.R b/R/verma.constraints.R index 5d33d74..ab986a1 100644 --- a/R/verma.constraints.R +++ b/R/verma.constraints.R @@ -1,34 +1,22 @@ verma.constraints <- function(G) { - if (length(edge.attributes(G)) == 0) { - G <- set.edge.attribute(G, "description", 1:length(E(G)), NA) + if (length(igraph::edge.attributes(G)) == 0) { + G <- igraph::set.edge.attribute(G, "description", 1:length(igraph::E(G)), NA) } G.obs <- observed.graph(G) G.unobs <- unobserved.graph(G) - G.adj <- as.matrix(get.adjacency(G.unobs)) - topo <- topological.sort(G.obs) - topo.u <- topological.sort(G.unobs) - v <- get.vertex.attribute(G, "name")[topo] - v.unobs <- get.vertex.attribute(G.unobs, "name")[topo.u] + G.adj <- as.matrix(igraph::get.adjacency(G.unobs)) + topo <- igraph::topological.sort(G.obs) + topo.u <- igraph::topological.sort(G.unobs) + v <- igraph::get.vertex.attribute(G, "name")[topo] + v.unobs <- igraph::get.vertex.attribute(G.unobs, "name")[topo.u] constraints <- list() s <- NULL for (i in 1:length(v)) { vi <- v[1:i] - G.vi <- induced.subgraph(G, vi) + G.vi <- igraph::induced.subgraph(G, vi) cc <- c.components(G.vi, v) if (length(cc) > 1) s <- Find(function(x) v[i] %in% x, cc) else s <- cc[[1]] - # pa.s <- parents(s, G.obs, v) - # pred <- setdiff(vi, pa.s) - # if (length(pred) > 0) { - # eff <- setdiff(pa.s, v[i]) - # if (dSep(G.adj, v[i], pred, eff)) { - # constraints <- c(constraints, list(c( - # "X" = v[i], - # "Y" = paste0(pred, collapse = ","), - # "Z" = paste0(eff, collapse = ",") - # ))) - # } - # } constraints <- c(constraints, q.constraints(s, v[i], G, G.obs, G.unobs, v, v.unobs, list())) } return(constraints) diff --git a/R/wrap.dSep.R b/R/wrap.dSep.R index 38f8182..dadb95d 100644 --- a/R/wrap.dSep.R +++ b/R/wrap.dSep.R @@ -1,7 +1,7 @@ -wrap.dSep <- function(amat, first, second, cond) { - if (identical(first, second)) return(FALSE) - if (length(first) == 0 | length(second) == 0) { +wrap.dSep <- function(G, x, y, z) { + if (identical(x, y)) return(FALSE) + if (length(x) == 0 || length(y) == 0) { return(TRUE) } - else return(dSep(amat, first, second, cond)) + else return(dSep(G, x, y, z)) } \ No newline at end of file diff --git a/R/zzaux.effect.R b/R/zzaux.effect.R index a21f96c..71e31bd 100644 --- a/R/zzaux.effect.R +++ b/R/zzaux.effect.R @@ -1,6 +1,6 @@ aux.effect <- function(y, x, z, G, expr = TRUE, simp = TRUE, steps = FALSE, primes = FALSE, stop_on_nonid = TRUE) { - if (length(edge.attributes(G)) == 0) { - G <- set.edge.attribute(G, "description", 1:length(E(G)), NA) + if (length(igraph::edge.attributes(G)) == 0) { + G <- igraph::set.edge.attribute(G, "description", 1:length(igraph::E(G)), NA) } res <- generalize(y = y, x = x, Z = list(z), D = list(G), expr = FALSE, simp = simp, steps = TRUE, primes = primes, stop_on_nonid = stop_on_nonid) res.prob <- res$P