Workflow of the MoTBFs package

The MoTBFs package is designed using S3 objects. The package implements functions for learning univariate, multidimensional, and conditional distributions, and provides support for parameter learning in hybrid Bayesian networks. In addition, it includes functions for incorporating prior knowledge when there is lack of data and for carrying out probabilistic inference. Moreover, two classes are incorporated in the package, motbf for defining univariate mixtures of truncated basis functions and jointmotbf for specifying multidimensional MoTBFs.

The functionality of the MoTBFs package is illustrated through an analysis carried out on a real world dataset. More precisely, we use the ecoli dataset [1], which is provided along with the package. The dataset contains information about Escherichia coli and consists of n=336 records, 8 input variables, and 1 output variable (the class). It is a bacterium of the genus Escherichia that is commonly found in the lower intestine of warm-blooded organisms. This dataset can be downloaded from http://archive.ics.uci.edu/ml/datasets/Ecoli.

How to install MoTBFs

The MoTBFs can be installed from CRAN, using the usual install.packages() function.

## Install and load the MoTBFs package
 install.packages("MoTBFs")
 library("MoTBFs")
## Load ecoli dataset
 data("ecoli", package = "MoTBFs")

The example dataset

The ecoli dataset is a data frame with 336 rows corresponding to proteins and 9 columns corresponding to variables. The dataset contains 4 discrete variables, stored as characters, and 5 continuous variables. The variables provide measurements of the cells used for predicting the localization site of proteins. The first variable, Sequence.Name, which is the accession number for the SWISS-PROT database, and the output variable class will not be used in this running example, and we will therefore remove them from the data frame. The discrete variables lip and chg are binary attributes, where character numbers are used as states; "0.48" and "1", and "0.5" and "1", respectively.

## Drop the first and last variables of the ecoli dataset 
  data <- ecoli[,-c(1,9)]

Split the dataset into train and test

For validation purposes, the dataset is split into a training and a test set using the TrainingandTestData() function.

## Split the dataset into train and test subsets
  set.seed(2)
  dataTT <- TrainingandTestData(data, percentage_test = 0.2)
  trainingData <- dataTT$Training
  testData <- dataTT$Test

The seed value determines the partitioning of the data into training and test, and is therefore key to reproducing the experiments. From now on, we will carry out all the analyses on the training data, leaving the test dataset for estimating the predictive capabilities of the learned models.

Learn the directed acyclic graph

Our illustrative example basically consists of fitting MoTBF densities to a previously learned Bayesian network structure over the variables in the dataset. The structure can, for instance, be obtained, using the function hc() from the bnlearn package. This function returns a directed acyclic graph obtained from the dataset using a local search method. For the sake of simplicity, we have included the function LearningHC() in our package, which automatically converts into factors those columns that are non-numeric, before calling the function hc() in bnlearn. LearningHC() can also be used to discretize the dataset before calling hc(), but we are not using this functionality in the running example.

## Learn the structure of the Bayesian network using the training data
  dag <- LearningHC(trainingData)
  dag
## 
##   Bayesian network learned via Score-based methods
## 
##   model:
##    [lip][alm1][mcg|lip:alm1][chg|lip][aac|alm1][gvh|mcg][alm2|gvh:lip:alm1] 
##   nodes:                                 7 
##   arcs:                                  8 
##     undirected arcs:                     0 
##     directed arcs:                       8 
##   average markov blanket size:           3.14 
##   average neighbourhood size:            2.29 
##   average branching factor:              1.14 
## 
##   learning algorithm:                    Hill-Climbing 
##   score:                                 BIC (cond. Gauss.) 
##   penalization coefficient:              2.797356 
##   tests used in the learning procedure:  102 
##   optimized:                             TRUE

We can visualize the network structure using the plot() generic function or the graphviz.plot() function of bnlearn package.

# Visualize the network structure in topological order
bnlearn::graphviz.plot(dag)

Before describing how to learn the MoTBF distributions associated with the network structure, we first present the basic functionality for learning different types of MoTBF representations, i.e., univariate, conditional, and joint MoTBF densities.

Univariate MoTBFs densities

We illustrate the learning of a univariate MoTBF density by considering the continuous variable mcg.

## Learn the density of variable mcg, using MTEs or MOP as basis functions
 f1 <- univMoTBF(trainingData[,1], POTENTIAL_TYPE = "MTE", nparam = 13)
 f2 <- univMoTBF(trainingData[,1], POTENTIAL_TYPE = "MOP", nparam = 11)

The function takes two mandatory arguments, data and POTENTIAL_TYPE, where the latter can either be "MOP" or "MTE" if polynomial or exponential basis functions should be used, respectively. univMoTBF() also accepts optional arguments: it is possible to specify the domain over which the model will be fitted, evalRange, the exact number of basis functions to be used, nparam, and the maximum number of parameters in the function, maxParam, which selects the best fit using the log-likelihood score. If nparam or maxParam are not given, then the Bayesian information criterion (BIC) [2] is used for scoring and function selection: it evaluates the two next functions and if the BIC value does not improve then the function with the best BIC score so far is returned.

The mathematical expression of the univariate density is shown via print(), whereas other hidden elements related to the learning task can be obtained using summary().

print(f2)
## [1] 0.00528468532660471+30.7730282255632*x-1125.05786889453*x^2+16034.7596391062*x^3-118127.908988612*x^4+516328.233596937*x^5-1402763.10738281*x^6+2378343.85773264*x^7-2437579.93461867*x^8+1377857.84608826*x^9-329159.020877368*x^10
summary(f2)
## 
##  MoTBFs FOR UNIVARIATE DISTRIBUTIONS 
## 
##  Model:
##  0.00528468532660471+30.7730282255632*x-1125.05786889453*x^2+16034.7596391062*x^3-118127.908988612*x^4+516328.233596937*x^5-1402763.10738281*x^6+2378343.85773264*x^7-2437579.93461867*x^8+1377857.84608826*x^9-329159.020877368*x^10 
## 
##  Class: univmotbf mop motbf
##  Subclass: mop 
## 
##  Coefficients:
##  0.005284685 30.77303 -1125.058 16034.76 -118127.9 516328.2 -1402763 2378344 -2437580 1377858 -329159 
## 
##  Domain:
##  (0, 0.89)

The object returned by univMoTBF() is a list of classes "univmotbf", "motbf", and either "mop" or "mte", depending on the basis functions used. The object returned contais several elements, including its mathematical expression and other hidden elements related to the learning task.

The learned densities can be plotted using the generic function plot(). The next figure shows the the model fits, provided by univMoTBF(), with blue dashed line for MOPs and red solid line for MTEs overlaying the histogram of the training data of the mcg variable.

## Plot the densities f1 and f2 over the histogram of variable mcg
  hist(trainingData[,1], prob = TRUE , main = "", xlab = "X")
  plot(f1, xlim = range(trainingData[,1]), col = "red", add = TRUE, lwd = 3)
  plot(f2, xlim = range(trainingData[,1]), col = "blue", add = TRUE, lwd = 3, lty = 2)
  legend("topleft", legend = c("MTE", "MOP"),col = c("red", "blue"), lty = 1:2, lwd = 3, inset = c(0, -0.5), xpd = TRUE)

To evaluate the predictive ability of the models we use the generic method as.function() developed for the "motbf" class to get the log-likelihood as well as BICMoTBF() to obtain the BIC score.

## Compute log-likelihood of the fitted densities
  sum(log(as.function(f1)(testData[,1])))
## [1] 12.32783
  sum(log(as.function(f2)(testData[,1])))
## [1] 10.51058
## Compute BIC score of the fitted densities  
  BICMoTBF(f1,testData[,1])
## [1] -17.10502
  BICMoTBF(f2,testData[,1])
## [1] -14.71757

An alternative way to visually check the goodness of fit of the estimated models is to simulate a data sample from the learned functions and compare it with the training data. For doing this, we use the inverse transform method, a technique for generating random samples from a specific probability distribution based on evaluating the inverse of the CDF on a uniform random number, yielding a value for the random variable being sampled. This is done by function rMoTBF(). For the sake of reproducibility, we fix the seed for the random numbers to be used by the rMoTBF() function, which is set to 5 in this example. In the next code snippet, the previous function fitted with a polynomial basis, f2, will be used.

## Simulate data from the estimated density f2
  set.seed(5)
  X <- rMoTBF(size = 400, fx = f2)
  
## Test whether or not the simulated sample and the observed data come from the same distribution
  ks.test(trainingData[,1], X)
## Warning in ks.test.default(trainingData[, 1], X): p-value will be approximate
## in the presence of ties
## 
##  Asymptotic two-sample Kolmogorov-Smirnov test
## 
## data:  trainingData[, 1] and X
## D = 0.065167, p-value = 0.5018
## alternative hypothesis: two-sided

In this example the two-sample Kolmogorov-Smirnov test is used. The p-value is notably above 0.05, so there is no evidence to reject the null hypothesis that both samples are drawn from the same population.

We can plot the histogram and the empirical cumulative distribution of both the training data of variable mcg and the sample simulated from the distribution learned from the same data, in order to compare both distributions.

## Compare the histogram of both distributions
 hist(trainingData[,1], prob = TRUE, col = "#FFC107", density = 10, angle = -45, main = "")
 hist(X, prob = TRUE, col = "#0C7BDC", main = "", ylim = c(0,2.2), cex.lab=1.5, cex.axis=1.5, density = 10,add = T)
 legend('topleft', legend = c("Training data", "Simulated data"), fill = c("#FFC107", "#0C7BDC"), density = 20, angle = c(-45, 45), inset = c(0, -0.5), xpd = TRUE)

## Compare the the CDF of both distributions
 plot(ecdf(trainingData[,1]), cex = 0, lwd = 3 , cex.lab = 1.5, cex.axis = 1.5, main = "")
 plot(integrate.motbf(f2), xlim = range(trainingData[,1]), col ="red", lwd = 3, add = TRUE)
 legend('topleft', legend = c("Training data", "Simulated data"), col = c("black", "red"), lwd = 3, inset = c(0, -0.5), xpd = TRUE)

We can also manipulate the distributions with a collection of methods for objects of class "motbf". Here is an example of the use of three of them, coef(), integrate.motbf(), and derivMoTBF().

## Compute the derivative and integral of the fitted density

  # Coefficients of the fitted density
  coef(f2)
##  [1]  5.284685e-03  3.077303e+01 -1.125058e+03  1.603476e+04 -1.181279e+05
##  [6]  5.163282e+05 -1.402763e+06  2.378344e+06 -2.437580e+06  1.377858e+06
## [11] -3.291590e+05
  # Indefinite integral of the fitted density
  integrate.motbf(f2)
## [1] 0.00528468532660471*x+15.3865141127816*x^2-375.01928963151*x^3+4008.68990977655*x^4-23625.5817977224*x^5+86054.7055994895*x^6-200394.729626116*x^7+297292.98221658*x^8-270842.21495763*x^9+137785.784608826*x^10-29923.547352488*x^11
  # Definite integral of the fitted density
  integrate.motbf(f2, lower = min(trainingData[,1]), upper = max(trainingData[,1]))
## [1] 1
  # Derivarive of the fited density
  derivMoTBF(f2)
## [1] 30.7730282255632-2250.11573778906*x+48104.2789173186*x^2-472511.635954448*x^3+2581641.16798468*x^4-8416578.64429686*x^5+16648407.0041285*x^6-19500639.4769494*x^7+12400720.6147943*x^8-3291590.20877368*x^9

Joint MoTBFs densities

The learning process for multidimensional variables is similar to the previous one. The function jointmotbf.fit() is used to solve the quadratic optimization problem and returns the analytical expression of the joint density. The returned object is of class "jointmotbf" and "motbf". The expression is the only visible element, while the others can be retrieved using attributes(). In this example only two variables are used, mcg and alm1, in order to be able to plot the results.

## Learn joint distributions
  P = jointmotbf.fit(X = trainingData[,c("mcg", "alm1")], dimensions = c(5,5))
  
  attributes(P)
## $names
## [1] "Function"   "Domain"     "Iterations" "Time"      
## 
## $class
## [1] "jointmotbf" "motbf"

The function print() can be used to obtain an expression of the learned joint density, while summary() yields a more thorough excerpt of the "jointmotbf" object.

print(P)
## [1] 1.00000000044876e-05-1.87890751245811e-13*alm1+1.34925514946196e-12*alm1^2-2.45341144148194e-12*alm1^3+1.31231637960609e-12*alm1^4-2.7069010790699*mcg+103.459764181587*mcg*alm1-461.084363497177*mcg*alm1^2+679.358798389865*mcg*alm1^3-319.267637943894*mcg*alm1^4-7.35251619834014*mcg^2+244.974957258816*mcg^2*alm1+38.6208995140072*mcg^2*alm1^2-1193.96736630721*mcg^2*alm1^3+920.72827509131*mcg^2*alm1^4+30.9216352421374*mcg^3-1102.42439690391*mcg^3*alm1+2410.7579267411*mcg^3*alm1^2-668.029611891877*mcg^3*alm1^3-677.36983290076*mcg^3*alm1^4-21.621357506081*mcg^4+782.648440776334*mcg^4*alm1-2103.42524400634*mcg^4*alm1^2+1294.26742623164*mcg^4*alm1^3+51.5825770400239*mcg^4*alm1^4
summary(P)
## 
##  MoTBFs FOR MULTIVARIATE DISTRIBUTIONS 
## 
##  Model:
##  1.00000000044876e-05-1.87890751245811e-13*alm1+1.34925514946196e-12*alm1^2-2.45341144148194e-12*alm1^3+1.31231637960609e-12*alm1^4-2.7069010790699*mcg+103.459764181587*mcg*alm1-461.084363497177*mcg*alm1^2+679.358798389865*mcg*alm1^3-319.267637943894*mcg*alm1^4-7.35251619834014*mcg^2+244.974957258816*mcg^2*alm1+38.6208995140072*mcg^2*alm1^2-1193.96736630721*mcg^2*alm1^3+920.72827509131*mcg^2*alm1^4+30.9216352421374*mcg^3-1102.42439690391*mcg^3*alm1+2410.7579267411*mcg^3*alm1^2-668.029611891877*mcg^3*alm1^3-677.36983290076*mcg^3*alm1^4-21.621357506081*mcg^4+782.648440776334*mcg^4*alm1-2103.42524400634*mcg^4*alm1^2+1294.26742623164*mcg^4*alm1^3+51.5825770400239*mcg^4*alm1^4 
## 
##  Class: jointmotbf motbf 
## 
##  Coefficients:
##  1e-05 -1.878908e-13 1.349255e-12 -2.453411e-12 1.312316e-12 -2.706901 103.4598 -461.0844 679.3588 -319.2676 -7.352516 244.975 38.6209 -1193.967 920.7283 30.92164 -1102.424 2410.758 -668.0296 -677.3698 -21.62136 782.6484 -2103.425 1294.267 51.58258 
## 
##  Domain mcg:
##  (0, 0.89)
##  Domain alm1:
##  (0.03, 1)
## 
##  Number of Iterations: 51 
## 
##  Processing Time: 0.002694845 secs

The processing time, P$Time, can vary depending on the CPU, but the learning outcome will always be the same for a specific data sample.

The generic function plot() can be used for objects of class "jointmotbf". This function accepts optional arguments such as type, where one can choose between "perspective" and "contour", ranges, used to specify the plotting range, orientation, which indicates the orientation of the perspective graph, and filled for getting a filled contour plot.

## Plot the joint distribution of 2 variables
 par(mar=c(2,3,2,2))
 
 # Filled contour
 plot(P, data = trainingData[,c(1,6)])

 # Simple contour
 plot(P, data = trainingData[,c(1,6)], filled = FALSE, cex.lab = 2, cex.axis = 1.85, lwd = 1.5) 

 # Perspective
 plot(P, type = "perspective", data = trainingData[,c(1,6)], orientation=c(25,25), cex.lab = 2,xaxs = "i")

The marginal.jointmotbf() function computes the marginals of joint densities. In this example we have two variables, so there are two marginal densities. The marginal variable can be specify by the index or name.

marginal.jointmotbf(P, var = "mcg")
## [1] 9.70000000934267e-06+1.3530290672647*mcg+13.772635514603*mcg^2-19.6390735145559*mcg^3+2.75960426007185*mcg^4
marginal.jointmotbf(P, var = 2)
## [1] -0.364291790424928+13.0274885842294*alm1-30.3090949604406*alm1^2+28.2511964413879*alm1^3-10.5727672738866*alm1^4

Conditional MoTBFs densities

The next step in our analysis is learning conditional densities, which is implemented by the function conditionalMethod(). Five of its arguments are compulsory: data, the dataset; nameParents, a character vector indicating the name of the parents; nameChild, a character string containing the name of the child; numIntervals, the maximum number of intervals for splitting the domain of the parent variables; POTENTIAL_TYPE, the type of basis function. Other arguments are optional, like maxParam, indicating the maximum number of parameters for each function, and s, the expert’s relative confidence in any prior knowledge, and priorData if prior knowledge is incorporated in the analysis.

We will do the conditional analysis for only two variables in order to be able to make a 2-dimensional plot of the obtained results. For example, taking into account the relationship found by the dag, we consider the child variable gvh with parent variable mcg.

## Learn conditional distributions
P <- conditionalMethod(trainingData, nameParents = "mcg", nameChild = "gvh",
                         numIntervals = 5, POTENTIAL_TYPE ="MOP", scale = FALSE)
printConditional(P)
## Parent: mcg       Range: 0 < mcg < 0.44 
## [1] 115.870454614414-1783.89438852081*x+10688.6518417108*x^2-32342.1748450932*x^3+54497.1685860695*x^4-52033.2595768611*x^5+26393.6468514034*x^6-5536.00792332333*x^7
## Parent: mcg       Range: 0.44 < mcg < 0.89 
## [1] -37.3007586075146+733.084458130356*x-5676.35686301634*x^2+22283.2019747918*x^3-47834.417884372*x^4+56907.8390342608*x^5-35242.6288447815*x^6+8867.52830777248*x^7

It can be noticed that the learning algorithm decides to split the domain of the parent into two intervals even though we have set the argument numIntervals to five. This is because the BIC score is not improved any further by splitting the domain into more than two intervals.

The resulting conditional density (a MOP in this case) can be plotted using plotConditional(). The sample points can be overlaid by setting the argument points to TRUE.

par(mar=c(2,3,2,2))
## Plot the conditional density of gvh given mcg
plotConditional(P, data = trainingData, nameChild = "gvh", points = TRUE)

MoTBF distributions associated with the network structure

The last step is to learn the distributions tied to the Bayesian network learned previously. For doing this task, the motbf.fit() function of the MoTBFs package is used. The graph is a mandatory argument, that can be of class "bn", "graphNEL" or "network". Other mandatory arguments are the data, the maximum number of intervals for splitting the domain of the parents (numIntervals), and the type of basis function (POTENTIAL_TYPE). The function also accepts additional arguments, but they are not listed here.

In the example, the DAG was obtained using the bnlearn package and therefore it is an object of class "bn". As an example, we will use a maximum of 4 intervals and "MTE" potentials when learning the densities (i.e. exponential basis functions).

## Learn the distributions of a Bayesian network
bn <- motbf.fit(dag, data = trainingData, numIntervals = 4, POTENTIAL_TYPE = "MTE")

The returned object is of class "motbf_fit" and "motbf". The results are reported using the generic print() function for objects of class "motbf_fit".

print(bn)
## Potential(mcg | alm1, lip)
## Parent: alm1      Range: 0.03 < alm1 < 0.45 
## Parent: lip       Value =  "0.48" 
## [1] 1.12359550561798+0*exp(5.28468530956588*mcg)
## Parent: alm1      Range: 0.03 < alm1 < 0.45 
## Parent: lip       Value =  "1" 
## [1] 1.12359550561798+0*exp(1.05693706191318*mcg)
## Parent: alm1      Range: 0.45 < alm1 < 1 
## Parent: lip       Value =  "0.48" 
## [1] 0.258507227216329+0.125038108071494*exp(5.28468530956588*mcg)-0.320433056518257*exp(-5.28468530956588*mcg)-0.00221027286466645*exp(10.5693706191318*mcg)-2.61527685677068*exp(-10.5693706191318*mcg)+9.86301836901382e-06*exp(15.8540559286976*mcg)+3.17849091704381*exp(-15.8540559286976*mcg)
## Parent: alm1      Range: 0.45 < alm1 < 1 
## Parent: lip       Value =  "1" 
## [1] 0.556915717393733+0.0796377643379347*exp(5.28468530956588*mcg)-2.98142905617842*exp(-5.28468530956588*mcg)-0.00076501787356927*exp(10.5693706191318*mcg)+3.13754385763789*exp(-10.5693706191318*mcg)
## 
## Potential(gvh | mcg)
## Parent: mcg       Range: 0 < mcg < 0.89 
## [1] 1.19047619047619+0*exp(6.81739017415969*gvh)
## 
## Potential(lip)
##       0.48          1
##  0.9630996 0.03690037
## 
## Potential(chg | lip)
##      lip
## chg          0.48         1
##   0.5 0.996183206 0.8181818
##   1   0.003816794 0.1818182
## 
## Potential(aac | alm1)
## Parent: alm1      Range: 0.03 < alm1 < 1 
## [1] 1.13636363636363+0*exp(8.10721368061671*aac)
## 
## Potential(alm1)
## [1] 1.03092783505155+0*exp(4.57129025985515*alm1)
## 
## Potential(alm2 | alm1, gvh, lip)
## Parent: alm1      Range: 0.03 < alm1 < 0.71 
## Parent: gvh       Range: 0.16 < gvh < 1 
## Parent: lip       Value =  "0.48" 
## [1] 1.01010101010101+0*exp(4.68957842101153*alm2)
## Parent: alm1      Range: 0.03 < alm1 < 0.71 
## Parent: gvh       Range: 0.16 < gvh < 1 
## Parent: lip       Value =  "1" 
## [1] 1.01010101010101+0*exp(4.68957842101153*alm2)
## Parent: alm1      Range: 0.71 < alm1 < 1 
## Parent: gvh       Range: 0.16 < gvh < 1 
## Parent: lip       Value =  "0.48" 
## [1] -1.66122336611897+0.0694678523633136*exp(4.68957842101153*alm2)+9.76498160863297*exp(-4.68957842101153*alm2)+0.00526173195933795*exp(9.37915684202306*alm2)-23.1570734636348*exp(-9.37915684202306*alm2)-0.000122314620146376*exp(14.0687352630346*alm2)+24.023608781198*exp(-14.0687352630346*alm2)+6.43439431827203e-07*exp(18.7583136840461*alm2)-9.04021189479816*exp(-18.7583136840461*alm2)
## Parent: alm1      Range: 0.71 < alm1 < 1 
## Parent: gvh       Range: 0.16 < gvh < 1 
## Parent: lip       Value =  "1" 
## [1] 0.0625235874946527+0.0437639725805153*exp(4.68957842101153*alm2)-0.101597981654157*exp(-4.68957842101153*alm2)

Notice how nodes in the DAG with only discrete parents contain as many functions as configurations of the parents, whereas nodes that have continuous parents have at most 4 functions for each parent and, finally, nodes that have mixed parents contain as many functions as configurations of the discrete parents times the number of regions into which the domain of the continuous parents is split.

The BIC criterion is used to decide the number of splitting points of the domain of the continuous parent nodes and to choose the number of basis functions used. The function BiC.MoTBFBN() can be used to compute the log-likelihood and the BIC score of a dataset given the Bayesian network.

BiC.MoTBFBN(bn, data = testData)
## $LogLikelihood
## [1] 32.09624
## 
## $BIC
## [1] -83.5328

Learn MoTBF distributions in a full hybrid network using prior knowledge

We will now exemplify the use of prior knowledge in the learning process. In order to illustrate the approach, we first select a small subset of the Ecoli dataset using TrainingandTestData(). In the next example the percentage of the test data is 99%, which means the training data is only 1% of the full dataset.

## Obtain small training subset
set.seed(4)
dataTT <- TrainingandTestData(data, percentage_test = 0.99)

trainingData <- dataTT$Training
testData <- dataTT$Test
nrow(trainingData)
## [1] 13

There are 13 entries in the training dataset. We are going to fit MoTBFs with and without prior information.

Learning univariate and conditional distributions and Bayesian networks can be done using the functions learnMoTBFpriorInformation() and motbf_fit(). The arguments for these functions are the same as previously explained and, in addition, it is necessary to specify the expert confidence in the prior knowledge, s, and the prior dataset priorData. On the one hand, to generate an artificial prior dataset the generateNormalPriorData() function can be used.

## Generate artificial prior dataset
means <- sapply(data, function(x){ifelse(is.numeric(x), mean(x),NA)})
set.seed(4)
priorData <- generateNormalPriorData(dag, data = trainingData, size = 5000, means = means)

On the other hand, argument s takes values on the interval [0,N], where N is the sample size, and is used to synchronize the support of the prior knowledge and the sample. We refer the reader to [3] for the details. In this example we will use the aac variable from the data set, have s = 5 as confidence level, and set "MOP" as potential type.

## Learn univariate distribution using prior information
f <- learnMoTBFpriorInformation(priorData$aac, trainingData$aac, s = 5, POTENTIAL_TYPE = "MOP", returnAll = TRUE)

print(f)  
## $coeffs
## [1] 0.7427962 0.2572038
## 
## $posteriorFunction
## [1] 0.312389351880304+1.09772806577837*x+0.157157947725107*x^2+3.54129103556391*x^3+227.08274180439*x^4-945.748457762723*x^5+1417.37276145449*x^6-931.480057505286*x^7+227.99099681937*x^8
## 
## $priorFunction
## [1] 0.111457607131215+1.47783208318457*x+0.211576131207797*x^2+4.76751362327549*x^3+305.713383704084*x^4-1273.22736575295*x^5+1908.15831899506*x^6-1254.01832816566*x^7+306.936134987115*x^8
## 
## $dataFunction
## [1] 0.892673709362328+0*x
## 
## $domain
## [1] -0.1397527  0.9804775

learnMoTBFpriorInformation() returns the fitted model using the training data only ($dataFunction), the model fit using the prior data only ($priorFunction) and the model fit that combines the training data and the prior data ($posteriorFunction). The three univariate densities can be plotted using the generic method plot().

plot(f$posteriorFunction, xlim = f$domain, ylim = c(0,2.1), lwd = 3)
plot(f$dataFunction, xlim = f$domain, add = TRUE, col = 2, lwd = 3, lty = 2)
plot(f$priorFunction, xlim = f$domain, add = TRUE, col = 4, lwd = 3, lty = 3)
legend("topleft", legend = c("Posterior", "Data", "Prior"), col = c(1,2,4), lty = 1:3, lwd = 3, inset = c(0, -0.7), xpd = TRUE)

## Log-likelihood of the model that uses prior information
  sum(log(as.function(f$posteriorFunction)(testData$aac)))
## [1] 131.5898
## Log-likelihood of the model that does not use prior information
  sum(log(as.function(f$dataFunction)(testData$aac)))
## [1] -36.67153

The best model, taking into account the log-likelihood, is the MoTBF which uses the prior data, f$posteriorFunction.

The last step is to incorporate the prior knowledge in the full Bayesian network. For this analysis we are not going to print out the results (which could be done using the generic function print()), because the structure is similar to the previous Bayesian network representations. As an example, we will use numIntervals = 2, POTENTIAL_TYPE = "MOP", and s = 5.

## Fit Bayesian network using prior information
priorBN <- motbf.fit(dag, trainingData, numIntervals = 2,
                           POTENTIAL_TYPE = "MOP", s = 5, priorData = priorData)


## Fit BN without using prior information
BN <- motbf.fit(dag, trainingData, numIntervals = 2,
                           POTENTIAL_TYPE = "MOP")

# Compute log-likelihood 
logLikelihood.MoTBFBN(priorBN, data = testData)
## [1] -35.64878
logLikelihood.MoTBFBN(BN, data = testData)
## [1] -22.20795

Looking at the log-likelihood corresponding to the network with and without prior data, we can see that, in this example, incorporating prior knowledge is better when data is scarce.

Inference

After a Bayesian network has been constructed, the MoTBFs package can be used to obtain the conditional density of any variable in the network given that some other variables have been observed. The conditional distribution is obtained by forward sampling. As an example, consider a network estimated from the ecoli dataset:

## Learn a Bayesian network
dag <- LearningHC(data)

bn <- motbf.fit(dag, data = data, numIntervals = 4,
                         POTENTIAL_TYPE = "MOP")

The observed values are specified using a data frame. We can obtain an approximation of the posterior probability distribution of a target variable given a set of observed variables using the get_approx_posterior() function, which runs the forward sampling algorithm if the evidence parameter is NULL, or the likelihood weighting algorithm otherwise.

In the example, we are assuming that we want to compute the conditional density of alm2 given that lip="0.48", alm1 = 0.55 and gvh = 0.9. We set the number of random samples to generate to size = 100.

# Specify the evidence set and target variable
obs <- data.frame(lip = "0.48", alm1 = 0.55, gvh = 0.9, stringsAsFactors=FALSE)
node <- "alm2"

# Get the conditional distribution of 'node' and the generated sample
set.seed(4)
ap_post = get_approx_posterior(bn, target = node, 
                     evidence = obs, size = 100, maxParam = 8)
## Processing Time: 1.73847889900208secs
ap_post$fx
## [1] 6.81076602642233-146.932848695937*x+1154.17777619494*x^2-4092.61854429109*x^3+7295.63939509001*x^4-6344.64924756399*x^5+2133.58529450247*x^6

The output consists of the posterior density and the sample from which the density parameters were estimated.

The posterior can also be computed using exact inference. The variable elimination algorithm is implemented for MOP distributions only. We will use the same query, applying the exact solution.

ex_post = variableElimination(bn, target = node, evidence = obs)
# Plot the posterior distribution obtained with each solution
plot(ap_post$fx, col = "red") 
plot(ex_post, add = T, col = "blue")
legend("topleft", legend = c("Approximate", "Exact"), col = c("red", "blue"), lwd = 1,  inset = c(0, -0.5), xpd = TRUE)

Structural learning

The MoTBFs package implements the Tree Augmented Naive Bayes model for MOP distributions, based on the Chow-Liu algorithm. The main function is fit_tan(), whose mandatory arguments are "target" (the class variable) and "data" (the dataset). Note that "target" might be either discrete or continuos and the variables in "data" can also be of either type.

Some optional arguments are "root" (if not specified, the function chooses the variable with highest mutual information with the target), and "mutualInfoCond" (a matrix containing the conditional mutual information, can be computed using the mutual_information_tan() function; if not given, fit_tan() computes it internally). Moreover, the optional argument "fit.args" allows to specify optional arguments accepted by function motbf.fit(). The returned object of fit_tan() is of class "motbf_fit".

# Build a TAN model for classification
  bn_tan_cl = fit_tan("lip",data)
  bnlearn::graphviz.plot(getDAG(bn_tan_cl))

# Compute mutual information of each variable with mcg. 
# Firstly, discrete variables must be coerced to factor
  data[sapply(data, is.character)] = lapply(data[sapply(data, is.character)], as.factor)
  MI = mutual_information_tan(data,"mcg")
# Build a TAN model for regression
  bn_tan_reg = fit_tan("mcg",data, mutualInfoCond = MI)
  bnlearn::graphviz.plot(getDAG(bn_tan_reg))

On the other hand, the package also contains a wrapper to the TAN implementation (function tree.bayes()) of the bnlearn package, which requires to discretize the continuous variables, and to the hill-climbing algoritm (function hc()). The returned object is of class "bn", i.e., it is just a DAG, which can be used as an argument of the motbf.fit() function to learn the full model. This wrapper is in function getStructure(), whose arguments are "data", "method" (a character string matching either "NB", for naive Bayes model; "TAN"; or "HC", for the hill-climbing algorithm), and "target" (needed only for NB and TAN).

# Naive Bayes structure
nb = getStructure(data, "NB", "mcg")
bnlearn::graphviz.plot(nb)

# TAN wrapper (continuous variables are internally discretized)
tan_disc = getStructure(data, "TAN", "mcg")
bnlearn::graphviz.plot(tan_disc)

# HC wrapper
hc_bnlearn = getStructure(data, "HC")
bnlearn::graphviz.plot(hc_bnlearn)

References

1.
Lichman, M. (2013). UCI machine learning repository. University of California, Irvine, School of Information; Computer Sciences. Retrieved from http://archive.ics.uci.edu/ml
2.
Schwarz, G. (1978). Estimating the dimension of a model. Annals of Statistics, 6, 461–464.
3.
Pérez-Bernabé, I., Fernández, A., Rumí, R., & Salmerón, A. (2016). Parameter learning in hybrid Bayesian networks using prior knowledge. Data Mining and Knowledge Discovery, 30, 576–604.