MCMC.GBN: MCMC algorithm to estimate causal effetcts of a GBN object
In andreamrau/GBNcausal: Causal estimation for Gaussian Bayesian networks

Description Usage Arguments Details Value See Also Examples

MCMC.GBN is used to infer causal effects (direct or undirect) of a GB network. It is a random walk MCMC Metropolis Hasting algorithm with a Mallows proposal distribution. The algorithm explores the space of nodes order and find the best estimation of the causal effects. In each iteration, it computes the maximum likelihood estimation of the GBN law, based on the data.

1
2

MCMC.GBN(data, firstGBN, nbSimulation = 20000, burnIn = 5000, seq = 25, verbose = FALSE, verbose.index = 500, 
 alpha = 0.05, alpha2 = 0.05, lambda = 0, listblocks = list(), str = matrix(1, length(firstGBN@resSigma), length(firstGBN@resSigma)), type = "")

`data`	data - Can be obtained by the function `dataFormat`.
`firstGBN`	GBN - An object of type GBN which is the initial value of the algorithm. data and firstGBN must have the same rownames and colnames for each elements.
`nbSimulation`	integer - The number of iterations the function has to do.
`burnIn`	integer - The number of iterations to do before starting to save the trajectory of the algorithm.
`seq`	integer - Every seq, if the number of iterations done is superior to the burnIn, the GBN found is saved into full.run.
`verbose`	logical - If TRUE, every verbose.index interation a message is written to give the number of iterations made so far by the algorithm.
`verbose.index`	integer - Indicates the number of iteration to do before writting a message if verbose is TRUE.
`alpha`	double - First temperature of the mallows function.
`alpha2`	double - Second temperature of the mallows function. If listblocks is empty, this parameter will not be used.
`lambda`	logarithmic - Coefficient of the penalty Ridge. Can be 0.
`listblocks`	list - A list of nodes of the form (c("N1","N2"),c("N3","N4","N5")) where "N1","N2","N3","N4" and "N5" are elements of rownames and colnames of firstGBN elements and data elements.
`str`	matrix - To improve the efficiency of the algorithm, a structure can be add. Colnames and rownames are not needed. Once the structure is set, the algorithm can't add any interaction that's not in the structure. It can remove some of them.
`type`	string - Type of observations. If type = obsOnly, it means the data are observational only, and the function will use a uniform sampling each step instead of a Mallows distribution.

This function only works with acyclic directed networks (DAG).

It allows to insert an a priori information, which can be a structure or a list of blocks.

A structure is in general a skeleton of a graph, or a CPDAG. It can be obtained by a PC algorithm (PCAlg) or a Lasso (glasso) method.

The structure is a higher constraint than the list of blocks, in the way that the algorithm can't add any direct causal effect that is not in the structure. It can remove some of them. However, this is a good way of improving the efficiency of the algorithm : it's faster and give better results if the structure corresponds to the network to estimate.

The list of blocks is a list composed of vector of nodes that we know they are close in the graph. If the list is not empty, the algorithm uses the rmallowsBlocks function to propose a new order of nodes. The elements of a block can't be separated.

This method can take observational and interventional data (only knock-out are implemented). Infering on observational data only can't give any good results.

`full.run`	The full trajectory of the algorithm : after the burnIn iteration, every seq iteration the GBN is saved into the full.run. Inference is done on full.run.
`accept`	The acceptance rate of the algorithm.

MCMC.step, newMallowsProposal, rmallows, rmallowsBlocks

# Data creation

seed = 1990
n = 3000
p <- 10
m<-rep(0,10)
sigma<-rep(0.1,10) 

W <- 1*upper.tri(matrix(0,p,p))

data <- dataCreate(nbData = 2*p, p = 10,KO = list(1,9), nbKO = c(p,p), W = W , m = m,sigma = sigma, seed = seed)$data

# Initial Value

W1=1*upper.tri(matrix(0,p,p)) 
m1=rep(0,p)
s1=rep(10e-4,p)
colnames(W1)=names(m1)=names(s1)=rownames(W1)=paste("N",1:p,sep="")

firstGBN = new("GBNetwork",WeightMatrix=W1,resMean=m1,resSigma=s1)
firstGBN = GBNmle(firstGBN,data,lambda= 0,sigmapre=s1)$GBN

# Algorithm
# Without constraint
results=MCMC.GBN(data, firstGBN, nbSimulation=2*n, burnIn=n, seq=1, verbose=TRUE,
             verbose.index=n/10, alpha=0.7,lambda=0.1)

# With constraint
listBlocks <- list(c("N1","N2","N3"),c("N4","N5"))

results_Blocks=MCMC.GBN(data, firstGBN, nbSimulation=2*n, burnIn=n, seq=1, verbose=TRUE,
                 verbose.index=n/10, alpha=0.7, alpha2 = 0.5, lambda=0.1, listblocks = listBlocks)
				 
# Observational data only :
data <- dataCreate(nbData = 2*p, p = 10,KO = list(), nbKO = c(), W = W , m = m,sigma = sigma, seed = seed)$data

firstGBN = new("GBNetwork",WeightMatrix=W1,resMean=m1,resSigma=s1)
firstGBN = GBNmle(firstGBN,data,lambda= 0,sigmapre=s1)$GBN

results_Blocks=MCMC.GBN(data, firstGBN, nbSimulation=100, burnIn=0, seq=1, verbose=TRUE,
                 verbose.index=n/10, alpha=0.7, alpha2 = 0.5, lambda=0.1, listblocks = listBlocks, type = "obsOnly")