Commit c76d7bcd authored by Gilles Kratzer's avatar Gilles Kratzer
Browse files

typos

parent c49ed4f3
Pipeline #1640 passed with stage
in 2 seconds
......@@ -7,17 +7,14 @@ Authors@R: c(person("Gilles", "Kratzer", role = c("aut", "cre"),
person("Reinhard", "Furrer", role = c("ctb"),
email = "reinhard.furrer@math.uzh.ch",
comment = c(ORCID = "0000-0002-6319-2332")))
Author: Gilles Kratzer [aut, cre] (<https://orcid.org/0000-0002-5929-8935>),
Reinhard Furrer [ctb] (<https://orcid.org/0000-0002-6319-2332>)
Maintainer: Gilles Kratzer <gilles.kratzer@math.uzh.ch>
Description: Flexible implementation of a structural MCMC sampler for Directed Acyclic Graphs (DAGs). It supports the new edge reversal move from Grzegorczyk and Husmeier (2008) <doi.10.1007/s10994-008-5057-7> and the Markov blanket resampling from Su and Borsuk (2016) <http://jmlr.org/papers/v17/su16a.html>. It supports three priors: a prior controlling for structure complexity from Koivisto and Sood (2004) <http://dl.acm.org/citation.cfm?id=1005332.1005352>, an uninformative prior and a user defined prior. The three main problems that can be addressed by this R package are selecting the most probable structure based on a cache of pre-computed scores, controlling for overfitting and sampling the landscape of high scoring structures. It allows to quantify the marginal impact of relationships of interest by marginalising out over structures or nuisance dependencies. Structural MCMC seems a very elegant and natural way to estimate the true marginal impact, so one can determine if it's magnitude is big enough to consider as a worthwhile intervention.
Depends: R (>= 3.5.0)
Depends: R (>= 3.0.0)
License: GPL-3
Encoding: UTF-8
LazyData: true
Imports: gRbase, abn, coda, ggplot2, ggpubr, cowplot
Imports: gRbase, abn, coda, ggplot2, cowplot, ggpubr
Suggests: bnlearn, knitr, rmarkdown, ggdag, testthat
RoxygenNote: 6.1.1
VignetteBuilder: knitr
URL: https://www.math.uzh.ch/pages/mcmcabn/
BugReports: https://git.math.uzh.ch/gkratz/mcmcabn/issues
export(mcmcabn,
query,
plot.mcmcabn
plot.mcmcabn,
summary.mcmcabn,
print.mcmcabn
)
importFrom(gRbase, topoSortMAT)
importFrom(stats,rbinom, acf)
importFrom(ggplot2,ggplot)
import(ggplot2)
importFrom(coda,effectiveSize)
import(abn)
importFrom(ggplot2,aes,geom_line,geom_hline,geom_text,geom_point,labs,scale_color_manual,geom_density,scale_fill_manual,coord_flip,scale_colour_manual, aes_string)
importFrom(ggpubr, theme_pubr)
importFrom(cowplot, ggdraw, axis_canvas, insert_yaxis_grob)
importFrom("graphics", "plot")
importFrom("stats", "as.formula", "dist", "quantile", "sd")
importFrom(graphics, plot)
importFrom(stats, as.formula, dist, quantile, sd)
#methods for class mcmcabn
S3method(plot, mcmcabn)
......
MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
## scores
score.G <- 0
score.G.prime <- 0
rejection <- 1
A <- NULL
# store number of edges
n.edges <- sum(dag.tmp)
# Randomly select a node Xi in current graph G, withprobability N^−1
selected.node <- sample(x = 1:n.var, size = 1)
# store current parent set
parent.set <- dag.tmp[selected.node, ]
# remove parent set + co parent
dag.G.0 <- dag.tmp
tmp <- which(x = dag.G.0[, selected.node] == 1, arr.ind = TRUE)
if (!is.integer0(tmp)) {
dag.G.0[tmp[1], -selected.node] <- 0
}
dag.G.0[selected.node, ] <- 0
# descendant of node i
descendents.i <- descendents(nodes = selected.node, dag = dag.G.0)
# from score take out parent + descendant
sc.tmp <- sc[score.cache$children == selected.node, ]
if (sum(parent.set) == 1) {
......@@ -41,7 +41,7 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
if (sum(parent.set) > 1) {
sc.tmp <- sc.tmp[rowSums(sc.tmp[, c(as.logical(parent.set), FALSE)]) == sum(parent.set), ]
}
# remove descendents of i
if (!is.null(descendents.i)) {
if (length(descendents.i) == 1) {
......@@ -50,9 +50,9 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
} else {
sc.tmp <- sc.tmp[sc.tmp[, descendents.i] == 0, ]
}
} else {
if (is.null(dim(sc.tmp))) {
sc.tmp <- sc.tmp[sum(sc.tmp[descendents.i] == 0) == length(descendents.i)]
} else {
......@@ -61,9 +61,9 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
}
}
# sample new parent set
if (nrow(sc.tmp) > 0 || is.vector(sc.tmp))
if (nrow(sc.tmp) > 0 || is.vector(sc.tmp))
{
if (is.vector(sc.tmp)) {
new.parent.i <- sc.tmp
new.parent.i <- new.parent.i[-length(new.parent.i)]
......@@ -71,35 +71,35 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
z.i.G.0 <- sc.tmp[length(sc.tmp)]
} else {
# sample a new parent set
new.parent.i <- sc.tmp[sample(x = 1:nrow(sc.tmp), size = 1, prob = range01(sc.tmp[, ncol(sc.tmp)] -
new.parent.i <- sc.tmp[sample(x = 1:nrow(sc.tmp), size = 1, prob = range01(sc.tmp[, ncol(sc.tmp)] -
sum(sc.tmp[, ncol(sc.tmp)]))), ]
new.parent.i <- new.parent.i[-length(new.parent.i)]
z.i.G.0 <- (sum(sc.tmp[, ncol(sc.tmp)]))
}
# dag G1
dag.G.1 <- dag.G.0
dag.G.1[selected.node, ] <- new.parent.i
## new co parent set
if (sum(dag.G.1[, selected.node]) != 0) {
score.G <- numeric(sum(dag.G.1[, selected.node]))
# randomization of children
order.child <- sample(1:sum(dag.G.1[, selected.node]))
for (j in order.child) {
# descendant of node j
descendents.j <- descendents(nodes = j, dag = dag.G.1)
if (!is.character(descendents.j)) {
# from score take out parent + descendant
sc.tmp <- sc[score.cache$children == j, ]
sc.tmp <- sc.tmp[sc.tmp[, selected.node] == 1, ]
# print(class(sc.tmp)) remove descendents of i
if (is.null(ncol(sc.tmp))) {
if (!is.null(descendents.j)) {
if (length(descendents.j) == 1) {
......@@ -117,54 +117,54 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
}
}
}
if (nrow(sc.tmp) > 0 || is.vector(sc.tmp)) {
if (!identical(unname(sc.tmp[length(sc.tmp)]), numeric(0))) {
if (is.vector(sc.tmp)) {
new.parent.j <- sc.tmp
new.parent.j <- new.parent.j[-length(new.parent.j)]
# store partition function print(sc.tmp[length(sc.tmp)])
score.G[j] <- sc.tmp[length(sc.tmp)]
} else {
# sample a new parent set
new.parent.j <- sc.tmp[sample(x = 1:nrow(sc.tmp), size = 1, prob = range01(sc.tmp[, ncol(sc.tmp)] -
new.parent.j <- sc.tmp[sample(x = 1:nrow(sc.tmp), size = 1, prob = range01(sc.tmp[, ncol(sc.tmp)] -
sum(sc.tmp[, ncol(sc.tmp)]))), ]
new.parent.j <- new.parent.j[-length(new.parent.j)]
score.G[j] <- sum(sc.tmp[, ncol(sc.tmp)])
}
dag.G.1[j, ] <- new.parent.j
}
}
}
}
}
dag.MBR <- dag.G.1
############################ complementary inverse move
# parent set
parent.set.G.1 <- dag.G.1[selected.node, ]
# delete co-parent
tmp <- which(x = dag.G.1[, selected.node] == 1, arr.ind = TRUE)
if (!is.integer0(tmp)) {
dag.G.1[tmp[1], -selected.node] <- 0
}
# delete parents
dag.G.1[selected.node, ] <- 0
# from score take out parent + descendant
sc.tmp <- sc[score.cache$children == selected.node, ]
if (sum(parent.set.G.1) != 0) {
......@@ -173,7 +173,7 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
} else {
sc.tmp <- sc.tmp[rowSums(sc.tmp[, c(as.logical(parent.set.G.1), FALSE)]) == sum(parent.set.G.1), ]
}
# remove descendents of i
if (!is.null(descendents.i)) {
if (length(descendents.i) == 1) {
......@@ -182,9 +182,9 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
} else {
sc.tmp <- sc.tmp[sc.tmp[, descendents.i] == 0, ]
}
} else {
if (is.null(dim(sc.tmp))) {
sc.tmp <- sc.tmp[sum(sc.tmp[descendents.i] == 0) == length(descendents.i)]
} else {
......@@ -192,31 +192,31 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
}
}
}
# store score
if (is.null(ncol(sc.tmp))) {
z.star.i.G.prime.0 <- sc.tmp[length(sc.tmp)]
} else {
z.star.i.G.prime.0 <- sum(sc.tmp[, ncol(sc.tmp)])
}
# equivalent to G1
dag.G.1[selected.node, ] <- parent.set
if (sum(dag.G.1[, selected.node]) != 0) {
score.G.prime <- numeric(sum(dag.G.1[, selected.node]))
for (j in 1:sum(dag.G.1[, selected.node])) {
# descendant of node j
descendents.j <- descendents(nodes = j, dag = dag.G.1)
if (!is.character(descendents.j)) {
# from score take out parent + descendant
sc.tmp <- sc[score.cache$children == j, ]
sc.tmp <- sc.tmp[sc.tmp[, selected.node] == 1, ]
# remove descendents of j
if (is.null(ncol(sc.tmp))) {
if (!is.null(descendents.j)) {
......@@ -226,7 +226,7 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
sc.tmp <- sc.tmp[sum(sc.tmp[descendents.j] == 0) == length(descendents.j)]
}
}
} else {
if (!is.null(descendents.j)) {
if (length(descendents.j) == 1) {
......@@ -236,34 +236,34 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
}
}
}
dag.G.1[j, ] <- dag.tmp[j, ]
if (is.null(ncol(sc.tmp))) {
score.G.prime[j] <- sum(sc.tmp[length(sc.tmp)])
} else {
score.G.prime[j] <- sum(sc.tmp[, ncol(sc.tmp)])
}
}
}
}
############################## Acceptance probability
score.A <- exp(z.i.G.0 - z.star.i.G.prime.0 + sum(score.G) - sum(score.G.prime))
A <- min(1, score.A)
if (rbinom(n = 1, size = 1, prob = A) == 1) {
rejection <- 1
dag.tmp <- dag.MBR
score.MBR <- 0
for (a in 1:n.var) {
sc.tmp <- sc[score.cache$children == a, ]
score.MBR <- sum(min(sc.tmp[(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag.tmp[a,
score.MBR <- sum(min(sc.tmp[(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag.tmp[a,
])))), n.var + 1]), score.MBR)
}
score <- score.MBR
......@@ -276,6 +276,3 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
############################## Return
return(list(dag.tmp = dag.tmp, score = score, alpha = A, rejection = rejection))
} #EOF
TRUE
TRUE
TRUE
REV <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
rejection <- 1
A <- 0
# stor number of edges
n.edges <- sum(dag.tmp)
# randomly select one
if (sum(dag.tmp) != 0)
if (sum(dag.tmp) != 0)
{
selected.edge <- which(x = dag.tmp == 1, arr.ind = TRUE)[sample(x = 1:n.edges, size = 1), ]
# store current parent set
parent.set <- dag.tmp[selected.edge[1], ]
# remove parent set
dag.M.dot <- dag.tmp
dag.M.dot[selected.edge[1], ] <- 0
# store descendents j row i col
descendents.M.dot.i <- descendents(nodes = unname(selected.edge[2]), dag = dag.M.dot)
descendents.M.dot.j <- descendents(nodes = unname(selected.edge[1]), dag = dag.M.dot)
# mark all parent sets of i score that do not include j
sc.tmp <- sc[score.cache$children == selected.edge[2], ]
sc.tmp <- sc.tmp[sc.tmp[, selected.edge[1]] == 1, ]
# remove descendents of i
if (is.null(ncol(sc.tmp))) {
if (!is.null(descendents.M.dot.i)) {
......@@ -50,13 +50,13 @@ REV <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
}
}
}
# sample new parent set
if (nrow(sc.tmp) > 0 || is.vector(sc.tmp)) {
if (is.vector(sc.tmp)) {
new.parent.i <- sc.tmp
new.parent.i <- new.parent.i[-length(new.parent.i)]
......@@ -64,37 +64,37 @@ REV <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
z.star.x.i.M.dot <- sc.tmp[length(sc.tmp)]
} else {
# !!!exp missing
new.parent.i <- sc.tmp[sample(x = 1:nrow(sc.tmp), size = 1, prob = range01(sc.tmp[, ncol(sc.tmp)] -
new.parent.i <- sc.tmp[sample(x = 1:nrow(sc.tmp), size = 1, prob = range01(sc.tmp[, ncol(sc.tmp)] -
sum(sc.tmp[, ncol(sc.tmp)]))), ]
new.parent.i <- new.parent.i[-length(new.parent.i)]
# store partition function
z.star.x.i.M.dot <- (sum(sc.tmp[, ncol(sc.tmp)]))
}
# M direct sum
dag.M.cross <- dag.M.dot
dag.M.cross[selected.edge[2], ] <- new.parent.i
# descendant of i
descendents.M.cross.i <- descendents(nodes = unname(selected.edge[2]), dag = dag.M.cross)
descendents.M.cross.j <- descendents(nodes = unname(selected.edge[1]), dag = dag.M.cross)
# score node j not descendant of M
sc.tmp <- sc[score.cache$children == selected.edge[1], ]
# remove descendents of i
if (!is.null(descendents.M.cross.j)) {
if (length(descendents.M.cross.j) == 1) {
sc.tmp <- sc.tmp[sc.tmp[, descendents.M.cross.j] == 0, ]
} else {
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.cross.j] == 0) == length(descendents.M.cross.j),
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.cross.j] == 0) == length(descendents.M.cross.j),
]
}
}
# sample new parent set
if (nrow(sc.tmp) > 0 || is.vector(sc.tmp)) {
if (is.vector(sc.tmp)) {
new.parent.j <- sc.tmp
# print(new.parent.i)
......@@ -103,22 +103,22 @@ REV <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
z.x.j.M.cross <- sc.tmp[length(sc.tmp)]
} else {
## !!! exp missing
new.parent.j <- sc.tmp[sample(x = 1:nrow(sc.tmp), size = 1, prob = range01(sc.tmp[, ncol(sc.tmp)] -
new.parent.j <- sc.tmp[sample(x = 1:nrow(sc.tmp), size = 1, prob = range01(sc.tmp[, ncol(sc.tmp)] -
sum(sc.tmp[, ncol(sc.tmp)]))), ]
new.parent.j <- new.parent.j[-length(new.parent.j)]
# store partition function
z.x.j.M.cross <- (sum(sc.tmp[, ncol(sc.tmp)]))
}
# M tilde
dag.M.tilde <- dag.M.cross
dag.M.tilde[selected.edge[1], ] <- new.parent.j
n.edges.tilde <- sum(dag.M.tilde)
############################## computing acceptance proba
# score node j that do not include i
sc.tmp <- sc[score.cache$children == selected.edge[1], ]
sc.tmp <- sc.tmp[sc.tmp[, selected.edge[2]] == 1, ]
......@@ -126,7 +126,7 @@ REV <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
if (is.vector(sc.tmp)) {
if (!is.null(descendents.M.dot.j)) {
if (length(descendents.M.dot.j) == 1) {
sc.tmp <- sc.tmp[sc.tmp[descendents.M.dot.j] == 0] #
} else {
# print(class(sc.tmp)) print(descendents.M.dot.j)
......@@ -136,59 +136,59 @@ REV <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
} else {
if (!is.null(descendents.M.dot.j)) {
if (length(descendents.M.dot.j) == 1) {
sc.tmp <- sc.tmp[sc.tmp[, descendents.M.dot.j] == 0, ]
} else {
# print(class(sc.tmp)) print(descendents.M.dot.j)
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.dot.j] == 0) == length(descendents.M.dot.j),
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.dot.j] == 0) == length(descendents.M.dot.j),
]
}
}
}
if (is.vector(sc.tmp)) {
z.star.x.j.M.dot <- sc.tmp[length(sc.tmp)]
} else {
z.star.x.j.M.dot <- (sum(sc.tmp[, ncol(sc.tmp)]))
}
dag.M.tilde.cross <- dag.M.dot
dag.M.tilde.cross[selected.edge[1], ] <- parent.set
# descendant
descendents.M.tilde.cross.i <- descendents(nodes = unname(selected.edge[2]), dag = dag.M.tilde.cross)
sc.tmp <- sc[score.cache$children == selected.edge[2], ]
if (!is.null(descendents.M.tilde.cross.i)) {
if (length(descendents.M.tilde.cross.i) == 1) {
sc.tmp <- sc.tmp[sc.tmp[, descendents.M.tilde.cross.i] == 0, ]
} else {
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.tilde.cross.i] == 0) == length(descendents.M.tilde.cross.i),
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.tilde.cross.i] == 0) == length(descendents.M.tilde.cross.i),
]
}
}
if (is.vector(sc.tmp)) {
z.x.i.M.tilde.cross <- sc.tmp[length(sc.tmp)]
} else {
z.x.i.M.tilde.cross <- (sum(sc.tmp[, ncol(sc.tmp)]))
}
############################## Acceptance probability
score.A <- n.edges/n.edges.tilde * exp(z.star.x.i.M.dot - z.star.x.j.M.dot + z.x.j.M.cross - z.x.j.M.cross)
A <- min(1, score.A)
if (rbinom(n = 1, size = 1, prob = A) == 1) {
rejection <- 0
dag.tmp <- dag.M.tilde
score.M.tilde <- 0
for (a in 1:n.var) {
sc.tmp <- sc[score.cache$children == a, ]
score.M.tilde <- sum(min(sc.tmp[(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag.tmp[a,
score.M.tilde <- sum(min(sc.tmp[(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag.tmp[a,
])))), n.var + 1]), score.M.tilde)
}
score <- score.M.tilde
......@@ -196,12 +196,9 @@ REV <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
}
}
} #EOFif
############################## Return
return(list(dag.tmp = dag.tmp, score = score, alpha = A, rejection = rejection))
} #EOF
TRUE
TRUE
TRUE
mc3 <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, prior.choice, prior.lambda, prior.dag, verbose) {
## construction of neighbours list