Commit 7faaf5ff authored by Gilles Kratzer's avatar Gilles Kratzer
Browse files

plot, vignette, s3 ...

parent f048fafb
Package: mcmcabn
Title: Flexible implementation of an MCMC sampler for DAGs.
Title: Flexible implementation of a structural MCMC sampler for DAGs
Version: 0.0.0.9000
Authors@R: c(person("Gilles", "Kratzer", role = c("aut", "cre"),
email = "gilles.kratzer@math.uzh.ch",
......@@ -10,10 +10,10 @@ Authors@R: c(person("Gilles", "Kratzer", role = c("aut", "cre"),
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 an MCMC sampler for 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) <doi:>
Description: Flexible implementation of a structural MCMC sampler for 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) <doi:>. The two main problems that can be addressed by this package is selecting the most probable structure based on a cache of pre-computed scores and the sampling of the landscape of high scoring networks which can be used as an alternative to parametric bootstrapping.
Depends: R (>= 3.5.1)
License: GPL-3
Encoding: UTF-8
LazyData: true
Imports: gRbase, abn
RoxygenNote: 6.1.0
Imports: gRbase, abn, bnlearn, coda, ggplot2, ggpubr, cowplot
RoxygenNote: 6.1.1
export(mcmc,
REV,
mc3)
query
)
importFrom(gRbase, topoSortMAT)
importFrom(stats,rbinom)
\ No newline at end of file
importFrom(stats,rbinom)
importFrom(ggplot2,ggplot)
import(abn)
import(ggpubr)
import(cowplot)
import(bnlearn)
#methods for class mcmcabn
S3method(plot, mcmcabn)
\ No newline at end of file
......@@ -3,6 +3,7 @@ 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
......@@ -29,7 +30,7 @@ MBR<-function(n.var, dag.tmp, max.parents, sc,score.cache, score, verbose){
#descendant of node i
descendants.i <- descendants(nodes = selected.node,dag = dag.G.0)
descendents.i <- descendents(nodes = selected.node,dag = dag.G.0)
#from score take out parent + descendant
......@@ -41,21 +42,21 @@ MBR<-function(n.var, dag.tmp, max.parents, sc,score.cache, score, verbose){
sc.tmp <- sc.tmp[rowSums(sc.tmp[,c(as.logical(parent.set),FALSE)])==sum(parent.set),]
}
# remove descendants of i
if(!is.null(descendants.i)){
if(length(descendants.i)==1){
# remove descendents of i
if(!is.null(descendents.i)){
if(length(descendents.i)==1){
if(is.null(dim(sc.tmp))){
sc.tmp <- sc.tmp[sc.tmp[descendants.i]==0]
sc.tmp <- sc.tmp[sc.tmp[descendents.i]==0]
}else{
sc.tmp <- sc.tmp[sc.tmp[,descendants.i]==0,]
sc.tmp <- sc.tmp[sc.tmp[,descendents.i]==0,]
}
}else{
if(is.null(dim(sc.tmp))){
sc.tmp <- sc.tmp[sum(sc.tmp[descendants.i]==0)==length(descendants.i)]
sc.tmp <- sc.tmp[sum(sc.tmp[descendents.i]==0)==length(descendents.i)]
}else{
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendants.i]==0)==length(descendants.i),]}
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendents.i]==0)==length(descendents.i),]}
}
}
#sample new parent set
......@@ -87,20 +88,20 @@ if(sum(dag.G.1[,selected.node])!=0){
for (j in order.child) {
#descendant of node j
descendants.j <- descendants(nodes = j,dag = dag.G.1)
if(!is.character(descendants.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 descendants of i
if(!is.null(descendants.j)){
if(length(descendants.j)==1){
sc.tmp <- sc.tmp[sc.tmp[,descendants.j]==0,]
# remove descendents of i
if(!is.null(descendents.j)){
if(length(descendents.j)==1){
sc.tmp <- sc.tmp[sc.tmp[,descendents.j]==0,]
}else{
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendants.j]==0)==length(descendants.j),]}
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendents.j]==0)==length(descendents.j),]}
}
if(nrow(sc.tmp)>0 || is.vector(sc.tmp)){
......@@ -154,21 +155,21 @@ for (j in order.child) {
sc.tmp <- sc.tmp[rowSums(sc.tmp[,c(as.logical(parent.set.G.1),FALSE)])==sum(parent.set.G.1),]
}
# remove descendants of i
if(!is.null(descendants.i)){
if(length(descendants.i)==1){
# remove descendents of i
if(!is.null(descendents.i)){
if(length(descendents.i)==1){
if(is.null(dim(sc.tmp))){
sc.tmp <- sc.tmp[sc.tmp[descendants.i]==0]
sc.tmp <- sc.tmp[sc.tmp[descendents.i]==0]
}else{
sc.tmp <- sc.tmp[sc.tmp[,descendants.i]==0,]
sc.tmp <- sc.tmp[sc.tmp[,descendents.i]==0,]
}
}else{
if(is.null(dim(sc.tmp))){
sc.tmp <- sc.tmp[sum(sc.tmp[descendants.i]==0)==length(descendants.i)]
sc.tmp <- sc.tmp[sum(sc.tmp[descendents.i]==0)==length(descendents.i)]
}else{
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendants.i]==0)==length(descendants.i),]}
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendents.i]==0)==length(descendents.i),]}
}
}
......@@ -190,18 +191,18 @@ for (j in order.child) {
for (j in 1:sum(dag.G.1[,selected.node])) {
#descendant of node j
descendants.j <- descendants(nodes = j,dag = dag.G.1)
if(!is.character(descendants.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 descendants of j
if(!is.null(descendants.j)){
if(length(descendants.j)==1){
sc.tmp <- sc.tmp[sc.tmp[,descendants.j]==0,]
# remove descendents of j
if(!is.null(descendents.j)){
if(length(descendents.j)==1){
sc.tmp <- sc.tmp[sc.tmp[,descendents.j]==0,]
}else{
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendants.j]==0)==length(descendants.j),]}
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendents.j]==0)==length(descendents.j),]}
}
dag.G.1[j,] <- dag.tmp[j,]
......@@ -218,8 +219,9 @@ if(!is.character(descendants.j)){
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){
if(rbinom(n = 1,size = 1,prob = A)==1){
rejection <- 1
dag.tmp <- dag.MBR
score.MBR <- 0
for(a in 1:n.var){
......@@ -233,5 +235,5 @@ if(is.null(A)){A<-0}
##############################
## Return
##############################
return(list(dag.tmp = dag.tmp, score = score, alpha=A))
return(list(dag.tmp = dag.tmp, score = score, alpha=A, rejection = rejection))
}#EOF
REV<-function(n.var, dag.tmp, max.parents, sc,score.cache, score, verbose){
#stor number of edges
rejection <- 1
#stor number of edges
n.edges <- sum(dag.tmp)
......@@ -16,24 +18,24 @@ parent.set <- dag.tmp[selected.edge[1],]
dag.M.dot <- dag.tmp
dag.M.dot[selected.edge[1],] <- 0
#store descendant
#store descendents
# j row
# i col
descendants.M.dot.i <- descendants(nodes = unname(selected.edge[2]),dag = dag.M.dot)
descendants.M.dot.j <- descendants(nodes = unname(selected.edge[1]),dag = dag.M.dot)
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 descendants of i
if(!is.null(descendants.M.dot.i)){
if(length(descendants.M.dot.i)==1){
sc.tmp <- sc.tmp[sc.tmp[,descendants.M.dot.i]==0,]
# remove descendents of i
if(!is.null(descendents.M.dot.i)){
if(length(descendents.M.dot.i)==1){
sc.tmp <- sc.tmp[sc.tmp[,descendents.M.dot.i]==0,]
}else{
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendants.M.dot.i]==0)==length(descendants.M.dot.i),]}
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendents.M.dot.i]==0)==length(descendents.M.dot.i),]}
}
#sample new parent set
......@@ -58,18 +60,18 @@ dag.M.cross <- dag.M.dot
dag.M.cross[selected.edge[2],] <- new.parent.i
#descendant of i
descendants.M.cross.i <- descendants(nodes = unname(selected.edge[2]),dag = dag.M.cross)
descendants.M.cross.j <- descendants(nodes = unname(selected.edge[1]),dag = dag.M.cross)
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 descendants of i
if(!is.null(descendants.M.cross.j)){
if(length(descendants.M.cross.j)==1){
sc.tmp <- sc.tmp[sc.tmp[,descendants.M.cross.j]==0,]
# 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[,descendants.M.cross.j]==0)==length(descendants.M.cross.j),]}
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendents.M.cross.j]==0)==length(descendents.M.cross.j),]}
}
#sample new parent set
......@@ -102,12 +104,12 @@ n.edges.tilde <- sum(dag.M.tilde)
# 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,]
# remove descendants of i
if(!is.null(descendants.M.dot.j)){
if(length(descendants.M.dot.j)==1){
sc.tmp <- sc.tmp[sc.tmp[,descendants.M.dot.j]==0,]
# remove descendents of i
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{
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendants.M.dot.j]==0)==length(descendants.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)){
......@@ -121,15 +123,15 @@ dag.M.tilde.cross <- dag.M.dot
dag.M.tilde.cross[selected.edge[1],] <- parent.set
#descendant
descendants.M.tilde.cross.i <- descendants(nodes = unname(selected.edge[2]),dag = dag.M.tilde.cross)
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(descendants.M.tilde.cross.i)){
if(length(descendants.M.tilde.cross.i)==1){
sc.tmp <- sc.tmp[sc.tmp[,descendants.M.tilde.cross.i]==0,]
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[,descendants.M.tilde.cross.i]==0)==length(descendants.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)){
......@@ -141,11 +143,12 @@ if(is.vector(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)
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){
......@@ -156,9 +159,11 @@ if(rbinom(n = 1,size = 1,prob = A)==1){
score <- score.M.tilde
}
}}#EOFif
if(is.null(A)){A<-0}
if(is.null(A)){A <- 0}
##############################
## Return
##############################
return(list(dag.tmp = dag.tmp, score = score, alpha=A))
return(list(dag.tmp = dag.tmp, score = score, alpha=A, rejection = rejection))
}#EOF
################################################################################
## ancestor.R ---
## Author : Gilles Kratzer
## Document created: 11/01/2018
## -
################################################################################
##-------------------------------------------------------------------------
## ancestors finding procedures
##
##-------------------------------------------------------------------------
mc3<-function(n.var, dag.tmp, max.parents, sc ,score.cache, score, prior, verbose){
mc3<-function(n.var, dag.tmp, max.parents, sc ,score.cache, score, prior.choice,prior.lambda, verbose){
##construction of neighbours list
neighbours.list <- NULL
alpha <- 0
rejection <- 1
##choose mc move
if(sum(dag.tmp)>0){
......@@ -116,11 +117,16 @@ mc3<-function(n.var, dag.tmp, max.parents, sc ,score.cache, score, prior, verbos
for(a in 1:n.var){
#Koivisto priors
if(prior==2){
if(prior.choice==2){
prior.G <- sum(1/choose(n = (n.var-1),k = sum(dag.tmp[a,])),prior.G)
prior.Gprime <- sum(1/choose(n = (n.var-1),k = sum(dag.gprime[a,])),prior.Gprime)
}
if(prior.choice==3){
prior.G <- sum(prior.lambda * exp(-prior.lambda*abs(dag.tmp[a,]-prior.lambda[a,])),prior.G)
prior.Gprime <-sum(1/choose(n = (n.var-1),k = sum(dag.gprime[a,])),prior.Gprime)
}
sc.tmp <- sc[score.cache$children==a,]
score.G <- sum(min(sc.tmp[(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag.tmp[a,])))),n.var+1]),score.G)
score.Gprime <- sum(min(sc.tmp[(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag.gprime[a,])))),n.var+1]),score.Gprime)
......@@ -134,8 +140,11 @@ mc3<-function(n.var, dag.tmp, max.parents, sc ,score.cache, score, prior, verbos
if(!is.null(dag.gprime) & rbinom(n = 1,size = 1,prob = alpha)==1){
dag.tmp <- dag.gprime
score <- score.Gprime
rejection <- 0
}}
return(list(dag.tmp=dag.tmp, score=score, alpha=alpha))
}
return(list(dag.tmp = dag.tmp,
score = score,
alpha = alpha,
rejection = rejection))
}#EOF
################################################################################
## heuristic.search.R ---
## mcmc.R ---
## Author : Gilles Kratzer
## Document created: 22/10/2018
## -
......@@ -18,11 +18,11 @@ mcmc <- function(score.cache = NULL,
seed = 42,
verbose = FALSE,
start.dag = NULL,
dag.retained = NULL,
dag.banned = NULL,
prior.dag = NULL,
prior.lambda = NULL,
prob.rev = 0.05,
prob.mbr = 0.05,
prior = 2){
prior.choice = 2){
####################################################
##Tests
......@@ -34,7 +34,14 @@ mcmc <- function(score.cache = NULL,
if(score %!in% c("bic", "aic", "mdl", "mlik"))stop("A method should be provided and be one of: bic, aic, mdl or mlik.")
if(max(rowSums(score.cache$node.defn))<max.parents)stop("Check max.parents. It should be the same as the one used in abn::buildscorecache() R function")
if(length(mcmc.scheme)!=3)stop("An MCMC scheme have to be provided. It should be such that c(returned,thinned,burned).")
if(max.parents<1 || max.parents>length(data.dists))stop("max.parents makes no sense.")
if(is.numeric(prob.rev) && prob.rev>1 && prob.rev<0)stop("prob.rev should be a probability.")
if(is.numeric(prob.mbr) && prob.mbr>1 && prob.mbr<0)stop("prob.mbr should be a probability.")
if(is.matrix(start.dag) && dim(start.dag)!=c(length(data.dists),length(data.dists)))stop("start.dag should be a squared matrix with dimension equal to the number of variables.")
if(length(data.dists)!=max(score.cache$children))stop("data.dists should be a named list of all variables used to build the cache of precomputed score.")
if(is.matrix(prior.dag)){prior.choice <- 3}
if(is.matrix(prior.dag) && is.null(prior.lambda)){prior.lambda <- 1}
if(is.matrix(prior.dag) && dim(prior.dag)==c(length(data.dists),length(data.dists)))stop("prior.dag should be a squared matrix with dimension equal to the number of variables.")
##end of tests
n.var <- length(data.dists)
......@@ -45,6 +52,8 @@ mcmc <- function(score.cache = NULL,
out.dags <- array(data = NA,dim = c(n.var,n.var,mcmc.scheme[1]+1))
out.scores <- rep(NA,length = mcmc.scheme[1]+1)
out.alpha <- rep(NA,length = mcmc.scheme[1]+1)
out.method <- rep(NA,length = mcmc.scheme[1]+1)
out.rejection <- rep(NA,length = mcmc.scheme[1]+1)
##seeding
set.seed(seed = seed)
......@@ -99,6 +108,10 @@ mcmc <- function(score.cache = NULL,
out.alpha[1] <- alpha <- 0
out.method[1] <- "MC3"
out.rejection[1] <- 0
##start mcmc search
if(verbose) cat("Start MCMC burn in \n");
......@@ -108,10 +121,10 @@ mcmc <- function(score.cache = NULL,
##method choice:
method.choice <- sample(x = c("mc3","REV","MBR"),size = 1,prob = c(prob.mc3,prob.rev,prob.mbr))
method.choice <- sample(x = c("MC3","REV","MBR"),size = 1,prob = c(prob.mc3,prob.rev,prob.mbr))
switch (method.choice, "mc3" = {
out <- mc3(n.var,(dag.tmp),max.parents,sc,score.cache,score,prior,verbose)
switch (method.choice, "MC3" = {
out <- mc3(n.var,(dag.tmp),max.parents,sc,score.cache,score,prior.choice,prior.lambda,verbose)
dag.tmp <- out$dag.tmp
score <- out$score
},
......@@ -143,13 +156,14 @@ switch (method.choice, "mc3" = {
j <- 0
while(j<=mcmc.scheme[2]){
method.choice <- sample(x = c("mc3","REV","MBR"),size = 1,prob = c(prob.mc3,prob.rev,prob.mbr))
method.choice <- sample(x = c("MC3","REV","MBR"),size = 1,prob = c(prob.mc3,prob.rev,prob.mbr))
switch (method.choice, "mc3" = {
out <- mc3(n.var,(dag.tmp),max.parents,sc,score.cache,score,prior,verbose)
switch (method.choice, "MC3" = {
out <- mc3(n.var,(dag.tmp),max.parents,sc,score.cache,score,prior.choice,prior.lambda,verbose)
dag.tmp <- out$dag.tmp
score <- out$score
alpha <- out$alpha
rejection <- out$rejection
},
"REV" = {
if(verbose){print("REV move")}
......@@ -157,6 +171,7 @@ switch (method.choice, "mc3" = {
dag.tmp <- out$dag.tmp
score <- out$score
alpha <- out$alpha
rejection <- out$rejection
},
"MBR" = {
if(verbose){print("MBR move")}
......@@ -164,6 +179,7 @@ switch (method.choice, "mc3" = {
dag.tmp <- out$dag.tmp
score <- out$score
alpha <- out$alpha
rejection <- out$rejection
})
j<-j+1
......@@ -175,9 +191,22 @@ switch (method.choice, "mc3" = {
out.alpha[mcmc.search+1] <- alpha
}#EOF mcmc search
out.method[mcmc.search+1] <- method.choice
out.rejection[mcmc.search+1] <- rejection
return(list(dags = out.dags,
}#EOF mcmc search
out<-list(dags = out.dags,
scores = out.scores,
alpha = out.alpha))
alpha = out.alpha,
method = out.method,
rejection = out.rejection,
iterations = mcmc.scheme[1]*mcmc.scheme[2],
thinning = mcmc.scheme[2]
)
class(out) <- "mcmcabn"
return(out)
}#eof
###############################################################################
## abn-internal.R ---
## mcmcabn-utiles.R ---
## Author : Gilles Kratzer
## Document created : 13/02/2018
## Last modification :
......@@ -127,14 +127,6 @@ formula.mcmcabn<-function(f, name){
return(out)
}
##-------------------------------------------------------------------------
##Ancestor function
##-------------------------------------------------------------------------
is.integer0 <- function(x)
{
is.integer(x) && length(x) == 0L
}
##-------------------------------------------------------------------------
##Ancestor function
......@@ -153,17 +145,17 @@ is.integer0 <- function(x)
# }
##-------------------------------------------------------------------------
##Descendant function
##descendent function
##-------------------------------------------------------------------------
descendants<-function(nodes, dag){
descendents<-function(nodes, dag){
if(!identical(topoSortMAT(amat = dag,index = FALSE),character(0))){
if(!is.integer0(which(x = dag[,nodes]==1,arr.ind = TRUE)) && diag(dag)==0){
if(length(nodes)==1){
tmp <- unname(which(x = dag[,nodes]==1,arr.ind = TRUE))}else{
tmp <- t(unname(which(x = dag[,nodes]==1,arr.ind = TRUE)))[1,]
}
return(unique(unname(c(tmp,descendants(nodes = (tmp),dag = dag)))))
return(unique(unname(c(tmp,descendents(nodes = (tmp),dag = dag)))))
}else{
return(NULL)
}
......@@ -171,6 +163,23 @@ descendants<-function(nodes, dag){
return("Not a DAG")
}
##-------------------------------------------------------------------------
##range function
##-------------------------------------------------------------------------
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
##-------------------------------------------------------------------------
##not in function
##-------------------------------------------------------------------------
'%!in%' <- function(x,y)!('%in%'(x,y))
##-------------------------------------------------------------------------
##is.integer0
##-------------------------------------------------------------------------
is.integer0 <- function(x)
{
is.integer(x) && length(x) == 0L
}
###############################################################################
## plot.mcmcabn.R ---
## Author : Gilles Kratzer
## Last modified : 19/02/2019
## :
###############################################################################
plot.mcmcabn <- function(x,
max.score = FALSE,
...){
dta <- data.frame(x[-1])
dta$X <- 1:length(x$scores)
max.score. <- max(x$scores)
if(!max.score){
original_plot <- ggplot(data = dta, aes(x = X,y = scores)) +
geom_line(alpha = 0.25) +
geom_hline(yintercept = max.score., linetype="dashed", color = "red")+
geom_point(aes(color=as.factor(method)))+
labs(x = "DAG index", y = "DAG scores",colour="Methods") + theme_pubr() +
scale_color_manual(values = c("#F2C500","#37454B", "#56B4E9"))
y_density <- axis_canvas(original_plot, axis = "y", coord_flip = TRUE) +
geom_density(data = dta, aes(x=scores,fill = factor(method)), color = NA, alpha = 0.5) +
scale_fill_manual(values = c( "#F2C500","#37454B", "#56B4E9")) +
coord_flip()
# create the combined plot
return(ggdraw(insert_yaxis_grob(plot = original_plot, grob = y_density, position = "right")))
invisible()
}else{
dta$cummax[1] <- dta$scores[1]