Bayesian Inference for Latent Chain Graphs

Deng Lu; Maria De Iorio; Ajay Jasra; Gary L. Rosner

arXiv:1908.04002·stat.ME·August 13, 2019

Bayesian Inference for Latent Chain Graphs

Deng Lu, Maria De Iorio, Ajay Jasra, Gary L. Rosner

PDF

Open Access

TL;DR

This paper develops a Bayesian inference framework using sequential Monte Carlo methods for AMP Gaussian chain graph models, which are useful in biological and financial applications but pose computational challenges.

Contribution

It introduces a novel SMC-based approach to handle prior modeling and inference in complex AMP Gaussian chain graph models.

Findings

01

Effective inference demonstrated on simulated data

02

Successful application to real-world case studies

03

Addresses computational challenges in chain graph models

Abstract

In this article we consider Bayesian inference for partially observed Andersson-Madigan-Perlman (AMP) Gaussian chain graph (CG) models. Such models are of particular interest in applications such as biological networks and financial time series. The model itself features a variety of constraints which make both prior modeling and computational inference challenging. We develop a framework for the aforementioned challenges, using a sequential Monte Carlo (SMC) method for statistical inference. Our approach is illustrated on both simulated data as well as real case studies from university graduation rates and a pharmacokinetics study.

Tables6

Table 1. Table 1: Posterior probability ℙ ( a i j = 0 | y 1 : m , α ) , 1 ≤ i < j ≤ p ℙ subscript 𝑎 𝑖 𝑗 conditional 0 subscript 𝑦 : 1 𝑚 𝛼 1 𝑖 𝑗 𝑝 \mathbb{P}(a_{ij}=0|y_{1:m},\alpha),1\leq i<j\leq p .

Nodes	$2$	$3$	$4$	$5$	$6$	$7$	$8$	$9$	$10$
$1$	0.882	0.892	0.920	0.932	0.920	0.870	0.874	0.908	0.934
$2$		0.902	0.906	0.828	0.898	0.934	0.784	0.804	0.906
$3$			0.914	0.890	0.880	0.890	0.908	0.900	0.890
$4$				0.900	0.866	0.882	0.770	0.918	0.900
$5$					0.788	0.952	0.912	0.830	0.918
$6$						0.806	0.932	0.906	0.908
$7$							0.904	0.738	0.916
$8$								0.918	0.740
$9$									0.952

$s p e n d$	average spending per student,
$s t r a t$	student-teacher ratio,
$s a l a r$	faculty salary,
$r e j r$	rejection rate,
$p a c c$	percentage of admitted students who accept university’s offer,
$t s t s c$	average test scores of incoming students,
$t o p 10$	class standing of incoming freshmen, and
$a p g r a$	average percentage of graduation.

Table 3. Table 4: Summaries of different chain graphs using package SEM.

Base chain graph				Chain graph selected by SIN
Edge	p-value	Edge	p-value	Edge	p-value
strat — spend	1.630e-14	pacc $\to$ salar	1.137e-06	strat — spend	1.630e-14
strat $\to$ salar	1.382e-06	pacc $\to$ rejr	1.470e-03	strat — salar	1.082e-05
spend — salar	2.629e-11	strat $\to$ apgra	8.237e-02	strat $\to$ top10	1.935e-09
strat $\to$ top10	4.743e-07	spend $\to$ apgra	6.067e-02	spend — salar	7.156e-13
spend $\to$ top10	2.822e-28	salar $\to$ apgra	4.794e-03	spend $\to$ top10	5.979e-34
top10 — salar	8.931e-03	top10 $\to$ apgra	4.253e-01	spend $\to$ tstsc	3.995e-12
strat $\to$ tstsc	3.140e-03	tstsc $\to$ apgra	1.096e-10	spend $\to$ rejr	2.909e-03
spend $\to$ tstsc	6.634e-01	pacc $\to$ apgra	2.711e-03	salar $\to$ tstsc	2.350e-05
salar $\to$ tstsc	2.008e-04			salar $\to$ rejr	1.323e-03
top10 $\to$ tstsc	6.831e-19			salar $\to$ pacc	1.827e-14
strat $\to$ rejr	1.954e-01			salar $\to$ apgra	1.570e-02
spend $\to$ rejr	3.621e-03			top10 — tstsc	1.256e-09
salar $\to$ rejr	2.575e-04			top10 — pacc	5.020e-01
top10 — rejr	1.816e-04			tstsc — rejr	8.297e-03
tstsc — rejr	1.003e-02			tstsc $\to$ apgra	8.352e-19
strat $\to$ pacc	2.585e-02			rejr — pacc	5.617e-03
spend $\to$ pacc	4.109e-07			pacc $\to$ apgra	5.481e-03
AIC		BIC		AIC	BIC
67.887		-13.319		80.838	-24.919
Chain graph selected by algorithm		Chain graph selected by algorithm		Chain graph selected by algorithm

( $α = (0.39, 0.25, 0.36, 0.05)$ )		( $α = (1, 1, 1, 1)$ )		( $α = (1, 3, 3, 3)$ )
Edge	p-value	Edge	p-value	Edge	p-value
strat — spend	1.630e-14	spend $\to$ strat	2.150e-52	spend $\to$ strat	1.536e-42
strat $\to$ salar	1.727e-06	strat $\to$ salar	9.127e-07	salar $\to$ strat	4.952e-06
strat $\to$ top10	3.597e-09	spend $\to$ salar	2.484e-24	strat $\to$ top10	1.636e-07
spend $\to$ salar	4.956e-33	top10 $\to$ spend	2.068e-25	tstsc $\to$ strat	8.237e-01
spend $\to$ top10	1.086e-33	salar $\to$ pacc	7.304e-10	salar $\to$ spend	9.659e-11
spend $\to$ tstsc	4.059e-12	apgra $\to$ salar	3.751e-11	spend $\to$ top10	3.881e-10
salar $\to$ tstsc	5.524e-05	top10 — tstsc	2.163e-14	tstsc $\to$ spend	2.886e-08
salar $\to$ pacc	1.012e-18	top10 $\to$ rejr	2.504e-03	tstsc $\to$ salar	3.654e-27
salar $\to$ apgra	3.201e-05	tstsc $\to$ rejr	2.847e-04	salar $\to$ rejr	8.844e-02
top10 — tstsc	2.295e-10	tstsc $\to$ apgra	1.428e-43	salar $\to$ pacc	1.062e-08
top10 $\to$ rejr	1.951e-03	rejr $\to$ pacc	1.480e-04	salar $\to$ apgra	1.036e-04
tstsc $\to$ rejr	2.348e-04	apgra $\to$ pacc	3.587e-03	tstsc $\to$ top10	1.033e-20
tstsc $\to$ apgra	1.367e-18			top10 $\to$ rejr	9.160e-03
rejr $\to$ pacc	1.032e-03			tstsc $\to$ rejr	5.443e-03
pacc — apgra	5.315e-03			tstsc $\to$ apgra	2.644e-17
				rejr $\to$ pacc	3.048e-04
				apgra $\to$ pacc	5.759e-03
AIC	BIC	AIC	BIC	AIC	BIC
61.372	-50.522	102.93	-18.169	58.401	-47.356

Table 4. Table 5: Tissues and cell types examined in the PK studies

Compound	Compartment	Notation
TFV	Blood plasma	TFV_plasma
TFV	Rectal biopsy tissue	TFV_tissue
TFV	Rectal fluid	TFV_rectal
TFVdp	Rectal biopsy tissue	TFVdp_tissue
TFVdp	Total mononuclear cells in rectal tissue	Total $_{MMC}$
TFVdp	CD4⁺ lymphocytes from MMC	CD4 ${}_{MMC}^{+}$
TFVdp	CD4^- lymphocytes from MMC	CD4 ${}_{MMC}^{-}$

Table 5. Table 6: Summaries of two different chain graphs using package SEM.

Chain graph selected by algorithm		Chain graph selected by algorithm

( $α = (1, 1, 1, 1)$ )		( $α = (1, 3, 3, 3)$ )
Edge	p-value	Edge	p-value
CD4 ${}_{MMC}^{-}$ $\to$ Total $_{MMC}$	9.881e-19	CD4 ${}_{MMC}^{+}$ $\to$ TFVdp_tissue	5.687e-01
CD4 ${}_{MMC}^{-}$ $\to$ CD4 ${}_{MMC}^{+}$	4.091e-81	CD4 ${}_{MMC}^{+}$ $\to$ CD4 ${}_{MMC}^{-}$	2.941e-46
CD4 ${}_{MMC}^{-}$ $\to$ TFV_rectal	1.496e-01	CD4 ${}_{MMC}^{+}$ $\to$ TFV_plasma	1.210e-02
Total $_{MMC}$ — TFVdp_tissue	1.589e-02	CD4 ${}_{MMC}^{+}$ $\to$ Total $_{MMC}$	7.028e-10
TFVdp_tissue — CD4 ${}_{MMC}^{+}$	1.812e-02	CD4 ${}_{MMC}^{+}$ $\to$ TFV_rectal	1.874e-03
CD4 ${}_{MMC}^{+}$ — TFV_rectal	8.583e-03	TFVdp_tissue — CD4 ${}_{MMC}^{-}$	8.477e-02
Total $_{MMC}$ $\to$ TFV_plasma	1.815e-03	TFVdp_tissue $\to$ TFV_rectal	2.991e-01
CD4 ${}_{MMC}^{+}$ $\to$ TFV_tissue	7.043e-01	TFVdp_tissue $\to$ TFV_plasma	2.584e-01
		TFVdp_tissue $\to$ Total $_{MMC}$	1.162e-13
		CD4 ${}_{MMC}^{-}$ $\to$ Total $_{MMC}$	2.352e-22
		TFV_plasma $\to$ TFV_tissue	4.259e-01
		TFV_plasma — Total $_{MMC}$	5.719e-02
		TFV_tissue $\to$ TFV_rectal	7.550e-01
		Total $_{MMC}$ $\to$ TFV_rectal	2.337e-02
AIC	BIC	AIC	BIC
Inf	Inf	46.875	-11.909

Table 6. Table 7: Summaries of modification indices for the model corresponding to the chain graph obtained under prior with α = ( 1 , 3 , 3 , 3 ) 𝛼 1 3 3 3 \alpha=(1,3,3,3) .

5 largest modification indices, $A$ matrix		5 largest modification indices, $P$ matrix

(regression coefficients)		(variances/covariances)
TFV_rectal $\to$ Total $_{MMC}$	3.281	TFV_rectal — Total $_{MMC}$	3.394
TFV_rectal $\to$ TFV_plasma	2.219	TFV_rectal — TFV_plasma	2.881
CD4 ${}_{MMC}^{-}$ $\to$ TFV_rectal	0.709	TFV_rectal — CD4 ${}_{MMC}^{-}$	0.709
TFV_rectal $\to$ TFVdp_tissue	0.654	TFV_rectal — TFVdp_tissue	0.709
TFV_rectal $\to$ CD4 ${}_{MMC}^{-}$	0.555	TFV_tissue — TFVdp_tissue	0.389

Equations35

a_{ij}=\left\{\begin{array}[]{ll}0&\textrm{ no edge between}\leavevmode\nobreak\ i\leavevmode\nobreak\ \textrm{and}\leavevmode\nobreak\ j\\ 1&\textrm{ undirected edge between}\leavevmode\nobreak\ i\leavevmode\nobreak\ \textrm{and}\leavevmode\nobreak\ j\\ 2&\textrm{ directed edge from}\leavevmode\nobreak\ i\leavevmode\nobreak\ \textrm{to}\leavevmode\nobreak\ j\\ 3&\textrm{ directed edge from}\leavevmode\nobreak\ j\leavevmode\nobreak\ \textrm{to}\leavevmode\nobreak\ i\end{array}\right.

a_{ij}=\left\{\begin{array}[]{ll}0&\textrm{ no edge between}\leavevmode\nobreak\ i\leavevmode\nobreak\ \textrm{and}\leavevmode\nobreak\ j\\ 1&\textrm{ undirected edge between}\leavevmode\nobreak\ i\leavevmode\nobreak\ \textrm{and}\leavevmode\nobreak\ j\\ 2&\textrm{ directed edge from}\leavevmode\nobreak\ i\leavevmode\nobreak\ \textrm{to}\leavevmode\nobreak\ j\\ 3&\textrm{ directed edge from}\leavevmode\nobreak\ j\leavevmode\nobreak\ \textrm{to}\leavevmode\nobreak\ i\end{array}\right.

a_{ji}=\left\{\begin{array}[]{ll}a_{ij}&\textrm{if}\leavevmode\nobreak\ a_{ij}\in\{0,1\}\\ 2&\textrm{if}\leavevmode\nobreak\ a_{ij}=3\\ 3&\textrm{if}\leavevmode\nobreak\ a_{ij}=2.\end{array}\right.

a_{ji}=\left\{\begin{array}[]{ll}a_{ij}&\textrm{if}\leavevmode\nobreak\ a_{ij}\in\{0,1\}\\ 2&\textrm{if}\leavevmode\nobreak\ a_{ij}=3\\ 3&\textrm{if}\leavevmode\nobreak\ a_{ij}=2.\end{array}\right.

Σ (B, Ω) = (I - B)^{- 1} Ω^{- 1} (I - B^{'})^{- 1}

Σ (B, Ω) = (I - B)^{- 1} Ω^{- 1} (I - B^{'})^{- 1}

p (y ∣ B, Ω, A) = τ \in T (G (A)) \prod p (y_{τ} ∣ y_{pa (τ)}, B, Ω) .

p (y ∣ B, Ω, A) = τ \in T (G (A)) \prod p (y_{τ} ∣ y_{pa (τ)}, B, Ω) .

p (y_{1 : m} ∣ B, Ω, A) = i = 1 \prod m p (y_{i} ∣ B, Ω, A)

p (y_{1 : m} ∣ B, Ω, A) = i = 1 \prod m p (y_{i} ∣ B, Ω, A)

P (a_{ij} = l ∣ π_{ij}) = π_{ij} (l), l = 0, \dots, r

P (a_{ij} = l ∣ π_{ij}) = π_{ij} (l), l = 0, \dots, r

P (a_{ij} = l ∣ α)

P (a_{ij} = l ∣ α)

= \frac{Γ ( 1 + α _{a_{ij}} )}{Γ ( α _{a_{ij}} )} \frac{Γ ( \sum _{k = 0}^{r} α _{k} )}{Γ ( 1 + \sum _{k = 0}^{r} α _{k} )}, l = 0, \dots, r

B (α) = \frac{\prod _{k = 0}^{r} Γ ( α _{k} )}{Γ ( \sum _{k = 0}^{r} α _{k} )}

B (α) = \frac{\prod _{k = 0}^{r} Γ ( α _{k} )}{Γ ( \sum _{k = 0}^{r} α _{k} )}

p\big{(}(a_{ij})_{i<j}\mid\alpha\big{)}\propto\mathbb{I}_{\mathcal{C}}(G(A))\Big{\{}\prod_{i<j}\mathbb{P}(a_{ij}|\alpha)\Big{\}}.

p\big{(}(a_{ij})_{i<j}\mid\alpha\big{)}\propto\mathbb{I}_{\mathcal{C}}(G(A))\Big{\{}\prod_{i<j}\mathbb{P}(a_{ij}|\alpha)\Big{\}}.

p (Ω, B ∣ A) = I_{P} (Σ (B, Ω)) p (Ω∣ A) p (B ∣ A)

p (Ω, B ∣ A) = I_{P} (Σ (B, Ω)) p (Ω∣ A) p (B ∣ A)

\pi(B,\Omega,(a_{ij})_{i<j}\mid y_{1:m},\alpha)\propto p(y_{1:m}|B,\Omega,A)p(\Omega,B|A)p\big{(}(a_{ij})_{i<j}|\alpha\big{)}.

\pi(B,\Omega,(a_{ij})_{i<j}\mid y_{1:m},\alpha)\propto p(y_{1:m}|B,\Omega,A)p(\Omega,B|A)p\big{(}(a_{ij})_{i<j}|\alpha\big{)}.

π_{t} (B, Ω, (a_{ij})_{i < j} ∣ y_{1 : m}, α) \propto ν_{t} (B, Ω, (a_{ij})_{i < j} ∣ α)

π_{t} (B, Ω, (a_{ij})_{i < j} ∣ y_{1 : m}, α) \propto ν_{t} (B, Ω, (a_{ij})_{i < j} ∣ α)

\nu_{t}(B,\Omega,(a_{ij})_{i<j}\mid\alpha)=\Big{[}p(y_{1:m}\mid B,\Omega,A)\Big{]}^{\phi_{t}}p(\Omega,B\mid A)p\big{(}(a_{ij})_{i<j}\mid\alpha\big{)}

\nu_{t}(B,\Omega,(a_{ij})_{i<j}\mid\alpha)=\Big{[}p(y_{1:m}\mid B,\Omega,A)\Big{]}^{\phi_{t}}p(\Omega,B\mid A)p\big{(}(a_{ij})_{i<j}\mid\alpha\big{)}

w_{t}^{(i)}=\Big{(}\frac{\nu_{t}(u_{t-1}^{(i)})}{\nu_{t-1}(u_{t-1}^{(i)})}\Big{)}\Big{(}\sum_{j=1}^{N}\frac{\nu_{t}(u_{t-1}^{(j)})}{\nu_{t-1}(u_{t-1}^{(j)})}\Big{)}^{-1}

w_{t}^{(i)}=\Big{(}\frac{\nu_{t}(u_{t-1}^{(i)})}{\nu_{t-1}(u_{t-1}^{(i)})}\Big{)}\Big{(}\sum_{j=1}^{N}\frac{\nu_{t}(u_{t-1}^{(j)})}{\nu_{t-1}(u_{t-1}^{(j)})}\Big{)}^{-1}

P (a_{ij} = 0 ∣ y_{1 : m}, α) = p = 1 \sum N W_{T}^{(p)} I_{{0}} (a_{T, ij}^{(p)}),

P (a_{ij} = 0 ∣ y_{1 : m}, α) = p = 1 \sum N W_{T}^{(p)} I_{{0}} (a_{T, ij}^{(p)}),

P (a_{ij} = k ∣ y_{1 : m}, α) = p = 1 \sum N W_{T}^{(p)} I_{{k}} (a_{T, ij}^{(p)}), k = 0, 1, 2, 3.

P (a_{ij} = k ∣ y_{1 : m}, α) = p = 1 \sum N W_{T}^{(p)} I_{{k}} (a_{T, ij}^{(p)}), k = 0, 1, 2, 3.

a_{ij} = ar g max_{k = 0, 1, 2, 3} P (a_{ij} = k ∣ y_{1 : m}, α) .

a_{ij} = ar g max_{k = 0, 1, 2, 3} P (a_{ij} = k ∣ y_{1 : m}, α) .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Taxonomy

TopicsBayesian Modeling and Causal Inference · Bayesian Methods and Mixture Models · Statistical Methods and Bayesian Inference

Full text

Abstract

In this article we consider Bayesian inference for partially observed Andersson-Madigan-Perlman (AMP) Gaussian chain graph (CG) models. Such models are of particular interest in applications such as biological networks and financial time series. The model itself features a variety of constraints which make both prior modeling and computational inference challenging. We develop a framework for the aforementioned challenges, using a sequential Monte Carlo (SMC) method for statistical inference. Our approach is illustrated on both simulated data as well as real case studies from university graduation rates and a pharmacokinetics study.

Key words: Chain Graphs, Bayesian Inference, Sequential Monte Carlo.

Bayesian Inference for Latent Chain Graphs

BY DENG LU1, MARIA DE IORIO2, AJAY JASRA3 & GARY L. ROSNER4

1Department of Statistics & Applied Probability, National University of Singapore, Singapore, 117546, SG. E-Mail: [email protected]

2Yale-NUS, Singapore, 138527, SG. & Department of Statistical Science, University College London, UK E-Mail: [email protected]

3Computer, Electrical and Mathematical Sciences and Engineering Division, King Abdullah University of Science and Technology, Thuwal, 23955, KSA. E-Mail: [email protected]

4Oncology Biostatistics and Bioinformatics, Sidney Kimmel Comprehensive Cancer Center, Johns Hopkins University, Baltimore, MD 21205, USA. E-Mail: [email protected]

1 Introduction

Two common approaches to probabilistically describe conditional dependence structures are based on undirected graphs (Markov networks) and directed acyclic graphs (DAG) (Bayesian networks), e.g. [25]. The vertices of the graph represent variables, while the presence (absence) of an edge between two vertices indicates possible dependence (independence) between the two corresponding variables. Applications of Markov networks include biological networks, e.g., to study the dependence structure among genes from expression data [9, 11] and financial time series for forecasting and predictive portfolio analysis [8, 33]; while Bayesian networks occur in expert system research for providing rapid absorption and propagation of evidence [18, 24], path analysis for computing implied correlations [34], and psychometrics for causal pathways [6, 16].

Chain graphs [17, 19, 35] provide an elegant unifying point of view on Markov and Bayesian networks. These models allow for edges that are both directed and undirected and do not contain any semi-directed cycles. Although chain graphs were first introduced in the late eighties, most research has focused on Bayesian networks and Gaussian Graphical Models. Recently they are receiving more attention as they have proved to be very useful in applications due to their ability to represent symmetric and non-symmetric relationships between the random variables of interest. However, there exist in the literature several different ways of interpreting chain graphs and what conditional independencies they encode, giving rise to different so-called chain graph interpretations. Whilst a directed edge works as in DAGs when it comes to representing independence, the undirected edge can be understood in different ways, giving rise to the different interpretations. This implies that chain graphs can represent every independence model achievable by any DAG, whereas the opposite does not hold. The most common interpretations are the Lauritzen-Wermuth-Frydenberg (LWF) interpretation, the Andersson-Madigan-Perlman interpretation and the multivariate regression (MVR) interpretation. Each interpretation has its own way of determining conditional independences in a CG and each interpretation subsumes another in terms of representable independence models; see [30]. Moreover, [20] discuss causal interpretation of chain graphs.

This work focuses on chain graphs of the AMP type [2, 22] and develops a statistical framework to infer the chain graph structure from data measured on a set of variables of interest with the goal of understanding their association pattern. Many inferential procedures have concentrated on the estimation of the parameters characterizing the graph, given the chain graph topology. For example, [13] consider likelihood-based parameter inference for chain graph models that are directly observed. Nevertheless, in the machine learning literature, methods based on optimization have been proposed to perform inference on the graph structure itself in the context of chain graphs, e.g. [23, 26]. We consider the case where not only are the parameters unknown, but the chain graph itself is unobserved. The main objective is to perform Bayesian inference in such a context.

There are several contributions that propose Bayesian methods in related contexts. In [28], the author proposes a method to perform full Bayesian inference for acyclic directed mixed graphs using Markov chain Monte Carlo (MCMC) methods. There are several similarities between the model there (based upon structural equation models (SEM) e.g . [7]) and the one developed here. However, the main difference is given by the structure of the latent graph and the associated constraints in our case. These latter constraints are more stringent and lead to more complexity in the computational strategy. In [32] a similar model to the one presented in [28] is developed, which uses a spike-and-slab prior to induce sparsity in the graph. The computational scheme is expected to outperform that in [28] for a large number of nodes. The prior structure in [32] could in principle be extended to accommodate chain graph structures at the cost of introducing constraints that could make computations infeasible.

The contributions of this article are as follows. Based upon the likelihood method in [13] we develop a new Bayesian model for latent AMP chain graphs. This model can also incorporate covariate information, if available. We introduce a sequential Monte Carlo (SMC) method as in [10] that improves upon MCMC-based methods in this context. Our approach is applied to real case studies from university graduation rates and a pharmacokinetics study. We find the performance of our algorithm to be stable and robust.

This article is structured as follows. In Section 2 we describe our model and prior specifications. In Section 2.3 we introduce the SMC algorithm. In Section 3 we present a simulation study and two real-life applications. In the appendix we detail further elements of our algorithm in Section 2.3.

2 Model

2.1 Likelihood

Let $Y=(Y_{1},\ldots,Y_{p})\in\mathbb{R}^{p}$ be a random vector whose elements correspond to $p\in\mathbb{N}$ nodes of a graph. We assume we have $m\in\mathbb{N}$ observations on $Y$ : $y_{1:m}$ , $y_{i}\in\mathbb{R}^{p}$ .

A graph $G=(V,E)$ is described by a set of nodes $V$ and edges $E$ , with variables $Y_{1},\ldots,Y_{p}$ placed at the nodes. The edges define the global conditional independence structure of the distribution. An AMP chain graph is a graph whose every edge is directed or undirected such that it does not contain any semi-directed cycles, that is, it contains no path from a node $v$ to itself with at least one directed edge such that all directed edges have the same orientation. Each graph can be identified by an adjacency matrix. Let $A$ be a $p\times p$ matrix with entries $(a_{ij})_{1\leq i,j\leq p}$ where $a_{ij}\in\{0,\dots,r\}$ (let $r=3$ , which is a notation used later on) with $a_{ii}=0$ , for $i\neq j$

[TABLE]

and, given the upper-triangular part of $A$ , for $j<i$

[TABLE]

Given $p$ , a labelling of the nodes and the adjacency matrix $A$ define a graph $G(A)$ . Let $\mathcal{C}$ denote the set of possible chain graphs for a set of $p$ vertices. Given $A$ , let $\Omega$ be a positive definite $p\times p$ (real-valued) matrix, such that if $a_{ij}\neq 1$ then $\omega_{ij}=0$ . In addition, given $A$ , for an arbitrary real-valued $p\times p$ matrix $B$ , then, if $a_{ij}\neq 2$ $b_{ij}=0$ . Finally set

[TABLE]

where $I$ is the $p\times p$ identity matrix. To inherit the AMP chain graph property, $\Sigma(B,\Omega)$ should be a $p\times p$ positive definite matrix.

The absence of semi-directed cycles implies that the vertex set of a chain graph can be partitioned into so-called chain components such that edges within a chain component are undirected whereas the edges between two chain components are directed and point in the same direction. More precisely, the vertex set of a chain graph $G(A)$ can be partitioned into subsets $\mathscr{T}(G(A))$ such that all edges within each subset $\tau$ are un-directed and edges between two different subsets $\tau\neq\tau^{\prime}$ are directed. In the following, we assume that the partition $\tau\in\mathscr{T}(G(A))$ is maximal, that is, any two vertices in a subset $\tau$ are connected by an un-directed path. Recall that the parents of a set of nodes $\mathcal{X}$ of $G$ is the set $\textrm{pa}(\mathcal{X})=\{V_{j}\text{ s.t. }V_{j}\rightarrow V_{i}\in G,V_{j}\not\in\mathcal{X}\quad\&\quad V_{i}\in\mathcal{X}\}$ . Let $\textrm{pa}(\tau)$ be parents of $\tau$ , $B_{\tau}=(B_{ij})_{i\in\tau,j\in\textrm{pa}(\tau)}$ and $\Omega_{\tau}=(\omega_{ij})_{i,j\in\tau}$ sub-matrices of $B$ and $\Omega$ , respectively. For a single observation $y\in\mathbb{R}^{p}$ , $y_{\tau}\mid y_{\textrm{pa}(\tau)},B,\Omega\sim\mathcal{N}_{|\tau|}(B_{\tau}y_{\textrm{pa}(\tau)},\Omega_{\tau}^{-1})$ , where, for $d\in\mathbb{N}$ , $\mathcal{N}_{d}(\mu,\Sigma)$ denotes the $d-$ dimensional Gaussian distribution with mean $\mu$ and covariance matrix $\Sigma$ . Then the joint density of $y\in\mathbb{R}^{p}$ given $B,\Omega,A$ is

[TABLE]

The likelihood is then given by

[TABLE]

that is, the observations are i.i.d. given $B,\Omega,A$ . We remark that this structure is similar to the structural equation model proposed by [29], where also a Gaussian distribution is assumed for the variables corresponding to the nodes of the graph. Moreover, these above assumptions on the precision matrix ensure that the corresponding graph satisfies the AMP Markox property as discussed in [2].

2.2 Prior Distribution on the Chain Graph Space

We specify the prior on the set of possible chain graphs, by specifying a prior on the elements of the corresponding adjacency matrix. For $1\leq i<j\leq p$ , let

[TABLE]

where $\pi_{ij}=(\pi_{ij}(0),\dots,\pi_{ij}(r))$ . We assume each vector $\pi_{ij}$ to follow a Dirichlet distribution with parameter $\alpha=(\alpha_{0},\dots,\alpha_{r})$ and the $\pi_{ij}$ to be conditionally independent given $\alpha$ . Note that, if available, covariate information could be incorporated to model $\pi_{ij}$ , as in [31] for example. Marginally, we have

[TABLE]

where

[TABLE]

The prior for a graph $G(A)$ then becomes:

[TABLE]

Note that the prior does not integrate to 1, but has a finite mass.

We assume that given $G(A)$ , $\Omega$ has a $G-$ Wishart prior distribution ([21]) of parameters $\delta,D$ , and given $A$ the non-zero entries of $B$ are i.i.d. univariate Gaussian random variables with mean $\xi$ and variance $\kappa$ . Let $\mathcal{P}$ be the cone of $p\times p$ real-valued positive definite matrices. Then the joint prior for $\Omega$ and $B$ is given by:

[TABLE]

which does not integrate to 1.

2.3 Posterior Inference

The posterior distribution $\pi$ of all the unknown of interest is given by:

[TABLE]

The structure of the model bears some resemblance to the models considered in [28, 32]. Besides the graph structure being different (namely, we do not allow bi-directional edges), we also impose the AMP constraint, mathematically represented by the term $\mathbb{I}_{\mathcal{P}}(\Sigma(B,\Omega))$ in $p(\Omega,B|A)$ , which leads to a much more complex structure of the graph space than the aforementioned references. As a result, posterior exploration of graph space is more challenging and computationally demanding and requires carefully devised algorithms, in a sense more sophisticated, than the ones considered in [28, 32].

Our approach is to design an adaptive SMC sampler as in [15] (see [10] for the original algorithm and [5] for convergence results). SMC based algorithms offer the advantage of being easily parallelisable, often reducing computation times over serial methods in some scenarios. The strategy involves simulating $N$ samples (particles) to approximate the sequence of densities

[TABLE]

where

[TABLE]

and $0=\phi_{0}<\cdots<\phi_{T}=1$ . The motivation for this algorithm is well-documented (see the aforementioned references for details), and SMC algorithms have been successfully employed in may contexts to sample from high dimensional posterior distributions.

In the implementation of the algorithm, we require a Markov kernel (e.g., a MCMC kernel) $K_{t}$ , $t\in\{1,2,\dots,T\}$ that admits $\pi_{t}$ as an invariant distribution; this step is detailed in the appendix. The algorithm is summarised in Algorithm 1. For notational convenience, we set $u_{t}=(B_{t},\Omega_{t},(a_{t,ij})_{i<j})$ . In the initialization stage, one simulates exactly from $\pi_{0}$ using rejection sampling. The sequence of $\{\phi_{t}\}$ is set as proposed in [37], and the choice of parameters for the MCMC kernel follows [15].

3 Simulation and Real Data Study

In the following numerical experiments, we set $\delta=3$ and $D=I_{p}$ for the G-Wishart prior, and $\xi=0,\kappa=1$ for the distribution of the non-zero element of $B$ . The number of particles for the SMC is $N=500$ . The MCMC steps are Metropolis-Hastings kernels.

3.1 Simulated Example

In the simulated example, we assume that the $p$ random variables corresponding to the nodes of the graph are independent, i.e., $a_{ij}=0$ for all $1\leq i,j\leq p$ . This is a benchmark example to evaluate the ability of our algorithm to recover the structure as this is a very simple graph structure. We consider $p=10$ vertices and $m=100$ observations. To generate the data, we set the precision matrix $\Omega$ equal to identity matrix $I_{p}$ , and the entries of matrix $B$ are set to be [math]. For the Dirichlet prior, we consider $\alpha=(3,1,1,1)$ . This prior implies a higher probability of no connection and assumes the same prior probability for any type of edge between two nodes. The reason we prefer this choice of hyper-parameters to $\alpha=(1,1,1,1)$ (which corresponds to a uniform prior) is due to the computational time required by the initialization step to generate the chain graph.

The simulation results are summarized in Figure 1. The effective sample size (ESS) drops very fast as the algorithm begins and goes into a steady lower state after several resampling procedures. The acceptance rate is acceptable as it does not fall below a very small value. Using the weighted samples $\{W_{T}^{(n)},A_{T}^{(n)}\}_{n=1}^{N}$ obtained at the last step of SMC, we estimate the posterior probability $\mathbb{P}(a_{ij}=0\mid y_{1:m},\alpha),1\leq i<j\leq p$ by

[TABLE]

and summarize the results in Table 1. From the table, we can see that all the estimated posterior probabilities are greater than 0.7, and most of them are above 0.9. This suggests that our algorithm is able to recover the structure of the graph (independence) used to generate the data.

3.2 University Graduation Rates

We investigate the performance of our algorithm when the latent graph has a more complex structure. We consider data first presented in [14] that stem from a study for college ranking carried out in 1993. After initial analysis, Druzdzel and Glymour focus on $p=8$ variables:

Based on $m=159$ universities, the correlation matrix of these eight variables is estimated in [14]. Conditional on this correlation matrix, [13] obtain a chain graph through the SIN model selection with significance level 0.15. To get a chain graph with different AMP and LWF Markov properties, [13] further deleted the undirected edge between $top10$ and $rejr$ and the undirected edge between $salar$ and $top10$ , and introduced an undirected edge between $top10$ and $rejr$ . The resulting graph is shown in Figure 2, and the corresponding adjacency matrix is shown in Table 2. This graph has three chain components $\tau_{1}=\{spend,strat,salar\},\tau_{2}=\{top10,tstsc,rejr,pacc\}$ and $\tau_{3}=\{apgra\}$ . In the original article [14] provide some insight about the causal relationship between some of the variables. The average spending per student ( $spend$ ), student-teacher ratio ( $strat$ ) and faculty salary ( $salar$ ) are determined based on budget considerations and are not influenced by any of the five remaining variables. Rejection rate ( $rejr$ ) and percentage of students who accept the university’s offer from among those who are offered admission ( $pacc$ ) precede the average test scores ( $tstsc$ ) and class standing ( $top10$ ) of incoming freshmen. The latter two are measures taken over matriculating students. Finally, graduation rate ( $apgra$ ) does not influence any of the other variables.

To test the performance of the proposed Bayesian chain graph model, we build an empirical baseline graph through the following procedure. Starting with a graph with only one node denoted by $G$ , each time we add a new node to the existing graph. For every node in the subgraph not containing the new node, we consider the candidate graphs that are generated by adding one of four types of correlation (no edge, undirected edge and directed edges in either direction) between the new node and the node in the subgraph. We choose the graph with smallest BIC value obtained by fitting a structural equation model (SEM). SEM is a multivariate statistical analysis technique that is commonly used to analyse structural relationships. In general, the formulation of the SEM is given by the basic equation $v=Av+u$ , where $v$ and $u$ are vectors of random variables. The parameter matrix $A$ contains regression coefficients and the matrix $P=\mathbb{E}(uu^{\prime})$ gives covariances among the elements of $u$ . The chain graph model described in this work can be written as $Y=BY+U$ , where $U$ is a multivariate Gaussian random vector with covariance matrix $\Omega^{-1}$ . This is consistent with the structural equation model formulation. Figure 3 shows the derived graph and Table 3 shows the corresponding adjacency matrix.

We compare posterior inference from our model with the chain graphs in Figure 2 and Figure 3. For the SMC sampler, we set the number of samples $N=5000$ . For the Dirichlet prior, we first consider a prior based on the analysis of [13], i.e., choosing $\alpha=(0.39,0.25,0.36,0.05)$ by matching the probabilities of each type of edges in Figure 2. More precisely, the number of no-edges, un-directed edges and directed edges from $i$ to $j$ are 11, 7 and 10 according to the graphs in Figure 2. So the corresponding proportions of these three types of edge are the first three components of $\alpha$ . The last component of $\alpha$ is chosen to be 0.05 to ensure the occurrence of directed edge from $j$ to $i$ and the probability of this type of edge is not large compared to the other three types. We also perform posterior inference with $\alpha=(1,1,1,1)$ , which is a uniform prior, and with $\alpha=(1,3,3,3)$ , which favours more connections. The remaining parameters are specified as in the previous section.

Based on the weighted samples $\{W_{T}^{(n)},A_{T}^{(n)}\}_{n=1}^{N}$ , we first estimate the posterior probability of occurrence of each edge $\widehat{\mathbb{P}}(a_{ij}=k\mid y_{1:m},\alpha),$ $1\leq i<j\leq p$ by

[TABLE]

Then the entries of the estimated adjacency matrix $A$ are given by

[TABLE]

Figure 4 shows the estimated chain graphs obtained under each prior. The structure of the chain graph obtained from setting $\alpha=(0.39,0.25,0.36,0.05)$ is obviously more similar to Figure 2 than that obtained using $\alpha=(1,1,1,1)$ . This is not surprising since the prior with $\alpha=(0.39,0.25,0.36,0.05)$ is very informative and forces the posterior distribution of the edges to be similar to Figure 2. The estimated chain graph obtained setting the hyper-parameters equal to $\alpha=(1,3,3,3)$ favours more connections compared with other two priors. If we compare our results with the conclusions in [14], the estimated chain graph obtained under the informative prior indeed shows that $spend$ , $strat$ and $salar$ are* parents* of other variables, and $apgra$ is a child of the others. However, it does not show that $rejr$ and $pacc$ precede $tstsc$ and $top10$ (the graph just shows the opposite relation). Similarly, the graph obtained under the prior with $\alpha=(1,3,3,3)$ shows that $tstsc$ is a parent of $salar$ , $strat$ and $spend$ , which is a contrast to the available prior information. This can be improved by modifying the ordering of nodes and choosing a small value of $\alpha_{3}$ .

Finally, we compare our results with those obtained using SEM. We use package sem in R to fit a SEM to these data. Note that when fitting the SEM model the graph topology is fixed. Table 4 presents a brief summary of three different chain graphs using the sem function. From the table, we can see that the chain graph selected by the Bayesian chain graph model with prior $\alpha=(0.39,0.25,0.36,0.05)$ has smaller AIC and BIC values compared with the graph selected by the SIN model selection procedure and the empirical graph, which suggests that our algorithm can take advantage of available prior information and perform better.

3.3 Tenofovir study

In this section, we illustrate our Bayesian model for AMP chain graphs on a real data application from the RMP-02/MTN-006 study ([3]). Tenofovir (TFV) is a medication used to treat chronic hepatitis B and to prevent and treat HIV. TFV $1\%$ gel demonstrated $39\%$ protective efficacy in women using the gel within 12 hours before and after sexual activity in the Centre for the AIDs Programme of Research in South Africa 004 study [1]. Daily dosing of tenofovir disoproxil fumarate (TDF)/emtricitabine provides 62 to $73\%$ protection against HIV transmission in serodiscordant men and women enrolled in the Partners PrEP study [4]. The RMP-02 study was designed to evaluate the systemic safety and biologic effects of oral TFV combined with a gel formulation of the drug for application rectally and vaginally. The study enrolled 18 patients, all of whom received a single oral dose of TFV, and randomized each patient to receive either the gel formulation or a placebo several weeks later. Details about the phase 1 study are given in [3]. [27] present analyses of the ancillary studies into the treatment’s biologic effects.

The biologic effects we examine here concern the pharmacokinetics (PK) of TFV and its active metabolite tenofovir diphosphate (TFVdp), as well as the pharmacodynamics of the drug. The PK studies evaluate subjects’ TFV and TFVdp concentrations in multiple physiologic compartments (i.e., tissues and cells) across multiple time points during the study. Table 5 lists the compartments.

[27] demonstrate that tissue HIV infectibility (cumulative p24) is correlated with in vivo concentrations of both TFV and TFVdp. Statistically significant, non-linear dose-response relationships with reduced tissue infectibility are found for one TFV compartment and four TFVdp compartments; the dose-response relationships are highly significant for TFVdp in whole rectal tissue, CD4 ${}^{+}_{\text{MMC}}$ , CD4 ${}^{-}_{\text{MMC}}$ and Total ${}_{\text{MMC}}$ compartments. Furthermore, [36] conduct a comprehensive pharmacokinetic study of rectally administered TFV gel that describes the distribution of TFV and TFVdp into various tissue compartments relevant to HIV infection. They argue that TFV rectal fluid concentrations may be reasonable bio-indicators of plasma and rectal tissue concentrations, making it easier to estimate adherence and TFV concentrations in the target tissue. Therefore, the correlations between the TFV and TFVdp in compartments can be helpful in studying HIV suppression procedure and providing a measure of drug efficacy, enabling more advanced population pharmacokinetic modelling methods. In the following analysis, we investigate the correlation structure of p = 7 concentration levels in Table 5 collected at visit 12. From clinical knowledge, we would expect the following associations:

•

TFVplasma is associated with TFVtissue (blood levels and tissue levels),

•

TFVtissue is associated with TFVrectal (rectal tissue and rectal fluid),

•

TFVtissue is associated with Total ${}_{\text{MMC}}$ (rectal tissue and mononuclear cells in rectal tissue),

•

Total ${}_{\text{MMC}}$ is associated with CD4 ${}^{+}_{\text{MMC}}$ and CD4 ${}^{-}_{\text{MMC}}$ (total and constituents).

We normalise the observations of each variable to have zero mean and standard deviation of one. The number of observation is $m=11$ , and the number of particles is set to $N=5000$ . We fit the model using as hyper-parameters in the Dirichlet prior both $\alpha=(1,1,1,1)$ , which corresponds to the uniform prior, and $\alpha=(1,3,3,3)$ , which favours the presence of connections. This latter prior choice is more suitable for small sample sizes. The remaining parameters are set as in the previous section.

Based on weighted samples $\{W_{T}^{(n)},(A,B,\Omega)_{T}^{(n)}\}_{n=1}^{N}$ , we calculate the posterior probabilities $\mathbb{P}\big{(}(A,B,\Omega)_{T}^{(n)}\mid y_{1:n},\alpha\big{)}$ . As posterior estimate of the resulting chain graph, we report the one obtained from the adjacency matrix $A_{T}^{(n^{*})}$ where $n^{*}=\mathop{\arg\max}_{n}\,\mathbb{P}\big{(}(A,B,\Omega)_{T}^{(n)}\mid y_{1:n},\alpha\big{)}$ , and is shown in Figure 5 (for both prior settings). Obviously, the chain graph obtained under the prior with $\alpha=(1,3,3,3)$ has more connections. Moreover, this graph also shows that the TFVtissue is related to TFVplasma and TFVrectal, the TFVtissue is associated with Total ${}_{\text{MMC}}$ through TFVplasma, and CD4 ${}^{+}_{\text{MMC}}$ causes Total ${}_{\text{MMC}}$ and CD4 ${}^{-}_{\text{MMC}}$ . This is consistent with the clinical knowledge mentioned before. However, some of these associations are missing in the chain graph obtained under the uniform prior, e.g. the edge between TFVplasma and TFVtissue (they are independent given Total ${}_{\text{MMC}}$ and CD4 ${}^{+}_{\text{MMC}}$ under the AMP chain graph property).

We compare our results to those obtained by fitting a SEM model using the sem function in the R package sem. We summarise the results in Table 7. The AIC and BIC values of the model corresponding to the chain graph obtained under the uniform prior are missing, which may due to the singular Hessian matrix when estimating the covariance matrix of the parameter. From the table, it is evident that the p-values of some edges are large. For example, the edge CD4 ${}^{+}_{\text{MMC}}$ $\rightarrow$ TFVtissue in the chain graph obtained under uniform prior has a p-value 0.704. These results suggest that these relationships are not significant in a SEM model.

When inadequate fit of a structural equation model is observed, model modification is often conducted followed by retesting of the modified model. The most popular statistic is the modification index, which is a chi-square score test statistic with degree of freedom one. The modification index provides an estimated value in which the model’s chi-square test statistic would decrease if the corresponding parameter is added to the model and respecified as a free parameter. We perform model modification using the modIndices function in sem package, which calculates modification indices and estimates parameter changes for the fixed and constrained parameters in a structural equation model. Table 5 shows the five largest modification indices for both the $A$ matrix and $P$ matrix for the model corresponding to the chain graph obtained under the Dirichlet prior with $\alpha=(1,3,3,3)$ . The modification indices suggest that a better fit to the data would be achieved by adding association between TFVrectal and TFVplasma to the model. The small sample size makes it obviously challenging to estimate associations. Furthermore, the level in rectal tissue may depend on whether or not the patient received the TFV gel or the placebo, which may make the correlation of TFVtissue with other variables unclear.

4 Conclusions

In this article we propose a novel Bayesian model for latent AMP chain graphs, for which observations are available only on the nodes of the graph. Posterior inference is performed through a specially devised SMC algorithm. We investigate the ability of the model to recover a range of structures, also when prior knowledge is available. The performance of the SMC sampler is stable and consistent in our numerical study. However, the sampler is not suitable when the number of nodes $p$ is large, as the initializing step is difficult. Moreover, the proposed algorithm does not scale well with respect to $p$ . The computational cost of computing the probabilities in the adjacency matrix is $\mathcal{O}(p^{2})$ , and the computation of the normalizing constant in the $G$ -Wishart distribution is quite expensive when $p$ is large (approximately $\mathcal{O}(p^{3})$ ).

Several extensions of this work are possible. First, the algorithm can be extended to large $p$ by choosing a more efficient initial proposal $q$ and by actually exploiting parallel computing techniques. Second, the model can be extended to accommodate multiple groups of observations, allowing borrowing information across groups or time periods. Third, we could avoid sampling $\Omega$ and $B$ by using a Laplace approximation when calculating the posterior probability.

Appendix A MCMC Kernel

The updates are taken in the order $A,\Omega,B$ together, then $\Omega$ and $B$ . We use Metropolis-Hastings steps. The updates for $\Omega$ and $B$ are preformed element-wise. For each element of $\Omega$ and $B$ a Gaussian random walk with variance $\sigma_{1}^{2}$ is used as proposal. The acceptance probability is calculated as usual and note that if any of the constraints on $\Omega$ or $B$ are violated, then the acceptance probability is zero.

The update for $A,\Omega,B$ together is much more complicated. We denote with $(A^{(n)},\Omega^{(n)},B^{(n)})$ the current values of $A,\Omega,B$ and the proposed states with $(A_{c}^{(n)},\Omega_{c}^{(n)},B_{c}^{(n)})$ . The move proceeds by picking an $(a_{ij}^{(n)})_{i<j}$ uniformly at random and then proposing one of the other possible 3 values with uniform probability. Note that at this stage, we will check if the proposed graph is a chain graph and if not, the move is rejected instantly. Based upon this proposal, we propose new values of $\Omega$ and $B$ as detailed below. Suppose the selected edge is $(i,j),i<j$ , consider the following scheme:

(i)

If the transformation of the selected edge is $``0\rightarrow 1"$ , we take $\Omega_{c}^{(n)}=\Omega^{(n)}$ except for elements $\Omega_{c}^{(n)}[i,j]$ and $\Omega_{c}^{(n)}[j,i]$ . For these elements, we set $\Omega_{c}^{(n)}[i,j]\sim\mathcal{N}(0,\sigma^{2}_{2})$ and $\Omega_{c}^{(n)}[j,i]=\Omega_{c}^{(n)}[i,j]$ . Finally, we set $B_{c}^{(n)}=B^{(n)}$ . 2. (ii)

If the transformation of the selected edge is $``0\rightarrow 2"$ , we take $B_{c}^{(n)}=B^{(n)}$ except for element $B_{c}^{(n)}[j,i]$ . For this element, we set $B_{c}^{(n)}[j,i]\sim\mathcal{N}(\xi,\kappa)$ . Finally, we set $\Omega_{c}^{(n)}=\Omega^{(n)}$ . 3. (iii)

If the transformation of the selected edge is $``0\rightarrow 3"$ , we take $B_{c}^{(n)}=B^{(n)}$ except for element $B_{c}^{(n)}[i,j]$ . For this element, we set $B_{c}^{(n)}[i,j]\sim\mathcal{N}(\xi,\kappa)$ . Finally, we set $\Omega_{c}^{(n)}=\Omega^{(n)}$ . 4. (iv)

If the transformation of the selected edge is $``1\rightarrow 0"$ , we take $\Omega_{c}^{(n)}=\Omega^{(n)}$ except for elements $\Omega_{c}^{(n)}[i,j]$ and $\Omega_{c}^{(n)}[j,i]$ . For these elements, we set $\Omega_{c}^{(n)}[i,j]=\Omega_{c}^{(n)}[j,i]=0$ . Finally, we set $B_{c}^{(n)}=B^{(n)}$ . 5. (v)

If the transformation of the selected edge is $``1\rightarrow 2"$ , we take $\Omega_{c}^{(n)}=\Omega^{(n)}$ except for elements $\Omega_{c}^{(n)}[i,j]$ and $\Omega_{c}^{(n)}[j,i]$ . For these elements, we set $\Omega_{c}^{(n)}[i,j]=\Omega_{c}^{(n)}[j,i]=0$ . Finally, we take $B_{c}^{(n)}=B^{(n)}$ except for element $B_{c}^{(n)}[j,i]$ . For this element, we set $B_{c}^{(n)}[j,i]\sim\mathcal{N}(\xi,\kappa)$ . 6. (vi)

If the transformation of the selected edge is $``1\rightarrow 3"$ , we take $\Omega_{c}^{(n)}=\Omega^{(n)}$ except for elements $\Omega_{c}^{(n)}[i,j]$ and $\Omega_{c}^{(n)}[j,i]$ . For these elements, we set $\Omega_{c}^{(n)}[i,j]=\Omega_{c}^{(n)}[j,i]=0$ . Finally, we take $B_{c}^{(n)}=B^{(n)}$ except for element $B_{c}^{(n)}[i,j]$ . For this element, we set $B_{c}^{(n)}[i,j]\sim\mathcal{N}(\xi,\kappa)$ . 7. (vii)

If the transformation of the selected edge is $``2\rightarrow 0"$ , we take $B_{c}^{(n)}=B^{(n)}$ except for element $B_{c}^{(n)}[j,i]$ . For this element, we set $B_{c}^{(n)}[j,i]=0$ . Finally, we set $\Omega_{c}^{(n)}=\Omega^{(n)}$ . 8. (viii)

If the transformation of the selected edge is $``2\rightarrow 1"$ , we take $B_{c}^{(n)}=B^{(n)}$ except for element $B_{c}^{(n)}[j,i]$ . For this element, we set $B_{c}^{(n)}[j,i]=0$ . Finally, we take $\Omega_{c}^{(n)}=\Omega^{(n)}$ except for elements $\Omega_{c}^{(n)}[i,j]$ and $\Omega_{c}^{(n)}[j,i]$ . For these elements, we set $\Omega_{c}^{(n)}[i,j]\sim\mathcal{N}(0,\sigma^{2}_{2})$ and $\Omega_{c}^{(n)}[j,i]=\Omega_{c}^{(n)}[i,j]$ . 9. (ix)

If the transformation of the selected edge is $``2\rightarrow 3"$ , we take $B_{c}^{(n)}=B^{(n)}$ except for elements $B_{c}^{(n)}[j,i]$ and $B_{c}^{(n)}[i,j]$ . For these elements, we set $B_{c}^{(n)}[j,i]=0$ , and $B_{c}^{(n)}[i,j]\sim\mathcal{N}(\xi,\kappa)$ . Finally, we set $\Omega_{c}^{(n)}=\Omega^{(n)}$ . 10. (x)

If the transformation of the selected edge is $``3\rightarrow 0"$ , we take $B_{c}^{(n)}=B^{(n)}$ except for element $B_{c}^{(n)}[i,j]$ . For this element, we set $B_{c}^{(n)}[i,j]=0$ . Finally, we set $\Omega_{c}^{(n)}=\Omega^{(n)}$ . 11. (xi)

If the transformation of the selected edge is $``3\rightarrow 1"$ , we take $B_{c}^{(n)}=B^{(n)}$ except for element $B_{c}^{(n)}[i,j]$ . For this element, we set $B_{c}^{(n)}[i,j]=0$ . Finally, we take $\Omega_{c}^{(n)}=\Omega^{(n)}$ except for elements $\Omega_{c}^{(n)}[i,j]$ and $\Omega_{c}^{(n)}[j,i]$ . For these elements, we set $\Omega_{c}^{(n)}[i,j]\sim\mathcal{N}(0,\sigma^{2}_{2})$ and $\Omega_{c}^{(n)}[j,i]=\Omega_{c}^{(n)}[i,j]$ . 12. (xii)

If the transformation of the selected edge is $``3\rightarrow 2"$ , we take $B_{c}^{(n)}=B^{(n)}$ except for elements $B_{c}^{(n)}[i,j]$ and $B_{c}^{(n)}[j,i]$ . For these elements, we set $B_{c}^{(n)}[i,j]=0$ , and $B_{c}^{(n)}[j,i]\sim\mathcal{N}(\xi,\kappa)$ . Finally, we set $\Omega_{c}^{(n)}=\Omega^{(n)}$ .

The acceptance probabiity for this move is easily calculated.

Bibliography37

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] Abdool , K. Q., Abdool , K. S. S., Frohlich , J. A., et al . (2010). Effectiveness and safety of tenofovir gel, an antiretroviral microbicide, for the prevention of HIV infection in women. Science , 329 , 1168-1174.
2[2] Andersson , S. A., Madigan , D. & Perlman , M. D. (2001). Alternative Markov properties for chain graphs. Scand. J. Statist. 28 , 33-85.
3[3] Anton , P. A., Cranston , R. D., Kashuba , A., Hendrix , C. W., Bumpus , N. N., Harman , N. R., Elliott , J., Janocko , L., Khanukhova , E., Dennis , R., Cumberland , W. G., Ju , C., Dieguez , A. C., Mauck , C., & Mc Gowan , I. (2012). RMP-02/MTN-006: A phase rectal safety, acceptability, pharmacokinetic, and pharmacodynamic study of tenofovir 1 % percent 1 1\% gel compared with oral tenofovir disoproxil fumarate. AIDS Res Hum Retroviruses 28 (11), 1412-1421.
4[4] Baeten , J. M., Donnell , D., Ndase , P., et al . (2012). Antiretroviral prophylaxis for HIV prevention in heterosexual men and women. N Engl J Med , 367 , 399-410.
5[5] Beskos , A., Jasra , A., Kantas , N. & Thiery , A. (2016). On the convergence of adaptive sequential Monte Carlo. Ann. Appl. Probab. , 26 , 1111-1146.
6[6] Boerebach , B. C., Lombarts , K. M., keijzer , C., Heineman , M. J., & Arah . O. A. (2012). The teacher, the physician and the person: how faculty’s teaching performance influences their role modeling. P Lo S One 7:e 32089. doi:10.1371/journal.pone.0032089.
7[7] Bollen , K. (1989). Structural Equation Models with Latent Variables . Wiley: New York.
8[8] Carvalho , C. M. & West , M. (2007). Dynamic Matrix-Variate Graphical Models. Bayesian Anal. , 2 , 69-98.