Commit 5106081c authored by Gilles Kratzer's avatar Gilles Kratzer
Browse files

CRAN checks: 1w 4n

parent 07948be3
Pipeline #1551 passed with stage
in 2 seconds
Package: mcmcabn
Title: Flexible implementation of a structural MCMC sampler for DAGs
Version: 0.0.0.9000
Version: 0.1
Authors@R: c(person("Gilles", "Kratzer", role = c("aut", "cre"),
email = "gilles.kratzer@math.uzh.ch",
comment = c(ORCID = "0000-0002-5929-8935")),
......@@ -10,10 +10,12 @@ 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 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)
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) and the Markov blanket resampling from Su and Borsuk (2016). It supports three priors: the one from Koivisto and Sood (2004), an uninformative prior and a user defined prior. The three main problems that can be addressed by this R package is selecting the most probable structure based on a cache of pre-computed scores, controlling for overfitting and sampling the landscape of high scoring structures which can be used as an alternative to parametric bootstrapping and allow the end user to make model averaging with DAGs.
Depends: R (>= 3.5.0)
License: GPL-3
Encoding: UTF-8
LazyData: true
Imports: gRbase, abn, bnlearn, coda, ggplot2, ggpubr, cowplot
Imports: gRbase, abn, coda, ggplot2, ggpubr, cowplot
Suggests: bnlearn, knitr, rmarkdown
RoxygenNote: 6.1.1
VignetteBuilder: knitr
......@@ -8,12 +8,13 @@ plot.mcmcabn
importFrom(ggplot2,ggplot)
importFrom(coda,effectiveSize)
import(abn)
import(ggpubr)
import(cowplot)
import(bnlearn)
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)
importFrom(ggpubr, theme_pubr)
importFrom(cowplot, ggdraw, axis_canvas, insert_yaxis_grob)
importFrom("graphics", "plot")
importFrom("stats", "as.formula", "dist", "quantile", "sd")
#methods for class mcmcabn
S3method(plot, mcmcabn)
S3method(summary, mcmcabn)
S3method(plot, mcmcabn)
\ No newline at end of file
S3method(print, mcmcabn)
\ No newline at end of file
## mcmcabn 0.0.900:
## mcmcabn 0.1:
* first release on CRAN
* mcmcabn(), query(), summary(), print() and plot() functions
......@@ -122,6 +122,7 @@ mc3<-function(n.var, dag.tmp, max.parents, sc ,score.cache, score, prior.choice,
prior.Gprime <- sum(1/choose(n = (n.var-1),k = sum(dag.gprime[a,])),prior.Gprime)
}
#user prior
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)
......
......@@ -145,6 +145,7 @@ switch (method.choice, "MC3" = {
}#EOF
out.scores[1] <- score
out.dags[,,1] <- dag.tmp
###EOF: Burnin phase!
......
......@@ -7,7 +7,7 @@
plot.mcmcabn <- function(x,
max.score = FALSE,
...){
...){
dta <- data.frame(x[-1])
dta$X <- 1:length(x$scores)
max.score. <- max(x$scores)
......
......@@ -13,17 +13,42 @@
query <- function(mcmcabn = NULL,
formula = NULL){
if(!(!is.matrix(formula) || !is.null(formula) || !class(formula)=="formula")){stop("Formula statment should be given using either a matrix or a formula. Alternativelly, the argument could be null to return all individual arc support")}
if(is.matrix(formula)){
n <- length(mcmcabn$scores)
n.var <- length(mcmcabn$dist)
m.array <- array(data = 0,dim = c(n.var,n.var,n))
for(i in 1:n.var){
for(j in 1:n.var){
if(formula[i,j]==1){m.array[i,j,]<-1}
if(formula[i,j]==-1){m.array[i,j,]<- -1}
}
}
prob.1 <- which(x = m.array==1,arr.ind = TRUE)
prob.2 <- which(x = m.array==-1,arr.ind = TRUE)
res.1 <- apply(matrix(data = mcmcabn$dags[prob.1],ncol = sum(m.1),byrow = TRUE), 1, prod)
res.2 <- apply(matrix(data = mcmcabn$dags[prob.2],ncol = sum(m.2),byrow = TRUE), 1, prod)
out <- sum(res.1*res.2)/n
return(out)
}
if(is.null(formula)){
out <- apply(mcmcabn$dags,1:2, mean)
colnames(out)<-rownames(out)<-names(mcmcabn$dist)
return(out)
}else{
}
if(class(formula)=="formula"){
f<-as.character(formula)
##start as a formula
if(!grepl('~',f[1],fixed = TRUE)){stop("DAG specifications should start with a ~")}
if(!grepl('~',f[1],fixed = TRUE)){stop("Formula specifications should start with a ~")}
if(grepl('-',f[2],fixed = TRUE)){
......@@ -77,6 +102,5 @@ query <- function(mcmcabn = NULL,
return(out)
}
}
}}
}#eof
---
output: github_document
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, echo = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "README-"
)
```
[![](https://www.r-pkg.org/badges/version-ago/mcmcabn)](http://cran.rstudio.com/web/packages/mcmcabn/index.html)
[![Downloads](http://cranlogs.r-pkg.org/badges/grand-total/mcmcabn)](http://cran.rstudio.com/web/packages/mcmcabn/index.html)
[![Downloads](http://cranlogs.r-pkg.org/badges/mcmcabn)](http://cran.rstudio.com/web/packages/mcmcabn/index.html)
[![License: GPL v3](https://img.shields.io/badge/License-GPL%20v3-blue.svg)](http://www.gnu.org/licenses/gpl-3.0)
___
## mcmcabn: Flexible implementation of a structural MCMC sampler for DAGs
Flexible implementation of a structural MCMC sampler for Directed Acyclic Graphs (DAGs). It supports the new edge reversal move from Grzegorczyk and Husmeier (2008) and the Markov blanket resampling from Su and Borsuk (2016). It supports three priors: the one from Koivisto and Sood (2004), an uninformative prior and a user defined prior. The three main problems that can be addressed by this R package is selecting the most probable structure based on a cache of pre-computed scores, controlling for overfitting and sampling the landscape of high scoring structures which can be used as an alternative to parametric bootstrapping and allow the end user to make model averaging with DAGs.
___
## What's New
- 23/04/2019 - mcmcabn is available on CRAN (v 0.1)
- 18/02/2019 - new pre-print [Is a single unique Bayesian network enough to accurately represent your data?](https://arxiv.org/pdf/1902.06641.pdf) on arXiv
____
**`mcmcabn` is developed and maintained by [Gilles Kratzer](https://gilleskratzer.netlify.com/) and [Prof. Dr. Reinhard Furrer](https://user.math.uzh.ch/furrer/) from
[Applied Statistics Group](https://www.math.uzh.ch/as/index.php?id=as) from the University of Zurich.**
___
# mcmcabn (PUBLIC)
!!! UNSTABLE VERSION !!!
mcmcabn is a one man show (me!) and made of more than 10000 lines of code which are not bug free! So use it with caution and awareness.
## Installation
install.packages("https://git.math.uzh.ch/gkratz/mcmcabn/raw/master/mcmcabn_0.1.tar.gz", repo=NULL, type="source")
## 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) and the Markov blanket resampling from Su and Borsuk (2016). It supports three priors: the one from Koivisto and Sood (2004), an uninformative prior and a user defined prior. The three main problems that can be addressed by this R package is selecting the most probable structure based on a cache of pre-computed scores, controlling for overfitting and sampling the landscape of high scoring structures which can be used as an alternative to parametric bootstrapping and allow the end user to make model averaging with DAGs.
## Future implementations (ordered by urgency)
No preview for this file type
......@@ -87,7 +87,7 @@
</a>
</li>
<li>
<a href="https://git.math.uzh.ch/gkratz/varrank">
<a href="https://git.math.uzh.ch/gkratz/mcmcabn">
<span class="fa fa-gitlab fa-lg"></span>
gitlab
......
......@@ -57,7 +57,7 @@
</a>
</li>
<li>
<a href="https://git.math.uzh.ch/gkratz/varrank">
<a href="https://git.math.uzh.ch/gkratz/mcmcabn">
<span class="fa fa-gitlab fa-lg"></span>
gitlab
......@@ -89,7 +89,13 @@
<p><strong>What is the package design for?</strong></p>
<p>The three main problems addressed by this R package are: 1. selecting the most probable structure based on a cache of pre-computed scores, 2. controlling for overfitting, and 3. sampling the landscape of high scoring networks which can be used as an alternative to parametric bootstrapping. The latter could be very useful in an applied perspective to avoid reducing the richeness of Bayesian network modelling to only <strong>one</strong> structure!</p>
<p>The three main problems addressed by this R package are:</p>
<ol style="list-style-type: decimal">
<li>selecting the most probable structure based on a cache of pre-computed scores</li>
<li>controlling for overfitting</li>
<li>sampling the landscape of high scoring networks which can be used as an alternative to parametric bootstrapping.</li>
</ol>
<p>The latter could be very useful in an applied perspective to avoid reducing the richeness of Bayesian network modelling to only <strong>one</strong> structure!</p>
<p><strong>How does it work?</strong></p>
<p>The <em>mcmcabn</em> R package takes a cache of pre-computed score as input (possibly computed using <strong>buildscorecache()</strong> function from <em>abn</em> R package). Then <em>mcmcabn</em> generates Monte Carlo Markoc Chain (MCMC) samples from the posterior distribution of the most supported Directed Acyclic Graphs (DAGs). To do so user should choose a structural prior and a learning scheme. <em>mcmcabn</em> is equipped with multiple structural priors: the so-called Koivisto prior, an uninformative prior, the Ellis prior, a modular flat prior and a possibly user define prior and three algorithms for inferring the structure of Bayesian networks: the classical Monte Carlo Markov Chain Model Choice (MC^3), the new edge reversal (REV) from Grzegorczyk and Husmeier (2008) and the Markov blanket resampling (MBR) from Su and Borsuk (2016).</p>
<p><strong>What are the R package functionalities?</strong></p>
......@@ -151,12 +157,12 @@ bsc.compute &lt;-<span class="st"> </span><span class="kw"><a href="https://www.
<span class="dt">score =</span> <span class="st">"mlik"</span>,
<span class="dt">data.dists =</span> dist,
<span class="dt">max.parents =</span> <span class="dv">2</span>,
<span class="dt">mcmc.scheme =</span> <span class="kw"><a href="https://www.rdocumentation.org/packages/base/topics/c">c</a></span>(<span class="dv">1000</span>,<span class="dv">9</span>,<span class="dv">500</span>),
<span class="dt">mcmc.scheme =</span> <span class="kw"><a href="https://www.rdocumentation.org/packages/base/topics/c">c</a></span>(<span class="dv">10000</span>,<span class="dv">9</span>,<span class="dv">1000</span>),
<span class="dt">seed =</span> <span class="dv">42</span>,
<span class="dt">verbose =</span> <span class="ot">FALSE</span>,
<span class="dt">start.dag =</span> <span class="st">"random"</span>,
<span class="dt">prob.rev =</span> <span class="fl">0.05</span>,
<span class="dt">prob.mbr =</span> <span class="fl">0.05</span>,
<span class="dt">prob.rev =</span> <span class="fl">0.025</span>,
<span class="dt">prob.mbr =</span> <span class="fl">0.025</span>,
<span class="dt">prior.choice =</span> <span class="dv">2</span>)</code></pre></div>
<p>The function <em>mcmcabn()</em> returns a list with multiple entries: “dags”, “scores”, “alpha” and “method”.</p>
<ul>
......@@ -166,7 +172,10 @@ bsc.compute &lt;-<span class="st"> </span><span class="kw"><a href="https://www.
<li>method is the algorithm chosen for the index MCMC step.</li>
</ul>
<p>The object <code>mcmc.out</code> is of class <code>mcmcabn</code>. We can check that we reach the maximum possible posterior DAG score even with a small number of MCMC steps.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co">#maximum score get from the MCMC sampler</span>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">##to speed up building of the vignette, we store the output
<span class="kw"><a href="https://www.rdocumentation.org/packages/base/topics/load">load</a></span>(<span class="dt">file =</span> <span class="st">"../data/mcmc_run_asia.rda"</span>)
<span class="co">#maximum score get from the MCMC sampler</span>
<span class="kw"><a href="https://www.rdocumentation.org/packages/base/topics/Extremes">max</a></span>(mcmc.out<span class="op">$</span>scores)
<span class="co">#&gt; [1] -11151</span>
......@@ -176,7 +185,7 @@ dag &lt;-<span class="st"> </span><span class="kw"><a href="https://www.rdocumen
<span class="co">#&gt; Total sets g(S) to be evaluated over: 256</span>
<span class="kw"><a href="https://www.rdocumentation.org/packages/abn/topics/fitabn">fitabn</a></span>(<span class="dt">dag.m =</span> dag,<span class="dt">data.df =</span> asia,<span class="dt">data.dists =</span> dist)<span class="op">$</span>mlik
<span class="co">#&gt; [1] -11151</span></code></pre></div>
<p>One can plot the output using the <em>plot()</em>. On the graph below, one can see the trace plot of the posterior structure score. The dashed red line is the maximum reached score (as expected = -11151). The colored dot on the trace plot indicate when different methods have been used. The densities on the left represent the relative occurence frequencies of the methods.</p>
<p>One can plot the output using the <em>plot()</em>. On the graph below, one can see the trace plot of the posterior structure score. The dashed red line is the maximum reached score (as expected = -11151). The colored dot on the trace plot indicate when different methods have been used. The densities on the left represent the relative occuring frequencies of the methods.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw"><a href="https://www.rdocumentation.org/packages/graphics/topics/plot">plot</a></span>(mcmc.out)</code></pre></div>
<p><img src="mcmcabn_files/figure-html/unnamed-chunk-7-1.png" width="576" style="display: block; margin: auto;"></p>
<p>One can also plot the cumulative best score over the MCMC steps. On this graph one can see the cumulated maximum posterior structure score (as expected = -11151). The dashed red line is the maximum reached score. The colored dot on the trace plot indicate when different methods have been used.</p>
......@@ -185,30 +194,61 @@ dag &lt;-<span class="st"> </span><span class="kw"><a href="https://www.rdocumen
<p>On can also plot the summary of the MCMC run:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw"><a href="https://www.rdocumentation.org/packages/base/topics/summary">summary</a></span>(mcmc.out)
<span class="co">#&gt; MCMC summary:</span>
<span class="co">#&gt; Number of Burn in steps: 500</span>
<span class="co">#&gt; Number of MCMC steps: 10000</span>
<span class="co">#&gt; Number of Burn in steps: 1000</span>
<span class="co">#&gt; Number of MCMC steps: 1e+05</span>
<span class="co">#&gt; Thinning: 9</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; Maximum score: -11151</span>
<span class="co">#&gt; Empirical mean: -11683</span>
<span class="co">#&gt; Empirical standard deviation: 524</span>
<span class="co">#&gt; Empirical mean: -11675</span>
<span class="co">#&gt; Empirical standard deviation: 531</span>
<span class="co">#&gt; Quantiles of the posterior network score:</span>
<span class="co">#&gt; 0.025 0.25 0.5 0.75 0.975</span>
<span class="co">#&gt; BN score -12827 -12187 -11499 -11271 -11159</span>
<span class="co">#&gt; BN score -13023 -12127 -11466 -11263 -11156</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; Acceptance rate: 0.303</span>
<span class="co">#&gt; Sample size adjusted for autocorrelation: 97</span></code></pre></div>
<p>The chosen number of MCMC steps is certainly not sufficient to draw any valid conclusion form the posterior sample. But one can still query it.</p>
<p>If one want to know <em>what is the probability of LungCancer node being </em></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">query</span>(<span class="dt">mcmcabn =</span> mcmc.out,<span class="dt">formula =</span> <span class="op">~</span>LungCancer<span class="op">|</span>Smoking,<span class="dt">dist =</span> dist)
<span class="co">#&gt; [1] 0.409</span></code></pre></div>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">query</span>(<span class="dt">mcmcabn =</span> mcmc.out,<span class="dt">formula =</span> <span class="op">~</span>LungCancer<span class="op">|</span>Smoking<span class="op">+</span>Smoking<span class="op">|</span>LungCancer,<span class="dt">dist =</span> dist)
<span class="co">#&gt; [1] 0.686</span></code></pre></div>
<span class="co">#&gt; Acceptance rate: 0.25</span>
<span class="co">#&gt; Sample size adjusted for autocorrelation: 439</span></code></pre></div>
<p>The chosen number of MCMC steps is certainly not sufficient to draw any valid conclusion form the posterior sample. But one can still query it. Below is displayed the individual frequencies of the arcs.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">query</span>(<span class="dt">mcmcabn =</span> mcmc.out)
<span class="co">#&gt; Asia Smoking Tuberculosis LungCancer Bronchitis Either</span>
<span class="co">#&gt; Asia 0.0000 0.0438 0.1008 0.0430 0.0406 0.0766</span>
<span class="co">#&gt; Smoking 0.0130 0.0000 0.0435 0.2718 0.3458 0.1531</span>
<span class="co">#&gt; Tuberculosis 0.0303 0.1824 0.0000 0.3291 0.0083 0.4575</span>
<span class="co">#&gt; LungCancer 0.0078 0.4355 0.1595 0.0000 0.0400 0.3393</span>
<span class="co">#&gt; Bronchitis 0.0162 0.5455 0.0112 0.0533 0.0000 0.1163</span>
<span class="co">#&gt; Either 0.0099 0.2340 0.2737 0.3036 0.0746 0.0000</span>
<span class="co">#&gt; XRay 0.0164 0.2201 0.2814 0.2405 0.0429 0.4167</span>
<span class="co">#&gt; Dyspnea 0.0102 0.1868 0.0916 0.1908 0.5922 0.4002</span>
<span class="co">#&gt; XRay Dyspnea</span>
<span class="co">#&gt; Asia 0.0519 0.0515</span>
<span class="co">#&gt; Smoking 0.1100 0.1712</span>
<span class="co">#&gt; Tuberculosis 0.1814 0.1053</span>
<span class="co">#&gt; LungCancer 0.2143 0.2009</span>
<span class="co">#&gt; Bronchitis 0.0559 0.3156</span>
<span class="co">#&gt; Either 0.3485 0.3162</span>
<span class="co">#&gt; XRay 0.0000 0.1158</span>
<span class="co">#&gt; Dyspnea 0.1415 0.0000</span></code></pre></div>
<p>The <code>query()</code> function recognize fomula statements. So it is possible to explicitelly query: <strong>what is the probability of LungCancer node being children of the Smoking node?</strong></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">query</span>(<span class="dt">mcmcabn =</span> mcmc.out,<span class="dt">formula =</span> <span class="op">~</span>LungCancer<span class="op">|</span>Smoking)
<span class="co">#&gt; [1] 0.435</span></code></pre></div>
<p>This means that an arrow from Smoking to LungCancer appears in 43.546% of the sampled DAGs.</p>
<p>One can also ask more complicated strctural requests. If want to know <strong>what is the probability of Smoking node being parent of both LungCancer and Bronchitis node?</strong></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">query</span>(<span class="dt">mcmcabn =</span> mcmc.out,<span class="dt">formula =</span> <span class="op">~</span><span class="st"> </span>LungCancer<span class="op">|</span>Smoking<span class="op">+</span>Bronchitis<span class="op">|</span>Smoking)
<span class="co">#&gt; [1] 0.229</span></code></pre></div>
<p>It means that an arrow from Smoking to LungCancer <strong>and</strong> an arrow from Smoking to Bronchitis appears in 22.868% of the sampled DAGs. This probability cannot be read in the first matrix above. Indeed, it is a joint statement.</p>
<p>If one want to query positive and negative statements such as <strong>What is the probability of previous statement when there is no arc from Smoking to Tuberculosis and from Bronchitis to XRay?</strong></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">query</span>(<span class="dt">mcmcabn =</span> mcmc.out,<span class="dt">formula =</span> <span class="op">~</span>LungCancer<span class="op">|</span>Smoking<span class="op">+</span>Bronchitis<span class="op">|</span>Smoking<span class="op">-</span>Tuberculosis<span class="op">|</span>Smoking<span class="op">-</span>XRay<span class="op">|</span>Bronchitis)
<span class="co">#&gt; [1] 9e-04</span></code></pre></div>
<p>Essentially zero!</p>
<div id="formula-statement-tutorial" class="section level3">
<h3 class="hasAnchor">
<a href="#formula-statement-tutorial" class="anchor"></a>Formula statement: tutorial</h3>
<p>The <strong>formula</strong> statement has been designed to ease querying over the MCMC sample. It without explicitly writing an adjacency matrix (which can be painful when the number of variables increases). The <code>formula</code> argument can be provided using typically a formula like: ~ node1|parent1:parent2 + node2:node3|parent3. The formula statement has to start with <code>~</code>. In this example, node1 has two parents (parent1 and parent2). node2 and node3 have the same parent3. The parents names have to exactly match those given in name. <code>:</code> is the separator between either children or parents, <code>|</code> separates children (left side) and parents (right side), <code>+</code> separates terms, <code>.</code> replaces all the variables in name. Additionnaly, when one want to exclude an arc simply put <code>-</code> in front of that statement. Then a formula like: ~ -node1|parent1 exclude all DAGs that have an arc between parent1 and node1. Alternativelly, one can query using an adjacency matrix. This matrix should have only: 0,1 and -1. the 1 indicates the requested arcs, the -1 the excludes and the 0 all other entries that are not subject to constrain. The rows indicated the set of parents of the index nodes. The order of rows and column should be the same as the one used in the <code><a href="../reference/mcmc.html">mcmcabn()</a></code> function in the <code>data.dist</code> argument. the The matrix should not be named but be squared.</p>
</div>
<div id="technical-fundations" class="section level2">
<h2 class="hasAnchor">
<a href="#technical-fundations" class="anchor"></a>Technical fundations</h2>
<p>This part aimed at describing the methodology used to</p>
<p>This part aimed at describing the implemented methodology in <em>mcmcabn</em>.</p>
<div id="algorithms" class="section level3">
<h3 class="hasAnchor">
<a href="#algorithms" class="anchor"></a>Algorithms</h3>
......
......@@ -87,7 +87,7 @@
</a>
</li>
<li>
<a href="https://git.math.uzh.ch/gkratz/varrank">
<a href="https://git.math.uzh.ch/gkratz/mcmcabn">
<span class="fa fa-gitlab fa-lg"></span>
gitlab
......
......@@ -57,7 +57,7 @@
</a>
</li>
<li>
<a href="https://git.math.uzh.ch/gkratz/varrank">
<a href="https://git.math.uzh.ch/gkratz/mcmcabn">
<span class="fa fa-gitlab fa-lg"></span>
gitlab
......
......@@ -87,7 +87,7 @@
</a>
</li>
<li>
<a href="https://git.math.uzh.ch/gkratz/varrank">
<a href="https://git.math.uzh.ch/gkratz/mcmcabn">
<span class="fa fa-gitlab fa-lg"></span>
gitlab
......@@ -113,17 +113,9 @@
</div>
<div id="varrank-0-2" class="section level2">
<div id="mcmcabn-0-0-900" class="section level2">
<h2 class="hasAnchor">
<a href="#varrank-0-2" class="anchor"></a>varrank 0.2:</h2>
<ul>
<li><p>pkgdown website</p></li>
<li><p>Bug in plot function</p></li>
</ul>
</div>
<div id="varrank-0-1" class="section level2">
<h2 class="hasAnchor">
<a href="#varrank-0-1" class="anchor"></a>varrank 0.1:</h2>
<a href="#mcmcabn-0-0-900" class="anchor"></a>mcmcabn 0.0.900:</h2>
<ul>
<li>first release on CRAN</li>
</ul>
......@@ -134,8 +126,7 @@
<div id="tocnav">
<h2>Contents</h2>
<ul class="nav nav-pills nav-stacked">
<li><a href="#varrank-0-2">0.2</a></li>
<li><a href="#varrank-0-1">0.1</a></li>
<li><a href="#mcmcabn-0-0-900">0.0.900</a></li>
</ul>
</div>
</div>
......
......@@ -87,7 +87,7 @@
</a>
</li>
<li>
<a href="https://git.math.uzh.ch/gkratz/varrank">
<a href="https://git.math.uzh.ch/gkratz/mcmcabn">
<span class="fa fa-gitlab fa-lg"></span>
gitlab
......
......@@ -90,7 +90,7 @@
</a>
</li>
<li>
<a href="https://git.math.uzh.ch/gkratz/varrank">
<a href="https://git.math.uzh.ch/gkratz/mcmcabn">
<span class="fa fa-gitlab fa-lg"></span>
gitlab
......
% mcmc.Rd ---
% mcmcabn.Rd ---
% Author : Gilles Kratzer
% Created on : 18.02.2019
% Last modification :
......@@ -6,36 +6,37 @@
\name{mcmcabn}
\alias{mcmcabn}
\title{Structural MCMC sampler for DAGs.}
\title{Structural MCMC sampler for DAGs}
\usage{
mcmcabn(score.cache = NULL,
score = "mlik",
data.dists = NULL,
max.parents = 1,
mcmc.scheme = c(100,1000,1000), #returned, #thinned, #burned
mcmc.scheme = c(100,1000,1000),
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.choice = 2
prior.choice = 2)
}
\arguments{
\item{score.cache}{output from \code{buildscorecache()} from the \code{abn} R package.}
\item{score}{character giving which network score should be used to sample the DAGs landscape.}
\item{data.dists}{a named list giving the distribution for each node in the network, see details}
\item{max.parents}{a constant giving the maximum number of parents allowed.}
\item{mcmc.scheme}{a sampling scheme. It is vector giving in that order: the number of returned DAGS, the number of thinned steps and length of the burn in phase.}
\item{seed}{a non-negative integer which sets the seed.}
\item{verbose}{extra output, see output for details.}
\item{start.dag}{a DAG given as a matrix, see details for format, which can be used to explicity provide a starting point for the structural search. Alternativelly charachter "random" will select a random DAG as strating point.}
\item{dag.banned}{a matrix or a formula statement defining which arcs are not permitted - banned - see details for format. Note that colnames and rownames must be set, otherwise same row/column names as \code{score.cache} will be assumed. If set as NULL an empty matrix is assumed.}
\item{dag.retained}{a matrix or a formula statement (see details for format) defining which arcs are must be retained in any model search, see details for format. Note that colnames and rownames must be set, otherwise same row/column names as \code{score.cache} will be assumed. If set as NULL an empty matrix is assumed.}
\item{start.dag}{a DAG given as a matrix, see details for format, which can be used to explicitly provide a starting point for the structural search. Alternatively character "random" will select a random DAG as starting point.}
\item{prior.dag}{user defined prior. It should be given as a matrix where entries range from zero to one. 0.5 is non-informative for the given arc.}
\item{prior.lambda}{hyper parameter representing the strength of belief in the user defined prior.}
\item{prob.rev}{probability of selecting a new edge reversal.}
\item{prob.mbr}{proability of selecting a Markov blanket resampling move.}
\item{prob.mbr}{probability of selecting a Markov blanket resampling move.}
\item{prior.choice}{an integer, 1 or 2, where 1 is a uniform structural prior and 2 uses a weighted prior, see details}
}
......@@ -43,21 +44,49 @@ mcmcabn(score.cache = NULL,
This function is a structural Monte Carlo Markov Chain Model Choice (MC^3) that is equipped with two radical possible MCMC moves that are purposed to accelerate chain mixing.
}
\details{The procedure runs a structural Monte Carlo Markov Chain Model Choice (MC^3) to find the most probable posterior network (DAG). The default algorithm is based on three MCMC move: edge addition, edge deletion and edge reversal. This algorithm is known as the MC^3 (Monte Carlo Markov Chain Model Choice). It is know to mix slowly and getting stuck in low probability regions. Indeed, changing of Markov equivalence region often requires multiple MCMC moves. Then two radical MCMC move are implemented and can be probability wise tunned by user. The new edge reversal move (REV) from Grzegorczyk and Husmeier (2008) and the Markov blanket resampling (MBR) from Su and Borsuk (2016). The classical reversal move depends on the global configuration of the parents and children, but then fails to propose MCMC jumps that produce valid but very different DAGs in a unique move. The MBR workaround applies the same idea but to the entire Markov blanket of a randomly chosen node.
\details{The procedure runs a structural Monte Carlo Markov Chain Model Choice (MC^3) to find the most probable posterior network (DAG). The default algorithm is based on three MCMC move: edge addition, edge deletion and edge reversal. This algorithm is known as the MC^3 (Monte Carlo Markov Chain Model Choice). It is know to mix slowly and getting stuck in low probability regions. Indeed, changing of Markov equivalence region often requires multiple MCMC moves. Then two radical MCMC move are implemented and can be probability wise tuned by user. The new edge reversal move (REV) from Grzegorczyk and Husmeier (2008) and the Markov blanket resampling (MBR) from Su and Borsuk (2016). The classical reversal move depends on the global configuration of the parents and children, but then fails to propose MCMC jumps that produce valid but very different DAGs in a unique move. The MBR workaround applies the same idea but to the entire Markov blanket of a randomly chosen node.
The classical MC^3 is unbiased but ineficient in mixing, the two radical MCMC alternative move are known to massivelly accelerate mixing without introducing biases. But those move are computationally expensive. Then low frequencies are advised.
The classical MC^3 is unbiased but inefficient in mixing, the two radical MCMC alternative move are known to massively accelerate mixing without introducing biases. But those move are computationally expensive. Then low frequencies are advised.
The parameter prior.choice determines the prior used within each individual node for a given choice of parent combination. In Koivisto and Sood (2004) p.554 a form of prior is used which assumes that the prior probability for parent combinations comprising of the same number of parents are all equal. Specifically, that the prior probability for parent set G with cardinality |G| is proportional to 1/[n-1 choose |G|] where there are n total nodes. Note that this favours parent combinations with either very low or very high cardinality which may not be appropriate. This prior is used when prior.choice=2. When prior.choice=1 an uninformative prior is used where parent combinations of all cardinalities are equally likely.
The parameter prior.choice determines the prior used within each individual node for a given choice of parent combination. In Koivisto and Sood (2004) p.554 a form of prior is used which assumes that the prior probability for parent combinations comprising of the same number of parents are all equal. Specifically, that the prior probability for parent set G with cardinality |G| is proportional to 1/[n-1 choose |G|] where there are n total nodes. Note that this favours parent combinations with either very low or very high cardinality which may not be appropriate. This prior is used when \code{prior.choice=2}. When prior.choice=1 an uninformative prior is used where parent combinations of all cardinalities are equally likely. When \code{prior.choice=3} a user defined prior is used, defined by \code{prior.dag}. It is given by an adjacency matrix (squared and same size as number of nodes) where entries ranging from zero to one give the user prior belief. An hyper parameter defining the global user belief in the prior is given by \code{prior.lambda}.
MCMC sampler came with asymptotic statistical guarantees. Therefore it is highly advised to run multiple long enough chains. The burn in phase length (i.e throwing away first MCMC iterations) should be adequatelly chosen.
MCMC sampler came with asymptotic statistical guarantees. Therefore it is highly advised to run multiple long enough chains. The burn in phase length (i.e throwing away first MCMC iterations) should be adequately chosen.
Binary variables must be declared as factors with two levels, and the argument data.dists must be a list with named arguments, one for each of the variables in \code{data.df} (except a grouping variable - if present), where each entry is either "poisson","binomial", or "gaussian"
}
\value{A list with an entry for the list of sampled DAGs, the list of scores and the acceptance probability.}
\value{A list with an entry for the list of sampled DAGs, the list of scores, the acceptance probability, the method used for each MCMC jump, the rejection status for each MCMC jump, the total number of iterations the thinning, the length of burn in phase and the named list of distribution per node. The returned object is of class mcmcabn.}
\author{Gilles Kratzer}
\references{Kratzer, G. Furrer, R. "Is a single unique Bayesian network enough to accurately represent your data?". arXiv preprint.}
\references{Kratzer, G. Furrer, R. "Is a single unique Bayesian network enough to accurately represent your data?". arXiv preprint arXiv:1902.06641.
\examples{
Grzegorczyk, M.Husmeier, D. "Improving the structure MCMC sampler for Bayesian networks by introducing a new edge reversal move", Machine Learning, vol. 71(2-3), pp. 265, 2008.
Su, C.Borsuk, M. E. "Improving structure MCMC for Bayesian networks through Markov blanket resampling", The Journal of Machine Learning Research, vol. 17(1), pp. 4042-4061, 2016.
Werhli, A. V., & Husmeier, D. (2007). "Reconstructing gene regulatory networks with Bayesian networks by combining expression data with multiple sources of prior knowledge". Statistical Applications in Genetics and Molecular Biology, 6 (Article 15).
Scutari, M. (2010). Learning Bayesian Networks with the bnlearn R Package. Journal of Statistical Software, 35(3), 1 - 22. doi:http://dx.doi.org/10.18637/jss.v035.i03.
}
\donotrun{
\examples{
## Example from the asia dataset from Lauritzen and Spiegelhalter (1988) provided by Scutari (2010)
data("mcmc_run_asia")
mcmc.out <- mcmcabn(score.cache = bsc.compute,
score = "mlik",
data.dists = dist,
max.parents = 2,
mcmc.scheme = c(10000,9,1000),
seed = 6,
verbose = FALSE,
start.dag = "random",
prob.rev = 0.005,
prob.mbr = 0.005,
prior.choice = 1)
summary(mcmc.out)
}}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment