I have matrix of differentially expressed genes and clinical parameters that are discrete values (like Outcome, Severity, timepoint). I would like to make a Bayesian Network to understand which Biological parameters are best explained by which set of genes/ connecting with which genes.

For this I am using bnlearn package and learning structure using hill-climbing (HC) method,

```
dag<-hc(combine)
```

Than I am learning Parameters for BN using bn.fit() with method = "mle-g", the maximum likelihood estimator for least squares regression models for hybrid networks.

```
fitted = bn.fit(dag, data= combine, method = "mle-g")
```

However, I feel like the BN is not getting properly learnt because most of the nodes are connected with each other. Rather than giving me specific modules.

can anyone guide me to what I can do to refine or improve on the methods that I am using to build the network and if there are any alternative that I can explore.

I am new to machine learning methods and would welcome any advice for further reading.

Thanks