Commit 498a746e authored by Gilles Kratzer's avatar Gilles Kratzer
Browse files

cran submission v 0.4 - change functions abn 2.5

parent aef5972c
Pipeline #5258 passed with stage
in 7 seconds
Package: mcmcabn
Title: Flexible Implementation of a Structural MCMC Sampler for DAGs
Version: 0.3
Version: 0.4
Authors@R: c(person("Gilles", "Kratzer", role = c("aut", "cre"),
email = "gilles.kratzer@math.uzh.ch",
comment = c(ORCID = "0000-0002-5929-8935")),
......@@ -8,7 +8,7 @@ Authors@R: c(person("Gilles", "Kratzer", role = c("aut", "cre"),
email = "reinhard.furrer@math.uzh.ch",
comment = c(ORCID = "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 us to quantify the marginal impact of relationships of interest by marginalizing out over structures or nuisance dependencies. Structural MCMC seems an 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.
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) <https://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 us to quantify the marginal impact of relationships of interest by marginalizing out over structures or nuisance dependencies. Structural MCMC seems an 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)
License: GPL-3
Encoding: UTF-8
......@@ -17,5 +17,5 @@ Imports: gRbase, abn, coda, ggplot2, cowplot, ggpubr
Suggests: bnlearn, knitr, rmarkdown, ggdag, testthat
VignetteBuilder: knitr
URL: https://www.math.uzh.ch/pages/mcmcabn/
BugReports: https://git.math.uzh.ch/gkratz/mcmcabn/issues
BugReports: https://git.math.uzh.ch/gkratz/mcmcabn/-/issues
RoxygenNote: 6.1.1
......@@ -11,3 +11,7 @@
* update mcmcabn to make it compatible with constraints imported from the cache of scores. Heating parameter to increase or decrease acceptance probability.
* new function CoupledHeatedmcmcabn() implementing parallel tempering
* new published article (preferred reference): Bayesian Network Modeling Applied to Feline Calicivirus Infection Among Cats in Switzerland in Front. Vet. Sci.
## mcmcmabn 0.4:
* update mcmcabn for interaction with abn 2.4
......@@ -11,7 +11,7 @@
stop("A cache of score should be provided. You can produce it using the R package abn.")
if (max(rowSums(score.cache$node.defn)) > (max.parents+1))
stop("Check max.parents. It should be the same as the one used in abn::buildscorecache() R function")
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) made of non negative integers.")
......
......@@ -87,7 +87,7 @@ for (i in 1:11) {
max.par <- i
# construction of the score cache
mycache <- buildscorecache(data.df = mice_output,
mycache <- buildScoreCache(data.df = mice_output,
data.dists = dists,
dag.banned = ~Sex|.+Age|.+Pedigree|.,
max.parents = max.par,method = "mle")
......@@ -145,7 +145,7 @@ dev.off()
## --------------------------
## adjustment for overfitting
mycache <- buildscorecache(data.df = mice_output,
mycache <- buildScoreCache(data.df = mice_output,
data.dists = dists,
dag.banned = ~Sex|.+Age|.+Pedigree|.,
max.parents = 7,method = "mle")
......@@ -466,7 +466,7 @@ tographviz(dag.m = dag.bic,data.df = mice_output,data.dists = dists_shame,outfil
##aic 10 parents
##bic 7 parents
mycache <- buildscorecache(data.df = mice_output,
mycache <- buildScoreCache(data.df = mice_output,
data.dists = dists,
dag.banned = ~Sex|.+Age|.+Pedigree|.,
max.parents = 10,method = "mle")
......
......@@ -27,7 +27,7 @@ CoupledHeatedmcmcabn(score.cache = NULL,
}
\arguments{
\item{score.cache}{output from \link[abn:build_score_cache]{buildscorecache} from the \code{abn} R package.}
\item{score.cache}{output from \link[abn:build_score_cache]{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.}
......
......@@ -9,7 +9,7 @@
\usage{bsc.compute.asia}
\format{
The data contains a cache of pre-computed scores with a maximum of two parents per node using \link[abn:build_score_cache]{buildscorecache} from the \code{abn} R package.
The data contains a cache of pre-computed scores with a maximum of two parents per node using \link[abn:build_score_cache]{buildScoreCache} from the \code{abn} R package.
\itemize{
\item \code{bsc.compute.asia}: cache of score with a maximum of two parents per node.
}}
......@@ -29,7 +29,7 @@ colnames(asia) <- c("Asia",
"XRay",
"Dyspnea")
bsc.compute.asia <- buildscorecache(data.df = asia,
bsc.compute.asia <- buildScoreCache(data.df = asia,
data.dists = dist.asia,
max.parents = 2)
}
......
......@@ -26,7 +26,7 @@ mcmcabn(score.cache = NULL,
}
\arguments{
\item{score.cache}{output from \link[abn:build_score_cache]{buildscorecache} from the \code{abn} R package.}
\item{score.cache}{output from \link[abn:build_score_cache]{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.}
......
......@@ -11,9 +11,9 @@
\format{
The data contains an object of class mcmcabn and a cache of score computed using
\link[abn:build_score_cache]{buildscorecache} from the \code{abn} R package.
\link[abn:build_score_cache]{buildScoreCache} from the \code{abn} R package.
\itemize{
\item \code{bsc.compute.asia}: cache of score with a maximum of two parents per node computed using \link[abn:build_score_cache]{buildscorecache} from the \code{abn} R package.
\item \code{bsc.compute.asia}: cache of score with a maximum of two parents per node computed using \link[abn:build_score_cache]{buildScoreCache} from the \code{abn} R package.
\item \code{dist.asia}: a named list giving the distribution for each node in the network.
\item \code{mcmc.out.asia}: an object of class mcmcabn.
}}
......@@ -44,7 +44,7 @@ dist.asia <- list(Asia = "binomial",
XRay = "binomial",
Dyspnea = "binomial")
bsc.compute.asia <- buildscorecache(data.df = asia,
bsc.compute.asia <- buildScoreCache(data.df = asia,
data.dists = dist.asia,
max.parents = 2)
......
......@@ -35,7 +35,7 @@ test_that("mc3",{
out.sim.0 <- simulateAbn(data.dists = dist,n.chains = 1,n.adapt = 1000,n.thin = 1,n.iter = 1000,data.param = 0.4*data.param, simulate = TRUE,seed = 132,verbose = FALSE)
bsc.compute.0 <- buildscorecache(data.df = out.sim.0, data.dists = dist, max.parents = 3)
bsc.compute.0 <- buildScoreCache(data.df = out.sim.0, data.dists = dist, max.parents = 3)
mc3.out <- mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
......@@ -49,8 +49,8 @@ test_that("mc3",{
prob.mbr = 0,
prior.choice = 1,heating = 0.5)
dag <- mostprobable(score.cache = bsc.compute.0,prior.choice = 1)
expect_equal(max(mc3.out$scores),fitabn(object = dag,data.df = out.sim.0,data.dists = dist)$mlik)
dag <- mostProbable(score.cache = bsc.compute.0,prior.choice = 1)
expect_equal(max(mc3.out$scores),fitAbn(object = dag,data.df = out.sim.0,data.dists = dist)$mlik)
mc3.out <- mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
......@@ -64,8 +64,8 @@ test_that("mc3",{
prob.mbr = 0,
prior.choice = 2)
dag <- mostprobable(score.cache = bsc.compute.0,prior.choice = 2)
expect_equal(max(mc3.out$scores),fitabn(object = dag,data.df = out.sim.0,data.dists = dist)$mlik, tol=0.0001)
dag <- mostProbable(score.cache = bsc.compute.0,prior.choice = 2)
expect_equal(max(mc3.out$scores),fitAbn(object = dag,data.df = out.sim.0,data.dists = dist)$mlik, tol=0.0001)
#test influence of user define prior
data.param.eq <- matrix(data = 0,nrow = 5,ncol = 5)
......@@ -119,7 +119,7 @@ test_that("REV",{
out.sim.0 <- invisible(simulateAbn(data.dists = dist,n.chains = 1,n.adapt = 1000,n.thin = 1,n.iter = 1000,data.param = 0.4*data.param, simulate = TRUE,seed = 132))
bsc.compute.0 <- buildscorecache(data.df = out.sim.0, data.dists = dist, max.parents = 3)
bsc.compute.0 <- buildScoreCache(data.df = out.sim.0, data.dists = dist, max.parents = 3)
mc3.out <- mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
......@@ -133,8 +133,8 @@ test_that("REV",{
prob.mbr = 0,
prior.choice = 1)
dag <- mostprobable(score.cache = bsc.compute.0,prior.choice = 1)
expect_equal(max(mc3.out$scores),fitabn(object = dag,data.df = out.sim.0,data.dists = dist)$mlik)
dag <- mostProbable(score.cache = bsc.compute.0,prior.choice = 1)
expect_equal(max(mc3.out$scores),fitAbn(object = dag,data.df = out.sim.0,data.dists = dist)$mlik)
expect_silent(a<-mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
......@@ -161,7 +161,7 @@ test_that("REV",{
out.sim.0 <- invisible(simulateAbn(data.dists = dist,n.chains = 1,n.adapt = 1000,n.thin = 1,n.iter = 10,data.param = 0.4*data.param, simulate = TRUE,seed = 132))
bsc.compute.0 <- buildscorecache(data.df = out.sim.0, data.dists = dist, max.parents = 3)
bsc.compute.0 <- buildScoreCache(data.df = out.sim.0, data.dists = dist, max.parents = 3)
#
start.m <- matrix(data = c(0,0,0,1,0,0,0,0,0),nrow = 3L,ncol = 3L,byrow = TRUE)
......@@ -198,7 +198,7 @@ test_that("MBR",{
out.sim.0 <- invisible(simulateAbn(data.dists = dist,n.chains = 1,n.adapt = 1000,n.thin = 1,n.iter = 1000,data.param = 0.4*data.param, simulate = TRUE,seed = 132))
bsc.compute.0 <- buildscorecache(data.df = out.sim.0, data.dists = dist, max.parents = 3)
bsc.compute.0 <- buildScoreCache(data.df = out.sim.0, data.dists = dist, max.parents = 3)
mc3.out <- mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
......@@ -212,8 +212,8 @@ test_that("MBR",{
prob.mbr = 0.1,
prior.choice = 2)
dag <- mostprobable(score.cache = bsc.compute.0,prior.choice = 1)
expect_equal(max(mc3.out$scores),fitabn(object = dag,data.df = out.sim.0,data.dists = dist)$mlik)
dag <- mostProbable(score.cache = bsc.compute.0,prior.choice = 1)
expect_equal(max(mc3.out$scores),fitAbn(object = dag,data.df = out.sim.0,data.dists = dist)$mlik)
expect_silent(mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
......@@ -242,7 +242,7 @@ test_that("mcmcabn",{
out.sim.0 <- invisible(simulateAbn(data.dists = dist,n.chains = 1,n.adapt = 20,n.thin = 1,n.iter = 100,data.param = data.param.0, simulate = TRUE,seed = 132))
bsc.compute.0 <- buildscorecache(data.df = out.sim.0, data.dists = dist, max.parents = 2)
bsc.compute.0 <- buildScoreCache(data.df = out.sim.0, data.dists = dist, max.parents = 2)
expect_error(mcmcabn(score.cache = bsc.compute.0,
score = "blabla",
......@@ -344,7 +344,7 @@ test_that("mcmcabn",{
colnames(asia) <- c("Asia","Smoking", "Tuberculosis", "LungCancer", "Bronchitis", "Either", "XRay", "Dyspnea")
bsc.compute.asia <- buildscorecache(data.df = asia,
bsc.compute.asia <- buildScoreCache(data.df = asia,
data.dists = dist.asia,
max.parents = 2)
......@@ -363,8 +363,8 @@ test_that("mcmcabn",{
#maximum scoring network using exact search (not MCMC based)
dag <- mostprobable(score.cache = bsc.compute.asia)
expect_equal(max(mcmc.out.asia$scores),fitabn(object = dag,data.df = asia,data.dists = dist.asia)$mlik)
dag <- mostProbable(score.cache = bsc.compute.asia)
expect_equal(max(mcmc.out.asia$scores),fitAbn(object = dag,data.df = asia,data.dists = dist.asia)$mlik)
mcmc.out.asia <- mcmcabn(score.cache = bsc.compute.asia,
score = "mlik",
......@@ -378,7 +378,7 @@ test_that("mcmcabn",{
prob.mbr = 0.2,
prior.choice = 1,heating = 0.7)
expect_equal(max(mcmc.out.asia$scores),fitabn(object = dag,data.df = asia,data.dists = dist.asia)$mlik)
expect_equal(max(mcmc.out.asia$scores),fitAbn(object = dag,data.df = asia,data.dists = dist.asia)$mlik)
## marks datasets
......@@ -393,7 +393,7 @@ test_that("mcmcabn",{
#colnames(asia) <- c("Asia","Smoking", "Tuberculosis", "LungCancer", "Bronchitis", "Either", "XRay", "Dyspnea")
bsc.compute.marks <- buildscorecache(data.df = marks,
bsc.compute.marks <- buildScoreCache(data.df = marks,
data.dists = dist.marks,
max.parents = 2)
......@@ -412,8 +412,8 @@ test_that("mcmcabn",{
#maximum scoring network using exact search (not MCMC based)
dag <- mostprobable(score.cache = bsc.compute.marks)
expect_equal(max(mcmc.out.marks$scores),fitabn(object = dag,data.df = marks,data.dists = dist.marks)$mlik)
dag <- mostProbable(score.cache = bsc.compute.marks)
expect_equal(max(mcmc.out.marks$scores),fitAbn(object = dag,data.df = marks,data.dists = dist.marks)$mlik)
##tests
#data(gaussian.test)
......@@ -428,7 +428,7 @@ test_that("mcmcabn",{
colnames(gaussian.test) <- c("A","B","C","D","E","G","H")
bsc.compute.gaussian.test <- buildscorecache(data.df = gaussian.test,
bsc.compute.gaussian.test <- buildScoreCache(data.df = gaussian.test,
data.dists = dist.gaussian.test,
max.parents = 2)
......@@ -448,8 +448,8 @@ test_that("mcmcabn",{
#maximum scoring network using exact search (not MCMC based)
dag <- mostprobable(score.cache = bsc.compute.gaussian.test)
expect_equal(max(mcmc.out.gaussian.test$scores),fitabn(object = dag,data.df = gaussian.test,data.dists = dist.gaussian.test)$mlik, tol=10)
dag <- mostProbable(score.cache = bsc.compute.gaussian.test)
expect_equal(max(mcmc.out.gaussian.test$scores),fitAbn(object = dag,data.df = gaussian.test,data.dists = dist.gaussian.test)$mlik, tol=10)
})
......@@ -467,7 +467,7 @@ test_that("query",{
colnames(gaussian.test) <- c("A","B","C","D","E","G","H")
bsc.compute.gaussian.test <- buildscorecache(data.df = gaussian.test,
bsc.compute.gaussian.test <- buildScoreCache(data.df = gaussian.test,
data.dists = dist.gaussian.test,
max.parents = 2)
......@@ -506,7 +506,7 @@ dist.asia <- list(Asia = "binomial",
colnames(asia) <- c("Asia","Smoking", "Tuberculosis", "LungCancer", "Bronchitis", "Either", "XRay", "Dyspnea")
bsc.compute.asia <- buildscorecache(data.df = asia,
bsc.compute.asia <- buildScoreCache(data.df = asia,
data.dists = dist.asia,
max.parents = 2)
......@@ -524,8 +524,8 @@ mcmc.out.asia <- CoupledHeatedmcmcabn(score.cache = bsc.compute.asia,
prior.choice = 1,heating = 5,n.chains = 4)
#maximum scoring network using exact search (not MCMC based)
dag <- mostprobable(score.cache = bsc.compute.asia)
expect_equal(max(mcmc.out.asia$score.coupled),fitabn(object = dag,data.df = asia,data.dists = dist.asia)$mlik, tol=5)
dag <- mostProbable(score.cache = bsc.compute.asia)
expect_equal(max(mcmc.out.asia$score.coupled),fitAbn(object = dag,data.df = asia,data.dists = dist.asia)$mlik, tol=5)
dist.marks <- list(MECH = "gaussian",
......@@ -537,7 +537,7 @@ dist.marks <- list(MECH = "gaussian",
#colnames(asia) <- c("Asia","Smoking", "Tuberculosis", "LungCancer", "Bronchitis", "Either", "XRay", "Dyspnea")
bsc.compute.marks <- buildscorecache(data.df = marks,
bsc.compute.marks <- buildScoreCache(data.df = marks,
data.dists = dist.marks,
max.parents = 2)
......@@ -555,8 +555,8 @@ mcmc.out.marks <- CoupledHeatedmcmcabn(score.cache = bsc.compute.marks,
prior.choice = 1)
#maximum scoring network using exact search (not MCMC based)
dag <- mostprobable(score.cache = bsc.compute.marks)
expect_equal(max(mcmc.out.marks$scores),fitabn(object = dag,data.df = marks,data.dists = dist.marks)$mlik, tol = 0.5)
dag <- mostProbable(score.cache = bsc.compute.marks)
expect_equal(max(mcmc.out.marks$scores),fitAbn(object = dag,data.df = marks,data.dists = dist.marks)$mlik, tol = 0.5)
})
......
......@@ -38,7 +38,7 @@ Note: in this vignette *structural* MCMC means that the MCMC moves are performed
**How does it work?**
The `mcmcabn` R package takes a cache of pre-computed network scores as input (possibly computed using `buildscorecache()` function from `abn` R package). Then, the `mcmcabn()` function generates MCMC samples from the posterior distribution of the most supported DAGs. To do so the user should choose a structural prior and a learning scheme. `mcmcabn` 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).
The `mcmcabn` R package takes a cache of pre-computed network scores as input (possibly computed using `buildScoreCache()` function from `abn` R package). Then, the `mcmcabn()` function generates MCMC samples from the posterior distribution of the most supported DAGs. To do so the user should choose a structural prior and a learning scheme. `mcmcabn` 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).
**What are the R package functionalities?**
......@@ -115,12 +115,12 @@ For the first search, it is advised to limit the number of parents per node to 2
```{r, eval=FALSE}
#loglikelihood scores
bsc.compute.asia <- buildscorecache(data.df = asia,
bsc.compute.asia <- buildScoreCache(data.df = asia,
data.dists = dist.asia,
max.parents = 2)
```
We use the `mcmcabn()` function on the cache of pre-computed networks scores. One needs to define the type of score used (here is the marginal likelihood `mlik` other choice are: `bic`, `aic` or `mdl`). The maximum number of parents per node (same as the one used in `buildscorecache()`). The MCMC learning scheme: number of MCMC samples, number of thinned sampled (to avoid autocorrelation), and the length of the burn-in phase. Here we choose: 1000 returned MCMC steps and 99 thinned steps in between. It means that for each returned MCMC move, 99 have been thinned. The burn-in phase is set to 10000 steps. We decide to initiate the algorithm with a random DAG `start.dag = "random"`, which means that the function randomly selects a valid DAG. In this context, valid means no cycle and an appropriate maximum number of parents. The frequencies of selecting a REV or the MBR jumps are set to 3%. It is usually a good idea to keep those frequencies low. Indeed, REV and MBR are efficient in producing high scoring structure but computationally costly compared to classical MCMC jumps. We select the Koivisto prior `prior.choice = 2`.
We use the `mcmcabn()` function on the cache of pre-computed networks scores. One needs to define the type of score used (here is the marginal likelihood `mlik` other choice are: `bic`, `aic` or `mdl`). The maximum number of parents per node (same as the one used in `buildScoreCache()`). The MCMC learning scheme: number of MCMC samples, number of thinned sampled (to avoid autocorrelation), and the length of the burn-in phase. Here we choose: 1000 returned MCMC steps and 99 thinned steps in between. It means that for each returned MCMC move, 99 have been thinned. The burn-in phase is set to 10000 steps. We decide to initiate the algorithm with a random DAG `start.dag = "random"`, which means that the function randomly selects a valid DAG. In this context, valid means no cycle and an appropriate maximum number of parents. The frequencies of selecting a REV or the MBR jumps are set to 3%. It is usually a good idea to keep those frequencies low. Indeed, REV and MBR are efficient in producing high scoring structure but computationally costly compared to classical MCMC jumps. We select the Koivisto prior `prior.choice = 2`.
```{r, eval=FALSE}
mcmc.out.asia <- mcmcabn(score.cache = bsc.compute.asia,
......
Markdown is supported
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