Package 'bnlearn'

Description

Bayesian network structure learning (via constraint-based, score-based and hybrid algorithms), parameter learning (via ML and Bayesian estimators) and inference (via approximate inference algorithms).

Details

bnlearn implements key algorithms covering all stages of Bayesian network modelling: data preprocessing, structure learning combining data and expert/prior knowledge, parameter learning, and inference (including causal inference via do-calculus). bnlearn aims to be a one-stop shop for Bayesian networks in R. It provides tools for discrete, Gaussian and conditional linear Gaussian, and zero-inflated count Bayesian networks on real-world data. It also supportsd incomplete data with missing values. Furthermore, the modular nature of bnlearn makes it easy to use it for simulation studies.

Implemented structure learning algorithms include:

• Constraint-based algorithms, which use conditional independence tests to learn conditional independence constraints from data. The constraints, in turn, are used to learn the structure of the Bayesian network under the assumption that conditional independence implies graphical separation (so, two variables that are independent cannot be connected by an arc).

• Score-based algorithms, which are general-purpose optimisation algorithms that rank network structures with respect to a goodness-of-fit score.

• Hybrid algorithms combine aspects of both constraint-based and score-based algorithms. They use conditional independence tests (to reduce the search space) and network scores (to find the optimal network in the reduced space) at the same time.

For more details about structure learning algorithms, see structure learning. Available conditional independence tests are described in independence tests and available network scores are described in network scores. Specialized algorithms for learning the structure of Bayesian network classifiers are described in network classifiers. All algorithms support the use of whitelists and blacklists to include or exclude arcs from the networks (see whitelists and blacklists). Many have parallel implementation built on the parallel package. Bayesian network scores also support graphical priors.

Parameter learning approaches include both frequentist and Bayesian estimators. Approximate inference is implemented particle filters approaches, such as likelihood weighting, and covers conditional probability queries, prediction and imputation. Exact inference is implemented for discrete and Gaussian networks.

Additional facilities include support for bootstrap and cross-validation; advanced plotting capabilities implemented on top of Rgraphviz and lattice; model averaging; random graphs and random samples generation; import/export functions to integrate bnlearn with software such as Hugin and GeNIe; an associated Bayesian network repository of golden-standard networks at https://www.bnlearn.com/bnrepository/.

Use citation("bnlearn") to find out how to cite bnlearn in publications and other materials; and visit https://www.bnlearn.com/ for more examples and code from publications using bnlearn.

Author(s)

Marco Scutari
Istituto Dalle Molle di Studi sull'Intelligenza Artificiale (IDSIA)

Maintainer: Marco Scutari [email protected]

References

reference books:

Koller D, Friedman N (2009). Probabilistic Graphical Models: Principles and Techniques. MIT Press.

Korb K, Nicholson AE (2010). Bayesian Artificial Intelligence. Chapman & Hall/CRC, 2nd edition.

Pearl J (1988). Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference. Morgan Kaufmann.

from the author:

Nagarajan R, Scutari M, Lebre S (2013). "Bayesian Networks in R with Applications in Systems Biology." Springer.

Scutari M (2010). "Learning Bayesian Networks with the bnlearn R Package." Journal of Statistical Software, 35(3):1–22.

Scutari M (20107). "Bayesian Network Constraint-Based Structure Learning Algorithms: Parallel and Optimized Implementations in the bnlearn R Package." Journal of Statistical Software, 77(2):1–20.

Examples

## the workflow of Bayesian network modelling in bnlearn:
# choose the data set to work on...
data(learning.test)
# ... choose an algorithm and learn the structure of the network from the data...
net = hc(learning.test)
# ... plot it...
## Not run: graphviz.plot(net)
# ... learn the parameters of the network...
bn = bn.fit(net, learning.test)
# ... explore the network with a classic barchart...
## Not run: graphviz.chart(bn)
# ... and perform inference to answer any question that interests you!
cpquery(bn, event = (A == "a"), evidence = (C == "a"))

Description

The ALARM ("A Logical Alarm Reduction Mechanism") is a Bayesian network designed to provide an alarm message system for patient monitoring.

Usage

data(alarm)

Format

The alarm data set contains the following 37 variables:

• CVP (central venous pressure): a three-level factor with levels LOW, NORMAL and HIGH.

• PCWP (pulmonary capillary wedge pressure): a three-level factor with levels LOW, NORMAL and HIGH.

• HIST (history): a two-level factor with levels TRUE and FALSE.

• TPR (total peripheral resistance): a three-level factor with levels LOW, NORMAL and HIGH.

• BP (blood pressure): a three-level factor with levels LOW, NORMAL and HIGH.

• CO (cardiac output): a three-level factor with levels LOW, NORMAL and HIGH.

• HRBP (heart rate / blood pressure): a three-level factor with levels LOW, NORMAL and HIGH.

• HREK (heart rate measured by an EKG monitor): a three-level factor with levels LOW, NORMAL and HIGH.

• HRSA (heart rate / oxygen saturation): a three-level factor with levels LOW, NORMAL and HIGH.

• PAP (pulmonary artery pressure): a three-level factor with levels LOW, NORMAL and HIGH.

• SAO2 (arterial oxygen saturation): a three-level factor with levels LOW, NORMAL and HIGH.

• FIO2 (fraction of inspired oxygen): a two-level factor with levels LOW and NORMAL.

• PRSS (breathing pressure): a four-level factor with levels ZERO, LOW, NORMAL and HIGH.

• ECO2 (expelled CO2): a four-level factor with levels ZERO, LOW, NORMAL and HIGH.

• MINV (minimum volume): a four-level factor with levels ZERO, LOW, NORMAL and HIGH.

• MVS (minimum volume set): a three-level factor with levels LOW, NORMAL and HIGH.

• HYP (hypovolemia): a two-level factor with levels TRUE and FALSE.

• LVF (left ventricular failure): a two-level factor with levels TRUE and FALSE.

• APL (anaphylaxis): a two-level factor with levels TRUE and FALSE.

• ANES (insufficient anesthesia/analgesia): a two-level factor with levels TRUE and FALSE.

• PMB (pulmonary embolus): a two-level factor with levels TRUE and FALSE.

• INT (intubation): a three-level factor with levels NORMAL, ESOPHAGEAL and ONESIDED.

• KINK (kinked tube): a two-level factor with levels TRUE and FALSE.

• DISC (disconnection): a two-level factor with levels TRUE and FALSE.

• LVV (left ventricular end-diastolic volume): a three-level factor with levels LOW, NORMAL and HIGH.

• STKV (stroke volume): a three-level factor with levels LOW, NORMAL and HIGH.

• CCHL (catecholamine): a two-level factor with levels NORMAL and HIGH.

• ERLO (error low output): a two-level factor with levels TRUE and FALSE.

• HR (heart rate): a three-level factor with levels LOW, NORMAL and HIGH.

• ERCA (electrocauter): a two-level factor with levels TRUE and FALSE.

• SHNT (shunt): a two-level factor with levels NORMAL and HIGH.

• PVS (pulmonary venous oxygen saturation): a three-level factor with levels LOW, NORMAL and HIGH.

• ACO2 (arterial CO2): a three-level factor with levels LOW, NORMAL and HIGH.

• VALV (pulmonary alveoli ventilation): a four-level factor with levels ZERO, LOW, NORMAL and HIGH.

• VLNG (lung ventilation): a four-level factor with levels ZERO, LOW, NORMAL and HIGH.

• VTUB (ventilation tube): a four-level factor with levels ZERO, LOW, NORMAL and HIGH.

• VMCH (ventilation machine): a four-level factor with levels ZERO, LOW, NORMAL and HIGH.

Note

The complete BN can be downloaded from https://www.bnlearn.com/bnrepository/.

Source

Beinlich I, Suermondt HJ, Chavez RM, Cooper GF (1989). "The ALARM Monitoring System: A Case Study with Two Probabilistic Inference Techniques for Belief Networks." Proceedings of the 2nd European Conference on Artificial Intelligence in Medicine, 247–256.

Examples

# load the data.
data(alarm)
# create and plot the network structure.
modelstring = paste0("[HIST|LVF][CVP|LVV][PCWP|LVV][HYP][LVV|HYP:LVF][LVF]",
  "[STKV|HYP:LVF][ERLO][HRBP|ERLO:HR][HREK|ERCA:HR][ERCA][HRSA|ERCA:HR][ANES]",
  "[APL][TPR|APL][ECO2|ACO2:VLNG][KINK][MINV|INT:VLNG][FIO2][PVS|FIO2:VALV]",
  "[SAO2|PVS:SHNT][PAP|PMB][PMB][SHNT|INT:PMB][INT][PRSS|INT:KINK:VTUB][DISC]",
  "[MVS][VMCH|MVS][VTUB|DISC:VMCH][VLNG|INT:KINK:VTUB][VALV|INT:VLNG]",
  "[ACO2|VALV][CCHL|ACO2:ANES:SAO2:TPR][HR|CCHL][CO|HR:STKV][BP|CO:TPR]")
dag = model2network(modelstring)
## Not run: graphviz.plot(dag)

Description

Estimate the optimal imaginary sample size for the BDe score, given a network structure and a data set, assuming a uniform prior.

Usage

alpha.star(x, data, debug = FALSE)

Arguments

Value

alpha.star() returns a positive number, the estimated optimal imaginary sample size value.

Author(s)

Marco Scutari

References

Steck H (2008). "Learning the Bayesian Network Structure: Dirichlet Prior versus Data." Proceedings of the 24th Conference on Uncertainty in Artificial Intelligence, 511–518.

Examples

data(learning.test)
dag = hc(learning.test, score = "bic")

for (i in 1:3) {

  a = alpha.star(dag, learning.test)
  dag = hc(learning.test, score = "bde", iss = a)

}#FOR

Description

Drop, add or set the direction of a directed or undirected arc (that is, an edge).

Usage

# arc operations.
set.arc(x, from, to, check.cycles = TRUE, check.illegal = TRUE, debug = FALSE)
drop.arc(x, from, to, debug = FALSE)
reverse.arc(x, from, to, check.cycles = TRUE, check.illegal = TRUE, debug = FALSE)

# edge (that is, undirected arc) operations
set.edge(x, from, to, check.cycles = TRUE, check.illegal = TRUE, debug = FALSE)
drop.edge(x, from, to, debug = FALSE)

Arguments

Details

The set.arc() function operates in the following way:

• If there is no arc between from and to, the arc from $\rightarrow$ to is added.

• If there is an undirected arc between from and to, its direction is set to from $\rightarrow$ to.

• If the arc to $\rightarrow$ from is present, it is reversed.

• If the arc from $\rightarrow$ to is present, no action is taken.

The drop.arc() function operates in the following way:

• If there is no arc between from and to, no action is taken.

• If there is a directed or an undirected arc between from and to, it is dropped regardless of its direction.

The reverse.arc() function operates in the following way:

• If there is no arc between from and to, it returns an error.

• If there is an undirected arc between from and to, it returns an error.

• If the arc to $\rightarrow$ from is present, it is reversed.

• If the arc from $\rightarrow$ to is present, it is reversed.

The set.edge() function operates in the following way:

• If there is no arc between from and to, the undirected arc from - to is added.

• If there is an undirected arc between from and to, no action is taken.

• If either the arc from $\rightarrow$ to or the arc to $\rightarrow$ from are present, they are replaced with the undirected arc from - to.

The drop.edge() function operates in the following way:

• If there is no undirected arc between from and to, no action is taken.

• If there is an undirected arc between from and to, it is removed.

• If there is a directed arc between from and to, no action is taken.

Value

All functions return an updated copy of x invisibly.

Author(s)

Marco Scutari

Examples

dag = cpdag(model2network("[A][C][F][B|A][D|A:C][E|B:F]"))
dag

## use debug = TRUE to get more information.
updated = set.arc(dag, "A", "B")
updated
updated = drop.arc(dag, "A", "B")
updated
updated = drop.edge(dag, "A", "B")
updated
updated = reverse.arc(dag, "A", "D")
updated

Description

Measure the strength of the probabilistic relationships expressed by the arcs of a Bayesian network, and use model averaging to build a network containing only the significant arcs.

Usage

# strength of the arcs present in x.
arc.strength(x, data, criterion = NULL, ..., debug = FALSE)
# strength of all possible arcs, as learned from bootstrapped data.
boot.strength(data, cluster, R = 200, m = nrow(data),
  algorithm, algorithm.args = list(), cpdag = TRUE, shuffle = TRUE,
  debug = FALSE)
# strength of all possible arcs, from a list of custom networks.
custom.strength(networks, nodes, weights = NULL, cpdag = TRUE, debug = FALSE)
# strength of all possible arcs, computed using Bayes factors.
bf.strength(x, data, score, ..., debug = FALSE)

# average arc strengths.
## S3 method for class 'bn.strength'
mean(x, ..., weights = NULL)

# averaged network structure.
averaged.network(strength, threshold)
# strength threshold for inclusion in the averaged network structure.
inclusion.threshold(strength)

Arguments

For bf.strength(), the label of the score used to compute the Bayes factors; see BF for details.

Details

arc.strength() computes a measure of confidence or strength for each arc, while keeping fixed the rest of the network structure.

If criterion is a conditional independence test, the strength is a p-value (so the lower the value, the stronger the relationship). The conditional independence test would be that used to drop the arc from the network. The only two possible additional arguments are alpha, which sets the significance threshold that is used in strength.plot(); and B, the number of permutations to be generated for each permutation test.

If criterion is the label of a score function, the strength is measured by the score gain/loss which would be caused by the arc's removal. In other words, it is the difference between the score of the network in which the arc is not present and the score of the network in which the arc is present. Negative values correspond to decreases in the network score, and positive values correspond to increases in the network score (the stronger the relationship, the more negative the difference). There may be additional arguments depending on choen score, see score for details. The significance threshold is set to 0.

boot.strength() estimates the strength of each arc as its empirical frequency over a set of networks learned from bootstrap samples. It computes the probability of each arc (modulo its direction) and the probabilities of each arc's directions conditional on the arc being present in the graph (in either direction). The significance threshold is computed automatically from the strength estimates.

bf.strength() estimates the strength of each arc using Bayes factors to overcome the fact that Bayesian posterior scores are not normalised, and uses the latter to estimate the probabilities of all possible states of an arc given the rest of the network. The significance threshold is set to 1.

custom.strength() takes a list of networks and estimates arc strength in the same way as
boot.strength().

Model averaging is supported for objects of class bn.strength returned by boot.strength, custom.strength and bf.strength. The returned network contains the arcs whose strength is greater than the threshold attribute of the bn.strength object passed to averaged.network().

Value

arc.strength(), boot.strength(), custom.strength(), bf.strength() and mean() return an object of class bn.strength; boot.strength() and custom.strength() also include information about the relative probabilities of arc directions.

averaged.network() returns an object of class bn.

See bn.strength class and bn-class for details.

Note

averaged.network() typically returns a completely directed graph; an arc can be undirected if and only if the probability of each of its directions is exactly 0.5. This may happen, for example, if the arc is undirected in all the networks being averaged.

Author(s)

Marco Scutari

References

for model averaging and bootstrap strength (confidence):

Friedman N, Goldszmidt M, Wyner A (1999). "Data Analysis with Bayesian Networks: A Bootstrap Approach." Proceedings of the 15th Annual Conference on Uncertainty in Artificial Intelligence, 196–201.

for the computation of the bootstrap strength (confidence) significance threshold:

Scutari M, Nagarajan R (2013). "On Identifying Significant Edges in Graphical Models of Molecular Networks." Artificial Intelligence in Medicine, 57(3):207–217.

See Also

strength.plot, score, ci.test.

Examples

data(learning.test)
dag = hc(learning.test)
arc.strength(dag, learning.test)

## Not run: 
arcs = boot.strength(learning.test, algorithm = "hc")
arcs[(arcs$strength > 0.85) & (arcs$direction >= 0.5), ]
averaged.network(arcs)

start = random.graph(nodes = names(learning.test), num = 50)
netlist = lapply(start, function(net) {
  hc(learning.test, score = "bde", iss = 10, start = net) })
arcs = custom.strength(netlist, nodes = names(learning.test),
         cpdag = FALSE)
arcs[(arcs$strength > 0.85) & (arcs$direction >= 0.5), ]
modelstring(averaged.network(arcs))

bf.strength(dag, learning.test, score = "bds", prior = "marginal")

## End(Not run)

Description

Synthetic data set on lung diseases and Asia travel from Lauritzen and Spiegelhalter (1988).

Usage

data(asia)

Format

The asia data set contains the following variables:

• D (dyspnoea), a two-level factor with levels yes and no.

• T (tuberculosis), a two-level factor with levels yes and no.

• L (lung cancer), a two-level factor with levels yes and no.

• B (bronchitis), a two-level factor with levels yes and no.

• A (visit to Asia), a two-level factor with levels yes and no.

• S (smoking), a two-level factor with levels yes and no.

• X (chest X-ray), a two-level factor with levels yes and no.

• E (tuberculosis versus lung cancer/bronchitis), a two-level factor with levels yes and no.

Note

Lauritzen and Spiegelhalter (1988) motivate this example as follows:

“Shortness-of-breath (dyspnoea) may be due to tuberculosis, lung cancer or bronchitis, or none of them, or more than one of them. A recent visit to Asia increases the chances of tuberculosis, while smoking is known to be a risk factor for both lung cancer and bronchitis. The results of a single chest X-ray do not discriminate between lung cancer and tuberculosis, as neither does the presence or absence of dyspnoea.”

Standard learning algorithms are not able to recover the true structure of the network because of the presence of a node (E) with conditional probabilities equal to both 0 and 1. Monte Carlo tests seem to behave better than their parametric counterparts.

The complete BN can be downloaded from https://www.bnlearn.com/bnrepository/.

Source

Lauritzen S, Spiegelhalter D (1988). "Local Computation with Probabilities on Graphical Structures and their Application to Expert Systems (with discussion)." Journal of the Royal Statistical Society: Series B, 50(2):157–224.

Examples

# load the data.
data(asia)
# create and plot the network structure.
dag = model2network("[A][S][T|A][L|S][B|S][D|B:E][E|T:L][X|E]")
## Not run: graphviz.plot(dag)

Description

Compute the Bayes factor between the structures of two Bayesian networks.

Usage

BF(num, den, data, score, ..., log = TRUE)

Arguments

Value

A single numeric value, the Bayes factor of the two network structures num and den.

Note

The Bayes factor for two network structures, by definition, is the ratio of the respective marginal likelihoods. If we assume the uniform prior over all possible DAGs, it is equivalent to the ratio of the corresponding posterior probabilities. However, note that it is possible to specify different priors using the “...” arguments of BF(). In that case, the value returned by the function will not be the standard Bayes factor.

Author(s)

Marco Scutari

See Also

score, compare, bf.strength.

Examples

data(learning.test)

dag1 = model2network("[A][B][F][C|B][E|B][D|A:B:C]")
dag2 = model2network("[A][C][B|A][D|A][E|D][F|A:C:E]")
BF(dag1, dag2, learning.test, score = "bds", iss = 1)

Description

The structure of an object of S3 class bn.

Details

An object of class bn is a list containing at least the following components:

• learning: a list containing some information about the results of the learning algorithm. It is never changed afterwards.

  • whitelist: a copy of the whitelist argument (a two-column matrix, whose columns are labelled from and to) as transformed by sanitisation functions.

  • blacklist: a copy of the blacklist argument (a two-column matrix, whose columns are labelled from and to) as transformed by sanitisation functions.

  • test: the label of the conditional independence test used by the learning algorithm (a character string); the label of the network score is used for score-based algorithms; the label of the network score used in the “Maximize” phase of hybrid algorithms; "none" for randomly generated graphs. For hybrid algorithms, test always has the same value as maxscore (see below).

  • ntests: the number of conditional independence tests or score comparisons used in the learning (an integer value).

  • algo: the label of the learning algorithm or the random generation algorithm used to generate the network (a character string).

  • args: a list. The values of the parameters of either the conditional tests or the scores used in the learning process. Only the relevant ones are stored, so this may be an empty list.

    • alpha: the target nominal type I error rate (a numeric value) of the conditional independence tests.

    • iss: a positive numeric value, the imaginary sample size used by the "bge" and "bde" scores.

    • k: a positive numeric value, the penalty coefficient used by the "aic", "aic-g", "aic-cg", "bic", "bic-g", "bic-cg", "pnal", "pnal-g" and "pnal-cg" scores.

    • prob: the probability of each arc to be present in a graph generated by the "ordered" graph generation algorithm.

    • burn.in: the number of iterations for the "ic-dag" graph generation algorithm to converge to a stationary (and uniform) probability distribution.

    • max.degree: the maximum degree for any node in a graph generated by the "ic-dag" graph generation algorithm.

    • max.in.degree: the maximum in-degree for any node in a graph generated by the "ic-dag" graph generation algorithm.

    • max.out.degree: the maximum out-degree for any node in a graph generated by the "ic-dag" graph generation algorithm.

    • training: a character string, the label of the training node in a Bayesian network classifier.

    • threshold: the threshold used to determine which arcs are significant when averaging network structures.

    • prior: the graphical prior used in combination with a Bayesian score such as "bde" or "bge".

    • beta: the parameters of the graphical prior.

• nodes: a list. Each element is named after a node and contains the following elements:

  • mb: the Markov blanket of the node (a vector of character strings).

  • nbr: the neighbourhood of the node (a vector of character strings).

  • parents: the parents of the node (a vector of character strings).

  • children: the children of the node (a vector of character strings).

• arcs: the arcs of the Bayesian network (a two-column matrix, whose columns are labelled from and to). Undirected arcs are stored as two directed arcs with opposite directions between the corresponding incident nodes.

Additional (optional) components under learning:

• optimized: whether additional optimisations have been used in the learning algorithm (a boolean value).

• illegal: arcs that are illegal according to the parametric assumptions used to learn the network structure (a two-column matrix, whose columns are labelled from and to).

• restrict: the label of the constraint-based algorithm used in the “Restrict” phase of a hybrid learning algorithm (a character string).

• rtest: the label of the conditional independence test used in the “Restrict” phase of a hybrid learning algorithm (a character string).

• maximize: the label of the score-based algorithm used in the “Maximize” phase of a hybrid learning algorithm (a character string).

• maxscore: the label of the network score used in the “Maximize” phase of a hybrid learning algorithm (a character string).

• max.sx: the maximum allowed size of the conditioning sets in the conditional independence tests used in constraint-based algorithms.

An object of class bn may also have auxiliary class labels, depending on how it has been created. In particular:

• bn.naive, identifying naive Bayes classifiers from naive.bayes();

• bn.tan, identifying tree-augmented naive Bayes (TAN) classifiers from tree.bayes();

• bn.twin, identifying twin networks from twin();

• bn.ctf, identifying counterfactual networks from twin().

Author(s)

Marco Scutari

Description

Apply a user-specified function to the Bayesian network structures learned from bootstrap samples of the original data.

Usage

bn.boot(data, statistic, R = 200, m = nrow(data), algorithm,
  algorithm.args = list(), statistic.args = list(), cluster,
  debug = FALSE)

Arguments

Details

The first argument of statistic is the bn object encoding the network structure learned from the bootstrap sample. The arguments specified in statistic.args are extracted from the list and passed to statistic as subsequent arguments.

Value

A list containing the results of the calls to statistic.

Author(s)

Marco Scutari

References

Friedman N, Goldszmidt M, Wyner A (1999). "Data Analysis with Bayesian Networks: A Bootstrap Approach." Proceedings of the 15th Annual Conference on Uncertainty in Artificial Intelligence, 196–201.

See Also

bn.cv, rbn.

Examples

## Not run: 
data(learning.test)
bn.boot(data = learning.test, R = 2, m = 500, algorithm = "gs",
  statistic = arcs)

## End(Not run)

Description

Perform a k-fold or hold-out cross-validation for a learning algorithm or a fixed network structure.

Usage

bn.cv(data, bn, loss = NULL, ..., algorithm.args = list(),
  loss.args = list(), fit, fit.args = list(), method = "k-fold",
  cluster, debug = FALSE)

## S3 method for class 'bn.kcv'
plot(x, ..., main, xlab, ylab, connect = FALSE)
## S3 method for class 'bn.kcv.list'
plot(x, ..., main, xlab, ylab, connect = FALSE)

loss(x)

Arguments

Value

bn.cv() returns an object of class bn.kcv.list if runs is at least 2, an object of class bn.kcv if runs is equal to 1.

loss() returns a numeric vector with a length equal to runs.

Cross-Validation Strategies

The following cross-validation methods are implemented:

• k-fold: the data are split in k subsets of equal size. For each subset in turn, bn is fitted (and possibly learned as well) on the other k - 1 subsets. The loss function is then computed using that subset. Loss estimates for each of the k subsets are then combined to give an overall loss for data.

• custom-folds: the data are manually partitioned by the user into subsets, which are then used as in k-fold cross-validation. Subsets are not constrained to have the same size, and every observation must be assigned to one subset.

• hold-out: k subsamples of size m are sampled independently without replacement from the data. For each subsample, bn is fitted (and possibly learned) on the remaining m - nrow(data) samples and the loss function is computed on the m observations in the subsample. The overall loss estimate is the average of the k loss estimates from the subsamples.

If cross-validation is used with multiple runs, the overall loss is the average of the loss estimates from the different runs.

To clarify, cross-validation methods accept the following optional arguments:

• k: a positive integer number, the number of groups into which the data will be split (in k-fold cross-validation) or the number of times the data will be split in training and test samples (in hold-out cross-validation).

• m: a positive integer number, the size of the test set in hold-out cross-validation.

• runs: a positive integer number, the number of times k-fold or hold-out cross-validation will be run.

• folds: a list in which each element corresponds to one fold and contains the indices for the observations that are included in that fold; or a list with an element for each run, in which each element is itself a list of the folds to be used for that run.

Loss Functions

The following loss functions are implemented:

• Log-Likelihood Loss ("logl"): also known as negative entropy or negentropy, it is the negated expected log-likelihood of the test set for the Bayesian network fitted from the training set. Lower values are better.

• Classification Error ("pred"): the prediction error for a single discrete node. Lower values are better.

• Exact Classification Error ("pred-exact"): closed-form exact posterior predictions are available for Bayesian network classifiers. Lower values are better.

• Predictive Correlation ("cor"): the correlation between the observed and the predicted values for a single continuous node. Higher values are better.

• Mean Squared Error ("mse"): the mean squared error between the observed and the predicted values for a single continuous node. Lower values are better.

• F1 score ("f1"): the F1 score between observed and predicted values for both binary and multiclass target variables.

• AUROC ("auroc"): the area under the ROC curve for both binary and multiclass target variables. The multiclass AUROC score is computed one-vs-rest by averaging the AUROC for each level of the target variable.

Optional arguments that can be specified in loss.args are:

• predict: a character string, the label of the method used to predict the observations in the test set. The default is "parents". Other possible values are the same as in predict().

• predict.args: a list containing the optional arguments for the prediction method. See the documentation for predict() for more details.

• target: a character string, the label of the target node for prediction in all loss functions but "logl".

Plotting Results from Cross-Validation

Both plot methods accept any combination of objects of class bn.kcv or bn.kcv.list (the first as the x argument, the remaining as the ... argument) and plot the respective expected loss values side by side. For a bn.kcv object, this mean a single point; for a bn.kcv.list object, this means a boxplot.

Author(s)

Marco Scutari

References

Koller D, Friedman N (2009). Probabilistic Graphical Models: Principles and Techniques. MIT Press.

See Also

bn.boot, rbn, bn.kcv-class.

Examples

bn.cv(learning.test, 'hc', loss = "pred",
  loss.args = list(predict = "bayes-lw", target = "F"))

folds = list(1:2000, 2001:3000, 3001:5000)
bn.cv(learning.test, 'hc', loss = "logl", method = "custom-folds",
  folds = folds)

xval = bn.cv(gaussian.test, 'mmhc', method = "hold-out",
         k = 5, m = 50, runs = 2)
xval
loss(xval)

## Not run: 
# comparing algorithms with multiple runs of cross-validation.
gaussian.subset = gaussian.test[1:50, ]
cv.gs = bn.cv(gaussian.subset, 'gs', runs = 10)
cv.iamb = bn.cv(gaussian.subset, 'iamb', runs = 10)
cv.inter = bn.cv(gaussian.subset, 'inter.iamb', runs = 10)
plot(cv.gs, cv.iamb, cv.inter,
  xlab = c("Grow-Shrink", "IAMB", "Inter-IAMB"), connect = TRUE)

# use custom folds.
folds = split(sample(nrow(gaussian.subset)), seq(5))
bn.cv(gaussian.subset, "hc", method = "custom-folds", folds = folds)

# multiple runs, with custom folds.
folds = replicate(5, split(sample(nrow(gaussian.subset)), seq(5)),
          simplify = FALSE)
bn.cv(gaussian.subset, "hc", method = "custom-folds", folds = folds)

## End(Not run)

Description

Fit, assign or replace the parameters of a Bayesian network conditional on its structure.

Usage

bn.fit(x, data, cluster, method, ..., keep.fitted = TRUE,
  debug = FALSE)
custom.fit(x, dist, ordinal, debug = FALSE)
bn.net(x)

Arguments

Details

bn.fit() fits the parameters of a Bayesian network given its structure and a data set; bn.net returns the structure underlying a fitted Bayesian network.

bn.fit() accepts data with missing values encoded as NA. If the parameter estimation method was not specifically designed to deal with incomplete data, bn.fit() uses locally complete observations to fit the parameters of each local distribution.

Available methods for discrete Bayesian networks are:

• "mle": the maximum likelihood estimator for conditional probabilities.

• "bayes": the classic Bayesian posterior estimator with a uniform prior matching that in the Bayesian Dirichlet equivalent ("bde") score.

• "hdir": the hierarchical Dirichlet posterior estimator for related data sets from Azzimonti, Corani and Zaffalon (2019).

• "hard-em": the Expectation-Maximimation implementation of the estimators above.

Available methods for Gaussian Bayesian networks are:

• "mle-g": the maximum likelihood estimator for least squares regression models.

• "hard-em-g": the Expectation-Maximisation implementation of the estimators above.

Available methods for conditional Gaussian Bayesian networks are:

• "mle-cg": a combination of the maximum likelihood estimators "mle" and "mle-g".

• "hard-em-cg": the Expectation-Maximisation implementation of the estimators above.

Available methods for zero-inflated Bayesian networks are:

• "mle-zihp": the maximum likelihood estimator for zero-inflated hyper-Poisson regression models, computed by expectation-maximization with the generalised linear model components fitted by iteratively reweighted least squares.

• "mle-zinb": the maximum likelihood estimator for zero-inflated negative binomial regression models, computed by expectation-maximization with the generalised linear model components fitted by iteratively reweighted least squares.

Additional arguments for the bn.fit() function:

• iss: a numeric value, the imaginary sample size used by the "bayes" method to estimate the conditional probability tables associated with discrete nodes (see score for details).

• replace.unidentifiable: a boolean value. If TRUE and method is a maximum likelihood estimator, unidentifiable parameters are replaced by zeroes (in the case of regression coefficients and standard errors in Gaussian and conditional Gaussian nodes) or by uniform conditional probabilities (in discrete nodes).

  If FALSE (the default), the conditional probabilities in the local distributions of discrete nodes have a maximum likelihood estimate of NaN for all parent configurations that are not observed in data. Similarly, regression coefficients are set to NA if the linear regressions correspoinding to the local distributions of continuous nodes are singular. Such missing values propagate to the results of functions such as predict().

• alpha0: a positive number, the amount of information pooling between the related data sets in the "hdir" estimator.

• group: a character string, the label of the node with the grouping of the observations into the related data sets in the "hdir" estimator.

• impute and impute.args: a character string, the label of the imputation method (and its arguments) used by "hard-em", "hard-em-g" and "hard-em-cg" to complete the data in the expectation step. The default method is the same as for impute().

• fit and fit.args: a character string, the label of the parameter estimation method used by "hard-em", "hard-em-g" and "hard-em-cg" to estimate the parameters in the maximisation step. The default method is the same as for bn.fit().

• loglik.threshold: a non-negative numeric value, the minimum improvement threshold to continue iterating in "hard-em", "hard-em-g" and "hard-em-cg". The threshold is defined as the relative likelihood improvement divided by the sample size of data, and defaults to 1e-3. Setting it to zero means that iterations stop only when the improvement is negative, which can happen due to stochastic noise if the imputation of missing data uses approximate inference.

• params.threshold: a non-negative numeric value, the minimum maximum relative change in the parameter values to continue iterating in "hard-em", "hard-em-g" and "hard-em-cg". The threshold is defined as the maximum difference between parameter values, divided by the parameter value in the latest iteration. The default value is 1e-3.

• max.iter: a positive integer value, the maximum number of iterations in "hard-em", "hard-em-g" and "hard-em-cg". The default value is 5.

• start: a bn.fit object, the fitted network used to initialise the "hard-em", "hard-em-g" and "hard-em-cg" estimators. The default is to use the bn.fit object obtained from x with the default parameter estimator for the data, which will use locally complete data to fit the local distributions.

• newdata: a data frame, a separate set of data used to assess the convergence of the "hard-em", "hard-em-g" and "hard-em-cg" estimators. The data in data are used by default for this purpose.

• em.max.iter: a positive integer value, the maximum number of expectation-maximization iterations in the "mle-zihp" and "mle-zinb" estimators. The default value is 100.

• em.tol: a non-negative numeric value, the relative convergence tolerance (on both the log-likelihood and the parameters) of the expectation-maximization iteration in the "mle-zihp" and "mle-zinb" estimators. The default value is 1e-8.

• m.step: a character string, either "full" (the default) to re-fit each generalised linear model component to convergence in every maximisation step, or "one-step" to take a single iteratively reweighted least squares step per iteration (a cheaper generalized-EM update) in the "mle-zihp" and "mle-zinb" estimators.

An in-place replacement method is available to change the parameters of each node in a bn.fit object; see the examples for discrete, continuous and conditional Gaussian networks below.

• For a discrete node (class bn.fit.dnode or bn.fit.onode), the new parameters must be in a table object.

• For a Gaussian node (class bn.fit.gnode), the new parameters can be defined either by an lm, glm or pensim object (the latter is from the penalized package) or in a list with elements named coef, sd and optionally fitted and resid.

• For a conditional Gaussian node (class bn.fit.cgnode), the new parameters can be defined by a list with elements named coef, sd and optionally fitted, resid and configs. In both cases coef should contain the new regression coefficients, sd the standard deviation of the residuals, fitted the fitted values and resid the residuals. configs should contain the configurations of the discrete parents of the conditional Gaussian node, stored as a factor.

• For a zero-inflated hyper-Poisson node (class bn.fit.zihpnode), the new parameters can be defined by a list with elements named inflation, intensity, dispersion and optionally fitted and resid.

• For a zero-inflated negative binomial node (class bn.fit.zinbnode), the new parameters can be defined by a list with elements named inflation, intensity, dispersion and optionally fitted and resid.

custom.fit() takes a set of user-specified distributions and their parameters and uses them to build a bn.fit object. Its purpose is to specify a Bayesian network (complete with parameters, not just the structure) using expert knowledge rather than learning it from a data set. The distributions must be passed to the function in a list, with elements named after the nodes of the network structure x. Each element of the list must be in one of the formats described above for in-place replacement.

Value

bn.fit() and custom.fit()returns an object of class bn.fit, bn.net() an object of class bn. See bn class and bn.fit class for details.

Note

Due to the way Bayesian networks are defined, it is possible to estimate their parameters only if the network structure is completely directed (that is, there are no undirected arcs). See set.arc and cextend for two ways of manually setting the direction of one or more arcs.

In the case of maximum likelihood estimators, bn.fit() may produce NA parameter estimates. For instance, for discrete and conditional Gaussian nodes, this happens when there are (discrete) parent configurations that are not observed in data. To avoid this, either set replace.unidentifiable to TRUE or, in the case of discrete networks, use method = "bayes".

Author(s)

Marco Scutari

References

Azzimonti L, Corani G, Zaffalon M (2019). "Hierarchical Estimation of Parameters in Bayesian Networks." Computational Statistics & Data Analysis, 137:67–91.

See Also

bn.fit utilities, bn.fit plots.

Examples

data(learning.test)

# learn the network structure.
cpdag = pc.stable(learning.test)
# set the direction of the only undirected arc, A - B.
dag = set.arc(cpdag, "A", "B")
# estimate the parameters of the Bayesian network.
fitted = bn.fit(dag, learning.test)
# replace the parameters of node B.
new.cpt = matrix(c(0.1, 0.2, 0.3, 0.2, 0.5, 0.6, 0.7, 0.3, 0.1),
            byrow = TRUE, ncol = 3,
            dimnames = list(B = c("a", "b", "c"), A = c("a", "b", "c")))
fitted$B = as.table(new.cpt)
# the network structure is still the same.
all.equal(dag, bn.net(fitted))

# learn the network structure.
dag = hc(gaussian.test)
# estimate the parameters of the Bayesian network.
fitted = bn.fit(dag, gaussian.test)
# replace the parameters of the node F.
fitted$F = list(coef = c(1, 2, 3, 4, 5), sd = 3)
# set the original parameters again.
fitted$F = lm(F ~ A + D + E + G, data = gaussian.test)

# discrete Bayesian network from expert knowledge.
dag = model2network("[A][B][C|A:B]")
cptA = matrix(c(0.4, 0.6), ncol = 2, dimnames = list(NULL, c("LOW", "HIGH")))
cptB = matrix(c(0.8, 0.2), ncol = 2, dimnames = list(NULL, c("GOOD", "BAD")))
cptC = c(0.5, 0.5, 0.4, 0.6, 0.3, 0.7, 0.2, 0.8)
dim(cptC) = c(2, 2, 2)
dimnames(cptC) = list("C" = c("TRUE", "FALSE"), "A" =  c("LOW", "HIGH"),
                   "B" = c("GOOD", "BAD"))
cfit = custom.fit(dag, dist = list(A = cptA, B = cptB, C = cptC))
# for ordinal nodes, it is nearly the same.
cfit = custom.fit(dag, dist = list(A = cptA, B = cptB, C = cptC),
         ordinal = c("A", "B"))

# Gaussian Bayesian network from expert knowledge.
distA = list(coef = c("(Intercept)" = 2), sd = 1)
distB = list(coef = c("(Intercept)" = 1), sd = 1.5)
distC = list(coef = c("(Intercept)" = 0.5, "A" = 0.75, "B" = 1.32), sd = 0.4)
cfit = custom.fit(dag, dist = list(A = distA, B = distB, C = distC))

# conditional Gaussian Bayesian network from expert knowledge.
cptA = matrix(c(0.4, 0.6), ncol = 2, dimnames = list(NULL, c("LOW", "HIGH")))
distB = list(coef = c("(Intercept)" = 1), sd = 1.5)
distC = list(coef = matrix(c(1.2, 2.3, 3.4, 4.5), ncol = 2,
               dimnames = list(c("(Intercept)", "B"), NULL)),
          sd = c(0.3, 0.6))
cgfit = custom.fit(dag, dist = list(A = cptA, B = distB, C = distC))

# zero-inflated hyper-Poisson nodes from expert knowledge.
dag = model2network("[A][B|A][C|A:B]")
ldistA = list(inflation = c("(Intercept)" = 0),
              intensity = c("(Intercept)" = 0),
              dispersion = 1)
ldistB = list(inflation = c("(Intercept)" = 0, A = 1),
              intensity = c("(Intercept)" = 0, A = 1),
              dispersion = 0.5)
ldistC = list(inflation = c("(Intercept)" = 0, A = 0.5, B = 0.5),
              intensity = c("(Intercept)" = 0, A = 0.4, B = 0.6),
              dispersion = 0.25)
zifit = custom.fit(dag, list(A = ldistA, B = ldistB, C = ldistC))

# zero-inflated negative binomial nodes from expert knowledge.
ldistA = list(inflation = c("(Intercept)" = 0),
              prsucc = c("(Intercept)" = 0),
              failures = 1)
ldistB = list(inflation = c("(Intercept)" = 0, A = 1),
              prsucc = c("(Intercept)" = 0, A = 1),
              failures = 0.5)
ldistC = list(inflation = c("(Intercept)" = 0, A = 0.5, B = 0.5),
              prsucc = c("(Intercept)" = 0, A = 0.4, B = 0.6),
              failures = 0.25)
zifit = custom.fit(dag, list(A = ldistA, B = ldistB, C = ldistC))

Description

The structure of an object of S3 class bn.fit.

Details

An object of class bn.fit is a list whose elements correspond to the nodes of the Bayesian network. If the network is discrete (the nodes are multinomial random variables), the object also has class bn.fit.dnet. Each node has class bn.fit.dnode and contains the following elements:

• node: a character string, the label of the node.

• parents: a vector of character strings, the labels of the parents of the node.

• children: a vector of character strings, the labels of the children of the node.

• prob: a (multi)dimensional numeric table, the conditional probability table of the node given its parents.

Nodes encoding ordinal variables (ordered factors) have class bn.fit.onode and contain the same elements as bn.fit.dnode nodes. Networks with only ordinal nodes also have class bn.fit.onet, while those containing both ordinal and multinomial nodes also have class bn.fit.donet.

If, on the other hand, the network is continuous (that is, the nodes are Gaussian random variables), the object also has class bn.fit.gnet; each node has class bn.fit.gnode and contains the following elements:

• node: a character string, the label of the node.

• parents: a vector of character strings, the labels of the parents of the node.

• children: a vector of character strings, the labels of the children of the node.

• coefficients: a numeric vector, the linear regression coefficients of the parents against the node.

• residuals: a numeric vector, the residuals of the linear regression.

• fitted.values: a numeric vector, the fitted mean values of the linear regression.

• sd: a numeric value, the standard deviation of the residuals (that is, the standard error).

Hybrid (conditional linear Gaussian) networks also have class bn.fit.gnet. Gaussian nodes have class bn.fit.gnode, discrete nodes have class bn.fit.dnode and conditional Gaussian nodes have class bn.fit.cgnode. Each node contains the following elements:

• node: a character string, the label of the node.

• parents: a vector of character strings, the labels of the parents of the node.

• children: a vector of character strings, the labels of the children of the node.

• dparents: an integer vector, the indexes of the discrete parents in parents.

• gparents: an integer vector, the indexes of the continuous parents in parents.

• dlevels: a list containing the levels of the discrete parents in parents.

• coefficients: a numeric matrix, the linear regression coefficients of the continuous parents. Each column corresponds to a configuration of the discrete parents.

• residuals: a numeric vector, the residuals of the linear regression.

• fitted.values: a numeric vector, the fitted mean values of the linear regression.

• configs: an integer vector, the indexes of the configurations of the discrete parents.

• sd: a numeric vector, the standard deviation of the residuals (the standard error) for each configuration of the discrete parents.

Zero-Inflated (non-integer count) networks also have class bn.fit.zihp (hyper-Poisson) or bn.fit.zinb (negative binomial). Zero-inflated hyper-Poisson nodes have class bn.fit.zihpnode and contain the following elements:

• node: a character string, the label of the node.

• parents: a vector of character strings, the labels of the parents of the node.

• children: a vector of character strings, the labels of the children of the node.

• inflation: a numeric vector, the coefficients of the zero-inflation (logistic) regression model associated with the parents.

• intensity: a numeric vector, the coefficients of the (log-linear) regression model for the Poisson intensity parameter associated with the parents.

• dispersion: a numeric value, the dispersion parameter that augments the Poisson (pseudo)counts into a hyper-Poisson distribution.

Zero-inflated negative binomial nodes have class bn.fit.zinbnode and contain the following elements:

• node: a character string, the label of the node.

• parents: a vector of character strings, the labels of the parents of the node.

• children: a vector of character strings, the labels of the children of the node.

• inflation: a numeric vector, the coefficients of the zero-inflation (logistic) regression model associated with the parents.

• prsucc: a numeric vector, the coefficients of the (logistic) regression model for the probability of success parameter associated with the parents.

• failures: a numeric value, the number of failures.

In addition, Bayesian network classifiers store the label of the training node in an attribute named training.

An object of class bn.fit may also have auxiliary class labels, depending on how it has been created. In particular:

• bn.naive, identifying naive Bayes classifiers from naive.bayes();

• bn.tan, identifying tree-augmented naive Bayes (TAN) classifiers from tree.bayes();

• bn.twin, identifying twin networks from twin();

• bn.ctf, identifying counterfactual networks from twin().

Author(s)

Marco Scutari

Description

Plot functions for the bn.fit, bn.fit.dnode and bn.fit.gnode classes, based on the lattice package.

Usage

## for Gaussian Bayesian networks.
bn.fit.qqplot(fitted, xlab = "Theoretical Quantiles",
  ylab = "Sample Quantiles", main, ...)
bn.fit.histogram(fitted, density = TRUE, xlab = "Residuals",
  ylab = ifelse(density, "Density", ""), main, ...)
bn.fit.xyplot(fitted, xlab = "Fitted values", ylab = "Residuals", main, ...)
## for discrete (multinomial and ordinal) Bayesian networks.
bn.fit.barchart(fitted, xlab = "Probabilities", ylab = "Levels", main, ...)
bn.fit.dotplot(fitted, xlab = "Probabilities", ylab = "Levels", main, ...)

Arguments

Details

bn.fit.qqplot() draws a quantile-quantile plot of the residuals.

bn.fit.histogram() draws a histogram of the residuals, using either absolute or relative frequencies.

bn.fit.xyplot() plots the residuals versus the fitted values.

bn.fit.barchart() and bn.fit.dotplot plot the probabilities in the conditional probability table associated with each node.

Value

The lattice plot objects. Note that if auto-printing is turned off (for example, when the code is loaded with the source function), the return value must be printed explicitly for the plot to be displayed.

Author(s)

Marco Scutari

See Also

bn.fit, bn.fit class.

Description

Assign, extract or compute various quantities of interest from objects of class bn.fit, bn.fit.dnode, bn.fit.gnode, bn.fit.cgnode or bn.fit.onode.

Usage

## methods available for "bn.fit"
## S3 method for class 'bn.fit'
fitted(object, ...)
## S3 method for class 'bn.fit'
coef(object, ...)
## S3 method for class 'bn.fit'
residuals(object, ...)
## S3 method for class 'bn.fit'
sigma(object, ...)
## S3 method for class 'bn.fit'
logLik(object, data, nodes, by.sample = FALSE, na.rm = FALSE, debug = FALSE, ...)
## S3 method for class 'bn.fit'
AIC(object, data, ..., k = 1)
## S3 method for class 'bn.fit'
BIC(object, data, ...)

## non-method functions for "bn.fit"
identifiable(x, by.node = FALSE)
singular(x, by.node = FALSE)

## methods available for "bn.fit.dnode"
## S3 method for class 'bn.fit.dnode'
coef(object, for.parents, ...)

## methods available for "bn.fit.onode"
## S3 method for class 'bn.fit.onode'
coef(object, for.parents, ...)

## methods available for "bn.fit.gnode"
## S3 method for class 'bn.fit.gnode'
fitted(object, ...)
## S3 method for class 'bn.fit.gnode'
coef(object, ...)
## S3 method for class 'bn.fit.gnode'
residuals(object, ...)
## S3 method for class 'bn.fit.gnode'
sigma(object, ...)

## methods available for "bn.fit.cgnode"
## S3 method for class 'bn.fit.cgnode'
fitted(object, ...)
## S3 method for class 'bn.fit.cgnode'
coef(object, for.parents, ...)
## S3 method for class 'bn.fit.cgnode'
residuals(object,  ...)
## S3 method for class 'bn.fit.cgnode'
sigma(object, for.parents, ...)

Arguments

Details

coef() (and its alias coefficients()) extracts model coefficients (that is, conditional probabilities for discrete nodes; linear regression coefficients for Gaussian and conditional Gaussian nodes; regression coefficients of the count and zero-inflation components for zero-inflated nodes).

residuals() (and its alias resid()) extracts model residuals, and fitted() (and its alias
fitted.values()) extracts fitted values from nodes other than discrete ones. If the bn.fit object does not include the residuals or the fitted values for the node of interest, both functions return NULL.

sigma() extracts the standard deviations of the residuals from Gaussian and conditional Gaussian networks and nodes.

logLik() returns the log-likelihood for the observations in data. If na.rm is set to TRUE, the log-likelihood will be NA if the data contain missing values. If na.rm is set to FALSE, missing values will be dropped, and the log-likelihood will be computed using only locally-complete observations (effectively returning the node-average log-likelihood times the sample size). Note that the log-likelihood may be NA even if na.rm = TRUE if the network contains NA parameters or is singular.

The for.parents argument in the methods for coef() and sigma() can be used to have both functions return the parameters associated with a specific configuration of the discrete parents of a node. If for.parents is not specified, all relevant parameters are returned.

Value

logLik() returns a numeric vector or a single numeric value, depending on the value of by.sample. AIC and BIC always return a single numeric value.

All the other functions return a list with an element for each node in the network (if object has class bn.fit) or a numeric vector or matrix (if object has class bn.fit.dnode, bn.fit.gnode, bn.fit.cgnode, bn.fit.onode, bn.fit.zihpnode, bn.fit.zinbnode.).

Author(s)

Marco Scutari

See Also

bn.fit, bn.fit-class.

Examples

data(gaussian.test)
dag = hc(gaussian.test)
fitted = bn.fit(dag, gaussian.test)
coefficients(fitted)
coefficients(fitted$C)
str(residuals(fitted))

data(learning.test)
dag2 = hc(learning.test)
fitted2 = bn.fit(dag2, learning.test)
coefficients(fitted2$E)
coefficients(fitted2$E, for.parents = list(F = "a", B = "b"))

Description

The structure of an object of S3 class bn.kcv or bn.kcv.list.

Details

An object of class bn.kcv.list is a list whose elements are objects of class bn.kcv.

An object of class bn.kcv is a list whose elements correspond to the iterations of a k-fold cross-validation. Each element contains the following objects:

• test: an integer vector, the indexes of the observations used as a test set.

• fitted: an object of class bn.fit, the Bayesian network fitted from the training set.

• learning: the learning element of the bn object that was used for parameter learning from the training set (either learned from the training set as well or specified by the user).

• loss: the value of the loss function (for method = "hold-out" or loss = "logl"), or NA (otherwise).

If the loss function requires predicting values from the test sets, each element also contains:

• predicted: a factor or a numeric vector, the predicted values for the target node for the observations in the test set.

• observed: a factor or a numeric vector, the observed values for the target node for the observations in the test set.

In addition, an object of class bn.kcv has the following attributes:

• loss: a character string, the label of the loss function.

• mean: the mean of the values of the loss function computed across all k iterations of the cross-validation, which is printed as the "expected loss" or averaged to compute the "average loss over the runs".

• bn: either a character string (the label of the learning algorithm to be applied to the training data in each iteration) or an object of class bn (a fixed network structure).

Author(s)

Marco Scutari

Description

The structure of an object of S3 class bn.strength.

Details

An object of class bn.strength is a data frame with the following columns and one row for each arc:

• from, to: the nodes incident on the arc.

• strength: the strength of the arc. See arc.strength, boot.strength, custom.strength and strength.plot for details.

and some additional attributes:

• nodes: a vector of character strings, the labels of the nodes of the network(s) the strength was computed from.

• method: a character string, the method used to compute the strength coefficients. It can be equal to "test", "score" or "bootstrap".

• threshold: a numeric value, the threshold for determining whether a strength coefficient is significant.

An optional column called direction may also be present, giving the probability of the arc's irection given its presence in the graph.

Only the plot() method is defined for this class; therefore, it can be manipulated as a standard data frame.

Author(s)

Marco Scutari

Description

Learn the structure of a (causal) Bayesian network with DirectLiNGAM.

Usage

direct.lingam(x, cluster, whitelist, blacklist, mi, maximize = "alasso",
  maximize.args = list(), debug = FALSE)

Arguments

Details

The "alasso" method has the following optional arguments:

• gamma: the coefficient of the weights used in the adaptive LASSO. The default value is 1.

• lambda.min.ratio and pmax: identical to the arguments of the same names in glmnet(), with the same default values and interpretation.

• k: the coefficient of the penalised likelihood score used to choose the optimal shrinkage in both ridge and LASSO models in the adaptive LASSO. It defaults to the BIC penalty.

Value

An object of class bn. See bn-class for details.

Note

See structure learning for a complete list of structure learning algorithms with the respective references.

Author(s)

Marco Scutari

See Also

local discovery algorithms, score-based algorithms, constraint-based algorithms, hybrid algorithms.

Description

Perform interventions and counterfactuals, using Bayesian networks as causal models.

Usage

# network transformations used in causal inference.
mutilated(x, evidence)
twin(x)

# causal inference types.
intervention(x, evidence)
counterfactual(x, evidence, merging = TRUE)

Arguments

Details

mutilated() constructs the mutilated network arising from an ideal intervention setting the nodes involved to the values specified by evidence.

twin() constructs the twin network required to compute counterfactuals, separating noise (exogenous) from functional (factual) nodes and adding counterfactual nodes.

intervention() is an alias of mutilated().

counterfactual() constructs the twin network and introduces a counterfactual (that is, an intervention on the counterfactual nodes) before applying node merging (optional). The intervention in evidence is specified using the labels of the nodes in x, and is then applied to the corresponding counterfactual nodes.

Both twin() and counterfactual() label nodes as follows:

• A counterfactual node has the label of the corresponding factual node with a period appended (as in paste0(label, ".")).

• A noise node has the label of the corresponding factual node with a "u" prepended (as in paste0("u", label)).

Evidence in intervention() can reference factual and counterfactual nodes. Evidence in counterfactual() can only reference counterfactual nodes; factual nodes will be replaced by the corresponding counterfactual nodes with a warning. Referencing exogenous nodes in either function will produce an error.

Value

A bn, bn.fit or scm object, depending on the class of x.

Objects returned by twin() will additionally have class bn.twin (if x is a bn or bn.fit object) or scm.twin (if x is an scm object).

Objects returned by counterfactual() will additionally have classes bn.twin and bn.ctf (if x is a bn or bn.fit object) or scm.twin (if x is an scm object).

Author(s)

Marco Scutari

References

Balke A, Pearl J (2022). "Probabilistic Evaluation of Counterfactual Queries." Probabilistic and Causal Inference: The Works of Judea Pearl, 237–254. The original is in the Proceedings of AAAI, 1994.

Graham D, Lee CM, Perov Y (2019). "Copy, Paste, Infer: A Robust Analysis of Twin Networks for Counterfactual Inference." NeurIPS19 CausalML workshop.

Koller D, Friedman N (2009). Probabilistic Graphical Models: Principles and Techniques. MIT Press.

Shpitser I, PearlJ (2007). "What counterfactuals can be tested." Proceedings of the 23rd Conference on Uncertainty in Artificial Intelligence (UAI), 352–359.

Examples

# example from Graham, Lee and Perov.
dag = model2network("[C][A|C][B|C][D|A:B]")
intervention(dag, evidence = list(A = 0))
twin(dag)
counterfactual(dag, evidence = list(A. = 0), merging = FALSE)
counterfactual(dag, evidence = list(A. = 0), merging = TRUE)

Description

Perform an independence or a conditional independence test.

Usage

ci.test(x, y, z, data, test, ..., debug = FALSE)

Arguments

Details

Additional arguments of the ci.test() function:

• B: a positive integer, the number of permutations used to compute the p-value of permutation tests. The default value is 5000 for nonparametric permutation tests and 100 for semiparametric permutation tests.

• fun: the function that computes the conditional independence test in the "custom-test" test. fun must have arguments x, y, z, data and args, in this order; in other words, it must have signature function(x, y, z, data, args). x and y will be the labels of the nodes to test for independence (one character string each); z will be the labels of the nodes in the conditioning set (a vector of character strings, possibily NULL for empty sets); data will contain the complete data set, with all the variables (a data frame); and args will be a list containing any additional arguments to the test.

• args: a list containing the optional arguments to fun, for tuning "custom-test" test functions.

Value

An object of class htest containing the following components:

Author(s)

Marco Scutari

See Also

independence tests, arc.strength.

Examples

data(gaussian.test)
data(learning.test)

# using a data frame and column labels.
ci.test(x = "F" , y = "B", z = c("C", "D"), data = gaussian.test)
# using a data frame.
ci.test(gaussian.test)
# using factor objects.
attach(learning.test)
ci.test(x = F , y = B, z = data.frame(C, D))

Description

This is a synthetic data set used as a test case in the bnlearn package.

Usage

data(clgaussian.test)

Format

The clgaussian.test data set contains one normal (Gaussian) variable, 4 discrete variables and 3 conditional Gaussian variables.

Note

The R script to generate data from this network is available from https://www.bnlearn.com/documentation/networks/.

Examples

# load the data.
data(clgaussian.test)
# create and plot the network structure.
dag = model2network("[A][B][C][H][D|A:H][F|B:C][E|B:D][G|A:D:E:F]")
## Not run: graphviz.plot(dag)

Description

Compare two different Bayesian networks; compute their Structural Hamming Distance (SHD) or the Hamming distance between their skeletons. Or graphically compare them by plotting them side by side,

Usage

compare(target, current, arcs = FALSE)
## S3 method for class 'bn'
all.equal(target, current, ...)

sid(learned, true, debug = FALSE)
shd(learned, true, wlbl = FALSE, cpdag = TRUE, debug = FALSE)
hamming(learned, true, debug = FALSE)

graphviz.compare(x, ..., groups, layout = "dot", shape = "rectangle",
  fontsize = 12, main = NULL, sub = NULL, diff = "from-first",
  diff.args = list())

Arguments

Details

graphviz.compare() can visualise differences between graphs in various ways depending on the value of the diff and diff.args arguments:

• "none": differences are not highlighted.

• "from-first": the first bn object, x, is taken as the reference network. All the other networks passed via the ... argument are compared to that first network and their true positive, false positive and false negative arcs are highlighted. Colours, line types and line widths for each category of arcs can be specified as the elements of a list via the diff.args argument, with names tp.col, tp.lty, tp.lwd, fp.col, fp.lty, fp.lwd, fn.col, fn.lty, fn.lwd. In addition, it is possible not to plot the reference network at all by setting show.first to FALSE.

Regardless of the visualisation, the nodes are arranged in the same positions across all networks to make comparison easier.

shd() computes the Structural Hamming Distance (SHD). Note that SHD, as defined in the reference, is defined on CPDAGs; therefore cpdag() is called on both learned and true before computing the distance.

sid() computes the Structural Interventional Distance (SID) between two DAGs, between a CPDAG and a DAG, between a DAG and a CPDAG, or between two CPDAGs.

Value

compare() returns a list containing the number of true positives (tp, the number of arcs in current also present in target), of false positives (fp, the number of arcs in current not present in target) and of false negatives (fn, the number of arcs not in current but present in target) if arcs is FALSE; or the corresponding arc sets if arcs is TRUE.

all.equal() returns either TRUE or a character string describing the differences between target and current.

shd() and hamming() return a non-negative integer number.

sid() returns a vector of non-negative integer numbers, with additional attributes describing the distribution over all the evaluated graphs.

graphviz.compare() plots one or more figures and invisibly returns a list containing the graph objects generated from the networks that were passed as arguments (in the same order). They can be further modified using the graph and Rgraphviz packages. Unlike the other functions, graphviz.compare() accepts networks with different node sets as long as they overlap.

Author(s)

Marco Scutari

References

Tsamardinos I, Brown LE, Aliferis CF (2006). "The Max-Min Hill-Climbing Bayesian Network Structure Learning Algorithm." Machine Learning, 65(1):31–78.

Peters J, Buhlmann P (2015). "Structural Intervention Distance for Evaluating Causal Graphs." Neural Computation, 27(3):771–779.

Examples

data(learning.test)

e1 = model2network("[A][B][C|A:B][D|B][E|C][F|A:E]")
e2 = model2network("[A][B][C|A:B][D|B][E|C:F][F|A]")
shd(e2, e1, debug = TRUE)
sid(e2, e1, debug = TRUE)
unlist(compare(e1,e2))
compare(target = e1, current = e2, arcs = TRUE)
## Not run: graphviz.compare(e1, e2, diff = "none")

Description

Create configurations of discrete variables, which are used in modelling conditional probability tables.

Usage

configs(data, all = TRUE)

Arguments

Value

A factor with one element for each row of data, and levels as specified by all.

Author(s)

Marco Scutari

Examples

data(learning.test)
configs(learning.test, all = TRUE)
configs(learning.test, all = FALSE)

Description

Learn the equivalence class of a directed acyclic graph (DAG) from data using the PC, Grow-Shrink (GS), Incremental Association (IAMB), Fast Incremental Association (Fast-IAMB), Interleaved Incremental Association (Inter-IAMB), Incremental Association with FDR (IAMB-FDR), Max-Min Parents and Children (MMPC), Semi-Interleaved HITON-PC or Hybrid Parents and Children (HPC) constraint-based algorithms.

Usage

pc.stable(x, cluster, whitelist, blacklist, test = NULL,
  alpha = 0.05, ..., max.sx = NULL, debug = FALSE, undirected = FALSE)
gs(x, cluster, whitelist, blacklist, test = NULL,
  alpha = 0.05, ..., max.sx = NULL, debug = FALSE, undirected = FALSE)
iamb(x, cluster, whitelist, blacklist, test = NULL,
  alpha = 0.05, ..., max.sx = NULL, debug = FALSE, undirected = FALSE)
fast.iamb(x, cluster, whitelist, blacklist, test = NULL,
  alpha = 0.05, ..., max.sx = NULL, debug = FALSE, undirected = FALSE)
inter.iamb(x, cluster, whitelist, blacklist, test = NULL,
  alpha = 0.05, ..., max.sx = NULL, debug = FALSE, undirected = FALSE)
iamb.fdr(x, cluster, whitelist, blacklist, test = NULL,
  alpha = 0.05, ..., max.sx = NULL, debug = FALSE, undirected = FALSE)
mmpc(x, cluster, whitelist, blacklist, test = NULL,
  alpha = 0.05, ..., max.sx = NULL, debug = FALSE, undirected = TRUE)
si.hiton.pc(x, cluster, whitelist, blacklist, test = NULL,
  alpha = 0.05, ..., max.sx = NULL, debug = FALSE, undirected = TRUE)
hpc(x, cluster, whitelist, blacklist, test = NULL,
  alpha = 0.05, ..., max.sx = NULL, debug = FALSE, undirected = TRUE)

Arguments

Value

An object of class bn. See bn-class for details.

Note

Note that even when undirected is set to FALSE, there is no guarantee that all arcs in the returned network will be directed; some arc directions are impossible to learn just from data due to score equivalence. cextend() provides a consistent extension of partially directed networks into directed acyclic graphs, which can then be used (for instance) for parameter learning.

See structure learning for a complete list of structure learning algorithms with the respective references. All algorithms accept incomplete data and handle it by computing individual conditional independence tests on locally complete observations.

Author(s)

Marco Scutari

See Also

independence tests, local discovery algorithms, score-based algorithms, hybrid algorithms, causal discovery algorithms, cextend.

Description

Probable risk factors for coronary thrombosis, comprising data from 1841 men.

Usage

data(coronary)

Format

The coronary data set contains the following 6 variables:

• Smoking (smoking): a two-level factor with levels no and yes.

• M. Work (strenuous mental work): a two-level factor with levels no and yes.

• P. Work (strenuous physical work): a two-level factor with levels no and yes.

• Pressure (systolic blood pressure): a two-level factor with levels <140 and >140.

• Proteins (ratio of beta and alpha lipoproteins): a two-level factor with levels <3 and >3.

• Family (family anamnesis of coronary heart disease): a two-level factor with levels neg and pos.

Source

Edwards DI (2000). Introduction to Graphical Modelling. Springer, 2nd edition.

Reinis Z, Pokorny J, Basika V, Tiserova J, Gorican K, Horakova D, Stuchlikova E, Havranek T, Hrabovsky F (1981). "Prognostic Significance of the Risk Profile in the Prevention of Coronary Heart Disease." Bratisl Lek Listy, 76:137–150. Published in Bratislava Medical Journal, in Czech.

Whittaker J (1990). Graphical Models in Applied Multivariate Statistics. Wiley.

Examples

# This is the undirected graphical model from Whittaker (1990).
data(coronary)
ug = empty.graph(names(coronary))
arcs(ug, check.cycles = FALSE) = matrix(
  c("Family", "M. Work", "M. Work", "Family",
    "M. Work", "P. Work", "P. Work", "M. Work",
    "M. Work", "Proteins", "Proteins", "M. Work",
    "M. Work", "Smoking", "Smoking", "M. Work",
    "P. Work", "Smoking", "Smoking", "P. Work",
    "P. Work", "Proteins", "Proteins", "P. Work",
    "Smoking", "Proteins", "Proteins", "Smoking",
    "Smoking", "Pressure", "Pressure", "Smoking",
    "Pressure", "Proteins", "Proteins", "Pressure"),
  ncol = 2, byrow = TRUE,
  dimnames = list(c(), c("from", "to")))
## Not run: graphviz.plot(ug, shape = "ellipse")

Description

Find the equivalence class and the v-structures of a Bayesian network, construct its moral graph, or create a consistent extension of an equivalent class.

Usage

cpdag(x, wlbl = FALSE, debug = FALSE)
cextend(x, strict = TRUE, debug = FALSE)
moral(x, debug = FALSE)

cextend.all(x, debug = FALSE)

colliders(x, arcs = FALSE, debug = FALSE)
shielded.colliders(x, arcs = FALSE, debug = FALSE)
unshielded.colliders(x, arcs = FALSE, debug = FALSE)
vstructs(x, arcs = FALSE, debug = FALSE)

Arguments

Details

Note that arcs whose directions are dictated by the parametric assumptions of the network are preserved as directed arcs in cpdag(). This means that, in a conditional Gaussian network, arcs from discrete nodes to continuous nodes will be treated as whitelisted in their only valid direction.

cextend.all() returns all possible consistent extensions of a CPDAG, whereas cextend() returns only one. cextend.all() also accepts a PDAG or an MPDAG (a partially directed graph carrying background knowledge, in which some directed arcs are not compelled), in which case it enumerates the consistent extensions that agree with the given arc directions; it stops with an error if no such extension exists.

Value

cpdag() returns an object of class bn, representing the equivalence class. moral on the other hand returns the moral graph. See bn-class for details.

cextend() returns an object of class bn, representing a DAG that is the consistent extension of x.

cextend.all() returns an object of class bn or a list of objects of class bn.

vstructs(), colliders(), shielded.colliders() and unshielded.colliders() return a matrix with either 2 or 3 columns, according to the value of the arcs argument.

Author(s)

Marco Scutari

References

Dor D (1992). A Simple Algorithm to Construct a Consistent Extension of a Partially Oriented Graph. UCLA, Cognitive Systems Laboratory.

Koller D, Friedman N (2009). Probabilistic Graphical Models: Principles and Techniques. MIT Press.

Pearl J (2009). Causality: Models, Reasoning and Inference. Cambridge University Press, 2nd edition.

Andersson SA, Madigan D, Perlman MD (1997). "A Characterization of Markov Equivalence Classes for Acyclic Digraphs." The Annals of Statistics, 25(2):505–541.

Wienobst M, Luttermann M, Bannach M, Liskiewicz M (2023). "Efficient Enumeration of Markov Equivalent DAGs." Proceedings of the Thirty-Seventh AAAI Conference on Artificial Intelligence, 12313–12320.

Examples

data(learning.test)
dag = hc(learning.test)
cpdag(dag)
vstructs(dag)

Description

Perform conditional probability queries (CPQs).

Usage

cpquery(fitted, event, evidence, cluster, method = "ls", ...,
  debug = FALSE)
cpdist(fitted, nodes, evidence, cluster, method = "ls", ...,
  debug = FALSE)

Arguments

Details

cpquery() estimates the conditional probability of event given evidence using the method specified in the method argument.

cpdist() generates random samples conditional on the evidence using the method specified in the method argument.

Both cpquery() and cpdist() are based on Monte Carlo particle filters, so they may return slightly different values on different runs due to simulation noise.

Value

cpquery() returns a numeric value, the conditional probability of event conditional on evidence.

cpdist() returns a data frame containing the samples generated from the conditional distribution of the nodes conditional on evidence. The data frame has class c("bn.cpdist", "data.frame"), and a method attribute storing the value of the method argument. In likelihood weighting, the weights are also attached as an attribute named weights.

Logic Sampling

Logic sampling is an approximate inference algorithm.

The event and evidence arguments must be two expressions describing the event of interest and the conditioning evidence in a format such that, if we denote with data the data set the network was learned from, data[evidence, ] and data[event, ] return the correct observations. If either event or evidence is set to TRUE, an unconditional probability query is performed with respect to that argument.

Three tuning parameters are available:

• n: a positive integer number, the number of random samples to generate from fitted. The default value is 5000 * log10(nparams(fitted)) for discrete and conditional Gaussian networks and 500 * nparams(fitted) for Gaussian networks.

• batch: a positive integer number, the number of random samples that are generated at one time. Defaults to 10^4. If the n is very large (say, 10^12), R would run out of memory if it tried to generate them all at once. Instead, random samples are generated in batches of size batch, and each batch is discarded before generating the next.

• query.nodes: a vector of character strings, the labels of the nodes involved in event and evidence. Simple queries do not require generating samples from all the nodes in the network, so cpquery() and cpdist() try to identify which nodes are used in event and evidence and reduce the network to their upper closure. query.nodes may be used to manually specify these nodes when automatic identification fails; there is no reason to use it otherwise.

Note that the number of samples returned by cpdist() is always smaller than n, because logic sampling is a form of rejection sampling. Therefore, only the observations matching evidence (out of the n that are generated) are returned, and their number depends on the probability of evidence. Furthermore, the probabilities returned by cpquery() are approximate estimates, and they will not sum up to 1 even when the corresponding underlying values do if they are computed in separate calls to cpquery().

Likelihood Weighting

Likelihood weighting is an approximate inference algorithm based on Monte Carlo sampling.

The event argument must be an expression describing the event of interest, as in logic sampling. The evidence argument must be a named list:

• Each element corresponds to one node in the network and must contain the value that node will be set to when sampling.

• In the case of a continuous node, two values can also be provided. In that case, the value for that node will be sampled from a uniform distribution on the interval delimited by the specified values.

• In the case of a discrete or ordinal node, two or more values can also be provided. In that case, the value for that node will be sampled with uniform probability from the set of specified values.

If either event or evidence is set to TRUE, an unconditional probability query is performed with respect to that argument.

Tuning parameters are the same as for logic sampling: n, batch and query.nodes.

Note that the samples returned by cpdist() are generated from the mutilated network and need to be appropriately weighted when computing summary statistics (for more details, see the references below). cpquery() does that automatically when computing the final conditional probability. Also note that the batch argument is ignored in cpdist() for speed and memory efficiency. Furthermore, the probabilities returned by cpquery() are approximate estimates, and they will not sum up to 1 even when the corresponding underlying values do if they are computed in separate calls to cpquery().

Author(s)

Marco Scutari

References

Koller D, Friedman N (2009). Probabilistic Graphical Models: Principles and Techniques. MIT Press.

Korb K, Nicholson AE (2010). Bayesian Artificial Intelligence. Chapman & Hall/CRC, 2nd edition.

Examples

## discrete Bayesian network (it is the same with ordinal nodes).
data(learning.test)
fitted = bn.fit(hc(learning.test), learning.test)
# the result should be around 0.025.
cpquery(fitted, (B == "b"), (A == "a"))
# programmatically build a conditional probability query...
var = names(learning.test)
obs = 2
str = paste("(", names(learning.test)[-3], " == '",
        sapply(learning.test[obs, -3], as.character), "')",
        sep = "", collapse = " & ")
str
str2 = paste("(", names(learning.test)[3], " == '",
         as.character(learning.test[obs, 3]), "')", sep = "")
str2

cmd = paste("cpquery(fitted, ", str2, ", ", str, ")", sep = "")
eval(parse(text = cmd))
# ... but note that predict works better in this particular case.
attr(predict(fitted, "C", learning.test[obs, -3], prob = TRUE), "prob")
# do the same with likelihood weighting.
cpquery(fitted, event = eval(parse(text = str2)),
  evidence = as.list(learning.test[2, -3]), method = "lw")
attr(predict(fitted, "C", learning.test[obs, -3],
               method = "bayes-lw", prob = TRUE), "prob")
# conditional distribution of A given C == "c".
table(cpdist(fitted, "A", (C == "c")))

## Gaussian Bayesian network.
data(gaussian.test)
fitted = bn.fit(hc(gaussian.test), gaussian.test)
# the result should be around 0.04.
cpquery(fitted,
  event = ((A >= 0) & (A <= 1)) & ((B >= 0) & (B <= 3)),
  evidence = (C + D < 10))

Description

Screen and transform the data to make them more suitable for structure and parameter learning.

Usage

# discretise continuous data into factors.
  discretize(data, method, breaks = 3, ordered = FALSE, ..., debug = FALSE)
  # screen continuous data for highly correlated pairs of variables.
  dedup(data, method, threshold, debug = FALSE)

Arguments

Details

discretize() takes a data frame as its first argument and returns a second data frame of discrete variables, transformed using one of three methods: "interval", "quantile" or "hartemink". Discrete variables are left unchanged.

The "hartemink" method has two additional tuning parameters:

• idisc: the method used for the initial marginal discretisation of the variables, either "interval" or "quantile".

• ibreaks: the number of levels the variables are initially discretised into, in the same format as in the breaks argument.

It is sometimes the case that the "quantile" method cannot discretise one or more variables in the data without generating zero-length intervals because the quantiles are not unique. If method = "quantile", discretize() will produce an error. If method = "quantile" and idisc = "quantile", discretize() will try to lower the number of breaks set by the ibreaks argument until quantiles are distinct. If this is not possible without making ibreaks smaller than breaks, discretize() will return an error.

dedup() screens the data for pairs of highly correlated variables, and discards one in each pair.

Both discretize() and dedup() accept data with missing values.

Value

discretize() returns a data frame with the same structure (number of columns, column names, etc.) as data, containing the discretised variables. The data frame has an attribute "cutpoints" containing the cutpoints used internally by discretize(), which are made available so that discretising another data set yields the same set of factors.

dedup() returns a data frame with a subset of the columns of data.

Author(s)

Marco Scutari

References

Hartemink A (2001). Principled Computational Methods for the Validation and Discovery of Genetic Regulatory Networks. Ph.D. thesis, School of Electrical Engineering and Computer Science, Massachusetts Institute of Technology.

Examples

data(gaussian.test)
d = discretize(gaussian.test, method = 'hartemink', breaks = 4, ibreaks = 10)
plot(hc(d))
d2 = dedup(gaussian.test)

Description

Check whether two nodes are d-separated.

Usage

dsep(bn, x, y, z)

Arguments

Value

dsep() returns TRUE if x and y are d-separated by z, and FALSE otherwise.

Author(s)

Marco Scutari

References

Koller D, Friedman N (2009). Probabilistic Graphical Models: Principles and Techniques. MIT Press.

Examples

bn = model2network("[A][C|A][B|C]")
dsep(bn, "A", "B", "C")
bn = model2network("[A][C][B|A:C]")
dsep(bn, "A", "B", "C")

Description

Read networks saved from other programs into bn.fit objects, and dump bn and bn.fit objects into files for other programs to read.

Usage

# Old (non-XML) Bayesian Interchange format.
read.bif(file, debug = FALSE)
write.bif(file, fitted)

# Microsoft Interchange format.
read.dsc(file, debug = FALSE)
write.dsc(file, fitted)

# HUGIN flat network format.
read.net(file, debug = FALSE)
write.net(file, fitted)

# Graphviz DOT format.
write.dot(file, graph)

Arguments

Value

read.bif(), read.dsc() and read.net() return an object of class bn.fit.

write.bif(), write.dsc(), write.net() and write.dot() return NULL invisibly.

Note

All the networks present in the Bayesian Network Repository have associated BIF, DSC and NET files that can be imported with read.bif(), read.dsc() and read.net().

HUGIN can import and export NET files; Netica can read (but not write) DSC files; and GeNIe can read and write both DSC and NET files.

DOT files can be read by Graphviz, Gephi and a variety of other programs.

Please note that these functions operate on a "best effort" basis, as the parsing of these formats was implemented by reverse-engineering the file formats from publicly available examples.

Author(s)

Marco Scutari

References

Bayesian Network Repository, https://www.bnlearn.com/bnrepository/.

Description

This is a synthetic data set used as a test case in the bnlearn package.

Usage

data(gaussian.test)

Format

The gaussian.test data set contains seven normal (Gaussian) variables.

Note

The R script to generate data from this network is available from https://www.bnlearn.com/documentation/networks/.

Examples

# load the data.
data(gaussian.test)
# create and plot the network structure.
dag = model2network("[A][B][E][G][C|A:B][D|B][F|A:D:E:G]")
## Not run: graphviz.plot(dag)

Description

Convert bn.fit objects to grain objects and vice versa.

Usage

## S3 method for class 'grain'
as.bn.fit(x, including.evidence = FALSE, ...)
## S3 method for class 'bn.fit'
as.grain(x)
## S3 method for class 'grain'
as.bn(x, ...)

Arguments

Value

An object of class grain (for as.grain), bn.fit (for as.bn.fit) or bn (for as.bn).

Note

Conditional probability tables in grain objects must be completely specified; on the other hand, bn.fit allows NaN values for unobserved parents' configurations. Such bn.fit objects will be converted to $m$ grain objects by replacing the missing conditional probability distributions with uniform distributions.

Another solution to this problem is to fit another bn.fit with method = "bayes" and a low iss value, using the same data and network structure.

Ordinal nodes will be treated as categorical by as.grain, disregarding the ordering of the levels.

Author(s)

Marco Scutari

Examples

## Not run: 
library(gRain)
a = bn.fit(hc(learning.test), learning.test)
b = as.grain(a)
c = as.bn.fit(b)
## End(Not run)

Description

Count directed acyclic graphs of various sizes and/or with specific characteristics.

Usage

count.graphs(type = "all.dags", ..., debug = FALSE)

Arguments

Details

The types of graphs, and the associated additional parameters, are:

• "all-dags": all directed acyclic graphs with nodes nodes.

• "dags-given-ordering": all directed acyclic graphs with a specific topological ordering and nodes nodes.

• "dags-with-k-roots": all directed acyclic graphs with k root nodes and nodes nodes.

• "dags-with-r-arcs": all directed acyclic graphs with r arcs and nodes nodes.

• "dags-in-equivalence-class": all directed acyclic arcs in the equivalence class described by eqclass.

Value

count.graphs() returns an object of class bigz from the gmp package, a vector with the graph counts.

Author(s)

Marco Scutari

References

Harary F, Palmer EM (1973). "Graphical Enumeration." Academic Press.

Rodionov VI (1992). "On the Number of Labeled Acyclic Digraphs." Discrete Mathematics, 105:319–321.

Liskovets VA (1976). "On the Number of Maximal Vertices of a Random Acyclic Digraph." Theory of Probability and its Applications, 20(2):401–409.

Wienobst M, Luttermann M, Bannach M, Liskiewicz (2023). "Efficient Enumeration of Markov Equivalent DAGs." Proceedings of the 37th AAAI Conference on Artificial Intelligence, 12313–12320.

Examples

## Not run: 
count.graphs("dags.with.r.arcs", nodes = 3:6, r = 2)

## End(Not run)

Description

Generate empty, complete or random directed acyclic graphs from a given set of nodes.

Usage

empty.graph(nodes, num = 1)
complete.graph(nodes, num = 1)
random.graph(nodes, num = 1, method = "ordered", ..., debug = FALSE)

Arguments

Details

Graph generation algorithms available in random.graph() are:

• full ordering based generation ("ordered"): generates graphs whose node ordering is given by the order of the labels in the nodes argument. The same algorithm is used in the randomDAG function in package pcalg.

• Ide's and Cozman's Generating Multi-connected DAGs algorithm ("ic-dag"): generates graphs with a uniform probability distribution over the set of multiconnected graphs.

• Melancon's and Philippe's Uniform Random Acyclic Digraphs algorithm ("melancon"): generates graphs with a uniform probability distribution over the set of all possible graphs.

Additional arguments for the random.graph() function are:

• prob: the probability of each arc to be present in a graph generated by the "ordered" algorithm. The default value is 2 / (length(nodes) - 1), which results in a sparse graph (the number of arcs should be of the same order as the number of nodes).

• burn.in: the number of iterations for the "ic-dag" and "melancon" algorithms to converge to a stationary (and uniform) probability distribution. The default value is 6 * length(nodes)^2.

• every: return only one graph every number of steps instead of all the graphs generated with "ic-dag" and "melancon". Since both algorithms are Markov Chain Monte Carlo methods, higher values of every yield more diverse set of networks. The default value is 1, that is, to return all generated networks.

• max.degree: the maximum degree for any node in a graph generated by the "ic-dag" and "melancon" algorithms. The default value is Inf.

• max.in.degree: the maximum in-degree for any node in a graph generated by the "ic-dag" and "melancon" algorithms. The default value is Inf.

• max.out.degree: the maximum out-degree for any node in a graph generated by the "ic-dag" and "melancon" algorithms. The default value is Inf.

empty.graph() generates num identical empty graphs, while complete.graph() generates num identical complete directed acyclic graphs.

Value

empty.graph(), complete.graph() and random[.graph() return an object of class bn (if num is equal to 1) or a list of objects of class bn (otherwise). If every is greated than 1, random.graph() always returns a list, regardless of the number of graphs it contains.

Author(s)

Marco Scutari

References

Ide JS, Cozman FG (2002). "Random Generation of Bayesian Networks." Proceedings of the 16th Brazilian Symposium on Artificial Intelligence, 366–375.

Melancon G, Dutour I, Bousquet-Melou M (2001). "Random Generation of Directed Acyclic Graphs." Electronic Notes in Discrete Mathematics, 10:202–207.

Melancon G, Philippe F (2004). "Generating Connected Acyclic Digraphs Uniformly at Random." Information Processing Letters, 90(4):209–213.

Examples

empty.graph(LETTERS[1:8])
random.graph(LETTERS[1:8])
plot(random.graph(LETTERS[1:8], method = "ic-dag", max.in.degree = 2))
plot(random.graph(LETTERS[1:8]))
plot(random.graph(LETTERS[1:8], prob = 0.2))

Description

Convert bn and bn.fit objects to graphNEL and graphAM objects and vice versa.

Usage

## S3 method for class 'graphNEL'
as.bn(x, ..., check.cycles = TRUE)
## S3 method for class 'graphAM'
as.bn(x, ..., check.cycles = TRUE)
## S3 method for class 'bn'
as.graphNEL(x)
## S3 method for class 'bn.fit'
as.graphNEL(x)
## S3 method for class 'bn'
as.graphAM(x)
## S3 method for class 'bn.fit'
as.graphAM(x)

Arguments

Value

An object of the relevant class.

Author(s)

Marco Scutari

Examples

## Not run: 
library(graph)
a = bn.fit(hc(learning.test), learning.test)
b = as.graphNEL(a)
c = as.bn(b)
## End(Not run)

Description

Check and manipulate graph-related properties of an object of class bn.

Usage

# check whether the graph is acyclic/completely directed.
acyclic(x, directed = FALSE, debug = FALSE)
directed(x)
valid.dag(x, debug = FALSE)
valid.cpdag(x, debug = FALSE)
valid.ug(x, debug = FALSE)
# check whether there is a path between two nodes.
path.exists(x, from, to, direct = TRUE, underlying.graph = FALSE, debug = FALSE)
# build the skeleton or a complete orientation of the graph.
skeleton(x)
pdag2dag(x, ordering)
# build a subgraph spanning a subset of nodes.
subgraph(x, nodes)
# perturb a network by random arc additions, removals and reversals.
perturb(x, nops, ops = c("set", "drop", "reverse"), maxp = Inf, debug = FALSE)

Arguments

Value

acyclic(), path() and directed() return a boolean value.
skeleton(), pdag2dag(), subgraph() and perturb() return an object of class bn.

valid.dag(), valid.cpdag() and valid.ug() return a boolean value.

Author(s)

Marco Scutari

References

Bang-Jensen J, Gutin G (2009). Digraphs: Theory, Algorithms and Applications. Springer, 2nd edition.

Andersson SA, Madigan D, Perlman MD (1997). "A Characterization of Markov Equivalence Classes for Acyclic Digraphs." The Annals of Statistics, 25(2):505–541.

Examples

data(learning.test)
cpdag = pc.stable(learning.test)

acyclic(cpdag)
directed(cpdag)
dag = pdag2dag(cpdag, ordering = LETTERS[1:6])
dag
directed(dag)
skeleton(dag)
perturb(dag, nops = 3, debug = TRUE)

Description

Plot a Bayesian network as a graph whose nodes are barplots representing the distributions of the corresponding variables. Requires the Rgraphviz package and, for discrete networks, also the gRain package.

Usage

graphviz.chart(x, type = "barchart", layout = "dot", draw.labels = TRUE,
  grid = FALSE, scale = c(0.75, 1.1), col = "black", bg = "transparent",
  text.col = "black", bar.col = "black", strip.bg = bg, main = NULL,
  sub = NULL)

Arguments

Value

graphviz.chart() invisibly returns NULL.

Author(s)

Marco Scutari

Examples

## Not run: 
modelstring = paste("[HIST|LVF][CVP|LVV][PCWP|LVV][HYP][LVV|HYP:LVF][LVF]",
  "[STKV|HYP:LVF][ERLO][HRBP|ERLO:HR][HREK|ERCA:HR][ERCA][HRSA|ERCA:HR][ANES]",
  "[APL][TPR|APL][ECO2|ACO2:VLNG][KINK][MINV|INT:VLNG][FIO2][PVS|FIO2:VALV]",
  "[SAO2|PVS:SHNT][PAP|PMB][PMB][SHNT|INT:PMB][INT][PRSS|INT:KINK:VTUB][DISC]",
  "[MVS][VMCH|MVS][VTUB|DISC:VMCH][VLNG|INT:KINK:VTUB][VALV|INT:VLNG][ACO2|VALV]",
  "[CCHL|ACO2:ANES:SAO2:TPR][HR|CCHL][CO|HR:STKV][BP|CO:TPR]", sep = "")
dag = model2network(modelstring)
fitted = bn.fit(dag, alarm)

# Netica style.
graphviz.chart(fitted, grid = TRUE, bg = "beige", bar.col = "black")
# Hugin style.
graphviz.chart(fitted, type = "barprob", grid = TRUE, bar.col = "green",
  strip.bg = "lightyellow")
# GeNIe style.
graphviz.chart(fitted, col = "darkblue", bg = "azure", bar.col = "darkblue")
# personal favourites.
graphviz.chart(fitted, type = "barprob", grid = TRUE, bar.col = "darkgreen",
  strip.bg = "lightskyblue")
graphviz.chart(fitted, type = "barprob", grid = TRUE, bar.col = "gold",
  strip.bg = "lightskyblue")
# dot-plot version.
graphviz.chart(fitted, type = "dotplot")

## End(Not run)

Description

Plot the graph associated with a Bayesian network using the Rgraphviz package.

Usage

graphviz.plot(x, highlight = NULL, groups, layout = "dot",
  shape = "rectangle", fontsize = 12, main = NULL, sub = NULL,
  render = TRUE)

Arguments

Details

The highlight argument is a list with at least one of the following elements:

• nodes: a character vector, the labels of the nodes to be highlighted.

• arcs: the arcs to be highlighted (a two-column matrix, whose columns are labelled from and to).

and optionally one or more of the following graphical parameters:

• col: an integer or character string (the highlight colour for the arcs and the node frames). The default value is red.

• textCol: an integer or character string (the highlight colour for the labels of the nodes). The default value is black.

• fill: an integer or character string (the colour used as a background colour for the nodes). The default value is transparent.

• lwd: a positive number (the line width of highlighted arcs). It overrides the line width settings in strength.plot(). The default value is to use the global settings of Rgraphviz.

• lty: the line type of highlighted arcs. Possible values are 0, 1, 2, 3, 4, 5, 6, "blank", "solid", "dashed", "dotted", "dotdash", "longdash" and "twodash". The default value is to use the global settings of Rgraphviz.

Note that all these parameters take a single value, which is then applied to all highlighted nodes and arcs.

Value

graphviz.plot() returns invisibly the graph object produced from the network passed as the x argument. It can be further modified using the graph and Rgraphviz packages.

Author(s)

Marco Scutari

See Also

plot.bn.

Description

Hailfinder is a Bayesian network designed to forecast severe summer hail in northeastern Colorado.

Usage

data(hailfinder)

Format

The hailfinder data set contains the following 56 variables:

• N07muVerMo (10.7mu vertical motion): a four-level factor with levels StrongUp, WeakUp, Neutral and Down.

• SubjVertMo (subjective judgment of vertical motion): a four-level factor with levels StrongUp, WeakUp, Neutral and Down.

• QGVertMotion (quasigeostrophic vertical motion): a four-level factor with levels StrongUp, WeakUp, Neutral and Down.

• CombVerMo (combined vertical motion): a four-level factor with levels StrongUp, WeakUp, Neutral and Down.

• AreaMesoALS (area of meso-alpha): a four-level factor with levels StrongUp, WeakUp, Neutral and Down.

• SatContMoist (satellite contribution to moisture): a four-level factor with levels VeryWet, Wet, Neutral and Dry.

• RaoContMoist (reading at the forecast center for moisture): a four-level factor with levels VeryWet, Wet, Neutral and Dry.

• CombMoisture (combined moisture): a four-level factor with levels VeryWet, Wet, Neutral and Dry.

• AreaMoDryAir (area of moisture and adry air): a four-level factor with levels VeryWet, Wet, Neutral and Dry.

• VISCloudCov (visible cloud cover): a three-level factor with levels Cloudy, PC and Clear.

• IRCloudCover (infrared cloud cover): a three-level factor with levels Cloudy, PC and Clear.

• CombClouds (combined cloud cover): a three-level factor with levels Cloudy, PC and Clear.

• CldShadeOth (cloud shading, other): a three-level factor with levels Cloudy, PC and Clear.

• AMInstabMt (AM instability in the mountains): a three-level factor with levels None, Weak and Strong.

• InsInMt (instability in the mountains): a three-level factor with levels None, Weak and Strong.

• WndHodograph (wind hodograph): a four-level factor with levels DCVZFavor, StrongWest, Westerly and Other.

• OutflowFrMt (outflow from mountains): a three-level factor with levels None, Weak and Strong.

• MorningBound (morning boundaries): a three-level factor with levels None, Weak and Strong.

• Boundaries (boundaries): a three-level factor with levels None, Weak and Strong.

• CldShadeConv (cloud shading, convection): a three-level factor with levels None, Some and Marked.

• CompPlFcst (composite plains forecast): a three-level factor with levels IncCapDecIns, LittleChange and DecCapIncIns.

• CapChange (capping change): a three-level factor with levels Decreasing, LittleChange and Increasing.

• LoLevMoistAd (low-level moisture advection): a four-level factor with levels StrongPos, WeakPos, Neutral and Negative.

• InsChange (instability change): three-level factor with levels Decreasing, LittleChange and Increasing.

• MountainFcst (mountains (region 1) forecast): a three-level factor with levels XNIL, SIG and SVR.

• Date (date): a six-level factor with levels May15_Jun14, Jun15_Jul1, Jul2_Jul15, Jul16_Aug10, Aug11_Aug20 and Aug20_Sep15.

• Scenario (scenario): an eleven-level factor with levels A, B, C, D, E, F, G, H, I, J and K.

• ScenRelAMCIN (scenario relevant to AM convective inhibition): a two-level factor with levels AB and CThruK.

• MorningCIN (morning convective inhibition): a four-level factor with levels None, PartInhibit, Stifling and TotalInhibit.

• AMCINInScen (AM convective inhibition in scenario): a three-level factor with levels LessThanAve, Average and MoreThanAve.

• CapInScen (capping withing scenario): a three-level factor with levels LessThanAve, Average and MoreThanAve.

• ScenRelAMIns (scenario relevant to AM instability): a six-level factor with levels ABI, CDEJ, F, G, H and K.

• LIfr12ZDENSd (LI from 12Z DEN sounding): a four-level factor with levels LIGt0, N1GtLIGt_4, N5GtLIGt_8 and LILt_8.

• AMDewptCalPl (AM dewpoint calculations, plains): a three-level factor with levels Instability, Neutral and Stability.

• AMInsWliScen (AM instability within scenario): a three-level factor with levels LessUnstable, Average and MoreUnstable.

• InsSclInScen (instability scaling within scenario): a three-level factor with levels LessUnstable, Average and MoreUnstable.

• ScenRel34 (scenario relevant to regions 2/3/4): a five-level factor with levels ACEFK, B, D, GJ and HI.

• LatestCIN (latest convective inhibition): a four-level factor with levels None, PartInhibit, Stifling and TotalInhibit.

• LLIW (LLIW severe weather index): a four-level factor with levels Unfavorable, Weak, Moderate and Strong.

• CurPropConv (current propensity to convection): a four-level factor with levels None, Slight, Moderate and Strong.

• ScnRelPlFcst (scenario relevant to plains forecast): an eleven-level factor with levels A, B, C, D, E, F, G, H, I, J and K.

• PlainsFcst (plains forecast): a three-level factor with levels XNIL, SIG and SVR.

• N34StarFcst (regions 2/3/4 forecast): a three-level factor with levels XNIL, SIG and SVR.

• R5Fcst (region 5 forecast): a three-level factor with levels XNIL, SIG and SVR.

• Dewpoints (dewpoints): a seven-level factor with levels LowEverywhere, LowAtStation, LowSHighN, LowNHighS, LowMtsHighPl, HighEverywher, Other.

• LowLLapse (low-level lapse rate): a four-level factor with levels CloseToDryAd, Steep, ModerateOrLe and Stable.

• MeanRH (mean relative humidity): a three-level factor with levels VeryMoist, Average and Dry.

• MidLLapse (mid-level lapse rate): a three-level factor with levels CloseToDryAd, Steep and ModerateOrLe.

• MvmtFeatures (movement of features): a four-level factor with levels StrongFront, MarkedUpper, OtherRapid and NoMajor.

• RHRatio (realtive humidity ratio): a three-level factor with levels MoistMDryL, DryMMoistL and other.

• SfcWndShfDis (surface wind shifts and discontinuities): a seven-level factor with levels DenvCyclone, E_W_N, E_W_S, MovigFtorOt, DryLine, None and Other.

• SynForcng (synoptic forcing): a five-level factor with levels SigNegative, NegToPos, SigPositive, PosToNeg and LittleChange.

• TempDis (temperature discontinuities): a four-level factor with levels QStationary, Moving, None, Other.

• WindAloft (wind aloft): a four-level factor with levels LV, SWQuad, NWQuad, AllElse.

• WindFieldMt (wind fields, mountains): a two-level factor with levels Westerly and LVorOther.

• WindFieldPln (wind fields, plains): a six-level factor with levels LV, DenvCyclone, LongAnticyc, E_NE, SEquad and WidespdDnsl.

Note

The complete BN can be downloaded from https://www.bnlearn.com/bnrepository/.

Source

Abramson B, Brown J, Edwards W, Murphy A, Winkler RL (1996). "Hailfinder: A Bayesian System for Forecasting Severe Weather." International Journal of Forecasting, 12(1):57–71.

Examples

# load the data.
data(hailfinder)
# create and plot the network structure.
modelstring = paste0("[N07muVerMo][SubjVertMo][QGVertMotion][SatContMoist][RaoContMoist]",
  "[VISCloudCov][IRCloudCover][AMInstabMt][WndHodograph][MorningBound][LoLevMoistAd][Date]",
  "[MorningCIN][LIfr12ZDENSd][AMDewptCalPl][LatestCIN][LLIW]",
  "[CombVerMo|N07muVerMo:SubjVertMo:QGVertMotion][CombMoisture|SatContMoist:RaoContMoist]",
  "[CombClouds|VISCloudCov:IRCloudCover][Scenario|Date][CurPropConv|LatestCIN:LLIW]",
  "[AreaMesoALS|CombVerMo][ScenRelAMCIN|Scenario][ScenRelAMIns|Scenario][ScenRel34|Scenario]",
  "[ScnRelPlFcst|Scenario][Dewpoints|Scenario][LowLLapse|Scenario][MeanRH|Scenario]",
  "[MidLLapse|Scenario][MvmtFeatures|Scenario][RHRatio|Scenario][SfcWndShfDis|Scenario]",
  "[SynForcng|Scenario][TempDis|Scenario][WindAloft|Scenario][WindFieldMt|Scenario]",
  "[WindFieldPln|Scenario][AreaMoDryAir|AreaMesoALS:CombMoisture]",
  "[AMCINInScen|ScenRelAMCIN:MorningCIN][AMInsWliScen|ScenRelAMIns:LIfr12ZDENSd:AMDewptCalPl]",
  "[CldShadeOth|AreaMesoALS:AreaMoDryAir:CombClouds][InsInMt|CldShadeOth:AMInstabMt]",
  "[OutflowFrMt|InsInMt:WndHodograph][CldShadeConv|InsInMt:WndHodograph][MountainFcst|InsInMt]",
  "[Boundaries|WndHodograph:OutflowFrMt:MorningBound][N34StarFcst|ScenRel34:PlainsFcst]",
  "[CompPlFcst|AreaMesoALS:CldShadeOth:Boundaries:CldShadeConv][CapChange|CompPlFcst]",
  "[InsChange|CompPlFcst:LoLevMoistAd][CapInScen|CapChange:AMCINInScen]",
  "[InsSclInScen|InsChange:AMInsWliScen][R5Fcst|MountainFcst:N34StarFcst]",
  "[PlainsFcst|CapInScen:InsSclInScen:CurPropConv:ScnRelPlFcst]")
dag = model2network(modelstring)
## Not run: graphviz.plot(dag, shape = "ellipse")

Description

Learn the structure of a Bayesian network with Max-Min Hill Climbing (MMHC), Hybrid HPC (H2PC), and the more general 2-phase Restricted Maximisation (RSMAX2) hybrid algorithms.

Usage

rsmax2(x, whitelist, blacklist, restrict = "si.hiton.pc",
  maximize = "hc", restrict.args = list(), maximize.args = list(), debug = FALSE)
mmhc(x, whitelist, blacklist, restrict.args = list(),
  maximize.args = list(), debug = FALSE)
h2pc(x, whitelist, blacklist, restrict.args = list(),
  maximize.args = list(), debug = FALSE)

Arguments

Value

An object of class bn. See bn-class for details.

Note

mmhc() is simply rsmax2() with restrict set to "mmpc" and maximize set to "hc". Similarly, h2pc is simply rsmax2() with restrict set to "hpc"and maximize set to "hc".

See structure learning for a complete list of structure learning algorithms with the respective references.

Author(s)

Marco Scutari

See Also

local discovery algorithms, score-based algorithms, constraint-based algorithms, causal discovery algorithms.

Description

Convert bn and bn.fit objects to igraph and vice versa.

Usage

## S3 method for class 'igraph'
as.bn(x, ..., check.cycles = TRUE)
## S3 method for class 'bn'
as.igraph(x, ...)
## S3 method for class 'bn.fit'
as.igraph(x, ...)

Arguments

Value

An object of the relevant class.

Author(s)

Marco Scutari

Examples

## Not run: 
a = bn.fit(hc(learning.test), learning.test)
b = as.igraph(a)
plot(b, edge.arrow.mode = 2L * !igraph::which_mutual(b))
c = as.bn(b)
## End(Not run)