mcmc_1sampleDAG0: Implementation of the parallel-tempered MCMC for one-sample...

In junsoukchoi/BayesDAG0: Bayesian Differential Zero-Inflated Negative Binomial Directed Acyclic Graphs

Implementation of the parallel-tempered MCMC for one-sample DAG0

Description

Implementation of the parallel-tempered MCMC for one-sample DAG0

Usage

mcmc_1sampleDAG0(
  data,
  starting,
  tuning,
  priors,
  n_sample = 5000,
  n_burnin = 2500,
  n_chain = 10,
  prob_swap = 0.1,
  temperature = NULL,
  do_par = FALSE,
  n_core = n_chain,
  verbose = TRUE,
  n_report = 100
)

Arguments

data	a matrix containing data.
starting	a list with each tag corresponding to a parameter name. Valid tags are 'E', 'alpha', 'beta', 'delta', 'gamma', 'psi', 'tau', and 'rho'. The value portion of each tag is the parameters' starting values for MCMC.
tuning	a list with each tag corresponding to a parameter name. Valid tags are 'E', 'E_rev', 'alpha', 'beta', 'delta', 'gamma', and 'psi'. The value portion of each tag defines the standard deviations of Normal proposal distributions for Metropolis sampler for each parameter. The tag 'E' corresponds to birth or death of an edge, while 'E_rev' indicates reversal of an edge. The value portion of both tags should consist of the standard deviations of Gaussian proposals for alpha, beta, delta, and gamma for an update of E.
priors	a list with each tag corresponding to a parameter name. Valid tags are 'psi', 'tau', and 'rho'. The value portion of each tag defines the hyperparameters for one-sample DAG0.
n_sample	the number of MCMC iterations.
n_burnin	the number of burn-in samples.
n_chain	the number of Markov chains for the parallel tempering.
prob_swap	the probability of a swapping step.
temperature	a sequence of increasing temperatures. Should have length of n_chain. Default value is T_l = (1.5)^{(l - 1) / (n_chain - 1)} for l = 1, ..., n_chain.
do_par	if TRUE, the parallel computation is used to update n_chain Markov chains in parallel. The parallel computation is only available on Unix-alike platforms as we create the worker process by forking.
n_core	the number of cores to be used for the parallel computation when do_par = TRUE.
verbose	if TRUE, progress of the sampler is printed to the screen. Otherwise, nothing is printed to the screen.
n_report	the interval to report MCMC progress.

Value

An object of class 1sampleDAG0, which is a list with the tags 'samples' and 'acceptance'. The value portion of the tag 'samples' gives MCMC samples from the posterior distribution of the parameters for one-sample DAG0. The value portion of the tag 'acceptance' shows the Metropolis sampling acceptance percents for alpha, beta, delta, gamma, and psi.

Examples

library(BayesDAG0)
library(igraph)
library(pscl)

# set random seed
set.seed(1)

# generate a random DAG E
p   = 10       # the number of variables 
n_e = p        # the number of edges
E_true = matrix(0, p, p)
while (sum(E_true == 1) < n_e)
{
   id_edge = matrix(sample(1 : p, 2), ncol = 2)
   E_true[id_edge] = 1
   g_true = graph_from_adjacency_matrix(t(E_true))
   if (!(is_dag(g_true)))
      E_true[id_edge] = 0
}

# generate model parameters for one-sample DAG0
alpha_true = matrix(0, p, p)
alpha_true[E_true == 1] = runif(n_e, 0.3, 1)
beta_true  = matrix(0, p, p)
beta_true[E_true == 1]  = runif(n_e, -1, -0.3)
delta_true = runif(p, -2, -1)
gamma_true = runif(p, 1, 2)
psi_true   = 1 / runif(p, 1, 5)

# generate data from one-sample DAG0 with given DAG and model parameters
dat = generate_data_DAG0(n = 100, E_true, alpha_true, beta_true, delta_true, gamma_true, psi_true)

### conduct posterior inference for one-sample DAG0
# starting values for MCMC 
starting = list()
starting$E     = matrix(0, ncol(dat), ncol(dat))
starting$alpha = matrix(0, ncol(dat), ncol(dat))
starting$beta  = matrix(0, ncol(dat), ncol(dat))
starting$delta = rep(0, ncol(dat))
starting$gamma = rep(0, ncol(dat))
starting$psi   = rep(0, ncol(dat))
for (j in 1 : ncol(dat))
{
   mle = zeroinfl(dat[ , j] ~ 1 | 1, dist = "negbin")
   starting$delta[j] = mle$coefficients$zero
   starting$gamma[j] = mle$coefficients$count
   starting$psi[j]   = 1 / mle$theta
}

starting$tau = c(1, 1, 1, 1)
starting$rho = 0.5

# sd's of the Metropolis sampler Normal proposal distribution
tuning = list()
tuning$E     = c(0.05, 0.05, 1.2, 0.3)
tuning$E_rev = c(0.5, 0.5, 1.2, 0.3) 
tuning$alpha = 1.2  
tuning$beta  = 0.4  
tuning$delta = 1.2  
tuning$gamma = 0.3 
tuning$psi   = 0.8  

# hyperparameters for one-sample DAG0
priors = list()
priors$psi = c(1, 1)
priors$tau = c(1, 1)
priors$rho = c(1, 1)

# run the parallel-tempered MCMC for one-sample DAG0
# the parallel computation (do_par = TRUE) is only available on Unix-alike platforms 
# Windows users should use do_par = FALSE not to allow the parallel computation
out = mcmc_1sampleDAG0(dat, starting, tuning, priors, n_sample = 2000, n_burnin = 1000, do_par = TRUE)

# report Metropolis sampling acceptance rates
out$acceptance

# graph estimate based on median probability model
E_est  = (apply(out$samples$E, c(1, 2), mean) > 0.5) + 0

# report contingency table for edge selection
table(E_est, E_true)

junsoukchoi/BayesDAG0 documentation built on Feb. 26, 2025, 2:37 p.m.