Hands-on exercise

5.1) Start with 1 parent as maximum limit

To make the search more efficient, we constrain the maximum number of parents allowed for each node. We will start from 1 and increase subsequently until the network score does not improve further even when more parents are allowed.

max.par <- 1 

5.2) Build a cache of all local computations

library(abn)

mycache <- buildscorecache(data.df = as.data.frame(abndata), data.dists = dist, 
                           dag.banned = banned, dag.retained = retain, 
                           max.parents = max.par)

5.3) Run the EXACT SEARCH for the specified parent limit

1. Find the optimal DAG given the parent limit

mydag <- mostprobable(score.cache = mycache)

Step1. completed max alpha_i(S) for all i and S
Total sets g(S) to be evaluated over: 256

2. Fit marginal densities

fabn <- fitabn(dag.m = mydag, data.df = as.data.frame(abndata), data.dists = dist)

3. check network score

fabn$mlik

[1] -2806.499

9.1) Extract parameters from best fitting model and save them for MCMC simulations

The parameters at each node are estimated as posterior probability density distributions. These marginal densities are the ones where JAGS needs to sample from in order to simulate data. In order to use these distributions with JAGS, the densities need to be approximated by a discrete distribution over a fine and equally spaced grid. Here is a mock example for a hypothetical parameter \(\beta\):

# Fit marginal densities over a fixed grid --> n.grid=1000
# --------------------------------------------------------
  marg.f <- fitabn(dag.m = mydag, data.df = as.data.frame(abndata),
                     data.dists = dist, compute.fixed=TRUE, n.grid=1000)

# Extract values 
# --------------
  m <- marg.f$marginals[[1]] 
  for(i in 2: length(marg.f$marginals))
  { m <- c(m, marg.f$marginals[[i]])}
  
# Bind all the marginals for the same node into a matrix
# ------------------------------------------------------
  AR.p <- cbind( m[[ "AR|(Intercept)"]], m[[ "AR|age"]])
  pneumS.p <- cbind( m[[ "pneumS|(Intercept)"]], m[[ "pneumS|age"]])
  female.p <- cbind( m[[ "female|(Intercept)"]])
  livdam.p <- cbind( m[[ "livdam|(Intercept)"]], m[[ "livdam|eggs"]])
  eggs.p <- cbind( m[[ "eggs|(Intercept)"]], m[[ "eggs|adg"]])
  wormCount.p <- cbind( m[[ "wormCount|(Intercept)"]], m[[ "wormCount|AR"]],
                        m[[ "wormCount|eggs"]], m[[ "wormCount|age"]], m[[ "wormCount|adg"]])
  age.p <- cbind( m[[ "age|(Intercept)"]], m[[ "age|female"]])
  prec.age.p <- cbind( m[[ "age|precision" ]])
  adg.p <- cbind( m[[ "adg|(Intercept)"]], m[[ "adg|age"]])
  prec.adg.p <- cbind( m[[ "adg|precision" ]])

# Save it to a file named PostParams to be read by JAGS
# -----------------------------------------------------
  dump(c("AR.p", "pneumS.p", "female.p", "livdam.p", "eggs.p", 
         "wormCount.p", "age.p", "prec.age.p", "adg.p", "prec.adg.p"),
       file="PostParams.R")

9.2) Write the BUG model

In order to simulate the variables of our dataset we need to provide a model for each of them, using the aforementioned parameters estimates. For instance, the binomial node AR in our DAG has one incoming arc coming from the node age. In a regression notation this would be translated into:

logit(AR) = \(\alpha\) + \(\beta\) x age + \(\epsilon\)

where \(\alpha\) is the intercept and \(\beta\) the regression coefficient for variable age and \(\epsilon\) is the error term modelled by a binomial distribution.

Given that we will simulate the data in a Bayesian framework, AR will be modelled as a probability distribution. Therefore it will look like:

AR ~ dbern(probAR);
logit(probAR)<- alpha + beta*age;

Then, the values of alpha and beta will be sampled (dcat) from our discrete distribution of parameters:

alpha.prob ~ dcat(AR.p[ ,2]); –> sample from the vector of density values f(x) (second column in matrix)
alpha ~ AR.p[alpha.prob,1]; –> corresponding x value for the sampled density (first column in matrix)
beta.prob ~ dcat(AR.p[ ,4]);
beta ~ AR.p[beta.prob,1];

The BUG file (model8vPois.bug) can be retrieved from the file directory.

9.3) Run JAGS and inspect the result of a simulated dataset

library(rjags)

# set inits
# ---------
init <- list(".RNG.name"="base::Mersenne-Twister", ".RNG.seed"=42)

# load data
# ---------
source("PostParams.R")

# run model once
# --------------
jj <- jags.model(file = "model8vPois.bug", 
                 data = list(  'AR.p'=AR.p , 'pneumS.p'=pneumS.p , 'female.p'=female.p, 
                               'livdam.p'=livdam.p , 'eggs.p'=eggs.p , 'wormCount.p'=wormCount.p , 
                               'age.p'=age.p ,'prec.age.p'=prec.age.p, 
                               'adg.p'=adg.p , 'prec.adg.p'=prec.adg.p),
                 inits = init,
                 n.chains = 1, 
                 n.adapt = 5000)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 0
   Unobserved stochastic nodes: 28
   Total graph size: 40091

Initializing model

# run more iterations
# -------------------
update(jj, 100000)

# set number of observation we want to extract for a dataset 
# ----------------------------------------------------------
n.obs=341

# sample data (same size as original: 341) with a sampling lag (20) to reduce avoid autocorrelation
# -------------------------------------------------------------------------------------------------
samp <- coda.samples(jj, c("AR", "pneumS", "female", "livdam", "eggs", 
                           "wormCount", "age", "prec.age", "adg", "prec.adg"),
                     n.iter= n.obs*20 , thin =20)

Now compare the simulated data with the original data. Observe that Gaussian nodes are by default centred when doing abn search, meaning that the simulated data for those nodes will be centred as well.

Simulated data looks fairly OK (perhaps wormCount is sub-optimal as the simulated data miss to represent the long right tail) so we can proceed with the bootstrapping.

9.4) Iterate dataset simulation + exact search over and over

What we will do is to create a loop to 1) simulate data, 2) do exact search on such data, and 3) store the best fitting DAG over and over for many times (e.g. 5000 iterations). Bootstrap data need to be saved in a folder to be further inspected.

vars <- colnames(abndata)

# load data for jags
source("PostParams.R")

# select nr. bootstrap samples to run
# -----------------------------------
set.seed(123)

# get 5000 random numbers to set different initial values
# -------------------------------------------------------
n <- sample(1:100000, 5000)


# specify max number of parents based on previous search
# ------------------------------------------------------
max.par <- 4


# Simulate data and run ABN on such dataset
# -----------------------------------------
for (i in 1:length(n)) {

  print(paste("/n Running simulation", i))

  # pick initials
  init <- list(".RNG.name"="base::Mersenne-Twister", ".RNG.seed"=n[i])

  # run model
  jj <- jags.model(file = "model8vPois.bug",
                   data = list(  'AR.p'=AR.p , 'pneumS.p'=pneumS.p , 'female.p'=female.p,
                                 'livdam.p'=livdam.p , 'eggs.p'=eggs.p , 'wormCount.p'=wormCount.p ,
                                 'age.p'=age.p ,'prec.age.p'=prec.age.p,
                                 'adg.p'=adg.p , 'prec.adg.p'=prec.adg.p),
                   inits = init,
                   n.chains = 1,
                   n.adapt = 5000)

  # run more iterations
  update(jj, 100000)

  # sample data (same size as original: 341) with a sampling lag (20) to reduce avoid autocorrelation
  samp <- coda.samples(jj, c("AR", "pneumS", "female", "livdam", "eggs",
                             "wormCount", "age", "prec.age", "adg", "prec.adg"),
                       n.iter= 6820 , thin =20)

  # build dataframe in the same shape as the original one
  dt.boot <- as.data.frame(as.matrix(samp)) # pay attention at order names

  dt.boot<- dt.boot[, vars]


  # now coerce to factors if need be and set levels - NOTES setting levels works as
  # "0" "1" is in the same order as "absent" "present" from original data
  abndata <- as.data.frame(abndata)

  for(j in 1:length(vars)) {
    if(is.factor(abndata[,j])) {
      dt.boot[,j]<-as.factor(dt.boot[,j]);
      levels(dt.boot[,j])<-levels(abndata[,j]);
    }
  }

# Build a cache of all local computations
# ---------------------------------------
 mycache <- buildscorecache(data.df = dt.boot, data.dists = dist, dag.banned = banned,
                             dag.retained = retain, max.parents = max.par)

# Run the EXACT SEARCH
# --------------------
 mp.dag <- mostprobable(score.cache = mycache)
 fabn <- fitabn(dag.m = mp.dag, data.df = dt.boot, data.dists = dist)

# Save the results obtained
# -------------------------
 save(mycache, mp.dag, fabn, dt.boot,  file = sprintf('BootData/dt.boot.%04d.RData',i))

}

9.5) Find the final pruned DAG

First we need to load all the bootstrap DAGs.

# set nr bootstrap samples
# ------------------------
n <- 5000

# load dags and boostrap data
# ------------------------------
dags <- list()
boot <- list()
for (i in 1:n) {
  load(file = sprintf('BootData/dt.boot.%04d.RData',i))
  dags[[i]] <- mp.dag
  boot[[i]] <- dt.boot
  rm(mp.dag, dt.boot, fabn, mycache)
} 

Then we will check what was the most frequent number of arcs:

# count total number of arcs in each dag
# --------------------------------------
arcs <- sapply(dags, sum)
barplot(table(arcs))

For comparison, in the original DAG there were 10 arcs, so it seems that there might have been some over-fitting.

Finally, we will count how many times each arc appears in the bootstrap data. The final pruned DAG will include only arcs present in at least 50% of bootstrap samples (e.g. 2500/5000 in our case).

# Count how many times each arc appear in the bootstrap data
# ----------------------------------------------------------
# function Reduce sums ("+") each element of each maxtrix in the list and store 
# results in a new matrix of same size
alldag <- Reduce("+", dags)  

# express it in percentage
perdag <- alldag/length(dags)


# Keep only arcs that appears in at least 50% of samples
# ------------------------------------------------------
trim.dag <- (alldag >=(n*0.5))*1
# sum(trim.dag)

# Send  final pruned DAG to Graphvis for visualization
# -----------------------------------------------------
tographviz(dag.m = trim.dag, data.df = abndata, data.dists = dist,
           outfile = "TrimDAG.dot")

This is the percentage of arcs retrieved within the bootstrap sample:

          AR pneumS female livdam eggs wormCount age adg
AR         0      0      0      0    0         0  76   5
pneumS     0      0      1      0    0         0  39   8
female     0      0      0      0    0         0   0   0
livdam     0      1      0      0   54         2   0   0
eggs       0      0      0     24    0         0   8  75
wormCount 71      1      0      0  100         0 100  61
age       18     17     56      0    4         0   0  25
adg        1      2     10      0    6         0  75   0

This is how the final pruned DAG looks like:

Now that we have a robust model (encoded by the pruned DAG), we can extract the parameters to (a) appreciate the magnitude (and precision) of the the associations and (b) further refine the DAG including whether the associations are positive or negative.

10.1) Visually inspect the marginal posterior distributions of the parameters

par(mfrow=c(1,4), mar=c(2,2,1.5,1))
for(i in 1:length(marg.f$marginals)){

# get the marginal for current node, which is a matrix [x, f(x)]
  cur.node <- marg.f$marginals[i]
  nom1 <- names(marg.f$marginals)[i]

# pick the first value (for models wothout random effects)
  cur.node <- cur.node[[1]]
  for(j in 1:length(cur.node) ) {
    nom2<-names(cur.node)[j]
    cur.param <- cur.node[[j]]
    plot(cur.param,type="l",main=paste(nom1, ":", nom2), cex=0.7)
  }
}

10.2) Calculate the area under the density curve

# extract marginals
marg.dens <- marg.f$marginals[[1]]
for (i in 2:length(marg.f$marginals)) {
  marg.dens <- c(marg.dens, marg.f$marginals[[i]])
}

# calculate AUC --> should be ~1
auc <- rep(NA, length(marg.dens))
names(auc) <- names(marg.dens)
for(i in 1:length(marg.dens)) {
  tmp <- spline(marg.dens[[i]])
  auc[i] <- sum(diff(tmp$x[-length(tmp$x)])*tmp$y[-1])
}

cbind(auc)

                            auc
AR|(Intercept)        0.9999997
AR|age                0.9999994
pneumS|(Intercept)    0.9999990
female|(Intercept)    0.9999994
livdam|(Intercept)    0.9999995
livdam|eggs           0.9999984
eggs|(Intercept)      0.9999993
eggs|adg              0.9999993
wormCount|(Intercept) 0.9999993
wormCount|AR          0.9999996
wormCount|eggs        0.9999997
wormCount|age         0.9999996
wormCount|adg         0.9999996
age|(Intercept)       0.9999995
age|female            0.9999992
age|precision         0.9999996
adg|(Intercept)       0.9999998
adg|age               0.9999998
adg|precision         0.9999985

par(las=2, mar=c(8.1,4.1,4.1,2.1))
barplot(auc, ylab="Area under Density", ylim=c(0,1.2))

Both the shape of the parameter density distributions and their AUC looks satisfactory, so we can proceed getting the estimates.

10.3) Get the table of quantiles for the marginals

mat <- matrix(rep(NA, length(marg.dens)*3), ncol=3)
rownames(mat) <- names(marg.dens)
colnames(mat) <- c("2.5%", "50%", "97.5%")
ignore.me <- union(grep("\\(Int", names(marg.dens)), grep("prec", names(marg.dens))) # take away background k and precision
comment <- rep("", length(marg.dens))
for (i in 1:length(marg.dens)) {
  tmp <- marg.dens[[i]]
  tmp2 <- cumsum(tmp[,2])/sum(tmp[,2])
  mat[i, ] <-c(tmp[which(tmp2>0.025)[1]-1,1],## -1 is so use value on the left of the 2.5%
               tmp[which(tmp2>0.5)[1],1],
               tmp[which(tmp2>0.975)[1],1])
  vec <- mat[i,]

  if( !(i%in%ignore.me) &&  (vec[1]<0 && vec[3]>0)){comment[i]<-"not sig. at 5%"}

  ## truncate for printing
  mat[i,]<-as.numeric(formatC(mat[i,],digits=3,format="f"))
}

cbind(mat)

                        2.5%    50%  97.5%
AR|(Intercept)         0.891  1.155  1.435
AR|age                 0.509  0.798  1.110
pneumS|(Intercept)    -2.134 -1.814 -1.519
female|(Intercept)    -0.196  0.018  0.230
livdam|(Intercept)     0.976  1.265  1.567
livdam|eggs            0.663  1.547  2.656
eggs|(Intercept)      -1.508 -1.232 -0.971
eggs|adg               0.346  0.614  0.895
wormCount|(Intercept) -1.646 -1.400 -1.170
wormCount|AR           0.152  0.273  0.395
wormCount|eggs         3.280  3.504  3.741
wormCount|age         -0.824 -0.685 -0.548
wormCount|adg         -0.374 -0.240 -0.107
age|(Intercept)       -0.349 -0.200 -0.052
age|female             0.188  0.397  0.606
age|precision          0.893  1.042  1.207
adg|(Intercept)       -0.053  0.000  0.053
adg|age               -0.921 -0.868 -0.816
adg|precision          3.485  4.070  4.712

As said, the marginals represent estimates of the parameters at each node (i.e. the arcs in the DAG). Being the variables of different nature, they have different meaning. Marginals represent correlation coefficients for Gaussian nodes (when the default centering of variables is kept), log rate ratios for Poisson nodes, and log odds ratios for binomial nodes. Therefore, the second and the latter need to be exponentiated, to get the odds ratios or rate ratios respectively.