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, providing the tools needed for learning and working with discrete Bayesian networks, Gaussian Bayesian networks and conditional linear Gaussian Bayesian networks on real-world data. Incomplete data with missing values are also supported. 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 optimization 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, as they use conditional independence tests (usually 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 to learn the structure of Bayesian network classifiers are described in network classifiers. All algorithms support the use of whitelists and blacklists to include and exclude arcs from the networks (see whitelists and blacklists); and many have parallel implementation built on the parallel package. Bayesian network scores support the use of graphical priors.

Parameter learning approaches include both frequentist and Bayesian estimators. Inference is implemented using approximate algorithms via particle filters approaches such as likelihood weighting, and covers conditional probability queries, prediction and imputation.

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 value of the imaginary sample size for the BDe score, assuming a uniform prior and given a network structure and a data set.

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 (also known as 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 (i.e. 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 invisibly an updated copy of x.

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 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 the choice of the 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 boostrap 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

Small synthetic data set from Lauritzen and Spiegelhalter (1988) about lung diseases (tuberculosis, lung cancer or bronchitis) and visits to Asia.

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 seems 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 which is equivalent to the ration of the corresponding posterior probabilities if we assume the uniform prior over all possible DAGs. 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 classic 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's never changed afterward.

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

  • blacklist: a copy of the blacklist argument (a two-column matrix, whose columns are labeled from and to) as transformed by sanitization 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 labeled 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 optimizations 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 labeled 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.

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 statistics.args are extracted from the list and passed to statistics as the 2nd, 3rd, etc. 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 and 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 averge 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 element corresponds to one fold and contains the indices for the observations that are included to 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 valuer are better.

• Gaussian Log-Likelihood Loss (logl-g): the negated expected log-likelihood for Gaussian Bayesian networks. Lower values are better.

• Classification Error (pred): the prediction error for a single node in a discrete network. Frequentist predictions are used, so the values of the target node are predicted using only the information present in its local distribution (from its parents). Lower values are better.

• Posterior Classification Error (pred-lw and pred-lw-cg): similar to the above, but predictions are computed from an arbitrary set of nodes using likelihood weighting to obtain Bayesian posterior estimates. pred-lw applies to discrete Bayesian networks, pred-lw-cg to (discrete nodes in) hybrid networks. 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 node in a Gaussian Bayesian network. Higher values are better.

• Posterior Predictive Correlation (cor-lw and cor-lw-cg): similar to the above, but predictions are computed from an arbitrary set of nodes using likelihood weighting to obtain Bayesian posterior estimates. cor-lw applies to Gaussian networks and cor-lw-cg to (continuous nodes in) hybrid networks. Higher values are better.

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

• Posterior Mean Squared Error (mse-lw and mse-lw-cg): similar to the above, but predictions are computed from an arbitrary set of nodes using likelihood weighting to obtain Bayesian posterior estimates. mse-lw applies to Gaussian networks and mse-lw-cg to (continuous nodes in) hybrid networks. Lower values are better.

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 target node for prediction in all loss functions but logl, logl-g and logl-cg.

• from: a vector of character strings, the labels of the nodes used to predict the target node in pred-lw, pred-lw-cg, cor-lw, cor-lw-cg, mse-lw and mse-lw-cg. The default is to use all the other nodes in the network. Loss functions pred, cor and mse implicitly predict only from the parents of the target node.

• n: a positive integer, the number of particles used by likelihood weighting for pred-lw, pred-lw-cg, cor-lw, cor-lw-cg, mse-lw and mse-lw-cg. The default value is 500.

Note that if bn is a Bayesian network classifier, pred and pred-lw both give exact posterior predictions computed using the closed-form formulas for naive Bayes and TAN.

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 respected 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(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-Maximization implementation of the estimators above.

Available methods for hybrid Bayesian networks are:

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

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

Available methods for discrete Bayesian networks are:

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

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

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 one of mle, mle-g or mle-cg, 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 mximum likelihood estimate of NaN for all parents configurations that are not observed in data. Similarly, regression coefficients are set to NA if the linear regressions correspoding 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 maximization 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 only stop in case of negative improvement, which can happen due to stochastic noise if the imputation of the missing data uses approximate inference.

• params.threshold: a non-negative numeric value, the minimum maximum elative 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 of the differences between parameter values divided scaled by the parameter value in the lastest 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 initialize 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.

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 hybrid 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 if the discrete parents of the conditional Gaussian node, stored as a factor.

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 the parameters, not only the structure) using knowledge from experts in the field instead of 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 (i.e. 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() produces NA parameter estimates for discrete and conditional Gaussian nodes when there are (discrete) parents 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 the 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 again the original parameters
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))

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 latter is discrete (i.e. 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 (i.e. ordered factors) have class bn.fit.onode and contain the same elements as bn.fit.dnode nodes. Networks containing 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 (i.e. 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 (i.e. the standard error).

Hybrid (i.e. 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 (i.e. the standard error) for each configuration of the discrete parents.

Furthermore, Bayesian network classifiers store the label of the training node in an additional attribute named training.

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 an object 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 (which are conditional probabilities for discrete nodes and linear regression coefficients for Gaussian and conditional Gaussian nodes).

residuals() (and its alias resid()) extracts model residuals and fitted() (and its alias
fitted.values()) extracts fitted values from Gaussian and conditional Gaussian nodes. 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 or bn.fit.onode).

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.

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

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

• observed: a factor or a numeric vector, the observed values for the target node 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 in the 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 (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 were 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 used to determine if a strength coefficient is significant.

An optional column called direction may also be present, giving the probability of the direction of an arc 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

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 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, ...)

shd(learned, true, wlbl = FALSE, 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 visualize differences between graphs in various way 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, false negative arcs relative to that first network 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, tp.lwd. In addition, it is possible not to plot the reference network at all by setting show.first to FALSE.

Regardless of the visualization, the nodes are arranged to be in the same position for all the networks to make it easier to compare them.

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.

graphviz.compare() plots one or more figures and returns invisibly 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.

Note

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.

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.

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)
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 can be 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 = NULL, blacklist = NULL, test = NULL,
  alpha = 0.05, ..., max.sx = NULL, debug = FALSE, undirected = FALSE)
gs(x, cluster, whitelist = NULL, blacklist = NULL, test = NULL,
  alpha = 0.05, ..., max.sx = NULL, debug = FALSE, undirected = FALSE)
iamb(x, cluster, whitelist = NULL, blacklist = NULL, test = NULL,
  alpha = 0.05, ..., max.sx = NULL, debug = FALSE, undirected = FALSE)
fast.iamb(x, cluster, whitelist = NULL, blacklist = NULL, test = NULL,
  alpha = 0.05, ..., max.sx = NULL, debug = FALSE, undirected = FALSE)
inter.iamb(x, cluster, whitelist = NULL, blacklist = NULL, test = NULL,
  alpha = 0.05, ..., max.sx = NULL, debug = FALSE, undirected = FALSE)
iamb.fdr(x, cluster, whitelist = NULL, blacklist = NULL, test = NULL,
  alpha = 0.05, ..., max.sx = NULL, debug = FALSE, undirected = FALSE)
mmpc(x, cluster, whitelist = NULL, blacklist = NULL, test = NULL,
  alpha = 0.05, ..., max.sx = NULL, debug = FALSE, undirected = TRUE)
si.hiton.pc(x, cluster, whitelist = NULL, blacklist = NULL, test = NULL,
  alpha = 0.05, ..., max.sx = NULL, debug = FALSE, undirected = TRUE)
hpc(x, cluster, whitelist = NULL, blacklist = NULL, 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, which they handle 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, 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 on 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)

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.

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.

vstructs() returns 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.

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)

mutilated(x, evidence)

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.

mutilated constructs the mutilated network arising from an ideal intervention setting the nodes involved to the values specified by evidence. In this case evidence must be provided as a list in the same format as for likelihood weighting (see below).

Note that both cpquery and cpdist are based on Monte Carlo particle filters, and therefore 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 meth, -8od attribute storing the value of the method argument. In the case of likelihood weighting, the weights are also attached as an attribute called weights.

mutilated returns a bn or bn.fit object, depending on the class of x.

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 (e.g. 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, discarding each batch 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 to generate 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 weighted appropriately 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))

## ideal interventions and mutilated networks.
mutilated(fitted, evidence = list(F = 42))

Description

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

Usage

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

Arguments

Details

discretize() takes a data frame as its first argument and returns a secdond data frame of discrete variables, transformed using 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 discretization of the variables, either interval or quantile.

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

It is sometimes the case that the quantile method cannot discretize 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 produce 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 discretized variables.

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 work on a "best effort" basis, as the parsing of these formats have been implemented by reverse engineering the file format from publicly available examples.

Author(s)

Marco Scutari

References

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

Description

This 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 with specific characteristics.

Usage

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

Arguments

Details

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

• all-dags: all directed acyclic graphs.

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

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

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

Value

count.graphs() returns an objects 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.

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 based on Markov Chain Monte Carlo approaches, high values of every result in a more diverse set of networks. The default value is 1, i.e. to return all the networks that are generated.

• 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)
# 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)

Arguments

Value

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

Author(s)

Marco Scutari

References

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

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)

Description

Plot a Bayesian network as a graph whose nodes are barplots representing the marginal distributions of the corresponding variables. Requires the Rgraphviz and gRain packages.

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 labeled 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 that is then applied to all nodes and arcs that will be highlighted.

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 Maximization (RSMAX2) hybrid algorithms.

Usage

rsmax2(x, whitelist = NULL, blacklist = NULL, restrict = "si.hiton.pc",
  maximize = "hc", restrict.args = list(), maximize.args = list(), debug = FALSE)
mmhc(x, whitelist = NULL, blacklist = NULL, restrict.args = list(),
  maximize.args = list(), debug = FALSE)
h2pc(x, whitelist = NULL, blacklist = NULL, 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 hpcand 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.

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)

Description

Overview of the conditional independence tests implemented in bnlearn, with the respective reference publications.

Details

Unless otherwise noted, the reference publication for conditional independence tests is:

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

Additionally for continuous permutation tests:

Legendre P (2000). "Comparison of Permutation Methods for the Partial Correlation and Partial Mantel Tests". Journal of Statistical Computation and Simulation, 67:37–73.

and for semiparametric discrete tests:

Tsamardinos I, Borboudakis G (2010). "Permutation Testing Improves Bayesian Network Learning". Machine Learning and Knowledge Discovery in Databases, 322–337.

Available conditional independence tests (and the respective labels) for discrete Bayesian networks (categorical variables) are:

• mutual information: an information-theoretic distance measure. It's proportional to the log-likelihood ratio (they differ by a $2n$ factor) and is related to the deviance of the tested models. The asymptotic $\chi^2$ test (mi and mi-adf, with adjusted degrees of freedom), the Monte Carlo permutation test (mc-mi), the sequential Monte Carlo permutation test (smc-mi), and the semiparametric test (sp-mi) are implemented.

• shrinkage estimator for the mutual information (mi-sh): an improved asymptotic $\chi^2$ test based on the James-Stein estimator for the mutual information.

  Hausser J, Strimmer K (2009). "Entropy inference and the James-Stein estimator, with application to nonlinear gene association networks". Statistical Applications in Genetics and Molecular Biology, 10:1469–1484.

• Pearson's $X^2$: the classical Pearson's $X^2$ test for contingency tables. The asymptotic $\chi^2$ test (x2 and x2-adf, with adjusted degrees of freedom), the Monte Carlo permutation test (mc-x2), the sequential Monte Carlo permutation test (smc-x2) and semiparametric test (sp-x2) are implemented.

Available conditional independence tests (and the respective labels) for discrete Bayesian networks (ordered factors) are:

• Jonckheere-Terpstra: a trend test for ordinal variables. The asymptotic normal test (jt), the Monte Carlo permutation test (mc-jt) and the sequential Monte Carlo permutation test (smc-jt) are implemented.

Available conditional independence tests (and the respective labels) for Gaussian Bayesian networks (normal variables) are:

• linear correlation: Pearson's linear correlation. The exact Student's t test (cor), the Monte Carlo permutation test (mc-cor) and the sequential Monte Carlo permutation test (smc-cor) are implemented.

  Hotelling H (1953). "New Light on the Correlation Coefficient and its Transforms". Journal of the Royal Statistical Society: Series B, 15(2):193–225.

• Fisher's Z: a transformation of the linear correlation with asymptotic normal distribution. The asymptotic normal test (zf), the Monte Carlo permutation test (mc-zf) and the sequential Monte Carlo permutation test (smc-zf) are implemented.

• mutual information: an information-theoretic distance measure. Again it is proportional to the log-likelihood ratio (they differ by a $2n$ factor). The asymptotic $\chi^2$ test (mi-g), the Monte Carlo permutation test (mc-mi-g) and the sequential Monte Carlo permutation test (smc-mi-g) are implemented.

• shrinkage estimator for the mutual information (mi-g-sh): an improved asymptotic $\chi^2$ test based on the James-Stein estimator for the mutual information.

  Ledoit O, Wolf M (2003). "Improved Estimation of the Covariance Matrix of Stock Returns with an Application to Portfolio Selection". Journal of Empirical Finance, 10:603–621.

Available conditional independence tests (and the respective labels) for hybrid Bayesian networks (mixed discrete and normal variables) are:

• mutual information: an information-theoretic distance measure. Again it is proportional to the log-likelihood ratio (they differ by a $2n$ factor). Only the asymptotic $\chi^2$ test (mi-cg) is implemented.

Description

Compute Shannon's entropy of a fitted Bayesian network and the Kullback-Leibler divergence between two fitted Bayesian networks.

Usage

H(P)
KL(P, Q)

Arguments

Value

H() and KL() return a single numeric value.

Note

Note that in the case of Gaussian and conditional Gaussian netwoks the divergence can be negative. Regardless of the type of network, if at least one of the two networks is singular the divergence can be infinite.

If any of the parameters of the two networks are NAs, the divergence will also be NA.

Author(s)

Marco Scutari

Examples

## Not run: 
# discrete networks
dag = model2network("[A][C][F][B|A][D|A:C][E|B:F]")
fitted1 = bn.fit(dag, learning.test, method = "mle")
fitted2 = bn.fit(dag, learning.test, method = "bayes", iss = 20)

H(fitted1)
H(fitted2)

KL(fitted1, fitted1)
KL(fitted2, fitted2)
KL(fitted1, fitted2)

## End(Not run)

# continuous, singular networks.
dag = model2network("[A][B][E][G][C|A:B][D|B][F|A:D:E:G]")
singular = fitted1 = bn.fit(dag, gaussian.test)
singular$A = list(coef = coef(fitted1[["A"]]) + runif(1), sd = 0)

H(singular)
H(fitted1)

KL(singular, fitted1)
KL(fitted1, singular)

Description

Insurance is a network for evaluating car insurance risks.

Usage

data(insurance)

Format

The insurance data set contains the following 27 variables:

• GoodStudent (good student): a two-level factor with levels False and True.

• Age (age): a three-level factor with levels Adolescent, Adult and Senior.

• SocioEcon (socio-economic status): a four-level factor with levels Prole, Middle, UpperMiddle and Wealthy.

• RiskAversion (risk aversion): a four-level factor with levels Psychopath, Adventurous, Normal and Cautious.

• VehicleYear (vehicle age): a two-level factor with levels Current and older.

• ThisCarDam (damage to this car): a four-level factor with levels None, Mild, Moderate and Severe.

• RuggedAuto (ruggedness of the car): a three-level factor with levels EggShell, Football and Tank.

• Accident (severity of the accident): a four-level factor with levels None, Mild, Moderate and Severe.

• MakeModel (car's model): a five-level factor with levels SportsCar, Economy, FamilySedan, Luxury and SuperLuxury.

• DrivQuality (driving quality): a three-level factor with levels Poor, Normal and Excellent.

• Mileage (mileage): a four-level factor with levels FiveThou, TwentyThou, FiftyThou and Domino.

• Antilock (ABS): a two-level factor with levels False and True.

• DrivingSkill (driving skill): a three-level factor with levels SubStandard, Normal and Expert.

• SeniorTrain (senior training): a two-level factor with levels False and True.

• ThisCarCost (costs for the insured car): a four-level factor with levels Thousand, TenThou, HundredThou and Million.

• Theft (theft): a two-level factor with levels False and True.

• CarValue (value of the car): a five-level factor with levels FiveThou, TenThou, TwentyThou, FiftyThou and Million.

• HomeBase (neighbourhood type): a four-level factor with levels Secure, City, Suburb and Rural.

• AntiTheft (anti-theft system): a two-level factor with levels False and True.

• PropCost (ratio of the cost for the two cars): a four-level factor with levels Thousand, TenThou, HundredThou and Million.

• OtherCarCost (costs for the other car): a four-level factor with levels Thousand, TenThou, HundredThou and Million.

• OtherCar (other cars involved in the accident): a two-level factor with levels False and True.

• MedCost (cost of the medical treatment): a four-level factor with levels Thousand, TenThou, HundredThou and Million.

• Cushioning (cushioning): a four-level factor with levels Poor, Fair, Good and Excellent.

• Airbag (airbag): a two-level factor with levels False and True.

• ILiCost (inspection cost): a four-level factor with levels Thousand, TenThou, HundredThou and Million.

• DrivHist (driving history): a three-level factor with levels Zero, One and Many.

Note

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

Source

Binder J, Koller D, Russell S, Kanazawa K (1997). "Adaptive Probabilistic Networks with Hidden Variables". Machine Learning, 29(2–3):213–244.

Examples

# load the data.
data(insurance)
# create and plot the network structure.
modelstring = paste0("[Age][Mileage][SocioEcon|Age][GoodStudent|Age:SocioEcon]",
  "[RiskAversion|Age:SocioEcon][OtherCar|SocioEcon][VehicleYear|SocioEcon:RiskAversion]",
  "[MakeModel|SocioEcon:RiskAversion][SeniorTrain|Age:RiskAversion]",
  "[HomeBase|SocioEcon:RiskAversion][AntiTheft|SocioEcon:RiskAversion]",
  "[RuggedAuto|VehicleYear:MakeModel][Antilock|VehicleYear:MakeModel]",
  "[DrivingSkill|Age:SeniorTrain][CarValue|VehicleYear:MakeModel:Mileage]",
  "[Airbag|VehicleYear:MakeModel][DrivQuality|RiskAversion:DrivingSkill]",
  "[Theft|CarValue:HomeBase:AntiTheft][Cushioning|RuggedAuto:Airbag]",
  "[DrivHist|RiskAversion:DrivingSkill][Accident|DrivQuality:Mileage:Antilock]",
  "[ThisCarDam|RuggedAuto:Accident][OtherCarCost|RuggedAuto:Accident]",
  "[MedCost|Age:Accident:Cushioning][ILiCost|Accident]",
  "[ThisCarCost|ThisCarDam:Theft:CarValue][PropCost|ThisCarCost:OtherCarCost]")
dag = model2network(modelstring)
## Not run: graphviz.plot(dag, shape = "ellipse")