# Improving Causal Gaussian Bayesian Network Inference using Parallel Tempering

Pascal Fieth, Gilles Monneret, Andrea Rau, Florence Jaffrezic, Alexander Hartmann, Grégory Nuel

► To cite this version:

Pascal Fieth, Gilles Monneret, Andrea Rau, Florence Jaffrezic, Alexander Hartmann, et al. Improving Causal Gaussian Bayesian Network Inference using Parallel Tempering. RSG DREAM 2016, Nov 2016, Phoenix, United States. hal-01608740

HAL Id: hal-01608740

https://hal.science/hal-01608740v1

Submitted on 5 Jun 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY-SA 4.0 - Attribution - ShareAlike - International
License
## Causality in Transcriptomic Data

<figure><img src="image_2.png" /><figcaption>Figure 1: Example of Gene Regulation Network (Zhou et al., 2007).</figcaption></figure>

### Gene Regulation Network:

* Biological tool to describe interaction between genes. It's very popular and easily readable (see example in Fig. 1).

* Links in these networks are supposed to have causal interpretation.

* Can be used to asked complex questions, such as what effect would have an intervention.

### Causality:

* Many high correlated variables do not have any effect on the variable of interest.

* Causality is usually inferred from interventional data (KO, siRNA, CRISPR-Cas9).

* Causal effects have much higher value than simple correlations.

### Already existing popular algorithms:

* Ges, which is a score based greedy algorithm (Chickering, 2003).

* Pcalg, that use individual partial correlation test to find causality (Kalisch and Bühlmann, 2007; Bühlmann et al., 2014);

* Glasso, do not find causality, but create Gaussian (unoriented) Graphical Model, by estimate inverse covariance matrix to do so (Friedman et al., 2008).

## Causal Gaussian Bayesian Network

* We use a linear Gaussian model (recursive structural equation modeling): each variable is disturbed with a noise $\epsilon_i$, which is a Gaussian random variable with mean $m_i$ and variance $\sigma_i^2$. All variables/genes $X_i$ of interest are linked with equations of the form:

$$X_i = \sum_k \alpha_k X_k + \epsilon_i \quad (1)$$

* These equations and random variables can be represented as a graphical network: in that case, each node $V_i$ correspond to a variable $X_i$ and each edge to a direct effect $\alpha_{i \to j}$. This gives a *Directed Acyclic Graph*. For example:

<figure><img src="image_3.png" /></figure>

$$ 
\begin{align*}
X_1 &= \epsilon_1 \\
X_2 &= \alpha_{1 \to 2} \times X_1 + \epsilon_2 \\
X_3 &= \alpha_{1 \to 3} \times X_1 + \epsilon_3 \\
X_4 &= \alpha_{1 \to 4} \times X_1 + \alpha_{2 \to 4} \times X_2 + \alpha_{3 \to 4} \times X_3 + \epsilon_4
\end{align*}
$$

* We use interventional data: on these, we have one or several *interventional* nodes. The equations which define these nodes are erased and the value of the node is fixed to 0. For the graph, we cut-off all edges coming to the corresponding nodes. For example, if we do an intervention on variables $X_2$ and $X_3$:

<figure><img src="image_4.png" /></figure>

$$ 
\begin{align*}
X_1 &= \epsilon_1 \\
X_2 &= 0 \\
X_3 &= 0 \\
X_4 &= \alpha_{1 \to 4} \times X_1 + \alpha_{2 \to 4} \times X_2 + \alpha_{3 \to 4} \times X_3 + \epsilon_4
\end{align*}
$$

* Using these two kind of data, and for a particular ordering, we can write the *likelihood*. We use the following notation: $N_j$ is the number of replicate where we do not have an intervention on variable $X_j$, $\mathcal{J}_k$ is the set of interventional node for the k-th replicate, $W$ is the adjacency matrix where we find all $\alpha$ and $e_j$ is a vector from the standard basis. We get:

$$ 
\ell(\mathbf{m}, \mathbf{W}, \boldsymbol{\sigma}^2) = -\frac{\log(2\pi)}{2} \sum_j N_j - \sum_j N_j \log(\sigma_j) -\frac{1}{2} \sum_k \sum_{j \in \mathcal{J}_k} \frac{1}{\sigma_j^2} (x_j^k - \mathbf{x}^k \mathbf{w} \mathbf{e}_j^T - m_j)^2. \quad (2)
$$

We can then maximize it, and we get maximum likelihood estimator for a fixed ordering (Rau et al., 2013)

## Exploring the Space of Causal Ordering

* The previous likelihood is true only for a *fixed* ordering of nodes. We do not know this one, we use an algorithm to explore this space. Note that the cardinality if this space can as huge as $N!$ (for $N$ nodes).

* We need a probability law on ordering: we use a Mallows law (Mallows, 1957). His law $q$ is:

$$q(v) \propto \exp^{-\beta d(u, v)}, \quad (3)$$

where $v$ is an ordering, $u$ a mode ordering, $\beta$ a dispersion parameter and $d$ is a distance for ordering (here the Kendall's tau). This is nice law with many good property, but in particular it's a *symmetric* law.

* We use an empirical Bayesian framework: we use an improper uniform prior on ordering and we try to find an a posterior law for these:

$$ 
\begin{align*}
\mathbb{P}(v|\text{data}) &\propto \int_{\theta} \mathbb{P}(\text{data}|\theta, v) \mathbb{P}(\theta|v) \mathbb{P}(v) d\theta \\
&\propto \mathbb{P}(\text{data}|\hat{\theta}_v) \mathbb{P}(v).
\end{align*}
$$

With our approximation, we see that the posterior law of ordering is proportional to the likelihood for a fixed ordering.

* We use a Metropolis-Hastings algorithm:

**Algorithm 1: Hastings-Metropolis**

for $i \in [1, n]$ **do**

Generate $y_n \sim \text{rmallow}(o_{n-1}, \beta)$;

$o_n = y_n$ with probability $\min \left( 1, \left( \frac{\mathbb{P}(\text{data}|\hat{\theta}_{y_n})}{\mathbb{P}(\text{data}|\hat{\theta}_{o_{n-1}})} \right)^{\frac{1}{T}} \right)$;

Else $o_n = o_{n-1}$

It's practicable for few nodes (~ 10), but we are very fast stuck in the dimensionality problem (space grow exponentially, very many local maxima).

## Using Temperature

* The Metropolis Hastings algorithm (Hukushima and Nemoto, 1996) enables the exploration of simple spaces (Fig. 2).

* Setting the temperature $T$ to value greater than 1 allows the system to leave local maxima.

* For different temperature $T$ after an equilibration time different typical likelihood values are sampled Fig. 3).

<figure><img src="image_5.png" /><figcaption>Figure 2: A possible simple likelihood landscape (the actual, high-dimensional, landscape can not be visualised).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 3: Typical time series of simulations for different scaling parameters $T$.</figcaption></figure>

* For more complex likelihood landscapes, the chance of leaving a local maximum gets comparatively small (Fig. 5).

## Parallel Tempering

* Idea: Use *parallel tempering* (Marinari and Parisi, 1992).

**Algorithm 2: Parallel tempering**

for $T_i \in [1, T_{\max}]$ **doparallel**

repeat

Simulate a system with temperature $T_i$;

Do one step of Metropolis Hastings;

Set $\Delta \ell_i = (\ell_i - \ell_{i+1})$ and $\Delta \beta_i = \frac{1}{T_i} - \frac{1}{T_{i+1}}$;

Make random exchange for neighboring systems

with probability: $P_{\text{swap}}(T_i) = \max(1, \exp(\Delta \ell_i \Delta \beta_i))$;

until convergence;

<figure><img src="image_7.png" /><figcaption>Figure 4: Swap neighbouring systems.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5: A possible complex likelihood landscape (the actual, high-dimensional, landscape can not be visualised).</figcaption></figure>

* Systems can leave local maxima with increasing $T$ and enter (different) maxima with decreasing $T$.

* Better exploration of space with the possibility to get out from local maxima.

* Improved calculation time thanks to the parallel programming.

## Simulation

* Networks were simulated according to Eq. (1);

* $m = 0.5$, $\sigma = 0.1$, and effects $a_{i \to j}$ uniformly drawn from $(-1, -0.25) \cup (0.25, 1)$;

* Number of nodes $N$ and number of experiments $E$ with random experiment schemes, controlled over an intervention probability per node $p_i$.

<figure><img src="image_9.png" /><figcaption>Table 1: Example of simulated adjacency matrix.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Table 2: Example of generated data, with different experiments.</figcaption></figure>

* Experiments data are obtained by removing the relevant coefficient in the adjacency matrix for observational data.

## Results

* Figure 6 shows a time series for the simulation of a realisation of a DAG with $N = 40$ nodes and $E = 80$ experiments and on average one knockout per experiment.

* Practically not typical but makes comparison possible and the maximum possible can be estimated as the one of the original ordering.

* Plot shows CPU time as if simulated on a single CPU → much faster on multi-core computers.

<figure><img src="image_11.png" /><figcaption>Figure 6: Time series for one realisation of DAG and experiment for the Metropolis Hastings algorithm and parallel tempering. The gray line shows the likelihood of the input ordering.</figcaption></figure>

* For complex landscapes the Metropolis Hastings algorithm often fails at reaching optimal likelihood values while parallel tempering does not.

* Effect increases with growing system size.

## Conclusion

* In case of complex likelihood landscapes parallel tempering can reach orderings with higher likelihood values more reliably.

* This becomes more evident with increasing system sizes.

* Due to the parallelization of system simulation, parallel tempering is faster in terms of Wallclocktime (for a sufficient amount of CPUs).

* While parallel tempering enables the exploration of ordering spaces for larger system sizes, a limitation is still given by the computation time of the likelihood, growing with $N^6$.

## References

* Peter Bühlmann, Markus Kalisch, and Lukas Meier. High-dimensional statistics with a view toward applications in biology. *Annual Review of Statistics and Its Application*, 1:255–278, 2014.

* David Maxwell Chickering. Optimal structure identification with greedy search. *The Journal of Machine Learning Research*, 3:507–554, 2003.

* Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Sparse inverse covariance estimation with the graphical lasso. *Biostatistics*, 9(3):432–441, 2008.

* Koji Hukushima and Koji Nemoto. Exchange monte carlo method and application to spin glass simulations. *Journal of the Physical Society of Japan*, 65(6):1604–1608, 1996.

* Markus Kalisch and Peter Bühlmann. Estimating high-dimensional directed acyclic graphs with the pc-algorithm. *The Journal of Machine Learning Research*, 8:613–636, 2007.

* Colin L Mallows. Non-null ranking models. i. Biometrika, pages 114–130, 1957.

* Enzo Marinari and Giorgio Parisi. Simulated tempering: a new monte carlo scheme. *EPL (Europhysics Letters)*, 19(6):451, 1992.

* Andrea Rau, Florence Jaffrézic, and Grégory Nuel. Joint estimation of causal effects from observational and intervention gene expression data. *BMC Systems Biology*, 7(1):111, 2013.

* Qing Zhou, Hiram Chipperfield, Douglas A Melton, and Wing Hung Wong. A gene regulatory network in mouse embryonic stem cells. *Proceedings of the National Academy of Sciences*, 104(42):16438–16443, 2007.

<footnote n="1">¹ Carl von Ossietzky Universitat Oldenburg, Oldenburg, Germany; ² LPMA, CNRS 7599, UPMC, Paris, France; ³ GABI, INRA, Jouy-en-Josas, Paris, France.</footnote>