1.1. Background and Motivation

Learning dependence relations between variables is a pervasive issue in many applied domains, such as biology, social sciences, and notably psychology (Briganti et al., Reference Briganti, Scutari and McNally2022; Isvoranu et al., Reference Isvoranu, Epskamp, Waldorp and Borsboom2022). In the latter context, the recent field of network psychometrics considers a network-based approach to represent psychological constructs and understand directed interactions between behavioral, cognitive and biological factors, possibly allowing for causal interpretations (Borsboom et al., Reference Borsboom, Deserno, Rhemtulla, Epskamp, Fried, McNally, Robinaugh, Perugini, Dalege, Costantini, Isvoranu, Wysocki, van Borkulo, van Bork and Waldorp2021). The 2022 Psychometrika special issue “Network Psychometrics in action” promoted the development of statistical methods for network modeling motivated by psychological problems, and collected several contributions to the field, covering both methodological and applied aspects (Marsman and Rhemtulla, Reference Marsman and Rhemtulla2022).

Early works in the network psychometrics area were conceived to support psychologists in providing insights on various psychological phenomena, such as those at the basis of psychopathology, and in particular the study of comorbidity and mental disorders (Borsboom, Reference Borsboom2008; Cramer et al., Reference Cramer, Waldorp, Van der Maas and Borsboom2010). Typical research questions thus relate to the identification of direct dependencies between manifest variables, differences in the underlying dependence structure across available groups of patients, or even to the design of clinical interventions based on an estimated network. Crucial to these purposes is the development of models that allow to infer a plausible network structure for the available data and to provide a coherent quantification of the uncertainty related to directed links, specific paths or the whole network structure. In this regard, methodologies that fully account for network uncertainty lead to parameter estimates that are more robust w.r.t. possible network-model misspecifications; see Haslbeck and Waldorp (Reference Haslbeck and Waldorp2018), Epskamp et al. (Reference Epskamp, Kruis and Marsman2017) and Marsman et al. (Reference Marsman, Huth, Waldorp and Ntzoufras2022) for recent contributions in this area inspired by psychological problems.

All of the issues introduced above have motivated the development of dedicated statistical methodologies, based both on a frequentist and on a Bayesian paradigm. In particular, probabilistic graphical models both based on undirected or directed networks provide an effective tool to infer conditional dependence relations from the data (Cowell et al., Reference Cowell, Dawid, Lauritzen and Spiegelhalter1999; Edwards, Reference Edwards2000). Additionally, directed acyclic graphs (DAGs) offer a powerful framework for causal reasoning, even from observational, namely non-experimental, studies and specifically to quantify effects of hypothetical interventions on target variables w.r.t. outcome responses of interest; see Pearl (Reference Pearl2000) for a general introduction to causal inference based on DAGs, Maathuis and Nandy (Reference Maathuis, Nandy, Bühlmann, Drineas, Kane and van der Laan2016) for a review. The next section offers an overview of the main recent contributions to graphical modeling.

1.2. Literature Review

From a statistical perspective, learning a network of dependencies from the data is a model selection problem also known as structure learning. Several related methodologies that can deal with Gaussian and categorical data separately have been proposed. Specifically, score-based methods implement score functions for network estimation, such as based on penalized maximum likelihood estimators (Friedman et al., Reference Friedman, Hastie and Tibshirani2008; Meinshausen and Bühlmann, Reference Meinshausen and Bühlmann2006), or marginal likelihoods for methodologies following a Bayesian perspective (Chickering, Reference Chickering2002; Heckerman et al., Reference Heckerman, Geiger and Chickering1995). Moreover, constraint-based methods implement conditional independence tests to learn the set of (in)dependence constraints characterizing the underlying DAG structure, as in the popular PC algorithm (Kalisch and Bühlmann, Reference Kalisch and Bühlmann2007; Spirtes et al., Reference Spirtes, Glymour and Scheines2000). On the other hand, Bayesian methodologies adopt Markov Chain Monte Carlo (MCMC) methods to approximate a posterior distribution over the space of network structures, or related features of interest; see, for instance, Castelletti et al. (Reference Castelletti, Consonni, Della Vedova and Peluso2018) and Castelletti and Peluso (Reference Castelletti and Peluso2021) for, respectively, Gaussian and categorical settings, Ni et al. (Reference Ni, Baladandayuthapani, Vannucci and Stingo2022) for a recent overview of Bayesian methods for structure learning with applications to biological problems.

Mixed-type data, i.e., observations from variables of different parametric families, are very common in many contexts and especially psychological studies, where ordinal, discrete and continuous measurements are simultaneously collected on subjects. In particular, dealing with ordinal measurements within multivariate models is a pervasive issue in psychometrics, as clinical psychological data typically include polytomous items measured on ordinal or Likert-type scales. In this context, item response theory (IRT) based on factor analysis assumes that ordered-categorical responses are discrete representations of continuous latent scores. From a modeling perspective observed categorical data are then assumed to generate from discretization of latent data, as in the classical probit model and its extensions to ordinal responses. See in particular Wirth and Edwards (Reference Wirth and Edwards2007) for a review of item factor analysis within a structural equation modeling framework.

A few methodologies for structure learning from mixed data have been proposed. Harris and Drton (Reference Harris and Drton2013) introduce the rank PC, an extension of the original PC algorithm to nonparanormal models, namely based on a semi-parametric latent Gaussian copula model, with purely continuous marginal distributions. Moreover, Cui et al. (Reference Cui, Groot, Heskes, Frasconi, Landwehr, Manco and Vreeken2016) propose the Copula PC, an adaptation of the PC algorithm to a mixture of discrete and continuous data assumed to be drawn from a Gaussian copula model. Cui et al. (Reference Cui, Groot and Heskes2018) extend the previous method to deal with data that are missing at random. Similar ideas, for the case of undirected graphs, are also considered by Müller and Czado (Reference Müller and Czado2019) and He et al. (Reference He, Zhang, Wang and Zhang2017). Still in the context of directed graphs, a more recent methodology for structure learning given both categorical and Gaussian data is proposed by Andrews et al. (Reference Andrews, Ramsey and Cooper2018). The authors introduce a mixed-variable polynomial score based on the notion of Conditional Gaussian (CG) distribution (Lauritzen and Wermuth, Reference Lauritzen and Wermuth1989), then extended to a highly scalable algorithm by Andrews et al. (Reference Andrews, Ramsey and Cooper2019). Conditional Gaussian distributions are also adopted for structure learning of undirected graphs by Lee and Hastie (Reference Lee, Hastie, Carvalho and Ravikumar2013) and Cheng et al. (Reference Cheng, Li, Levina and Zhu2017) who implement penalized likelihoods and regression models with weighted lasso penalties, respectively. In a Bayesian setting, Bhadra et al. (Reference Bhadra, Rao and Baladandayuthapani2018) propose a unified framework for both categorical and Gaussian data based on Gaussian scale mixtures.

One main difficulty in developing statistical models for general mixed-type data is related to the non-standard joint support of the available variables. Typically however, interest lies in estimating dependence parameters of the joint distribution, corresponding to a network structure or correlation-type measures, rather than parameters indexing the marginal distributions of the variables. In this context, copula models, which allow to model the two sets of parameters separately, can provide an effective solution for statistical inference of network models. In addition, semiparametric copula models lacks any parametric assumption on the marginal c.d.f.’s which are estimated through their empirical distributions (Hoff, Reference Hoff2007). Contributions to copula graphical modeling based on undirected graphs are provided by Dobra and Lenkoski (Reference Dobra and Lenkoski2011) and Mohammadi et al. (Reference Mohammadi, Abegaz, van den Heuvel and Wit2017).

1.3. Contribution and Structure of the Paper

We propose a novel methodology for structure learning of networks which applies to mixed data, i.e., comprising continuous, categorical as well as discrete and ordinal measurements. Specifically, we consider a Gaussian copula DAG model in which each observable is associated to a latent counterpart, and the dependence parameter (covariance matrix) among latent variables reflects the conditional independencies imposed by a directed acyclic graph. We consider a Bayesian framework and proceed by assigning suitable prior distributions to DAG structures and DAG-dependent parameters. Inference is carried out by implementing an MCMC scheme which approximates the posterior distribution over network structures and covariance matrices. The main contributions of the proposed method can be summarized as follows: (i) we introduce a Bayesian framework for the analysis of complex dependence relations in multivariate settings characterized by mixed data; (ii) we provide a coherent quantification of the uncertainty around the estimated network or features of interest such as directed links, and a full posterior distribution of the underlying dependence parameter (correlation matrix), possibly summarized by Bayesian Model Averaging (BMA) estimates; (iii) our model allows to incorporate prior knowledge of the underlying network in terms of a partial ordering of the variables or edge orientations that are known in advance, thus improving DAG identification and enhancing causal inference.

The rest of the paper is organized as follows. In Sect. 2, we introduce Gaussian graphical models based on DAGs and the copula DAG model that we adopt for the analysis of mixed data. Section 3 completes our Bayesian model formulation by assigning prior distributions to DAG structures and DAG-model parameters. We implement in Sect. 4 an MCMC scheme which approximates the posterior distribution of DAGs and parameters. Our method is evaluated through extensive simulation experiments in Sect. 5. Section 6 is devoted to empirical studies, including the analysis of well-being data from a social survey promoted by the United Nations and mental health data collected from a cohort of medical students. In Sect. 7, we finally provide a discussion together with possible extensions of the proposed method to heterogeneous settings and latent trait models. Additional simulation results, comparisons with alternative methods and examples of MCMC diagnostics of convergence are included in Appendix.

2.1. Directed Acyclic Graphs

A Directed Acyclic Graph (DAG) is a pair $D=(V,E)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}=(V,E)$$\end{document} consisting of a set of vertices (or nodes) $V=\{1,\dots ,q\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V=\{1,\ldots ,q\}$$\end{document} and a set of directed edges $E\subseteq V\times V$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E \subseteq V \times V$$\end{document} . For any two nodes $u,v\in V$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u,v \in V$$\end{document} , we denote an edge from u to v as (u, v) or $u\to v$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u\rightarrow v$$\end{document} indifferently; also, the set E is such that if $(u,v)\in E$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(u,v)\in E$$\end{document} then $(v,u)\notin E$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(v,u)\notin E$$\end{document} . A sequence of nodes $(v_1,v_2,\dots ,v_k)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(v_1,v_2, \ldots ,v_k)$$\end{document} is a path if there exists $v_1\to v_2\to \cdots \to v_k$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_1\rightarrow v_2\rightarrow \cdots \rightarrow v_k$$\end{document} in $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document} . We assume that $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document} does not contain cycles, that is paths such that $v_1\equiv v_k$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_1\equiv v_k$$\end{document} . For a given node $v\in V$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v\in V$$\end{document} we let ${\text{pa}}_D\left(v\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {pa}_{\mathcal {D}}(v)$$\end{document} be the set of parents of v in $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document} , i.e., the set of all nodes u such that $(u,v)\in E$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(u,v)\in E$$\end{document} . Moreover, we say that u is a descendant of v if there exists a path from v to u; by converse, v is an ancestor of u. The set of all descendants and ancestors of a node v in $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document} are ${\text{de}}_D\left(v\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {de}_{\mathcal {D}}(v)$$\end{document} and ${\text{an}}_D\left(v\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {an}_{\mathcal {D}}(v)$$\end{document} , respectively.

A DAG encodes a set of conditional independencies of the form $A\perp \perp B|C$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A\bot \hspace{-6pt}\bot B \,\vert \,C$$\end{document} , reading as “A and B are conditionally independent given C”, where A, B, C are disjoint subsets of the vertex set V. The set of all conditional independencies characterizing the DAG determines the DAG Markov property and can be read-off from the graph using graphical criteria such as d-separation (Pearl, Reference Pearl2000). In particular, each node is conditionally independent from its non-descendants given its parents. Simple examples are provided in Fig. 1, where using d-separation it is possible to show that $u\perp \perp z|v$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u \bot \hspace{-6pt}\bot z \,\vert \,v$$\end{document} in both $D_1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}_1$$\end{document} and $D_2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}_2$$\end{document} ; differently, $u\perp \perp z$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u \bot \hspace{-6pt}\bot z$$\end{document} in $D_3$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}_3$$\end{document} meaning that u and z are marginally independent. We refer the reader to Lauritzen (Reference Lauritzen1996) for further notions on graph theory.

Figure 1 Three DAGs on the set of nodes $V=\{u,v,z\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V=\{u,v,z\}$$\end{document} . $D_1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}_1$$\end{document} and $D_2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}_2$$\end{document} encode the conditional independence $u\perp \perp z|v$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u \bot \hspace{-5pt}\bot z \,\vert \,v$$\end{document} . In $D_3$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}_3$$\end{document} we instead have $u\perp \perp z$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u \bot \hspace{-5pt}\bot z$$\end{document} .

2.2. Gaussian DAG Models

Let $D=(V,E)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}=(V,E)$$\end{document} be a DAG and $z={(Z_1,\dots ,Z_q)}^{\top }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{z}=(Z_1,\ldots ,Z_q)^{\top }$$\end{document} a collection of q real-valued random variables, each associated with a node in $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document} , and with joint p.d.f. $f(·)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f(\cdot )$$\end{document} . We assume that

(1) $\begin{array}{c}{Z}_1,\dots ,Z_q|\Omega ,D\sim N_q(0,{\Omega }^{-1}),\Omega \in P_D,\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} Z_1,\ldots ,Z_q \,\vert \,\varvec{\Omega },\mathcal {D}\sim \mathcal {N}_q (\varvec{0}, \varvec{\Omega }^{-1}), \quad \varvec{\Omega }\in \mathcal {P}_{\mathcal {D}}, \end{aligned}$$\end{document}

where $\Omega$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Omega }$$\end{document} is the precision matrix (inverse of the covariance matrix $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} ) and $P_D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {P}_{\mathcal {D}}$$\end{document} denotes the set of all symmetric positive definite (s.p.d.) precision matrices Markov w.r.t. DAG $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document} . Accordingly, we impose to $\Omega$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Omega }$$\end{document} the conditional independencies encoded by $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document} that are deducible from d-separation (Sect. 2.1).

An equivalent representation of Model (1), useful for later developments, is given by the allied Structural Equation Model (SEM). To this end, let $D=\text{diag}(D_{11},\dots ,D_{\mathrm{qq}})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{D}=\text {diag}(\varvec{D}_{11},\ldots ,\varvec{D}_{qq})$$\end{document} and $L$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{L}$$\end{document} a (q, q) matrix of (regression) coefficients with diagonal elements equal to 1 and (u, v)-element $L_{u,v}\ne 0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{L}_{u,v}\ne 0$$\end{document} if and only if $u\to v$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u\rightarrow v$$\end{document} in $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document} . Nonzero elements of $L$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{L}$$\end{document} correspond to directed links between nodes, while zero entries to missing edges in the DAG. Accordingly, $L$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{L}$$\end{document} resembles the DAG structure and the set of all parent–child relations characterizing its Markov property. A Gaussian SEM can be written as

(2) $\begin{array}{c}{L}^{\top }z=\epsilon ,\epsilon \sim N_q(0,D),\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{L}^{\top }\varvec{z}= \varvec{\varepsilon }, \quad \varvec{\varepsilon } \sim \mathcal {N}_q(\varvec{0},\varvec{D}), \end{aligned}$$\end{document}

where $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{D}$$\end{document} corresponds to the covariance matrix of the error terms, which is assumed to be diagonal and collects the conditional variances of $(Z_1,\dots ,Z_q)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(Z_1,\ldots ,Z_q)$$\end{document} . Moreover, Eq. (2) implies $V\text{ar}\left(z\right)=\Sigma =L^{-\top }D{L}^{-1}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {V}\!\text {ar}(\varvec{z})=\varvec{\Sigma }=\varvec{L}^{-\top }\varvec{D}\varvec{L}^{-1}$$\end{document} ; equivalently, $\Omega =L{D}^{-1}{L}^{\top }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Omega }=\varvec{L}\varvec{D}^{-1}\varvec{L}^{\top }$$\end{document} . The latter decomposition provides a re-parameterization of $\Omega$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Omega }$$\end{document} in terms of $(D,L)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\varvec{D},\varvec{L})$$\end{document} . From (2) we have, for each $j=1,\dots ,q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1,\ldots , q$$\end{document} , $Z_j=-L_{\prec j]}^{\top }{z}_{{\text{pa}}_D\left(j\right)}+{\epsilon }_j,$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_j = -\varvec{L}_{\prec j\,]}^{\top }\varvec{z}_{\text {pa}_{\mathcal {D}}(j)} + \varepsilon _j, $$\end{document} with ${\epsilon }_j\sim N(0,D_{\mathrm{jj}})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon _j\sim \mathcal {N}(0,\varvec{D}_{jj})$$\end{document} , where $\prec j]={\text{pa}}_D\left(j\right)\times j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\prec j\,]\,=\text {pa}_{\mathcal {D}}(j)\times j$$\end{document} and $L_{A\times B}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{L}_{A\times B}$$\end{document} denotes the sub-matrix of $L$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{L}$$\end{document} with elements belonging to rows and columns indexed by A and B, respectively. Each equation above resembles the structure of a linear regression model for variable $Z_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_j$$\end{document} , with $-L_{\prec j]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\varvec{L}_{\prec j\,]}$$\end{document} corresponding to the regression coefficients associated with variables in $z_{{\text{pa}}_D\left(j\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{z}_{\text {pa}_{\mathcal {D}}(j)}$$\end{document} , namely the parents of node/variable $Z_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_j$$\end{document} ; see also Sect. 3.2. for a comparison with the Bayesian analysis of normal linear regression models. Accordingly, Model (1) can be equivalently written as

(3) $\begin{array}{c}f(z_1,\dots ,z_q|D,L,D)=\prod _{j=1}^{q}dN(z_j|-L_{\prec j]}^{\top }{z}_{{\text{pa}}_D\left(j\right)},D_{\mathrm{jj}}),\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} f(z_1,\ldots ,z_q\,\vert \,\varvec{D},\varvec{L},\mathcal {D})=\prod _{j=1}^{q}\,d\,\mathcal {N}\big (z_j\,\vert \,-\varvec{L}_{\prec j\,]}^{\top }\varvec{z}_{\text {pa}_{\mathcal {D}}(j)}, \varvec{D}_{jj}\big ), \end{aligned}$$\end{document}

where $dN(·|\mu ,{\sigma }^2)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d\,\mathcal {N}(\cdot \,\vert \,\mu ,\sigma ^2)$$\end{document} denotes the p.d.f. of a univariate $N(\mu ,{\sigma }^2)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {N}(\mu ,\sigma ^2)$$\end{document} . Finally, given n i.i.d. samples from (3), $z_i={(z_{i,1},\dots ,z_{i,q})}^{\top }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{z}_i=(z_{i,1},\ldots ,z_{i,q})^{\top }$$\end{document} , $i=1,\dots ,n$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,\ldots ,n$$\end{document} , collected in the (n, q) matrix $Z$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Z}$$\end{document} (row-binding of the $z_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{z}_i$$\end{document} ’s), the likelihood function can be written as

(4) $\begin{array}{c}f\left(Z\right|D,L,D)=\prod _{i=1}^n\left(\prod _{j=1}^q,d,,N,(,z_{i,j},|-,L_{\prec j]}^{\top },z_{i,{\text{pa}}_D\left(j\right)},,,D_{\mathrm{jj}},)\right).\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} f(\varvec{Z}\,\vert \,\varvec{D},\varvec{L}, \mathcal {D})= \prod _{i=1}^{n}\left\{ \prod _{j=1}^{q}d\,\mathcal {N}\big (z_{i,j}\,\vert \,-\varvec{L}_{\prec j\,]}^{\top }\varvec{z}_{i,\text {pa}_{\mathcal {D}}(j)}, \varvec{D}_{jj}\big )\right\} . \end{aligned}$$\end{document}

2.3. Copula DAG Models

Consider now a collection of q random variables, $X_1,\dots ,X_q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_1,\ldots ,X_q$$\end{document} , comprising binary, ordinal, continuous or count variables, each with marginal cumulative distribution function (c.d.f.)  $F_j(·)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_j(\cdot )$$\end{document} , $j=1,\dots ,q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1,\ldots ,q$$\end{document} . In what follows we will consider a collection of n q-dimensional observations from $X_1,\dots ,X_q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_1,\ldots ,X_q$$\end{document} . To model these mixed data we need to specify a joint distribution for $X_1,\dots ,X_q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_1,\ldots , X_q$$\end{document} that we define through a Gaussian copula DAG model. Specifically, let $Z_1,\dots ,Z_q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_1,\ldots ,Z_q$$\end{document} be a collection of q latent random variables with joint Gaussian distribution as in (1), We establish a link between each observed variable $X_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_j$$\end{document} and its latent counterpart $Z_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_j$$\end{document} by assuming that

(5) $\begin{array}{c}{X}_j=F_j^{-1}\left(\Phi ,\left(,Z_j,\right)\right),\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} X_j = F_j^{-1}\left\{ \Phi (Z_j)\right\} , \end{aligned}$$\end{document}

where $F_j^{-1}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F^{-1}_j$$\end{document} is the (pseudo) inverse c.d.f. of $X_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_j$$\end{document} and $\Phi \left(Z_j\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi (Z_j)$$\end{document} the c.d.f. of a standard Normal distribution. The joint c.d.f. of $X_1,\dots ,X_q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_1,\ldots ,X_q$$\end{document} can be written as

(6) $\begin{array}{cc}& P(X_1\le x_1,\dots ,X_q\le x_q|\Omega ,F_1,\dots ,F_q)\\ & ={\Phi }_q\left({\Phi }^{-1},\left(F_1\left(x_1\right)\right),,,\dots ,,,{\Phi }^{-1},\left(F_q\left(x_q\right)\right),|\Omega \right),\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&P(X_1\le x_1,\ldots , X_q \le x_q\,\vert \,\varvec{\Omega },F_1,\ldots ,F_q) \nonumber \\&\quad = \Phi _q\left( \Phi ^{-1}(F_1(x_1)),\ldots ,\Phi ^{-1}(F_q(x_q)) \,\vert \,\varvec{\Omega }\right) , \end{aligned}$$\end{document}

where ${\Phi }_q(·|\Omega )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _q(\cdot \,\vert \,\varvec{\Omega })$$\end{document} denotes the c.d.f. of $N_q(0,{\Omega }^{-1})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {N}_q(\varvec{0},\varvec{\Omega }^{-1})$$\end{document} in (1). Also notice that Model (6) depends on the marginal distributions $F_1,\dots ,F_q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_1,\ldots ,F_q$$\end{document} and (although not emphasized in the equation) their parameters which would need to be “estimated”. A semiparametric estimation strategy would replace $F_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_j$$\end{document} with the corresponding empirical estimates ${\hat{F}}_j\left(k_j\right)=n^{-1}{\sum }_{i=1}^{n}1(x_{i,j}<k_j),$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \widehat{F}_j(k_j)=n^{-1}\sum _{i=1}^{n}\mathbbm {1}(x_{i,j}<k_j), $$\end{document} where $k_j\in \text{unique}\{x_{1,j},\dots ,x_{n,j}\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_j\in \text {unique}\{x_{1,j},\ldots ,x_{n,j}\}$$\end{document} . A comparison with a parametric strategy based on probabilistic-model assumptions for the marginal c.d.f.’s is instead provided in Appendix.

As an alternative to the estimation procedures above, Hoff (Reference Hoff2007) proposes a rank-based nonparametric approach, that we also employ in our methodology. Specifically, let $x_i={(x_{i,1},\dots ,x_{i,q})}^{\top }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}_i=(x_{i,1},\ldots ,x_{i,q})^{\top }$$\end{document} , $i=1,\dots ,n$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,\ldots ,n$$\end{document} , be n i.i.d. samples from (6) and $X$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{X}$$\end{document} the (n, q) data matrix. Since the $F_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_j$$\end{document} ’s are non-decreasing, for each pair of distinct observations $x_{i,j}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{i,j}$$\end{document} and $x_{l,j}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{l,j}$$\end{document} , if $x_{i,j}<x_{l,j}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{i,j}<x_{l,j}$$\end{document} then $z_{i,j}<z_{l,j}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_{i,j}<z_{l,j}$$\end{document} . Therefore, observing $X$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{X}$$\end{document} implies that the latent data $Z$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Z}$$\end{document} must lie in the set

(7) $\begin{array}{c}A\left(X\right)=\left(Z,\in ,R^{n\times q},:,\text{max},\left(z_{k,j},:,x_{k,j},<,x_{i,j}\right),<,z_{i,j},<,\text{max},\left(z_{k,j},:,x_{i,j},<,x_{k,j}\right)\right)\\ \end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} A(\varvec{X}) = \left\{ \varvec{Z}\in \mathbb {R}^{n\times q}: \text {max}\left\{ z_{k,j}:x_{k,j}<x_{i,j}\right\}< z_{i,j}< \text {max}\left\{ z_{k,j}:x_{i,j}<x_{k,j} \right\} \right\} \nonumber \\ \end{aligned}$$\end{document}

and one can take the occurrence of such an event as the data. Thus, the extended rank likelihood (Hoff, Reference Hoff2007) can be written as

(8) $\begin{array}{c}p(Z\in A\left(X\right)|D,L,D)={\int }_{A\left(X\right)}f\left(Z\right|D,L,D)\text{d}Z,\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} p(\varvec{Z}\in A(\varvec{X})\,\vert \,\varvec{D}, \varvec{L}, \mathcal {D})=\int _{A(\varvec{X})}f(\varvec{Z}\,\vert \,\varvec{D}, \varvec{L}, \mathcal {D})\, {\text {d}}\varvec{Z}, \end{aligned}$$\end{document}

with $f\left(Z\right|D,L,D)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f(\varvec{Z}\,\vert \,\varvec{D}, \varvec{L}, \mathcal {D})$$\end{document} as in Eq. (4).

Our model formulation assumes that a DAG Markov property holds at a latent-variable stage, namely between $Z_1,\dots ,Z_q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_1,\ldots ,Z_q$$\end{document} , in force of factorization (3). Through the copula-transfer link (5), this translates into a Markov property for $X_1,\dots ,X_q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_1,\ldots ,X_q$$\end{document} , provided that all the marginal distributions are continuous (Liu et al., Reference Liu, Lafferty and Wasserman2009). The presence of non-continuous variables (e.g., ordinal, discrete) might induce additional dependencies among the observed variables (w.r.t. those included in the latent model); however, such dependencies can be regarded as having a secondary relevance since they emerge from the marginals, rather than from the joint distribution; see in particular Dobra & Lenkoski (Reference Dobra and Lenkoski2011, Section 4) for a related discussion.

3.1. Prior on $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document}

Let $S_q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {S}_q$$\end{document} be the set of all DAG structures on q nodes. In many applied problems, exact knowledge about the orientation of some edges, whenever present in the graph, is typically available and accordingly one would like to incorporate such information in the model. Without loss of generality, let $S_q^C$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {S}_q^C$$\end{document} be the set of all DAGs on q nodes satisfying a set of structural constraints C, corresponding to edge orientations that are known in advance (equivalently, a set of “reversed” edge orientations that are forbidden). As an example, suppose node u represents a response variable of interest, so that there are no outgoing edges from u (equivalently, u cannot have children). In such a case we have $C=\{u\nrightarrow v|v=1,\dots ,q,v\ne u\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C=\{u\not \rightarrow v \,\vert \,v=1,\ldots ,q, v \ne u\}$$\end{document} and accordingly an edge between u and v whenever present will be oriented as $v\to u$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v\rightarrow u$$\end{document} . See also Sect. 6 for examples on real-data. We assign a prior to DAGs belonging to $S_q^C$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {S}_q^C$$\end{document} as follows.

For a given DAG $D=(V,E)\in S_q^C$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}=(V,E)\in \mathcal {S}_q^C$$\end{document} , let $S^D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{S}^{\mathcal {D}}$$\end{document} be the $0-1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0-1$$\end{document} adjacency matrix of its skeleton, that is the underlying undirected graph obtained after removing the orientation of all its edges. For each (u, v)-element of $S^D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{S}^{\mathcal {D}}$$\end{document} , we have $S_{u,v}^D=1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{S}_{u,v}^{\mathcal {D}}=1$$\end{document} if and only if $(u,v)\in E$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(u,v)\in E$$\end{document} or $(v,u)\in E$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(v,u)\in E$$\end{document} , zero otherwise. Conditionally on a prior probability of inclusion $\pi \in (0,1)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi \in (0,1)$$\end{document} we assume for each $u>v$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u>v$$\end{document} , $S_{u,v}^D|\pi \stackrel{\text{iid}}{\sim }\text{Ber}\left(\pi \right),$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \varvec{S}_{u,v}^{\mathcal {D}} \,\vert \,\pi \overset{\text {iid}}{\sim }\text {Ber}(\pi ), $$\end{document} which implies

(9) $\begin{array}{c}p\left(S^D\right|\pi )={\pi }^{|S^D|}{(1-\pi )}^{\frac{q(q-1)}{2}-\left|S^D\right|},\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} p(\varvec{S}^{\mathcal {D}}\,\vert \,\pi )=\pi ^{|\varvec{S}^{\mathcal {D}}|} (1-\pi )^{\frac{q(q-1)}{2}-|\varvec{S}^{\mathcal {D}}|}, \end{aligned}$$\end{document}

where $|S^D|$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|\varvec{S}^{\mathcal {D}}|$$\end{document} is the number of edges in $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document} (equivalently in its skeleton) and $q(q-1)/2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q(q-1)/2$$\end{document} is the maximum number of edges in a DAG on q nodes. We then assume $\pi \sim \text{Beta}(c,d)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi \sim \text {Beta}(c,d)$$\end{document} , so that, by integrating out $\pi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi $$\end{document} , the resulting prior on $S^D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{S}^{\mathcal {D}}$$\end{document} is

(10) $\begin{array}{c}p\left(S^D\right)=\frac{\Gamma \left(|,S^D,|+c\right)\Gamma \left(\frac{q(q-1)}{2},-,\left|S^D\right|,+,d\right)}{\Gamma \left(\frac{q(q-1)}{2},+,c,+,d\right)}·\frac{\Gamma \left(c,+,d\right)}{\Gamma \left(c\right)\Gamma \left(d\right)}.\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} p(\varvec{S}^{\mathcal {D}}) = \frac{\Gamma \left( |\varvec{S}^{\mathcal {D}}| + c\right) \Gamma \left( \frac{q(q-1)}{2} - |\varvec{S}^{\mathcal {D}}| + d \right) }{\Gamma \left( \frac{q(q-1)}{2} + c + d\right) } \cdot \frac{\Gamma \left( c + d\right) }{\Gamma \left( c\right) \Gamma \left( d\right) }. \end{aligned}$$\end{document}

Finally, we set $p\left(D\right)\propto p\left(S^D\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ p(\mathcal {D})\propto p(\varvec{S}^{\mathcal {D}})$$\end{document} for each $D\in S_q^C$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}\in \mathcal {S}_q^C $$\end{document} . See also Scott and Berger (Reference Scott and Berger2010) for a comparison with multiplicity correction priors adopted in a linear model selection setting. Hyperparameters c and d can be chosen to reflect a prior knowledge of sparsity in the graph, whenever available; in particular any choice $c<d$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c < d$$\end{document} will imply $E\left(\pi \right)<0.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {E}(\pi )<0.5$$\end{document} , thus favoring sparse graphs. The default choice $c=d=1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c=d=1$$\end{document} , which corresponds to $\pi \sim \text{Unif}(0,1)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi \sim \text {Unif}(0,1)$$\end{document} , can be instead adopted in the absence of substantive prior information.

3.2. Prior on $(D,L)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({\varvec{D}},{\varvec{L}})$$\end{document}

Conditionally on DAG $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document} we assign a prior to $(D,L)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\varvec{D},\varvec{L})$$\end{document} through a DAG-Wishart distribution with position hyperparameter $U$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{U}$$\end{document} (a (q, q) s.p.d. matrix) and shape hyperparameter $a^D={(a_1^D,\dots ,a_q^D)}^{\top }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a^{\mathcal {D}}=(a_1^{\mathcal {D}},\ldots ,a_q^{\mathcal {D}})^{\top }$$\end{document} . The DAG-Wishart distribution (Cao et al., Reference Cao, Khare and Ghosh2019) provides a conjugate prior for the Gaussian DAG model (3). Accordingly, conditionally on the latent data $Z$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Z}$$\end{document} , the posterior distribution of $(D,L)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\varvec{D},\varvec{L})$$\end{document} as well as the marginal likelihood of the model are available in closed form; see Sect. 4 for more details. A feature of the DAG-Wishart distribution is also that node-parameters $\{(D_{\mathrm{jj}},L_{\prec j]}),j=1,\dots ,q\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\big \{ (\varvec{D}_{jj}, \varvec{L}_{\prec j\,]} ), \, j=1,\ldots ,q \big \}$$\end{document} are a priori independent with distribution

(11) $\begin{array}{cc}{D}_{\mathrm{jj}}|D& \sim \text{I-Ga}\left(\frac{1}{2},a_j^D,,,\frac{1}{2},U_{j|{\text{pa}}_D\left(j\right)}\right),\\ L_{\prec j]}|D_{\mathrm{jj}},D& \sim N_{|{\text{pa}}_D\left(j\right)|}\left(-,U_{\prec j\succ }^{-1},U_{\prec j]},,,D_{\mathrm{jj}},U_{\prec j\succ }^{-1}\right),\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{D}_{jj} \,\vert \,\mathcal {D}&\sim \,\,\, \text {I-Ga}\left( \frac{1}{2}a_j^{\mathcal {D}}, \frac{1}{2}\varvec{U}_{j\,\vert \,\text {pa}_{\mathcal {D}}(j)}\right) , \nonumber \\ \varvec{L}_{\prec j\,]}\,\vert \,\varvec{D}_{jj}, \mathcal {D}&\sim \,\,\, \mathcal {N}_{|\text {pa}_{\mathcal {D}}(j)|}\left( -\varvec{U}_{\prec j \succ }^{-1}\varvec{U}_{\prec j\,]},\varvec{D}_{jj}\varvec{U}_{\prec j \succ }^{-1}\right) , \end{aligned}$$\end{document}

where $\text{I-Ga}(\alpha ,\beta )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {I-Ga}(\alpha ,\beta )$$\end{document} stands for an Inverse-Gamma distribution with shape $\alpha >0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha >0$$\end{document} and rate $\beta >0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta >0$$\end{document} having expectation $\beta /(\alpha -1)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta /(\alpha -1)$$\end{document} ( $\alpha >1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha >1$$\end{document} ). Moreover, $U_{j|{\text{pa}}_D\left(j\right)}=U_{\mathrm{jj}}-U_{[j\succ }{U}_{\prec j\succ }^{-1}{U}_{\prec j]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{U}_{j\,\vert \,\text {pa}_{\mathcal {D}}(j)} = \varvec{U}_{jj}-\varvec{U}_{[\, j\succ } \varvec{U}^{-1}_{\prec j\succ } \varvec{U}_{\prec j\,]}$$\end{document} , with $\prec j]={\text{pa}}_D\left(j\right)\times j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\prec j\,] = \text {pa}_{\mathcal {D}}(j)\times j$$\end{document} , $[j\succ =j\times {\text{pa}}_D\left(j\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[\, j\succ \, = j \times \text {pa}_{\mathcal {D}}(j)$$\end{document} , $\prec j\succ ={\text{pa}}_D\left(j\right)\times {\text{pa}}_D\left(j\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\prec j\succ \, = \text {pa}_{\mathcal {D}}(j)\times \text {pa}_{\mathcal {D}}(j)$$\end{document} .

With regard to hyperparameters $a_1^D,\dots ,a_q^D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_1^{\mathcal {D}},\ldots ,a_q^{\mathcal {D}}$$\end{document} we also consider the default choice $a_j^D=a+\left|{\text{pa}}_D\left(j\right)\right|-q+1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_j^{\mathcal {D}}=a+|\text {pa}_{\mathcal {D}}(j)|-q+1$$\end{document} $(a>q-1)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(a>q-1)$$\end{document} which guarantees compatibility (same marginal likelihood) among prior distributions for Markov equivalent DAGs; see also Peluso and Consonni (Reference Peluso and Consonni2020). Finally, the prior on $(D,L)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\varvec{D},\varvec{L})$$\end{document} is given by

(12) $\begin{array}{c}p(D,L|D)=\prod _{j=1}^{q}p\left(L_{\prec j]}\right|D_{\mathrm{jj}})p\left(D_{\mathrm{jj}}\right).\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} p(\varvec{D},\varvec{L}\,\vert \,\mathcal {D})=\prod _{j=1}^{q}p(\varvec{L}_{\prec j\,]}\,\vert \,\varvec{D}_{jj})\,p(\varvec{D}_{jj}). \end{aligned}$$\end{document}

A DAG-Wishart prior on $(D,L)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\varvec{D},\varvec{L})$$\end{document} implicitly assigns (independent) Normal-Inverse-Gamma distributions to each pair of node-parameters $(D_{\mathrm{jj}},L_{\prec j]})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\varvec{D}_{jj},\varvec{L}_{ \prec j\,]})$$\end{document} , a conditional variance and vector-regression coefficient for the j-th term of the SEM, as in the standard conjugate Bayesian analysis of a normal linear regression model. Moreover, under the default choice $U=g{I}_q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{U}=g\varvec{I}_q$$\end{document} (Sect. 5), with $g>0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g>0$$\end{document} and $I_q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{I}_q$$\end{document} the (q, q) identity matrix, it is easy to show that $E\left(L_{\prec j]}\right|D_{\mathrm{jj}},D)=0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {E}(\varvec{L}_{ \prec j\,]}\,\vert \,\varvec{D}_{jj}, \mathcal {D})=\varvec{0}$$\end{document} and $V\text{ar}\left(L_{\prec j]}\right|D_{\mathrm{jj}},D)=D_{\mathrm{jj}}/g{I}_{|{\text{pa}}_D\left(j\right)|}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {V}\!\text {ar}(\varvec{L}_{ \prec j\,]}\,\vert \,\varvec{D}_{jj}, \mathcal {D})=\varvec{D}_{jj}/g\varvec{I}_{|\text {pa}_{\mathcal {D}}(j)|}$$\end{document} for each $j=1,\dots ,q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1,\ldots ,q$$\end{document} , so that priors on regression coefficients are centered at zero and with diagonal covariance matrix reflecting an assumption of prior independence across elements of $L_{\prec j]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{L}_{ \prec j\,]}$$\end{document} ; moreover, smaller values of g make such prior less informative; we refer to Sect. 5 for details about the choice of hyperparameters.

4.1. Update of $(D,L,D)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({\varvec{D}},{\varvec{L}},\mathcal {D})$$\end{document}

The joint full conditional distribution of $(D,L,D)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\varvec{D},\varvec{L},\mathcal {D})$$\end{document} is given by

(14) $\begin{array}{c}p(D,L,D|X,Z)=p(D,L,D\left|Z\right)\propto f\left(Z\right|D,L,D\left)p\right(D,L\left|D\right)p\left(D\right).\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} p(\varvec{D},\varvec{L},\mathcal {D}\,\vert \,\varvec{X},\varvec{Z}) = p(\varvec{D},\varvec{L},\mathcal {D}\,\vert \,\varvec{Z}) \propto f(\varvec{Z}\,\vert \,\varvec{D},\varvec{L},\mathcal {D})p(\varvec{D},\varvec{L}\,\vert \,\mathcal {D})p(\mathcal {D}). \end{aligned}$$\end{document}

Conditionally on the latent data $Z$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Z}$$\end{document} , we can sample from $p(D,L,D|Z)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(\varvec{D},\varvec{L},\mathcal {D}\,\vert \,\varvec{Z})$$\end{document} using the MCMC scheme for posterior inference of Gaussian DAGs presented in Castelletti and Consonni (Reference Castelletti and Consonni2021, Supplementary Material). The latter consists of a partial analytic structure (PAS) algorithm (Godsill, Reference Godsill2012) based on the two following steps.

Given DAG $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document} , parameters $(D,L)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\varvec{D},\varvec{L})$$\end{document} are first sampled from their full conditional distribution, which, because of conjugacy of the DAG-Wishart prior with the distribution of the latent data, is such that, for $j=1,\dots ,q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1,\ldots ,q$$\end{document} ,

(15) $\begin{array}{cc}{D}_{\mathrm{jj}}|D,Z& \sim \text{I-Ga}\left(\frac{1}{2},{\tilde{a}}_j^D,,,\frac{1}{2},{\tilde{U}}_{j|{\text{pa}}_D\left(j\right)}\right),\\ L_{\prec j]}|D_{\mathrm{jj}},D,Z& \sim N_{|{\text{pa}}_D\left(j\right)|}\left(-,{\tilde{U}}_{\prec j\succ }^{-1},{\tilde{U}}_{\prec j]},,,D_{\mathrm{jj}},{\tilde{U}}_{\prec j\succ }^{-1}\right),\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{D}_{jj} \,\vert \,\mathcal {D}, \varvec{Z}&\sim \,\,\, \text {I-Ga}\left( \frac{1}{2}\widetilde{a}_j^{\mathcal {D}}, \frac{1}{2}\widetilde{\varvec{U}}_{j\,\vert \,\text {pa}_{\mathcal {D}}(j)}\right) , \nonumber \\ \varvec{L}_{\prec j\,]}\,\vert \,\varvec{D}_{jj}, \mathcal {D}, \varvec{Z}&\sim \,\,\, \mathcal {N}_{|\text {pa}_{\mathcal {D}}(j)|}\left( -\widetilde{\varvec{U}}_{\prec j \succ }^{-1}\widetilde{\varvec{U}}_{\prec j\,]},\varvec{D}_{jj}\widetilde{\varvec{U}}_{\prec j \succ }^{-1}\right) , \end{aligned}$$\end{document}

with ${\tilde{a}}_j^D=\tilde{a}+\left|{\text{pa}}_D\left(j\right)\right|-q+1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{a}_j^{\mathcal {D}} = \widetilde{a}+|\text {pa}_{\mathcal {D}}(j)|-q+1$$\end{document} , $\tilde{a}=a+n$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{a} = a+n$$\end{document} and $\tilde{U}=U+Z^{\top }Z$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{\varvec{U}} = \, \varvec{U}+ \varvec{Z}^\top \varvec{Z}$$\end{document} .

Next, update of DAG $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document} is performed through a Metropolis Hastings step in which a DAG $D^{\ast }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}^*$$\end{document} is drawn from a proposal distribution $q\left(D^{\ast }\right|D)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q(\mathcal {D}^*\,\vert \,\mathcal {D})$$\end{document} . For a given DAG $D\in S_q^C$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}\in \mathcal {S}_q^C$$\end{document} , the adopted proposal distribution is built upon the set of all DAGs belonging to $S_q^C$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {S}_q^C$$\end{document} that can be reached from $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document} through insertion, deletion or reversal of an edge in $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document} . Specifically, we construct the set of all direct successors of $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document} , $O_D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O}_{\mathcal {D}}$$\end{document} , and then draw uniformly a DAG $D^{\ast }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}^*$$\end{document} from $O_D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O}_{\mathcal {D}}$$\end{document} . It follows that $q\left(D^{\ast }\right|D)=1/\left|O_D\right|$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q(\mathcal {D}^*\,\vert \,\mathcal {D})=1/|\mathcal {O}_{\mathcal {D}}|$$\end{document} . Each proposed DAG then differs locally by a single edge (u, v) which is inserted (a), deleted (b) or reversed (c) in $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document} . It can be shown that the acceptance probability for $D^{\ast }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}^*$$\end{document} under a PAS algorithm is given by ${\alpha }_{D^{\ast }}=min\{1;r_{D^{\ast }}\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _{\mathcal {D}^*}=\min \{1;r_{\mathcal {D}^*}\}$$\end{document} with

(16) $\begin{array}{c}{r}_{D^{\ast }}=\left(\begin{array}{cc}\frac{m_{}\left(Z_v\right|Z_{{\text{pa}}_{D^{\ast }}\left(v\right)},D^{\ast })}{m_{}\left(Z_v\right|Z_{{\text{pa}}_D\left(v\right)},D)}·\frac{p\left(D^{\ast }\right)}{p\left(D\right)}·\frac{q\left(D\right|D^{\ast })}{q\left(D^{\ast }\right|D)}& \\ \text{if}\left(a\right)\text{or}\left(b\right)\\ \frac{m_{}\left(Z_v\right|Z_{{\text{pa}}_{D^{\ast }}\left(v\right)},D^{\ast })}{m_{}\left(Z_v\right|Z_{{\text{pa}}_D\left(v\right)},D)}·\frac{m_{}\left(Z_u\right|Z_{{\text{pa}}_{D^{\ast }}\left(u\right)},D^{\ast })}{m_{}\left(Z_u\right|Z_{{\text{pa}}_D\left(u\right)},D)}·\frac{p\left(D^{\ast }\right)}{p\left(D\right)}·\frac{q\left(D\right|D^{\ast })}{q\left(D^{\ast }\right|D)}& \\ \text{if}\left(c\right),\end{array}\right)\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} r_{\mathcal {D}^*}\,=\, \left\{ \begin{array}{rl} \dfrac{m_{}(\varvec{Z}_v\,\vert \,\varvec{Z}_{\text {pa}_{\mathcal {D}^*}(v)}, \mathcal {D}^*)}{m_{}(\varvec{Z}_v\,\vert \,\varvec{Z}_{\text {pa}_{\mathcal {D}}(v)}, \mathcal {D})} \cdot \dfrac{p(\mathcal {D}^*)}{p(\mathcal {D})} \cdot \dfrac{q(\mathcal {D}\,\vert \,\mathcal {D}^*)}{q(\mathcal {D}^*\,\vert \,\mathcal {D})} &{} \quad \\ \text { if } (a) \text { or } (b) \\ \dfrac{m_{}(\varvec{Z}_v\,\vert \,\varvec{Z}_{\text {pa}_{\mathcal {D}^*}(v)}, \mathcal {D}^*)}{m_{}(\varvec{Z}_v\,\vert \,\varvec{Z}_{\text {pa}_{\mathcal {D}}(v)}, \mathcal {D})} \cdot \dfrac{m_{}(\varvec{Z}_u\,\vert \,\varvec{Z}_{\text {pa}_{\mathcal {D}^*}(u)}, \mathcal {D}^*)}{m_{}(\varvec{Z}_u\,\vert \,\varvec{Z}_{\text {pa}_{\mathcal {D}}(u)}, \mathcal {D})} \cdot \dfrac{p(\mathcal {D}^*)}{p(\mathcal {D})} \cdot \dfrac{q(\mathcal {D}\,\vert \,\mathcal {D}^*)}{q(\mathcal {D}^*\,\vert \,\mathcal {D})} &{} \quad \\ \text { if } (c), \end{array} \right\} \end{aligned}$$\end{document}

where

$\begin{array}{c}m\left(Z_v,|,Z_{{\text{pa}}_D\left(v\right)},,,D\right)={\left(2\pi \right)}^{-\frac{n}{2}}·\frac{|U_{\prec v\succ }{|}^{\frac{1}{2}}}{|{\tilde{U}}_{\prec v\succ }{|}^{\frac{1}{2}}}·\frac{\Gamma \left(\frac{1}{2},{\tilde{a}}_v^D\right)}{\Gamma \left(\frac{1}{2},a_v^D\right)}·\frac{(\frac{1}{2}{U}_{v|{\text{pa}}_D\left(v\right)}{)}^{\frac{1}{2}{a}_v^D}}{(\frac{1}{2}{\tilde{U}}_{v|{\text{pa}}_D\left(v\right)}{)}^{\frac{1}{2}{\tilde{a}}_v^D}}\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} m\left( \varvec{Z}_v\,\vert \,\varvec{Z}_{\text {pa}_{\mathcal {D}}(v)}, \mathcal {D}\right) = (2\pi )^{-\frac{n}{2}}\cdot \frac{\big |\varvec{U}_{\prec v \succ }\big |^{\frac{1}{2}}}{\big |\widetilde{\varvec{U}}_{\prec v \succ }\big |^{\frac{1}{2}}}\cdot \frac{\Gamma \left( \frac{1}{2}\widetilde{a}_v^{\mathcal {D}}\right) }{\Gamma \left( \frac{1}{2}a_v^{\mathcal {D}}\right) }\cdot \frac{\Big (\frac{1}{2}\varvec{U}_{v\,\vert \,\text {pa}_{\mathcal {D}}(v)}\Big )^{\frac{1}{2}a_v^{\mathcal {D}}}}{\Big (\frac{1}{2}\widetilde{\varvec{U}}_{v\,\vert \,\text {pa}_{\mathcal {D}}(v)}\Big )^{\frac{1}{2}\widetilde{a}_v^{\mathcal {D}}}} \end{aligned}$$\end{document}

is the marginal (i.e., integrated w.r.t.  $L_{\prec v]}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{L}_{\prec v\,]}$$\end{document} and $D_{\mathrm{vv}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{D}_{vv}$$\end{document} ) data distribution relative to the conditional density of $Z_v$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_v$$\end{document} given $z_{{\text{pa}}_D\left(v\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{z}_{\text {pa}_{\mathcal {D}}(v)}$$\end{document} appearing in (4). We refer the reader to Castelletti and Mascaro (Reference Castelletti and Mascaro2022) for further computational details.

4.2. Update of $Z$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{Z}}$$\end{document}

Since the latent data $Z$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Z}$$\end{document} are known only relative to the event $A\left(X\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A(\varvec{X})$$\end{document} (Eq. 7), we also need to sample them from their full conditional distribution. The latter is given by

(17) $\begin{array}{c}p\left(Z\right|D,L,D,X)\propto p(Z\in A\left(X\right)|D,L,D).\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} p(\varvec{Z}\,\vert \,\varvec{D},\varvec{L},\mathcal {D},\varvec{X}) \propto p(\varvec{Z}\in A(\varvec{X}) \,\vert \,\varvec{D},\varvec{L},\mathcal {D}). \end{aligned}$$\end{document}

Also, by exploiting the structure of the extended rank likelihood and the set $A\left(X\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A(\varvec{X})$$\end{document} in (8) and (7), respectively, the conditional distribution of $Z_{i,j}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_{i,j}$$\end{document} corresponds to a $N\left(z_{i,j}\right|-L_{\prec j]}^{\top }{z}_{i,\text{pa}\left(j\right)},D_{\mathrm{jj}})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {N}(z_{i,j}\,\vert \,-\varvec{L}_{ \prec j\,]}^\top \varvec{z}_{i,\text {pa}(j)}, \varvec{D}_{jj})$$\end{document} truncated at $(z_{i,j}^{\left(l\right)},z_{i,j}^{\left(u\right)})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\big (z_{i,j}^{(l)}, z_{i,j}^{(u)}\big )$$\end{document} , where $z_{i,j}^{\left(l\right)}=\text{max}\{z_{k,j}:x_{k,j}<x_{i,j}\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_{i,j}^{(l)}=\text {max}\{z_{k,j}:x_{k,j}<x_{i,j}\}$$\end{document} and $z_{i,j}^{\left(u\right)}=\text{max}\{z_{k,j}:x_{i,j}<x_{k,j}\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_{i,j}^{(u)}=\text {max}\{z_{k,j}:x_{i,j}<x_{k,j} \}$$\end{document} . Therefore, conditionally of the current values of $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{D}$$\end{document} and $L$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{L}$$\end{document} , update of the (n, q) data matrix $Z$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Z}$$\end{document} can be performed by drawing each latent observation $z_{i,j}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_{i,j}$$\end{document} , $j=1,\dots ,q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1,\ldots ,q$$\end{document} , $i=1,\dots ,n$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,\ldots ,n$$\end{document} , from its corresponding truncated-Normal conditional distribution.

4.3. Posterior Inference

The output of our MCMC scheme is a collection of DAGs $\{D^{\left(s\right)}{\}}_{s=1}^S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\big \{\mathcal {D}^{(s)}\big \}_{s=1}^S$$\end{document} and DAG-parameters $\{(D^{\left(s\right)},L^{\left(s\right)}){\}}_{s=1}^S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\big \{\big (\varvec{D}^{(s)}, \varvec{L}^{(s)}\big )\big \}_{s=1}^S$$\end{document} approximately sampled from (13), where S is the number of fixed MCMC iterations. An approximate posterior distribution over the space of DAGs can be obtained by computing, for each $D\in S_q^C$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}\in \mathcal {S}_q^C$$\end{document} ,

(18) $\begin{array}{c}\hat{p}\left(D\right|X)=\frac{1}{S}\sum _{s=1}^{S}1\left(D^{\left(s\right)},=,D\right),\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \widehat{p}(\mathcal {D}\,\vert \,\varvec{X}) =\frac{1}{S} \sum _{s=1}^{S}\mathbbm {1}\left\{ \mathcal {D}^{(s)}=\mathcal {D}\right\} , \end{aligned}$$\end{document}

which approximates the posterior probability of each DAG structure through its MCMC frequency of visits. As a summary of the previous output, we can also recover a (q, q) matrix of (marginal) posterior probabilities of edge inclusion, whose (u, v)-element corresponds to

(19) $\begin{array}{c}\hat{p}(u\to v|X)\equiv {\hat{p}}_{u\to v}=\frac{1}{S}\sum _{s=1}^S{1}_{u\to v}\{D^{\left(s\right)}\},\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \widehat{p} (u\rightarrow v \,\vert \,\varvec{X}) \equiv \widehat{p}_{u\rightarrow v} =\frac{1}{S}\sum _{s=1}^S\mathbbm {1}_{u\rightarrow v}\big \{\mathcal {D}^{(s)}\big \}, \end{aligned}$$\end{document}

an MCMC frequency-based estimated posterior probability of $u\to v$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u \rightarrow v$$\end{document} , where $1_{u\to v}\{D^{\left(s\right)}\}=1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbbm {1}_{u\rightarrow v}\big \{\mathcal {D}^{(s)}\big \}=1$$\end{document} if $D^{\left(s\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}^{(s)}$$\end{document} contains $u\to v$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u\rightarrow v$$\end{document} , 0 otherwise. From the previous quantities, a single graph estimate summarizing the entire MCMC output can be also recovered. Specifically, by fixing a threshold for edge inclusion $k\in (0,1)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k \in (0,1)$$\end{document} , a DAG estimate can be obtained by including those edges whose posterior probability of inclusion in (19) exceeds k. When $k=0.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=0.5$$\end{document} , we name the resulting estimate Median Probability DAG Model (MPM DAG), following the idea proposed by Barbieri and Berger (Reference Barbieri and Berger2004) in a linear regression context. Finally, we can recover a posterior summary of the covariance matrix $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} through the corresponding Bayesian model averaging (Hoeting et al., Reference Hoeting, Madigan, Raftery and Volinsky1999, BMA) estimate

(20) $\begin{array}{c}\hat{\Sigma }=\frac{1}{S}\sum _{s=1}^S{\Sigma }^{\left(s\right)},\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \widehat{\varvec{\Sigma }} =\frac{1}{S} \sum _{s=1}^{S}\varvec{\Sigma }^{(s)}, \end{aligned}$$\end{document}

where for each $s=1,\dots ,S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s=1,\ldots ,S$$\end{document} , ${\Sigma }^{\left(s\right)}=(L^{\left(s\right)}{)}^{-\top }{D}^{\left(s\right)}(L^{\left(s\right)}{)}^{-1}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }^{(s)}=\big (\varvec{L}^{(s)}\big )^{-\top }\varvec{D}^{(s)}\big (\varvec{L}^{(s)}\big )^{-1}$$\end{document} .

5.1. Scenarios

We fix the number of variables $q=20$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q=20$$\end{document} and consider different simulation settings where data are generated by varying the following features:

(a)

The class of DAG structure;

(b)

The type of variables;

(c)

The sample size.

With regard to (a), we consider three different classes of DAG structures:

(a.1)

Free: no structural constraints imposed to the DAG;

(a.2)

Regression: we fix each node $u\in \{1,2,3\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u\in \{1,2,3\}$$\end{document} as a response, so that edges $u\to v$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u\rightarrow v$$\end{document} , for $v\in V\backslash \{1,2,3\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v \in V {\setminus } \{1,2,3\}$$\end{document} are not allowed;

(a.3)

Block: nodes are partitioned into two sets of equal size A and B and only edges within each set or from set A to B are allowed.

Under each scenario defined by (a), $N=40$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=40$$\end{document} DAGs are generated by fixing a probability of edge inclusion $\pi =0.1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi =0.1$$\end{document} for all edges which are allowed in the corresponding class of DAG structures (a). For each given DAG $D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {D}$$\end{document} , we generate the parameters of the underlying Gaussian DAG model in (3) by fixing $D=I_q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{D}= \varvec{I}_q$$\end{document} , while drawing nonzero elements of $L$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{L}$$\end{document} uniformly in the interval $[-1,-0.1]\cup [0.1,1]$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[-1,-0.1] \cup [0.1,1]$$\end{document} . Latent observations collected in the (n, q) data matrix $Z$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Z}$$\end{document} are then generated according to (3) for different sample sizes $n\in \{100,200,500,1000,2000\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n \in \{100,200,500,1000,2000\}$$\end{document} (c). Starting from the latent data matrix $Z$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Z}$$\end{document} , observed data are finally generated as in (5) for different choices of the marginal c.d.f’s $F_1,\dots ,F_q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_1,\ldots ,F_q$$\end{document} (b). Specifically, we consider the following assumptions:

(b.1)

Binary: $X_j\sim \text{Ber}\left({\eta }_j\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_j\sim \text {Ber}(\eta _j)$$\end{document} , $j=1,\dots ,q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1,\ldots ,q$$\end{document} , with ${\eta }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _j$$\end{document} randomly drawn in [0.2, 0.8];

(b.2)

Ordinal: $X_j\sim \text{Binomial}({\theta }_j,5)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_j\sim \text {Binomial}(\theta _j,5)$$\end{document} , $j=1,\dots ,q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1,\ldots ,q$$\end{document} , with ${\theta }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _j$$\end{document} randomly drawn in [0.2, 0.8];

(b.3)

Count: $X_j\sim \text{Poisson}\left({\lambda }_j\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_j\sim \text {Poisson}(\lambda _j)$$\end{document} , $j=1,\dots ,q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1,\ldots ,q$$\end{document} , with ${\lambda }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _j$$\end{document} randomly drawn in [1, 10];

(b.4)

Mixed: $X_1\dots ,X_7$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_1\ldots ,X_7$$\end{document} as in Binary; $X_8,\dots ,X_{15}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_8,\ldots ,X_{15}$$\end{document} as in Ordinal; $X_{16},\dots ,X_{20}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_{16},\ldots ,X_{20}$$\end{document} as in Count.

Each combination of (a), (b) and (c) defines a simulation scenario which consists of $N=40$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=40$$\end{document} (true) DAGs and allied datasets.

5.2. Results

We run our MCMC scheme to approximate the posterior distribution of DAG structures and parameters. In particular, under each scenario among Free, Regression, Block, we limit the DAG space $S_q^C$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {S}_q^C$$\end{document} to those DAGs satisfying the constraints imposed by the corresponding class of DAG structures; see also Sect. 3. The number of iterations was fixed as $S=$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S=$$\end{document}  15,000 after some pilot simulations that were used to assess the convergence of the MCMC chain. We also set $U=g{I}_q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{U}=g\varvec{I}_q$$\end{document} with $g=1/n$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g = 1/n$$\end{document} , $a=q$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = q$$\end{document} in the DAG-Wishart prior (11), while $c=1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c = 1$$\end{document} , $d=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d = 5$$\end{document} in (10), which corresponds to a $\text{Beta}(1,5)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Beta}(1,5)$$\end{document} prior on the probability of edge inclusion $\pi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi $$\end{document} and implies $E\left(\pi \right)=0.1\overline{6}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {E}(\pi )=0.1\bar{6}$$\end{document} . While this specific choice of c and d can favor sparsity in the graphs, some further analyses not reported for brevity showed that results are quite insensitive to the values of the two hyperparameters.

We now evaluate the performance of our method in recovering the true DAG structure. To this end, we first estimate the posterior probabilities of edge inclusion in (19) for each pair of distinct nodes u, v. Given a threshold for edge inclusion $k\in [0,1]$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k\in [0,1]$$\end{document} , we construct a graph estimate by including those edges (u, v) such that ${\hat{p}}_{u\to v}\ge k$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widehat{p}_{u\rightarrow v} \ge k$$\end{document} . We compare the resulting graph with the true DAG by computing the sensitivity (SEN) and specificity (SPE) indexes, defined as

$\begin{array}{c}\text{SEN}=\frac{\text{TP}}{\text{TP}+\text{FN}},\text{SPE}=\frac{\text{TN}}{\text{TN}+\text{FP}},\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\text {SEN}} = \frac{{\text {TP}}}{{\text {TP}}+{\text {FN}}}, \quad {\text {SPE}} = \frac{{\text {TN}}}{{\text {TN}}+{\text {FP}}}, \end{aligned}$$\end{document}

where TP, TN, FP, FN are the numbers of true positives, true negatives, false positives and false negatives, based on the adjacency matrix of the estimated graph. By varying the threshold k within a grid in [0, 1] and computing SEN and SPE for each value of k, we obtain a receiver operating characteristic (ROC) curve where each point corresponds to the values of SEN and $(1-\text{SPE})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(1-{\text {SPE}})$$\end{document} computed for a given threshold k. Moreover, under each scenario, an average (w.r.t. the $N=40$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=40$$\end{document} simulations) curve is constructed as follows. For each threshold k, we compute SEN and $(1-\text{SPE})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(1-{\text {SPE}})$$\end{document} under each of the 40 simulated DAGs and compute the average values of the two indexes. By repeating this procedure for each k we obtain a collection of points which are joined by the average ROC curve. Results for Scenario Mixed (b.4), different sample sizes (c) and types of DAG structures (a) are summarized in Fig. 2. We also proceed similarly to compute the 5th and 95th percentiles and obtain the blue band which is included in each plot of the same figure. As expected, the performance of the method in recovering the true DAG structure improves as the sample size n increases under all settings Free, Regression, Block.

As a single graph estimate summarizing the MCMC output we also consider the median probability (MPM) DAG model $\hat{D}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widehat{\mathcal {D}}$$\end{document} which is obtained by fixing the threshold for edge inclusion as $k=0.5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=0.5$$\end{document} . We compare each DAG estimate $\hat{D}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widehat{\mathcal {D}}$$\end{document} with the corresponding true DAG by measuring the Structural Hamming Distance (SHD). The latter corresponds to the number of edge insertions, deletions or flips needed to transform the estimated DAG into the true DAG; accordingly, lower values of SHD correspond to better performances. Results for each simulation scenario are summarized in the box-plots of Fig. 3, where each plot reports the distribution (across the $N=40$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=40$$\end{document} simulated datasets) of the ratio between SHD and the maximum number of edges in the graphs (SHD/edges) for a combination of (a) and (c) and increasing sample sizes n.

The same behavior observed in Fig. 2 for Scenario Mixed is even more apparent, with SHD decreasing at zero as n increases under all settings. In addition, DAG learning is more difficult in the Binary scenario, where the collected categorical data provide a limited source of information to estimate dependence relations. By converse, structural learning is much easier under the Ordinal and Count scenarios. The performance in the case of mixed data is somewhat intermediate w.r.t. the previous scenarios.

Figure 2 Simulations. Receiver operating characteristic (ROC) curve obtained under varying thresholds for the posterior probabilities of edge inclusion for type of variables Mixed (b.4), sample size $n\in \{100,200,500,1000,2000\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\in \{100,200,500,1000,2000\}$$\end{document} (c) and type of DAG structure Free, Regression, Block (a). Dotted lines represent the (average over the 40 simulated DAGs) ROC curve, while the blue area represents the 5th–95th percentile band.

Figure 3 Simulations. Distribution (across $N=40$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=40$$\end{document} simulated datasets) of the relative structural Hamming distance (SHD/edges) between estimated and true DAG for type of variables Binary, Ordinal, Count, Mixed (b), sample size $n\in \{100,200,500,1000,2000\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\in \{100,200,500,1000,2000\}$$\end{document} (c) and type of DAG structure Free, Regression, Block (a).

We finally summarize in Tables 1, 2, 3 and 4 the sensitivity (SEN) and specificity (SPE) indexes computed under the median probability DAG model estimate. Each table refers to one scenario among Binary, Ordinal, Count, Mixed and reports the average (w.r.t. the $N=40$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=40$$\end{document} simulations) percentage value of the two indexes under settings Free, Regression, Block corresponding to the three different classes of DAG structures. While SPE attains high levels even for small sample sizes, SEN significantly increases as n grows. As a consequence, the improvement in the performance observed in Fig. 3 is mainly due to a reduction in the inclusion of “false negative” edges as the number of available observations grows. In addition, it appears that DAG identification is somewhat easier under Scenarios Regression and Block. Here, structural constraints limiting the DAG space can help identifying dependence relations between variables which otherwise may not be uniquely identifiable because of DAG Markov equivalence (Andersson et al., Reference Andersson, Madigan and Perlman1997).

Table 1 Simulations. Binary data. Specificity (SPE) and sensitivity (SEN) indexes (averaged over the 40 simulations) computed w.r.t. the median probability DAG estimate, for sample size $n\in \{100,200,500,1000,2000\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\in \{100,200,500,1000,2000\}$$\end{document} and under each scenario Free, Regression, Block corresponding to different classes of DAG structures.

Table 2 Simulations. Ordinal data. Specificity (SPE) and sensitivity (SEN) indexes (averaged over the 40 simulations) computed w.r.t. the median probability DAG estimate, for sample size $n\in \{100,200,500,1000,2000\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\in \{100,200,500,1000,2000\}$$\end{document} and under each scenario Free, Regression, Block corresponding to different classes of DAG structures.

Table 3 Simulations. Count data. Specificity (SPE) and sensitivity (SEN) indexes (averaged over the 40 simulations) computed w.r.t. the median probability DAG estimate, for sample size $n\in \{100,200,500,1000,2000\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\in \{100,200,500,1000,2000\}$$\end{document} and under each scenario Free, Regression, Block corresponding to different classes of DAG structures.

Table 4 Simulations. Mixed data. Specificity (SPE) and sensitivity (SEN) indexes (averaged over the 40 simulations) computed w.r.t. the median probability DAG estimate, for sample size $n\in \{100,200,500,1000,2000\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\in \{100,200,500,1000,2000\}$$\end{document} and under each scenario Free, Regression, Block corresponding to different classes of DAG structures.

5.3. Simulation Experiments with Unbalanced Correlation Structure

In the simulation scenarios considered before we randomly drew the nonzero elements of matrix $L$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{L}$$\end{document} , corresponding to regression coefficients in the latent linear SEM, uniformly in the symmetric interval $[-1,-0.1]\cup [0.1,1]$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[-1,-0.1]\cup [0.1,1]$$\end{document} . This implies an expected “balanced” correlation structure between variables, meaning that both positive and negative associations in the generating model will be present, and with same expected proportion. However, in many psychological applications variables are mostly positively correlated each other; see also our results in the data analyses of Sect. 6. An important example is represented by cognitive test scores, whose pattern of positive correlations was already advocated by Spearman (Reference Spearman1904) in his positive manifold theory. To assess the ability of our method to capture such “unbalanced” correlation structure, we now consider a random choice of the nonzero elements of $L$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{L}$$\end{document} in the interval [0.1, 1] and implement the same scenarios introduced in Sect. 5.1. Accordingly, all variables are now positively correlated at the latent level. Results, in terms of relative SHD between true and estimated DAG, are reported in Fig. 4, which summarizes the distribution of SHD across simulated datasets under all scenarios defined in Sect. 5.1. The performance of our method is in line with what we observed in Fig. 3 under a balanced setting, suggesting that our model specification is insensitive to the specific proportion of positive/negative correlations.

Figure 4 Simulations with unbalanced correlation structure. Distribution (across $N=40$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=40$$\end{document} simulated datasets) of the relative structural Hamming distance (SHD/edges) between estimated and true DAG for type of variables Binary, Ordinal, Count, Mixed (b), sample size $n\in \{100,200,500,1000,2000\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\in \{100,200,500,1000,2000\}$$\end{document} (c) and type of DAG structure Free, Regression, Block (a)

6.1. Well-Being Data

The Your Work-Life Balance is a project promoted by the United Nations (UN) and implemented through a public survey available at http://www.authentic-happiness.com/your-life-satisfaction-score. Scope of the survey is to evaluate how people thrive in both professional and personal lives based on several indicators that are related with life satisfaction. Variables considered in the study are classified into five dimensions:

• • Healthy body, features reflecting fitness and healthy habits;

• • Healthy mind, indicating how well subjects embrace positive emotions;

• • Expertise, measuring the ability to grow expertise and achieve something unique;

• • Connection, assessing the strength of social relationships;

• • Meaning, evaluating compassion, generosity and happiness.

The survey supports the UN Sustainable Development Goals (https://sdgs.un.org) and aims at providing insights on the determinants of human well-being. Accordingly, some questions of interest are the following:

• • “What are the strongest correlations between the various dimensions?”

• • “What are the best predictors of a balanced life?”

The complete dataset is publicly available at https://www.kaggle.com/datasets/ydalat/lifestyle-and-wellbeing-data. It includes observations collected across years 2015–2020 of 20 ordinal variables (with levels ranging in 1–5 or 1–10) each measuring closeness of a subject w.r.t. to one perceived dimension, besides gender (binary) and age (ordinal with four classes). We include in our analysis the $n=459$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n = 459$$\end{document} observations available for year 2020. We consider variable stress as the response and accordingly allow edges from each of the remaining variables to the response only. In addition, age and gender are considered as objective features and accordingly cannot have incoming edges. We do not impose further constraints among the remaining variables in terms of edge directions that are known in advance.

Following the scope of the original project, we are interested in understanding how the various dimensions relate each other and what are the direct/indirect determinants of the perceived level of stress. To this end, we implement our MCMC algorithm for a number of iterations $S=$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S=$$\end{document}  20,000, after an initial burnin period of 5000 runs that are discarded from posterior analysis. Diagnostic tools based on multiple chains and graphical inspections of the behavior across iterations of sampled parameters were also adopted to assess the convergence of the MCMC; see also Appendix for details.

We use the MCMC output to provide an estimate of the posterior probabilities of edge inclusion as well as a BMA estimate of the correlation matrix between (latent) variables. Results are reported in the two heat maps of Fig. 5. The upper map, which collects the estimated posterior probabilities of edge inclusion clearly suggests a sparse structure in the underlying network, with only a few edges whose posterior probability exceeds 0.5. This also emerges from the graph estimate reported in Fig. 6, which corresponds to the CPDAG representing the equivalence class of the MPM DAG estimate.

Figure 5 Well-being data. Upper panel: Heat map with estimated posterior probabilities of edge inclusion $\hat{p}(u\to v|X)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widehat{p} (u\rightarrow v \,\vert \,\varvec{X})$$\end{document} for each edge (u, v). Lower panel: Estimated correlation matrix.

Figure 6 Well-being data. Estimated CPDAG.

The estimated graph reveals that two variables directly affect the perceived level of stress, namely shouting (“How often do you shout or sulk at somebody?”) and to do list (“how well do you complete your weekly to-do list?”). Moreover, variable stress is conditionally independent from flow (“How many hours you experience flow, i.e., you fell fully immersed in performing an activity?”) given to do list. Also, flow and to do list are positively correlated, which suggests that people performing better in their activities also follow through with many more of their weekly goals. This in turn has a direct impact on the perceived level of stress. It also appears that shouting is positively correlated with the level of stress, while there is negative correlation with to do list; accordingly better results in completing to-do-lists, imply a reduction in the stress level perceived by individuals.

6.2. Student Mental Health Data

We consider a dataset from a cross-sectional study conducted on 886 medical students from the University of Lausanne (Switzerland) and presented by Carrard et al. (Reference Carrard, Bourquin, Berney, Schlegel, Gaume, Bart, Preisig, Mast and Berney2022). The data collection was based on a longitudinal open cohort study survey. In particular, all medical students enrolled in the university were eligible for participation and received a link to an online questionnaire; we refer to the original paper for details about the survey design and the collected data. Target of the study was to provide insights on students’ well-being, in order to implement policies aimed at improving their academic-life satisfaction and conditions. A number of dimensions related to empathy are measured through self-reported questionnaires based on the Questionnaire of Cognitive and Affective Empathy (QCAE) and the Jefferson Scale of Physician Empathy (JSPE); related variables are the QCAE affective empathy score (qcae aff), the QCAE cognitive empathy score (qcae cog) and the JSPE total empathy score (jspe). Burnout is a state of emotional, physical, and mental exhaustion which is caused by excessive exposure to stress. The burnout dimension is measured through the Maslach Burnout Inventory-Student Survey (MBI-SS); the latter is based on 15 items and provides three scores evaluating the following dimensions: emotional exhaustion (mbi ex), cynicism (mbi cy), and academic efficacy (mbi ea). In addition, students’ anxiety and depression is measured through the Center for Epidemiologic Studies Depression (CESD) score and the State-Trait Anxiety Inventory (STAI) score, both based on a questionnaire with self-report items on Likert scales (cesd and stai, respectively). Finally, the dataset contains information on demographic factors such as age and gender, besides variables measuring job satisfaction (job), partnership status (part) and self-reported health status (health), represented by an ordinal variable with 5 categories corresponding to increasing levels of perceived health satisfaction. We again refer to Carrard et al. (Reference Carrard, Bourquin, Berney, Schlegel, Gaume, Bart, Preisig, Mast and Berney2022) for a detailed description of the complete dataset, which is publicly available at https://zenodo.org/record/5702895. We emphasize that the structure of the analyzed dataset is quite heterogeneous, as it collects binary, ordinal (with different ranges of levels) as well as continuous measurements simultaneously.

One specific aim of the original study is to identify how variables included in the survey, in particular depressive symptoms, anxiety, and burnout, are related to empathy and mental health. In our analysis, we consider age, sex, part and job as exogenous variables, while we regard health as a response of interest.

Our MCMC algorithm is implemented for S = 20000 iterations, after a burnin period of 5000 runs that we adopt to assess the convergence of the chain. We use the MCMC output to provide an estimate of the (marginal) posterior probability of inclusion (PPI) for each possible directed edge in the DAG space; see Eq. (19). The resulting heat-map with the collection of estimated PPIs (Fig. 7, upper plot) suggests the existence of a few strong dependence relations among variables that correspond to directed links and paths in the graphs visited by the MCMC chain; this also appears from the MPM CPDAG estimate reported in Fig. 8 which contains a moderate number of edges. Furthermore, to investigate how variables correlate each other, we provide a BMA estimate of the correlation matrix between (latent) variables (Fig. 7, lower plot).

Notably, most variables belonging to the burnout dimension as well as to the anxiety/depression sphere are positively correlated with the exception of mbi ea, namely the score measuring academic efficacy, which as expected is negatively associated with the other variables, and in particular with the perceived level of anxiety (stai). mbi ea is also positively influenced by variable study (hours of study per week), which in turn has a positive, although less marked, effect on student’s anxiety. Importantly however, students’ depression, here quantified by the cesd score, has a negative effect on health, as it appears from the graph estimate of Fig. 8 and the correlation matrix in Fig. 7. Also, age has a positive impact on variable jspe which summarizes the empathy dimension, implying that medical students develop throughout their academic life a growing ability in understanding and sharing feelings that are experienced by others. Finally, an interesting set of structural dependencies is represented by the directed path study $\to$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rightarrow $$\end{document} stai $\to$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rightarrow $$\end{document} cesd $\to$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rightarrow $$\end{document} health. The latter structure suggests that students more involved in studying activities may incur in higher levels of anxiety and depression, which consequently affects personal health conditions.

Figure 7 Student mental health data. Upper panel: heat map with estimated posterior probabilities of edge inclusion $\hat{p}(u\to v|X)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widehat{p} (u\rightarrow v \,\vert \,\varvec{X})$$\end{document} for each edge (u, v). Lower panel: estimated correlation matrix

Figure 8 Student mental health data. Estimated CPDAG.