1

1

Generative Models, Interpretable AI, Identifiability, Latent Variables, Bayesian, Rank likelihood.

1 Introduction↩︎

The widespread popularity of AI tools are due to the impressive predictive capabilities of underlying deep generative models (DGMs). Many DGMs rely on black-box, deep learning architectures governed by extremely large numbers of parameters relative to the training sample size. Such DGMs are fundamentally unidentified, and the estimated parameters lack interpretation of the substantive relationships in the data. From a statistical perspective, this is problematic; indeed, statistical inference is a cornerstone for scientific discovery.

In this work, we contribute to a rapidly developing literature focusing on developing statistically identified and simultaneously interpretable DGMs. Specifically, we introduce the Deep Discrete Encoder (DDE) Copula, a generative model with multiple latent layers for multivariate data with arbitrary marginal distributions. The DDE copula features a hierarchical, pyramid-like directed graph where the latent variables in each layer are binary and conditionally independent given parents (See Figure [fig:DDE95right]). We then nest this structure – capable of capturing complex dependencies among diverse data types – within a copula to accommodate variables with arbitrary scales and marginal features.

We first provide identifiability conditions for the DDE copula and show that when these are satisfied, the DDE copula is strictly identifiable. This is vitally important for real data analyses, as the parameters governing each layer of the network provide directly interpretable and increasingly fine grained insight into multivariate relationships. As we show in a real data analysis (Section 6), these novel insights demonstrate how the DDE copula is a powerful exploratory model.

From a computational standpoint, we introduce several algorithmic novelties that enable scalable Bayesian estimation of the DDE Copula with unknown layer widths. First, we decouple marginal modeling from DDE posterior inference using an extended rank likelihood [1]. The extended rank likelihood leverages multivariate ranks – not the observed values – among the data to estimate the DDE parameters. As such, the univariate marginals are nuisance parameters, which greatly facilitates estimation. Due to this attractive property, the extended rank likelihood has been used for Bayesian estimation of Gaussian copulas and Gaussian mixture copulas [2]–[4]. We show how to leverage the rank likelihood for deep copula models.

We also clarify the large-sample target of rank likelihood inference for the DDE copula: for continuous margins, we prove quotient-space posterior consistency under the exact rank likelihood; for discrete or mixed margins, where ties make the extended rank likelihood [1] a generalized likelihood rather than an exact likelihood, we state posterior concentration under an explicit extended-rank contrast condition (Section 3.3). We also confirm excellent finite sample behavior in simulation (Section 5).

Next, we construct an efficient expectation-maximization algorithm for maximum a posteriori (MAP) estimation of the DDE parameters under the rank likelihood (Section 4). The algorithm conveniently accommodates the rank likelihood, as an intractable E-step is replaced by a simulation step for all latent variables [5], [6]. The stochastic E-step also further simplifies the M-step, enabling reasonable computation times for large networks and data sets. We also provide novel initialization strategies that facilitate rapid convergence (Section 4.3).

Finally, we demonstrate how state-of-the-art Bayesian priors for rank-selection in Gaussian factor models [7] may be extended to automate layer-specific width selection for the DDE copula architecture. Our formulation enables information sharing between layers and provides structured sparsity that improves interpretability. The ability to estimate layer-specific widths under the DDE copula provides analysts with an additional tool to discover hierarchical latent structures in their data.

The DDE copula is directly motivated by recent work on Bayesian Pyramids [8] and Deep Discrete Encoders (Figure [fig:DDE95left]; [9]), two deep and statistically identifiable generative models for multivariate data. However, the distinguishing features of the DDE copula are motivated by two primary limitations of these preceding models. First, the Deep Discrete Encoder features parametric models for each study variable at the observed data layer of the architecture, while the Bayesian Pyramid is only compatible with multivariate categorical data. The DDE copula is broadly compatible with mixed discrete and continuous variables without any marginal distributional assumptions. Second, both models require pre-specification of at least some of the unknown layer widths. We instead treat the dimensions of the DDE network as unknown model components and infer each layer-specific width within our MAP estimation algorithm.

More broadly, the DDE copula sits at the intersection of identifiable deep generative models and semiparametric copula models. Classical deep latent-variable architectures, including deep belief networks [10], deep Boltzmann machines [11], and variational autoencoders [12], provide flexible representation learning but are typically not identifiable without additional structural assumptions [13] and are difficult to interpret. In contrast, the recent DDE and Bayesian Pyramid frameworks use sparse discrete latent layers to obtain interpretable hierarchical representations in parametric models [8], [9]. On the other hand, semiparametric Gaussian copula models accommodate mixed margins through rank likelihoods [1], [2], [4], but typically rely on shallow continuous latent factors. The DDE copula combines these strengths by retaining a statistically identifiable deep discrete latent hierarchy while allowing arbitrary marginal distributions, and further addresses the crucial question of automatically estimating the latent layer widths in a statistically principled manner.

This paper is organized as follows: In Section 2, we review the Deep Discrete Encoder of [9]. In Section 3, we introduce the DDE copula, demonstrate how to leverage the rank likelihood for estimation, and establish key theoretical properties including identifiability and rank-likelihood large-sample guarantees. Section 4 includes details on Bayesian specification of the DDE copula, as well as estimation and initialization algorithms. We provide an extensive simulation study in Section 5, and conclude with an analysis of survey responses on human personality in Section 6. Concluding thoughts and ideas for future research are provided in Section 7. The Supplementary Material contains all proofs of the theoretical results and additional numerical results. Codes implementing the proposed approach are available https://github.com/jfeldman396/DDE-Copula/tree/main.

2 Deep Discrete Encoders↩︎

We first describe the Deep Discrete Encoder [9], which provides the basis for our proposed innovations. A \(D\)-layer DDE forms a directed graphical model with a pyramid-like shape. At the shallowest and widest layer (\(d=0\)) are the \(J\)-dimensional observed data \(\boldsymbol{Y} = (Y_{1}, \dots, Y_{J})\). Each subsequent layer \(d \in \{1, \dots, D\}\) is governed by \(K^{(d)}\) binary latent variables \(\boldsymbol{A}^{(d)} = (A^{(d)}_{1}, \dots, A^{(d)}_{K^{(d)}}) \in \{0,1\}^{K^{(d)}}\) with \(K^{D}< \dots < K^{1} < J\). For identifiability, the dimension of each layer \(d\) is constrained to be \(\lfloor K^{(d-1)}/3\rfloor\), with \(K^{(1)} = \lfloor J /3 \rfloor\) [8], [9]. In Figure [fig:DDE95left], we provide an illustration of this architecture.

To describe the data-generating model under a \(D\)-dimensional DDE, the top layer latent variables are assumed to be independent Bernoulli random variables with parameter \(\boldsymbol{\pi} = (\pi_{1}, \dots, \pi_{K^{(D)}})\). Then, latent variables in each subsequent middle layer are conditionally independent given the binary variables in the previous layer, modeled with additive logistic regressions. That is, conditional on any \(d\)-th level binary representation \(\boldsymbol{\alpha}^{(d)} \in \{0,1\}^{K^{(d)}}\), the joint probability for any \(\boldsymbol{\alpha}^{(d-1)} \in \{0,1\}^{K^{(d-1)}}\) factors as \[\begin{align} \tag{1} p(\boldsymbol{A}^{(d-1)} = \boldsymbol{\alpha}^{(d-1)} \mid \boldsymbol{\alpha}^{(d)}) &= \prod_{k=1}^{K^{(d-1)}} p(A^{(d-1)}_{k} = \alpha^{(d-1)}_{k}\mid \boldsymbol{\alpha}^{(d)})\\ &= \prod_{k=1}^{K^{(d-1)}} g_{\text{logistic}}(\beta^{(d)}_{k,0}+ \sum_{\ell =1}^{K^{(d)}} \beta^{(d)}_{k,l} \alpha^{(d)}_{l}\tag{2} ) \end{align}\] Here, \(g_{\text{logistic}}\) is the logistic link function \(g(x) = 1/1+e^{-x}\). Modifications are available for non-linear and non-additive relationships between the layers [14], [15]. This formulation implies a set of layer-specific weight matrices \(\boldsymbol{B}^{(d)} = \{\beta_{jk}\}, \; \forall d \in \{1, \dots, D\}\), where the \(d\)-th layer-specific dimension is \(K^{(d-1)} \times K^{(d)}\).

At the observed data layer, the likelihood is \[\label{outcome} p(\boldsymbol{Y} \mid \boldsymbol{\alpha}^{(1)}) = \prod_{j=1}^{J} p_{j}\{g_{j}(Y_{j};\beta^{(1)}_{j,0} + \sum_{k =1}^{K^{(1)}} \beta^{(1)}_{j,k} \alpha^{(1)}_{k}, \gamma_{j})\}\tag{3}\] where \(p_{j}\) is the probability density/mass function of a user-specified parametric family corresponding to the scale of \(Y_{j}\) (e.g., count, continuous, binary), \(\gamma_{j}\) is a dispersion parameter when required, and \(g_{j}\) is the appropriate link function. For example, for count-valued \(Y_{j}\), the user may specify \(p_{j}\) to be the Poisson probability mass function with link function \(g_{j}(x) = e^{x}\) and no dispersion parameter.

From the construction 1 –3 , the joint probability of the observed and latent variables is \[\label{eq::joint} p(\boldsymbol{Y}, \{\boldsymbol{A}^{(d)}\}_{d \in [D]}) = p(\boldsymbol{Y} \mid \boldsymbol{A}^{(1)})\prod_{d=2}^{D}p(\boldsymbol{A}^{(d-1)} \mid \boldsymbol{A}^{(d)})p(\boldsymbol{A}^{(D)})\tag{4}\] This distribution is parameterized by the set \(\boldsymbol{\theta} = (\boldsymbol{B}, \boldsymbol{\pi}, \boldsymbol{\gamma})\), where \(\boldsymbol{B}= \{\boldsymbol{B}^{(d)}\}_{d=1}^{D}\) are the layer-specific weight matrices and \(\boldsymbol{\gamma} = (\gamma_{1},\dots, \gamma_{J})\). Marginalization over the latent variables in 4 reveals that the distribution of \(\boldsymbol{Y}\) is a mixture with \(2^{\sum_{d=1}^{D} K^{(d)}}\) components. This yields a highly expressive model capable of capturing hierarchical dependence in increasingly complex and high-dimensional data, including text, images, and multi-modal data sources [9]. A primary statistical advantage lies in the observation that, despite the saturated marginal model, the DDE requires estimation of only \(J + K^{(D)} + \sum_{d=2}^{D} K^{(d)} \times K^{(d-1)}\) parameters, which comprises the exponential family dispersion parameters \(\boldsymbol{\gamma}\), the top-level Bernoulli probabilities \(\boldsymbol{\pi}\), and the layer-specific weight matrices \(\boldsymbol{B}\), respectively.

In practice, the specification of a DDE requires the user to fix the widths of the layers given the dimension of \(\boldsymbol{Y}\), along with the parametric family for each marginal model \(p_{j}\) [9]. For the former, the authors proposed several ad-hoc approaches based on score-based metrics and layer-specific spectral estimators. For the latter, depending on the structure of the data, the model is amenable to specification of different exponential families, including conditionally Normal, Poisson, and Bernoulli distributions. With these choices, the authors developed a stochastic approximation expectation-maximization algorithm [16] to infer \(\boldsymbol{\theta}\). Next, we propose methods that enable learning of the DDE layer widths, in addition to a relaxation of the parametric marginal assumption.

3 The DDE Copula and Its Theoretical Properties↩︎

The reliance of the DDE on parametric marginal models for each \(Y_j\) limits its applicability. For example, in survey data with ordinal responses (Section 6), each \(Y_j\) has bounded, discrete support, making it difficult to specify a suitable parametric model \(p_j\) in 3 . As demonstrated in our simulations (Section 5), misspecification of marginal features such as boundedness, skewness, or multimodality can degrade estimation of the latent dependence structure through \(\boldsymbol{B}\).

Copulas provide a natural remedy. By Sklar’s theorem [17], any joint distribution \(\boldsymbol{F}\) can be decomposed into marginal distribution functions \(\{F_{Y_j}\}_{j=1}^J\) and a copula \(\boldsymbol{C}\) parameterized by \(\boldsymbol{\Theta}\), \(\boldsymbol{F}(y_1,\dots,y_J) = \boldsymbol{C}_{\boldsymbol{\Theta}}\{F_{Y_1}(y_1),\dots,F_{Y_J}(y_J)\}\). We construct a DDE copula by combining the hierarchical dependence structure of a \(D\)-layer DDE with arbitrary marginals. This yields a flexible model capable of handling data with diverse marginal features, while preserving interpretable dependence through the DDE parameters.

We introduce a latent Gaussian layer \(\boldsymbol{Z}\) between \(\boldsymbol{A}^{(1)}\) and \(\boldsymbol{Y}\). Fixing the middle-layer structure 1 –2 , we define \[\begin{align} p(\boldsymbol{Z} \mid \boldsymbol{\alpha}^{(1)}) &= \prod_{j=1}^{J} p_j\{g_j(Z_j; \beta^{(1)}_{j,0} + \sum_{k=1}^{K^{(1)}} \beta^{(1)}_{j,k} \alpha_k^{(1)}, \gamma_j)\}, \tag{5} \\ Y_j &= F_{Y_j}^{-1}\{F_{Z_j}(Z_j)\}, \quad j=1,\dots,J. \tag{6} \end{align}\] The transformation 6 induces a copula in which dependence is governed by the DDE on \(\boldsymbol{Z}\), while the marginals \(\{F_{Y_j}\}\) remain unrestricted. The DDE copula architecture is visualized in Figure [fig:DDE95right]. For computational convenience, we take \(p_j\) to be Gaussian and \(g_j\) the identity, \(p(\boldsymbol{Z} \mid \boldsymbol{\alpha}^{(1)}) = \prod_{j=1}^{J} \phi\!\left(Z_j; \beta^{(1)}_{j,0} + \sum_{k=1}^{K^{(1)}} \beta_{j,k} \alpha_k^{(1)}, \gamma_j\right)\), which yields substantial computational gains while retaining flexible marginal modeling.

Marginally, \(\boldsymbol{Z}\) follows a Gaussian mixture induced by \(\boldsymbol{A}\). The DDE copula thus resembles a Gaussian mixture copula [4], but achieves substantially greater modeling capacity and statistical parsimony. For example, with \(J=50\), \(K^{(1)}=16\), and \(K^{(2)}=5\), the implied mixture has \(2^{K^{(1)} + K^{(2)}}=2^{21}\) components, yet requires only \(50 + 5 + 50\times16 + 16\times5\) parameters. This feature enables the DDE copula to be estimated on modestly-sized data sets.

3.1 Identifiability Theory of the Canonical DDE Copula↩︎

We next establish identifiability of the dependence parameters in the DDE copula. The result is stated for the copula law induced by the latent Gaussian vector \(Z\), because the univariate marginal distributions of the observed variables \(Y_1,\ldots,Y_J\) are left unrestricted and are therefore nuisance parameters.

Recall that the DDE copula introduces a latent Gaussian layer \(\boldsymbol{Z}=(Z_1,\ldots,Z_J)\) between the shallowest binary latent layer \(\boldsymbol{A}^{(1)}\) and the observed data \(\boldsymbol{Y}\). For \(\boldsymbol{\alpha}\in\{0,1\}^{K^{(1)}}\), write \[Z_j\mid \boldsymbol{A}^{(1)}=\boldsymbol{\alpha} \sim N(\mu_{j\boldsymbol{\alpha}},\gamma_j), \qquad \mu_{j\boldsymbol{\alpha}} = \beta^{(1)}_{j0} + \sum_{k=1}^{K^{(1)}}\beta^{(1)}_{jk}\alpha_k, \qquad \gamma_j>0 .\] The observed variable is generated by the monotone marginal transformation \(Y_j = F_{Y_j}^{-1}\{F_{Z_j}(Z_j)\}\), \(j=1,\ldots,J,\) where \(F_{Y_j}\) is unrestricted and \(F_{Z_j}\) is the marginal distribution function of \(Z_j\) under the DDE model. Let \(U_j = F_{Z_j}(Z_j)\), \(j=1,\ldots,J.\) Then \(\boldsymbol{U}=(U_1,\ldots,U_J)\) has uniform margins, and its joint distribution is the copula induced by the latent DDE model. We denote this copula by \(C_{\boldsymbol{\Theta}}\), where \[\boldsymbol{\Theta} = \Bigl( \boldsymbol{\pi},~\{\boldsymbol{B}^{(d)}\}_{d=1}^D,~\boldsymbol{\gamma},~\{\boldsymbol{G}^{(d)}\}_{d=1}^D \Bigr)\] collects the DDE parameters. Here \(\boldsymbol{\pi}=(\pi_1,\ldots,\pi_{K^{(D)}})\) are the independent Bernoulli probabilities in the top latent layer, \(\boldsymbol{B}^{(d)}\) is the coefficient matrix between layers \(d\) and \(d-1\), \(\boldsymbol{\gamma}=(\gamma_1,\ldots,\gamma_J)\), and \(g^{(d)}_{\ell k}=1\{\beta^{(d)}_{\ell k}\ne 0\}\). defines the graphical matrix \(\boldsymbol{G}^{(d)}\). We use the convention \(K^{(0)}=J\), so that rows of \(\boldsymbol{B}^{(1)}\) correspond to \(Z_1,\ldots,Z_J\), whereas rows of \(\boldsymbol{B}^{(d)}\), \(d\ge2\), correspond to the variables in \(\boldsymbol{A}^{(d-1)}\).

*Canonical parameterization. The copula is invariant under strictly increasing transformations of each coordinate \(Z_j\). Within the Gaussian first layer, positive affine transformations \(Z_j\mapsto a_j Z_j+b_j\), \(a_j>0\), preserve the copula while changing \((\beta^{(1)}_{j0},\beta^{(1)}_{j1},\ldots,\beta^{(1)}_{jK^{(1)}},\gamma_j)\). We therefore impose the canonical normalization \[E_{\boldsymbol{\Theta}}(Z_j)=0, \qquad \operatorname{Var}_{\boldsymbol{\Theta}}(Z_j)=1, \qquad j=1,\ldots,J . \label{eq:canonical-normalization}\tag{7}\] Equivalently, if \(\eta^{(1)}_{\boldsymbol{\alpha}} = P_{\boldsymbol{\Theta}}(\boldsymbol{A}^{(1)}={\boldsymbol{\alpha}}),\) \({\boldsymbol{\alpha}}\in\{0,1\}^{K^{(1)}},\) then \[\sum_{{\boldsymbol{\alpha}}}\eta^{(1)}_{\boldsymbol{\alpha}} \mu_{j{\boldsymbol{\alpha}}}=0, \qquad \gamma_j+ \sum_{{\boldsymbol{\alpha}}}\eta^{(1)}_{\boldsymbol{\alpha}} \mu_{j{\boldsymbol{\alpha}}}^2 = 1, \qquad j=1,\ldots,J . \label{eq:canonical-moment-equations}\tag{8}\] The first equation is \(E(Z_j)=0\); the second follows from \(E(Z_j^2)=E[\gamma_j+\mu_{j,A^{(1)}}^2]=1\), which equals \(\operatorname{Var}(Z_j)=1\) precisely because the mean is zero.

As in latent-variable models generally, the latent coordinates within each layer are only identifiable up to label permutations. We say that two parameters \(\boldsymbol{\Theta}\) and \(\widetilde{\boldsymbol{\Theta}}\) are equivalent, written \(\boldsymbol{\Theta}\sim_K\widetilde{\boldsymbol{\Theta}}\), if there exist permutations \(\sigma^{(d)}\in S_{K^{(d)}}\), \(d=1,\ldots,D\), such that all latent-layer coordinates, columns and rows of the corresponding coefficient matrices, graphical matrices, and top-layer probabilities are transformed according to these permutations, while the observed coordinates \(j=1,\ldots,J\) are fixed.

3.2 Extended Rank Likelihoods↩︎

Though the DDE copula is identified with arbitrary marginals, estimation remains challenging in part because the marginal distributions \(F_{Y_{j}}\) enter the likelihood. Thus, the modeler needs to specify \(F_{Y_{j}}\) for each variable, which, in high dimensions, is arduous. In addition, developing scalable estimation algorithms for both sets of parameters – the DDE dependence structure and the marginal distribution functions – is increasingly complicated.

We instead separate dependence from marginal modeling using the extended rank likelihood [1], which provides an approximation to the full data likelihood under the DDE copula with several advantages. The RL renders \(\{F_{Y_j}\}\) as nuisance parameters for estimation of the DDE copula. As a result, this avoids the need to specify each marginal model which improves scalability and extends the DDE copula to accommodate mixed data types [3], [4]. Section 3.3 states the corresponding large-sample guarantees for the exact continuous-margin rank likelihood and for the extended rank likelihood viewed as a generalized likelihood with ties. Under 6 , \(Y_j\) and \(Z_j\) are monotone transforms, implying \(y_{ij} < y_{\ell j} \;\Rightarrow\; z_{ij} < z_{\ell j}.\) Define the rank-consistent set \[\label{eq:rank} \mathcal{R}(\boldsymbol{Y}) = \{\boldsymbol{Z} \in \mathbb{R}^{n \times J}: y_{ij} < y_{\ell j} \Rightarrow z_{ij} < z_{\ell j}, i = 1, \dots, n, \;j = 1,\dots, J, \}.\tag{11}\] The likelihood for the DDE copula may be written \[\begin{align} &p(\boldsymbol{Y}, \boldsymbol{A}\mid \boldsymbol{\theta}, \{F_{Y_{j}}\}_{j=1}^{J}) = p(\boldsymbol{Y}, \boldsymbol{A}, \boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y}) \mid \boldsymbol{\theta}, \{F_{Y_{j}}\}_{j=1}^{J}) \tag{12}\\ &= \underbrace{p( \boldsymbol{A},\boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y}) \mid \boldsymbol{\theta})}_{\text{extended rank likelihood}}p(\boldsymbol{Y} \mid \boldsymbol{A},\boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y}),\boldsymbol{\theta}, \{F_{Y_{j}}\}_{j=1}^{J})\tag{13}. \end{align}\]

The equivalence in 12 is by construction: observing the data \(\boldsymbol{Y}\) implies that the corresponding latent data \(\boldsymbol{Z}\) must fall into the rank-consistent set 11 . The first term on the right-hand side of 13 is the extended rank likelihood, which depends only on \(\boldsymbol{\theta}\) and not on \(\{F_{Y_j}\}\). Thus, inference for \(\boldsymbol{\theta}\) under the extended rank likelihood \(p\{\boldsymbol{A}, \boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y}) \mid \boldsymbol{\theta}\}\) can proceed without specifying marginals.

Intuitively, \(\boldsymbol{Z}\) captures much of the information about the DDE dependence structure. Prior work shows that the RL is sufficient for estimating Gaussian copula parameters [1] and motivates Bayesian rank-likelihood inference for mixed data [2]. More recently, [4] studied posterior consistency for rank-likelihood inference in Gaussian mixture copulas.

3.3 Rank-likelihood Identifiability and Posterior Consistency↩︎

Section 3.1 establishes identifiability of the canonical DDE-copula parameter from the population copula law \(C_{\boldsymbol{\Theta}}\). We now connect this result to rank-likelihood inference. The key distinction is between continuous margins and margins with ties. With continuous margins, the rank likelihood is the exact likelihood of the observed rank data, so a Doob martingale argument can be applied on the quotient parameter space. With discrete or mixed margins, the extended rank likelihood [1] remains margin-free, but ties make it a generalized likelihood rather than the exact likelihood of the observed weak-rank/tie pattern. We therefore treat the continuous-margin and mixed-margin cases separately.

Let \(\boldsymbol{U}_i=(U_{i1},\ldots,U_{iJ})\), \(i=1,\ldots,n\), be i.i.d. from the DDE copula \(C_{\boldsymbol{\Theta}}\), and let \(Y_{ij}=F_{Y_j}^{-1}(U_{ij})\) for arbitrary nondegenerate margins \(F=(F_{Y_1},\ldots,F_{Y_J})\). For a realized data array \(\boldsymbol{Y}^{(n)}\), define the rank-consistent set on the copula scale by \[\mathcal{D}_n(\boldsymbol{Y}^{(n)}) = \Bigl\{ \boldsymbol{u}\in(0,1)^{n\times J}: y_{ij}<y_{\ell j}\Rightarrow u_{ij}<u_{\ell j}, \quad i,\ell=1,\ldots,n, \;j=1,\ldots,J \Bigr\}. \label{eq:Dn-rank-set}\tag{14}\] The marginal extended rank likelihood is \[L_n^R(\boldsymbol{\Theta};\boldsymbol{Y}^{(n)}) = P_{\boldsymbol{\Theta}}\{\boldsymbol{U}^{(n)}\in \mathcal{D}_n(\boldsymbol{Y}^{(n)})\}. \label{eq:extended-rank-likelihood}\tag{15}\] Because \(F_{Z_j}\) is continuous and strictly increasing under Assumption 1, this is equivalent to the latent-Gaussian rank likelihood based on \(\boldsymbol{Z}^{(n)}\in\mathcal{R}(\boldsymbol{Y}^{(n)})\). It depends on \(\boldsymbol{\Theta}\) only through \(C_{\boldsymbol{\Theta}}\), and not on the unknown margins.

Let \(\mathcal{K}=(K^{(1)},\ldots,K^{(D)})\) denote the active layer widths, and let \(\mathcal{T}_K\) be the canonical parameter space satisfying the normalization 7 . Let \(\Psi_K=\mathcal{T}_K/\!\sim_K\) be the quotient space under within-layer permutation equivalence. We write \(\psi=[\boldsymbol{\Theta}]\in\Psi_K\), and equip \(\Psi_K\) with a metric \(d_K\) generating the quotient topology. Let \(\bar\Pi\) denote the prior induced on \(\Psi_K\). Since \(L_n^R(\boldsymbol{\Theta};\boldsymbol{Y}^{(n)})\) is invariant under latent-label permutations, write \(L_n^R(\psi;\boldsymbol{Y}^{(n)})\) for its common value on \(\psi\). For \(A\in\mathcal{B}(\Psi_K)\), define the quotient-space rank-likelihood posterior by \[\bar\Pi_n^R(A\mid \boldsymbol{Y}^{(n)}) = \frac{ \int_A L_n^R(\psi;\boldsymbol{Y}^{(n)})\,\bar\Pi(d\psi) }{ \int_{\Psi_K}L_n^R(\psi;\boldsymbol{Y}^{(n)})\,\bar\Pi(d\psi) }. \label{eq:quotient-rank-posterior}\tag{16}\] The denominator is positive because \(\mathcal{D}_n(\boldsymbol{Y}^{(n)})\) is nonempty and canonical DDE copulas in \(\mathcal{T}_K\) have positive densities on \((0,1)^J\).

Let \(\mathscr R_n = \sigma\!\bigl( 1\{Y_{ij}<Y_{\ell j}\}, 1\{Y_{ij}=Y_{\ell j}\}: i,\ell=1,\ldots,n, \;j=1,\ldots,J \bigr)\) be the sigma-field generated by the coordinatewise weak ranks and ties, and define \(\mathscr R_\infty=\sigma\!\left(\bigcup_{n\ge1}\mathscr R_n\right).\) When all margins are continuous, ties occur with probability zero, and \(\mathscr R_n\) coincides almost surely with the ordinary coordinatewise rank sigma-field.

The compactness condition is a theoretical truncation to a large nondegenerate subset of the canonical parameter space; the computational prior in Section 4 is used as a practical approximation to this population target.

3.3.1 Continuous Margins↩︎

With continuous margins, ties occur with probability zero, \(\mathcal{D}_n(\boldsymbol{Y}^{(n)})\) is the observed-rank cell, and \(L_n^R\) is the exact rank-data likelihood.

*Interpretation. For continuous margins, Theorem 2 has the status of an ordinary posterior consistency result for an exact likelihood. Because ties occur with probability zero, \(\mathcal{D}_n(\boldsymbol{Y}^{(n)})\) is the observed-rank cell and \(L_n^R\) is the likelihood of the rank data. Thus the posterior is a regular conditional posterior on the quotient space, and the Doob argument gives consistency for \(\bar\Pi\)-almost every true equivalence class once the infinite ranks identify \(\psi_0\).

3.3.2 Discrete or Mixed Margins↩︎

When ties occur, the extended rank likelihood in [1] remains margin-free, but it is no longer an exact likelihood for the observed weak-rank/tie pattern. Tied observations impose no strict ordering constraints in \(\mathcal{D}_n(\boldsymbol{Y}^{(n)})\), so the extended-rank event need not have probability equal to the probability of the observed tie pattern. Thus, with ties, the extended rank likelihood should be treated as a generalized likelihood rather than as the probability mass function of the observed weak ranks. Concrete examples, including the \(n=2,J=1\) tied case and a binary-threshold calculation, are given in the Supplement.

For mixed margins, posterior concentration of the extended-rank posterior requires more than coarsened-rank identifiability. We use a high-level regularity condition, stated formally as Assumption S.2 on “extended-rank contrast consistency” in the Supplementary Material: under the true pair \((\psi_0,F)\), the normalized extended-rank log likelihood admits a deterministic large-sample contrast that is uniquely maximized at \(\psi_0\), and convergence to this contrast is uniform over the compact quotient space \(\Psi_K\). This is the generalized-posterior analogue of ordinary likelihood separation. We do not verify this condition for all mixed-margin DDE copulas in this manuscript; verification would require a separate analysis of high-dimensional order-constrained probabilities and is left for future work.

The proof is a generalized posterior concentration argument. Uniform convergence transfers the strict separation of the limiting contrast to the finite-sample extended-rank log likelihood outside any fixed neighborhood of \(\psi_0\). The posterior numerator outside the neighborhood is then exponentially smaller than a denominator lower bound obtained by integrating over a positive-prior-mass set near the contrast maximizer. The full proof is given in the Supplement.

*Interpretation. Assumption S.2 is stronger than coarsened-rank identifiability. Injectivity of \(\psi\mapsto Q_{\psi,F}\) says that infinite weak ranks contain enough information to distinguish \(\psi\); the contrast condition additionally requires the particular extended-rank generalized likelihood to exploit this information in a uniformly separating way. The simulation settings below are favorable for this condition, with rich discrete margins and identifiable DDE graph structure, and the empirical recovery patterns are consistent with the predicted concentration behavior.

4 Bayesian Specification and Estimation↩︎

We perform Bayesian inference for \(\boldsymbol{\theta}\) by targeting the extended rank likelihood posterior \[p\{\boldsymbol{\theta} \mid \boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y})\} \propto \int p\{\boldsymbol{A}, \boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y}) \mid \boldsymbol{\theta}\} \, p(\boldsymbol{\theta}) \, d \boldsymbol{A}. \label{eq:RLDDEpost}\tag{17}\]

Two key challenges arise. First, we seek to infer the layer-specific widths \(\{K^{(d)}\}_{d=1}^{D}\) rather than fixing them a priori, as in [9]. To this end, we introduce cumulative shrinkage process (CSP; [7]) priors on the elements of each \(\boldsymbol{B}^{(d)}\), which induce sparsity and enable automatic dimension selection within each layer. While CSP priors have been used to infer shallow-layer widths in related models [8], our formulation extends this mechanism to all layers of the DDE. In doing so, we enable information sharing using layer-specific activity indicators.

Second, direct evaluation of 17 is intractable. The rank likelihood involves the integral \(p\{\boldsymbol{A}, \boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y}) \mid \boldsymbol{\theta}\} = \int_{\mathcal{R}(\boldsymbol{Y})} p(\boldsymbol{A}, \boldsymbol{Z} \mid \boldsymbol{\theta}) \, d\boldsymbol{Z}\), while the posterior 17 also marginalizes over the high-dimensional latent \(\boldsymbol{A}\). Although data augmentation [18] is a natural approach for MCMC-based approaches, the joint latent space \((\boldsymbol{Z}, \boldsymbol{A})\) is high-dimensional and strongly dependent, making such methods computationally inefficient.

To address this, we develop a stochastic optimization procedure inspired by expectation–maximization [19] and its recent variants for Bayesian dimension selection in linear regression and Gaussian factor models [6], [20]–[22]. The resulting algorithm targets maximum a posteriori (MAP) estimation of 17 , yielding point estimates of the layer-specific weight matrices along with MAP estimates of the effective dimension in each layer.

4.1 Prior Elicitation↩︎

To complete the Bayesian specification of the DDE copula under extended rank likelihood, we assume prior independence, \(p(\boldsymbol{\theta}) = \prod_{d=1}^{D} p(\boldsymbol{B}^{(d)}) \prod_{j=1}^{J} p(\gamma_{j}) \prod_{k=1}^{K^{(D)}} p(\pi_{k}),\) with \(\gamma_{j} \sim InverseGamma(a,b)\) and flat priors \(\pi_{k} \propto 1\).

We initialize a maximal-width DDE for \(\boldsymbol{Z}\) and induce sparsity through the prior on \(\boldsymbol{B}^{(d)}\), which also enables layer-wise dimension selection. Let \(K^{(d)}_{\text{max}}\) denote the initialized width of layer \(d\), chosen according to the identifiability conditions in Section 3 as \(K^{(d)}_{\text{max}} = \lfloor K^{(d-1)}/3 \rfloor\). For example, in a two-layer DDE with \(J=50\), we take \(K^{(1)}_{\text{max}} = 16\) and \(K^{(2)}_{\text{max}} = 5\). The goal is then to shrink redundant dimensions through the prior.

For each element \(\beta_{jk}^{(d)}\) of \(\boldsymbol{B}^{(d)}\), we impose independent cumulative shrinkage process (CSP) priors. Specifically, for \(j = 1,\dots,K^{(d-1)}_{\text{max}}\) and \(k = 1,\dots,K^{(d)}_{\text{max}}\), \[\begin{align} \beta^{(d)}_{jk} \sim \{1- \mathbb{1}(c^{(d)}_{k} \leq k)\}\,\phi(\lambda^{(d)}_{0k}) &+ \mathbb{1}(c^{(d)}_{k} \leq k)\,\phi(\lambda^{(d)}_{1}), \tag{18}\\ p(c^{(d)}_{k} = \ell \mid \omega^{(d)}_\ell) &= \omega^{(d)}_{\ell}, \quad \ell = 1,\dots,K^{(d)}_{\text{max}},\\ \omega_{\ell}^{(d)} = v_{\ell}^{(d)} \prod_{m<\ell}(1-v_{m}^{(d)}), & \quad v_{m}^{(d)} \sim Beta(1,\alpha^{(d)}). \tag{19} \end{align}\] Here, \(\phi(\lambda)\) denotes the Laplace density with rate \(\lambda\), differing from the Gaussian formulation in [7]. This choice yields exact zeros under MAP estimation and facilitates dimension selection; See Section 4.2. We assign hyperpriors \(\lambda_{0k}^{(d)} \sim Gamma(a^{(d)}, b^{(d)})\) and fix \(\lambda_{1}^{(d)} \ll \mathbb{E}[\lambda_{0k}^{(d)}]\) to distinguish spike and slab components. The parameter \(\lambda_{1}^{(d)}\) may be layer-specific; see the supplement for tuning strategies.

Marginalizing over \(c_{k}^{(d)}\), the prior for \(\beta_{jk}^{(d)}\) is a spike-and-slab mixture with ordered shrinkage: the probability of assignment to the spike increases with \(k\). Consequently, the effective layer width is \(K^{(d)*} = \sum_{k=1}^{K^{(d)}_{\text{max}}} \mathbb{1}(c_{k}^{(d)} > k),\) since \(\mathbb{1}(c_{k}^{(d)} > k)\) indicates whether the \(k\)th column of \(\boldsymbol{B}^{(d)}\) is drawn from the slab. When \(c_{k}^{(d)} \leq k\), the corresponding column is shrunk toward zero and the associated latent unit is redundant. Thus, the CSP prior induces a prior for the unknown layer-specific dimension \(K^{(d)^{*}}\) through the spike and slab indicator variables \(c_{k}^{(d)}\). Finally, we fix \(\alpha^{(d)} = K^{(d)}_{\text{max}}\), motivated by that in the infinite-width limit the expected number of active components satisfies \(\sum_{\ell=1}^{\infty} \mathbb{E}[\mathbb{1}(c_{k}^{(d)} > k)] = \alpha^{(d)}\) [7].

4.2 MAP Estimation via Coordinate Ascent Monte Carlo EM↩︎

We proceed to outline our expectation–maximization (EM; [19]) algorithm targeting maximum a posteriori (MAP) estimation of 17 . Under the CSP prior 18 –19 , the parameter set \(\boldsymbol{\theta}\) is augmented to include the layer-specific dimensions \(\{K^{(d)*}\}_{d=1}^{D}\).

Let \(\boldsymbol{\lambda}_{0}^{(d)} = \{\lambda^{(d)}_{0k}\}_{k=1}^{K^{(d)}_{\text{max}}}\) and \(\boldsymbol{v}^{(d)} = \{v^{(d)}_{k}\}_{k=1}^{K^{(d)}_{\text{max}}}\). We partition the full parameter set into model parameters \(\boldsymbol{\theta} = (\{\boldsymbol{B}^{(d)}\}_{d=1}^{D}, \{K^{(d)*}\}_{d=1}^{D}, \{\gamma_j\}_{j=1}^{J}, \boldsymbol{\pi})\), latent variables \(\boldsymbol{\theta}_{\text{lat}} = (\{\boldsymbol{A}^{(d)}\}_{d=1}^{D}, \boldsymbol{Z})\), and hyperparameters \(\boldsymbol{\theta}_{\text{hyp}} = (\{\boldsymbol{v}^{(d)}\}_{d=1}^{D}, \{\boldsymbol{\lambda}_{0}^{(d)}\}_{d=1}^{D})\). Given an initialization \(\boldsymbol{\theta}^{0}\), the EM update at iteration \(t+1\) is \[\begin{align} \boldsymbol{\theta}^{t+1} &= \arg\max_{\boldsymbol{\theta}} \mathbb{E}_{\boldsymbol{\theta}_{\text{lat}}, \boldsymbol{\theta}_{\text{hyp}} \mid \boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y}), \boldsymbol{\theta}^t} \big[\log p\{\boldsymbol{\theta}, \boldsymbol{\theta}_{\text{lat}}, \boldsymbol{\theta}_{\text{hyp}} \mid \boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y})\}\big] \nonumber\\ &= \arg\max_{\boldsymbol{\theta}} Q(\boldsymbol{\theta} \mid \boldsymbol{\theta}^t). \label{eq:Qfn95tight95mid2} \end{align}\tag{20}\]

The hierarchical structure of the DDE and prior independence across layers yields a layer-wise decomposition of 20 \[\label{eq:Qfndecomp} Q(\boldsymbol{\theta} \mid \boldsymbol{\theta}^{t}) = \sum_{d=1}^{D} Q(\boldsymbol{\theta}^{(d)} \mid \boldsymbol{\theta}^{t}).\tag{21}\] This enables layer-wise updates within the optimization routine, reducing the M-step to a collection of independent optimization problems. In particular, for \(2 \leq d < D\), each \(\boldsymbol{B}^{(d)}\) is updated via \(K^{(d-1)}_{\text{max}}\) independent regression problems, while the first layer is optimized using an expectation–conditional maximization step [23], alternating updates of \(\boldsymbol{B}^{(1)}\) and \(\boldsymbol{\gamma}\). At the top layer, \(\boldsymbol{\theta}^{(D)} = \boldsymbol{\pi}\). The effective layer widths \(\{K^{(d)*}\}\) are determined as a function of the updated spike indicators \(\{c_{k}^{(d)}\}\) (Algorithm 3).

The main computational challenge is evaluating expectations under \(p\{\boldsymbol{\theta}_{\text{lat}}, \boldsymbol{\theta}_{\text{hyp}} \mid \boldsymbol{Z} \in \mathcal{D}(\boldsymbol{Y}), \boldsymbol{\theta}^t\}\), which involves (intractable) integration over the rank-constrained set \(\mathcal{D}(\boldsymbol{Y})\) and summation over exponentially many configurations of \(\{\boldsymbol{A}^{(d)}\}\).

We address this using Monte Carlo EM [5]. At each iteration, we draw \(C\) samples \((\boldsymbol{\theta}_{\text{lat}}^{(c)}, \boldsymbol{\theta}_{\text{hyp}}^{(c)})\) from the conditional posterior and approximate \(Q(\boldsymbol{\theta} \mid \boldsymbol{\theta}^t) \approx \frac{1}{C} \sum_{c=1}^{C} \log p\{\boldsymbol{\theta}, \boldsymbol{\theta}_{\text{lat}}^{(c)}, \boldsymbol{\theta}_{\text{hyp}}^{(c)} \mid \boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y})\} = \hat{Q}(\boldsymbol{\theta} \mid \boldsymbol{\theta}^t).\) Following [9], we set \(C=1\) and observe strong empirical performance.

Sampling from \(p\{\boldsymbol{\theta}_{\text{lat}}, \boldsymbol{\theta}_{\text{hyp}} \mid \boldsymbol{Z} \in \mathcal{D}(\boldsymbol{Y}), \boldsymbol{\theta}^t\}\) is carried out via Gibbs updates for \(\boldsymbol{Z}\), \(\{\boldsymbol{A}^{(d)}\}\), \(\{\boldsymbol{\lambda}_{0}^{(d)}\}\), and \(\{\boldsymbol{v}^{(d)}\}\). The latter two remain conditionally conjugate under the CSP prior and are deferred to the supplement, while we detail the updates for \(\boldsymbol{Z}\) and \(\{\boldsymbol{A}^{(d)}\}\) below.

*Extended Rank Likelihood Re-sampling For \(\boldsymbol{Z} \sim p\{\boldsymbol{Z} \mid \boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y}), \boldsymbol{\theta}^t, \text{--}\}\), the RL implies \[p\{\boldsymbol{Z}_i \mid \boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y}), \theta^{(t)},\text{--}\}\sim \prod_{j=1}^{J} N(\beta^{(1)}_{j,0} + \sum_{k=1}^{K^{(1)}_{\text{max}}} \beta_{j,k}^{(1)} \boldsymbol{\alpha}^{(1)}_{i,k}, \gamma_{j})\mathbb{1}\{\boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y})\}\label{eq:RLsamp}.\tag{22}\] That is, each component of the vector \(\boldsymbol{Z}_i\) is normally distributed subject to the rank consistency between the entire latent data matrix \(\boldsymbol{Z}\) and observed data matrix \(\boldsymbol{Y}\). This implies that each \(Z_{ij}\) may be sampled from a truncated normal distribution, with lower and upper bounds determined by the largest \(Z_{\ell j}\) such that \(Y_{\ell j} < Y_{ij}\) and smallest \(Z_{r j}\) such that \(Y_{r j} > Y_{ij}\), respectively. Due to conditional independence, the joint distribution of \(\boldsymbol{Z}_i\) is simply the product of similarly formed truncated normals. Thus, sampling from 22 is simple, requiring a sequence of univariate truncated normal sampling steps, which are parallelizable across \(J\). We summarize this procedure in Algorithm 2.

*Binary Latent Variables Sampling from \(p[\{\boldsymbol{A}^{(d)}\}_{d=1}^{D} \mid \boldsymbol{\theta}^t,\boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y}), \text{---}]\) is computationally challenging due to the exponential size of the joint state space. Rather than performing a sequential Gibbs sweep over individual coordinates, we employ an approximation based on a product-form (mean-field) representation of the full conditional. For each observation \(i\), we approximate the joint conditional distribution by independent Bernoulli updates across layers and coordinates, \[\label{eq:MFVI} p(\boldsymbol{A}_i\mid \boldsymbol{Z}^t_i, \boldsymbol{\theta}^t, \text{--}) \approx \prod_{d=1}^{D} \prod_{k=1}^{K^{(d)}} \text{Bernoulli}(\pi_{ik}^{(d)}),\tag{23}\]

where \(\boldsymbol{A}_i = (\boldsymbol{A}_i^{(1)}, \dots, \boldsymbol{A}_i^{(D)})\). When used within the EM algorithm, we find that sampling from this approximation enables efficient exploration of the posterior relative to a true Gibbs step, resulting in a significantly more accurate estimation of the layer-specific dimensions \(\{K^{(d)}\}_{d=1}^{D}\) (Section 5).

At iteration \(t\), the update proceeds by computing local log-odds for each latent variable while holding all other coordinates fixed at their current values. Specifically, for each \(d \in \{1,\dots, D\}\) and \(k \in \{1, \dots, K^{(d)}_{\text{max}}\}\) we evaluate \(\Delta_{ik}^{(d),t} = \log p(A_{ik}^{(d)} = 1 \mid \boldsymbol{A}^{-(d), t}_{i,-k}, \boldsymbol{Z}^t_i, \boldsymbol{\theta} ^t, \text{--}) - \log p(A_{ik}^{(d)} = 0 \mid \boldsymbol{A}^{-(d), t}_{i,-k}, \boldsymbol{Z}^t_i, \boldsymbol{\theta} ^t, \text{--}),\) and set \(\pi_{ik}^{(d), t} = \operatorname{logistic}\!\big(\Delta_{ik}^{(d),t}\big)\). Here, \(\boldsymbol{A}^{-(d), t}_{i,-k}\) are the binary latent variables without \(A_{ik}^{(d)}\). The quantities \(\Delta_{ik}^{(d),t}\) admit simple closed forms due to the conditional independence structure of the DDE, which may be found in the supplementary materials. All log-odds are computed using the same previous configuration \(\boldsymbol{A}_i^{t-1}\), allowing updates to be performed in parallel across coordinates and layers. Given \(\{\pi_{ik}^{(d),t}\}\), we then sample independently: \(A_{ik}^{(d),t} \sim \text{Bernoulli}(\pi_{ik}^{(d),t}), \;\text{for all}\;d,k\).

This update can be viewed as an approximate Gibbs step: rather than sequentially sampling from exact full conditionals, we construct a product-form approximation and draw all coordinates simultaneously. The resulting update resembles a mean-field variational step [24], but instead of optimizing the Bernoulli parameters governing each \(A_{ik}^{(d)}\), we sample directly from the approximate family. Although this procedure doesn’t produce exact samples from the target distribution, it yields substantial computational gains relative to a full Gibbs sampling step and performs well in high-dimensional settings; see Section 5 for detailed comparisons.

The full coordinate ascent Monte Carlo EM algorithm for MAP estimation of the \(D\)-layer DDE copula is summarized in Algorithm 3, with details provided in the supplement. The Monte Carlo E-step also simplifies the maximization step. When \(C=1\), the decomposition of \(Q(\boldsymbol{\theta} \mid \boldsymbol{\theta}^t)\) in 21 , together with conditional independence across nodes within each layer, reduces the M-step to a collection of independent, parallelizable weighted-\(\ell_1\)-regularized logistic or linear regressions. For the \(j\)th row in layer \(d\), the penalty is \(\sum_{k=1}^{K^{(d)}_{\text{max}}} \Big[ (1-\mathbb{1}\{c_{k}^{(d),t+1} \leq k\}) \lambda^{(d),t+1}_{0k} + \mathbb{1}\{c_{k}^{(d),t+1} \leq k\} \lambda^{(d)}_{1} \Big] \lvert \beta_{jk} \rvert,\) arising from the Laplace mixture specification of the CSP prior. The current estimate of \(c_{k}^{(d)}\) determines whether the \(k\)th coefficient is penalized by the spike rate \(\lambda^{(d)}_{1}\) or the slab rate \(\lambda^{(d),t+1}_{0k}\).

This thresholding scheme differs from [6], which employs a soft weighting of spike and slab penalties in factor models. This hard-thresholding also induces hierarchical gating: if \(c_{k}^{(d),t+1} \leq k\) in layer \(d\), we deterministically set the \(k\)th row of \(\boldsymbol{B}^{(d+1),t+1}\) to zero. Empirically, this improves convergence and estimation accuracy.

Despite the efficiency gains from the approximate Gibbs update 23 , the EM algorithm may exhibit sensitivity to local modes in the high-dimensional posterior. To mitigate this, we incorporate deterministic annealing [20] by introducing a temperature parameter \(\tau \in (0,1]\) and targeting a tempered posterior proportional to \(p\{\boldsymbol{\theta}, \boldsymbol{\theta}_{\text{lat}}, \boldsymbol{\theta}_{\text{hyp}} \mid \boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y})\}^{\tau}\). This modification induces simple, layer-wise adjustments to both the E-step and M-step; see the supplement for more details.

4.3 Spectral Initialization of Latent Variables and Parameters↩︎

Initialization of both \(\boldsymbol{\theta}\) and \(\boldsymbol{\theta}_{\text{hyp}}\) is paramount to accurate estimation of the DDE copula. Although the rank likelihood enforces marginal rank-consistency between \(\boldsymbol{Y}\) and \(\boldsymbol{Z}\), the marginal distributions of \(Y_j\) and \(Z_j\) may differ substantially. In particular, \(Y_j\) may be discrete, skewed, or have zero-inflation, while the corresponding latent variable \(Z_j\) is typically a factorial mixture of Gaussians induced by the binary latent layers. Moreover, dependence in \(\boldsymbol{Z}\) arises from a mixture of Gaussian components that is subsequently warped by the marginal transformations \(F_{Y_j}\), making direct initialization difficult.

To address this, we propose a spectral initialization procedure that exploits local dependence structure in \(\boldsymbol{Y}\) to recover latent clustering behavior. The key heuristic is that under the DDE copula, \(Z_j = G_{Z_j}^{-1}\{F_{Y_j}(Y_j)\}\), so that \(Z_j\) is a monotone transformation of \(Y_j\). We approximate this transformation by the identity and instead aim to recover latent structure through pairwise relationships in \(\boldsymbol{Y}\).

We begin by computing the singular value decomposition \(\boldsymbol{Y} = \boldsymbol{U} \boldsymbol{\Sigma} \boldsymbol{V}^\top\) and defining a distance matrix \(D_{ij} = \|\boldsymbol{V}_i - \boldsymbol{V}_j\|_2^2\), where \(\boldsymbol{V}_j\) denotes the \(j\)th row of \(\boldsymbol{V}\), the right singular vectors of \(\boldsymbol{Y}\). This representation captures dependence between variables through their spectral embeddings. To reduce noise and mitigate issues arising from discrete observations, we construct a low-rank approximation \(\boldsymbol{Y}^* = \boldsymbol{U}_k \boldsymbol{\Sigma}_k \boldsymbol{V}_k^\top\) using the leading \(k\) singular values. We select \(k\) as the smallest integer satisfying \(\frac{\sum_{\ell=1}^{k} \Sigma_{\ell\ell}^{2}}{\sum_{\ell=1}^{\min(N,J)} \Sigma_{\ell\ell}^{2}} \geq 0.8,\) where \(\Sigma_{\ell\ell}\) denotes the \(\ell\)th singular value in the diagonal matrix \(\boldsymbol{\Sigma}\).

For each variable \(j\), we identify its nearest neighbor \(i^*(j) = \arg\min_{i \neq j} D_{ij}\) and consider the bivariate data \((Y^*_{i^*(j)}, Y^*_j)\). Empirically, these pairs often exhibit clustering structure that reflects shared latent factors; see the bottom-left panel of Figure 4. We apply \(k\)-means clustering to this bivariate data, selecting the number of clusters via the elbow method, and use the resulting cluster centers and weights to define an isotropic Gaussian mixture model. Drawing \(n\) i.i.d.samples \((Z^{*}_{i^*(j)},Z^{*}_j)\) from this fitted model yields an approximation to the latent distribution (\(Z_{i(j)}, Z_{j})\). Finally, we enforce rank-consistency by matching order statistics, setting \(Z_{(r),j} = Z^{*}_{(r),j}\) for \(r=1,\dots,n\). This process is outlined in Algorithm 5.

This procedure leverages the fact that pairs \((Z_i, Z_j)\) are mixtures of low-dimensional Gaussian components induced by shared latent structure. The spectral embedding identifies strongly dependent variables, while the de-noised representation reveals clustering patterns that are otherwise obscured on the observed scale. Matching order statistics ensures compatibility with the rank likelihood, producing an initialization that respects the marginal ordering constraints.

Once we have initialized \(\boldsymbol{Z}\), we must then initialize \(\boldsymbol{\theta}\) and \(\{\boldsymbol{A}^{(d)}\}_{d=1}^{D}\). Because we are effectively estimating a Gaussian DDE on rank-consistent \(\boldsymbol{Z}\), we may proceed using the layer-wise Double-SVD initialization strategy in [9]. The basic idea is to denoise \(\boldsymbol{Z}\) with a singular value decomposition, which provides initial estimates of \((\hat{\boldsymbol{B}}^{(1)}, \hat{\boldsymbol{A}}^{(1)})\). Then, the algorithm proceeds recursively, treating initialized \(\hat{\boldsymbol{A}}^{(d)}\) as observed data for a denoising process that enables initialization of \((\boldsymbol{B}^{(d-1)}, \boldsymbol{A}^{(d-1)})\). Complete details on this process are provided in the Supplementary Material.

Importantly, the initialization produced by Algorithm 5 need not satisfy the normalization \(\mathbb{E}(Z_j)=0\) and \(\mathrm{Var}(Z_j)=1\) in Section 3. This is inconsequential for estimation, as the DDE copula depends only on the induced rank structure of \(\boldsymbol{Z}\) and is invariant to marginal affine transformations. In practice, we therefore fit the DDE model directly to the unnormalized \(\boldsymbol{Z}\) obtained from the spectral procedure. After fitting, we impose the identification constraint by re-centering and re-scaling each coordinate of \(\boldsymbol{Z}\), with a corresponding adjustment to the first-layer parameters \((\boldsymbol{B}^{(1)}, \boldsymbol{\gamma})\) to preserve the likelihood. This post hoc normalization yields an identified representation without affecting the estimated dependence structure.

5 Simulation: Dimension and Parameter Estimation↩︎

We conduct a simulation study to evaluate the performance of the proposed Bayesian extended rank likelihood DDE copula model under layer-wise CSP priors in recovering latent structure and estimating model parameters. The primary goals are to assess: (i) recovery of the latent dimensions \(\{K^{(d)}\}_{d=1}^{(D)}\), and (ii) accuracy of the estimated loading matrices \(\{\boldsymbol{B}^{(d)}\}_{d=1}^{D}\).

Data are generated from a two-layer DDE copula model with \(K^{(1)} = 10\) first-layer latent variables and \(K^{(2)} = 3\) second-layer latent variables. The observed dimension \(J\) of \(\boldsymbol{Y}\) varies across settings, with \(J \in \{50, 100, 150\}\).

The loading matrices are constructed to exhibit structured sparsity, see Figure 6 for the exact patterns. To mimic a plausible real data setting, groups of variables in \(\boldsymbol{Y}\) load strongly (both positively and negatively) on certain shallow-level latent variables, which also share deeper latent features. Crucially, \(K^{(d)}<< K_{\text{max}}^{(d)}\) for each \(J\), providing a challenging setting for dimension selection and weight matrix estimation.

Given these parameters, data are generated from the DDE copula using the Gaussian latent construction 5 –6 . To induce heterogeneous marginals, \(\{F_{Y_j}\}_{j=1}^{J}\) are assigned cyclically via \(\mathrm{mod}(j-1,3)+1\). Let \(Q_j = F_{Z_j}(Z_j)\). For type 1 variables, we generate zero-inflated Poisson marginals by setting \(X_j=0\) with probability \(\pi_0=0.3\), and otherwise applying \(X_j = F^{-1}_{\text{Pois}(r_j)}\!\big((Q_j-\pi_0)/(1-\pi_0)\big)\) for \(Q_j \ge \pi_0\). For type 2 variables, we use discretized Gamma marginals, \(X_j = \mathrm{round}\!\big(F^{-1}_{\Gamma(2,r_j)}(Q_j)\big)\), yielding over-dispersed counts with mean approximately \(2r_j\). For type 3 variables, we use standard Poisson marginals, \(X_j = F^{-1}_{\text{Pois}(r_j)}(Q_j)\). For each \(j\), \(r_j\) is drawn uniformly from \(\{1,\dots,10\}\). The resulting data exhibit dependence governed by the DDE copula, while retaining heterogeneous marginal scales and shapes.

Across 100 simulated datasets with \(n \in \{500,1000,2000,4000,8000,16000\}\), we initialize rank-consistent latent variables \(\boldsymbol{Z}\) via Algorithm 5, and initialize \(\boldsymbol{\theta}\) and \(\{\boldsymbol{A}^{(d)}\}_{d=1}^{2}\) using the Double SVD approach of [9] given \(\boldsymbol{Z}\). All weight matrices are initialized at their maximal admissible dimensions under the identifiability conditions of Section 3: \(\boldsymbol{B}^{(1)} \in \mathbb{R}^{J \times (\lfloor J/3 \rfloor + 1)}\) and \(\boldsymbol{B}^{(2)} \in \mathbb{R}^{(\lfloor J/3 \rfloor) \times (\lfloor J/9 \rfloor + 1)}\). As \(J\) increases, this induces increasing misspecification in the initial dimensionality.

We evaluate the ability of the Bayesian extended rank likelihood DDE copula with independent layer-wise CSP priors (RL CSP) to recover (i) the true layer dimensions \(\{K^{(d)}\}\), (ii) the sparsity structure \(\{G^{(d)}\}_{d=1}^{2}\) with \(g^{(d)}_{jk} = \mathbb{I}(\beta^{(d)}_{jk} \neq 0)\), and (iii) the weight matrices \(\{\boldsymbol{B}^{(d)}\}_{d=1}^{2}\). Hyperparameter choices are provided in the supplement; in practice, as \(J\) increases, we increase \(\lambda_{1}^{(d)}\) and decrease the initial temperature to encourage stronger shrinkage in increasingly sparse regimes.

We compare against three alternative approaches. First, RL Gibbs replaces the approximate Gibbs step in Algorithm 3 with exact Gibbs updates for \(\boldsymbol{A}^{(d)}\), holding other steps fixed. Second, RL SAEM replaces the Bayesian CSP prior with a maximum likelihood approach based on a stochastic approximation EM algorithm [9], augmented with rank-likelihood sampling (Algorithm 2) and a truncated lasso penalty [25] to induce sparsity and estimate latent dimension. Finally, we consider a Poisson DDE that omits the Gaussian copula layer 5 –6 , assuming conditionally Poisson outcomes in 3 that is compatible with the discrete marginals in \(\boldsymbol{Y}\). All methods are initialized using the same \(\boldsymbol{Z}\) (when applicable) and \(\boldsymbol{\theta}\).

We report results for \(J=100\) observed variables, and observe qualitatively similar patterns in other settings. Table 1 summarizes average layer-wise graph recovery and MAP estimates of the latent dimensions across methods. Only RL CSP and RL Gibbs estimate \(K^{(d)*}\), as latent dimension is treated as an unknown parameter and updated via Algorithm 3. RL CSP achieves consistently higher sparsity recovery and near-perfect dimension estimation, whereas competing approaches degrade substantially in deeper layers.

These patterns are reflected in Figures 7-8, which display the average entry-wise MSE between \(\boldsymbol{B}^{(d)}\) and estimates \(\hat{\boldsymbol{B}}^{(d)}\) and the average estimate \(\boldsymbol{\hat{B}}^{(2)}\) across sample sizes and methods, respectively. Relative to other rank-likelihood-based estimators, RL CSP yields uniformly lower error, highlighting the combined benefits of the approximate latent-variable sampling step and the CSP shrinkage prior. The Poisson DDE performs poorly across all settings, indicating that parametric DDE formulations—even when aligned with the data type—lack sufficient flexibility to capture complex marginal structure. Additional results in the supplement, including for deeper models (\(D > 2)\), show similarly strong performance for RL CSP in both dimension recovery and weight estimation.

Finally, we note that fitting was fast: for The RL CSP DDE Copula, average run times were 2.83 minutes for \(n = 1000\), 4.25 minutes for \(n = 4000\), and 26.89 minutes for \(n = 16000\) without parallel computation (we ran simulation iterations in parallel for each sample size). However, our code implementation is parallelized for the sampling of \(\boldsymbol{A}, \boldsymbol{Z}\) and optimization of \(\boldsymbol{B}^{(d)}\), which significantly decreases run time in large \(n\) and \(J\) regimes.

6 Application to Personality and Political Leanings↩︎

We apply the proposed method to analyze an online survey of 1500 respondents conducted by YouGov linking personality traits, demographic characteristics, and political ideology [26]. Respondents were matched to a sampling frame on gender, age, race, and education, which was constructed by stratified sampling from the full 2016 American Community Survey (ACS; [27]). We subsequently refer to these data as the Big5 survey.

Participants were asked to self-rate themselves across 98 items (we exclude two items measuring likelihood to vote for liberal or conservative political candidates due to our downstream analysis), consistent with the Big Five construct of personality [28]. Their answers were recorded on an ordinal scale ranging from 1 ("very inaccurate") to 5 ("very accurate"). An example item is displayed in Table 2.

Surveys of personality exhibit complex between- and within-trait dependencies [29]. In addition to the widely studied Big Five traits (Agreeableness, Openness, Neuroticism, Conscientiousness, and Extraversion), prior work suggests the presence of higher-order latent structure [30]. Estimating such structure is challenging, as few methods simultaneously (i) accommodate ordinal discrete data, (ii) capture nonlinear and hierarchical dependence, and (iii) yield interpretable summaries of these relationships.

We fit a two-layer RL CSP DDE copula to the Big5 survey, with \(K^{(1)}_{\max} = \lfloor 98/3 \rfloor\) and \(K^{(2)}_{\max} = \lfloor K^{(1)}_{\max}/3 \rfloor\). Model and hyperparameter selections were based on predictive evaluation. For each candidate configuration of hyperparameters, we estimated \(\hat{\boldsymbol{\theta}}\) and generated \(50\) synthetic datasets of size \(n = 1500\) via sequential sampling of \(\boldsymbol{\tilde{A}}_i^{(2), \ell}\), \(\boldsymbol{\tilde{A}}_i^{(1), \ell}\), and \(\boldsymbol{\tilde{Z}}_i^{\ell}\) for \(\ell = 1,\dots,50\), followed by \(\tilde{\boldsymbol{Y}}_i^{\ell} = \hat{F}_{Y_j}\{F_{Z_j}(\tilde{Z}_{ij}^{\ell})\}\), where \(\hat{F}_{Y_j}\) denotes the empirical distribution function of \(Y_j\).

We then combined the observed and synthetic data and trained a Random Forest classifier [31] using five-fold cross-validation. For each hold-out set, we computed the propensity score mean squared error (pMSE), \(\text{pMSE}_k^{\ell} = n_{\text{test}}^{-1} \sum_{i=1}^{n_{\text{test}}} (\hat{p}_i - 0.5)^2\) [32], and averaged across folds to obtain \(\overline{\text{pMSE}}^{\ell} = 5^{-1} \sum_{k=1}^{5} \text{pMSE}_k^{\ell}\). The pMSE quantifies joint distributional similarity between real and synthetic data; with equal proportions of each, a strong generative model produces synthetic data that is indistinguishable from the real observations and so the classifier assigns probabilities near 0.5. Thus, the DDE copula with the lowest average pMSE best approximates the Big5 joint distribution. Additional univariate and multivariate diagnostics further support the selected model. We selected the model minimizing \(50^{-1} \sum_{\ell=1}^{50} \overline{\text{pMSE}}^{\ell}\), yielding \(K^{(1)*} = 9\) and \(K^{(2)*} = 1\).

To aid interpretation, we identified, for each active shallow-layer latent variable \(k = 1,\dots,9\), the three items with strongest relative loadings using \(\max\{\min_{l \neq k} (\beta_{jk} - \beta_{jl}),\, 0\}\). The resulting items (Table 3) are non-overlapping and thematically coherent within each node, suggesting a structured latent representation beyond the classical Big Five. The supplement includes detailed comparisons to an identifiable variational autoencoder (iVAE; [13]), constructed to mirror the discovered latent dimension of the DDE copula. The neural network weights for both the encoder and decoder in iVAE are not nearly as interpretable as those discovered using our method.

2pt

At the deeper layer, \(A_1^{(2)}\) partitions the shallow-layer factors into two groups: \((A_1^{(1)}, A_4^{(1)}, A_6^{(1)}, A_7^{(1)}, A_9^{(1)})\) load negatively, while \((A_2^{(1)}, A_3^{(1)}, A_5^{(1)}, A_8^{(1)})\) load positively in \(\hat{\boldsymbol{B}}^{(2)}\). The associated items suggest a dichotomy between more withdrawn and reactive traits and those reflecting greater social engagement and self-regulation.

Finally, we examined the relationship between latent representations and political ideology. After convergence, we extracted \(\{(\boldsymbol{\hat{A}}_i^{(1)}, \boldsymbol{\hat{A}}_i^{(2)})\}_{i=1}^{1500}\) and focused on 389 individuals identifying as “Very Liberal” or “Very Conservative” (48.59% and 51.41%, respectively).

We computed group-specific averages \(\bar{A}_{k,\text{ideology}} = |n_{\text{ideology}}|^{-1} \sum_{i \in \text{ideology}} \hat{A}_{ik}^{(1)}, \;\text{ideology} \in \{\text{Very Conservative}, \text{Very Liberal}\}\) and compared these across groups (Figure 9). Several factors, notably \((A_4^{(1)}, A_5^{(1)}, A_7^{(1)}, A_9^{(1)})\), exhibit systematic differences by ideology. These correspond to variation in traits related to artistic inclination, reactivity, self-confidence, and work ethic. Additionally, 46% of Very Liberal individuals have \(\hat{A}_1^{(2)} = 1\), compared to 56% among Very Conservative individuals.

7 Discussion↩︎

We introduce the Bayesian Deep Discrete Encoder Copula, a generative model that accommodates arbitrary and mixed data types. By incorporating layer-wise cumulative shrinkage priors, the framework enables principled estimation of unknown layer widths within the graphical model architecture. We develop scalable initialization and MAP estimation procedures, leveraging a Monte Carlo EM algorithm with variationally motivated approximate sampling steps for latent variables, yielding tractable maximization for the DDE parameters.

Theoretically, we establish identifiability of the canonical DDE copula from the copula law under transparent layerwise pure-child and separation conditions. The argument combines the tensor-identifiability machinery for DDEs with a new recovery step for the normalized Gaussian copula layer. The rank-likelihood results then clarify the large-sample inferential target: exact-rank posterior consistency on the quotient space for continuous margins, and generalized-posterior concentration for tied margins under the supplementary extended-rank contrast condition.

In a comprehensive simulation study, the proposed method demonstrates improved recovery of layer-specific weight matrices and latent dimensions relative to existing parametric DDE approaches across a range of challenging settings. Gains are particularly pronounced in deeper layers, underscoring the importance of the copula construction, approximate Gibbs updates, and Bayesian shrinkage. We further apply the model to a moderately sized real-world dataset of personality survey responses measured on discrete ordinal scales. The analysis reveals fine-grained personality structure that departs from classical theory while also recovering coarser higher-order organization. These findings are supported by strong associations between learned binary embeddings and political ideology, with clear contrasts emerging between liberal and conservative respondents.

This work contributes to the growing literature on interpretable unsupervised learning by extending identifiable deep generative models to settings with arbitrary marginal distributions via rank likelihoods. This extension is more than a marginal modeling device: it preserves the interpretation of DDE parameters as features of multivariate dependence after arbitrary monotone marginal transformations. Consequently, the learned binary embeddings and directed edges can be viewed as summaries of dependence rather than artifacts of a chosen parametric model for each observed variable. The proposed framework also provides scalable estimation procedures compatible with modern Bayesian shrinkage priors.

Several directions for future research remain. To enable uncertainty quantification without full MCMC, one approach is to pre-train the binary latent variables \(\boldsymbol{A}\) via a MAP procedure [33], and then condition on \(\boldsymbol{\hat{A}}\) when estimating \(\boldsymbol{\theta}\). This suggests the approximation \(p\{\boldsymbol{\theta} \mid \boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y})\} \approx p\{\boldsymbol{\theta} \mid \boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y}), \boldsymbol{\hat{A}}\},\) which may be theoretically justified in high dimensions. More broadly, the flexibility of the DDE copula makes it well suited for extracting interpretable structure from multi-modal data. Finally, our analysis of Big Five survey data suggests that personality may exhibit richer structure than traditionally assumed; future work can investigate whether similar patterns arise across other psychological instruments.

*Data Availability The data that support the findings of this study are available from the corresponding author, J.F., upon reasonable request.

*Supplementary Material. The Supplementary Material contains proofs of the theoretical results, computational details, and additional numerical results.

*Funding Yuqi Gu’s research is partially supported by NSF grant DMS-2210796

In Section 8, we present the proofs of the theoretical results. In Section 9, we include detailed descriptions of our proposed coordinate ascent Monte Carlo EM algorithm, as well as details about the SAEM RL algorithm, hyperparameter selections, and methods for addressing latent variable permutations in our simulations. Section 10 contains additional simulation results, completing the details excluded from the Simulation studies in the main paper, and extending the DDE copula to deeper architectures. In Section 11, we include results and comparisons from fitting identifiable variational autoencoders to the Big5 survey data.

8 Proofs of Theoretical Results↩︎

8.1 Proof of Theorem 1.↩︎

The proof proceeds in four steps. First, on the copula scale, the distribution of \(\boldsymbol{U}=F_{\boldsymbol{Z}}(\boldsymbol{Z})\) is a one-layer saturated latent class model with latent variable \(\boldsymbol{A}^{(1)}\). Second, a tensor-decomposition argument identifies the first-layer latent mixing proportions \(\boldsymbol{\eta}^{(1)}\) and the conditional copula-scale distribution functions \[H_{j\boldsymbol{\alpha}}(u)=P(U_j\le u\mid \boldsymbol{A}^{(1)}=\boldsymbol{\alpha})\] up to the unavoidable relabeling of the binary latent coordinates. Third, the special Gaussian structure of the first layer identifies the normalized Gaussian parameters \((\boldsymbol{B}^{(1)}, \boldsymbol{\gamma})\) from the functions \(H_{j\boldsymbol{\alpha}}\). Fourth, once the distribution of \(\boldsymbol{A}^{(1)}\) is identified, the same one-layer argument is applied recursively to identify \((\boldsymbol{B}^{(2)},\boldsymbol{G}^{(2)}),\ldots,(\boldsymbol{B}^{(D)},\boldsymbol{G}^{(D)})\) and the top-layer Bernoulli probabilities \(\pi_1,\ldots,\pi_{K^{(D)}}\).

8.2 Proof of Theorem 2↩︎

The proof is organized through four lemmas.

8.3 Proofs and Additional Details for Section 3.3.2↩︎

This section gives the details omitted from Section 3.3.2: first, why the extended rank likelihood is not an exact likelihood when ties occur; second, what infinite weak-rank/tie information reveals; and third, the proof of Theorem 3.

8.3.1 The extended rank likelihood is not an exact likelihood with ties↩︎

Let \[\boldsymbol{U}_i=(U_{i1},\ldots,U_{iJ}), \qquad i=1,\ldots,n,\] be i.i.d. from the copula \(C_\psi\). For a realized data array \(\boldsymbol{Y}^{(n)}\), the extended rank set is \[\mathcal{D}_n(\boldsymbol{Y}^{(n)}) = \Bigl\{ \boldsymbol{u}\in(0,1)^{n\times J}: y_{ij}<y_{\ell j}\Rightarrow u_{ij}<u_{\ell j}, \quad i,\ell=1,\ldots,n, \;j=1,\ldots,J \Bigr\}.\] The marginal extended rank likelihood is \[L_n^R(\psi;\boldsymbol{Y}^{(n)}) = P_\psi\{\boldsymbol{U}^{(n)}\in\mathcal{D}_n(\boldsymbol{Y}^{(n)})\}.\] When all margins are continuous, ties occur with probability zero and the sets \(\mathcal{D}_n(\boldsymbol{Y}^{(n)})\) are exactly the rank cells corresponding to the observed coordinatewise ranks. These rank cells form a measurable partition of \((0,1)^{nJ}\), so \(L_n^R\) is the exact likelihood of the observed ranks.

When ties can occur, tied observations impose no strict ordering constraints. Consequently, \(\mathcal{D}_n(\boldsymbol{Y}^{(n)})\) is a rank-compatible superset rather than the event corresponding exactly to the observed weak-rank/tie pattern. The sets associated with different weak-rank/tie patterns are no longer disjoint and therefore cannot define a probability mass function for those patterns.

Example S.1: one margin and two observations. Take \(n=2\) and \(J=1\). There are three weak-rank patterns: \[Y_1<Y_2, \qquad Y_2<Y_1, \qquad Y_1=Y_2.\] The corresponding extended rank sets are \[\mathcal{D}_< = \{(u_1,u_2):u_1<u_2\}, \qquad \mathcal{D}_> = \{(u_1,u_2):u_2<u_1\},\] and \[\mathcal{D}_= = (0,1)^2,\] because a tie imposes no strict inequality constraint. Since \(U_1\) and \(U_2\) are i.i.d. uniform random variables, \[P(\mathcal{D}_<)=\frac{1}{2}, \qquad P(\mathcal{D}_>)=\frac{1}{2}, \qquad P(\mathcal{D}_=)=1.\] These three numbers cannot be probabilities of the three mutually exclusive weak-rank patterns, since they sum to \(2\), not \(1\).

The difference from the exact likelihood is already visible for a binary margin. Suppose \[Y_i=1\{U_i>\tau\}, \qquad 0<\tau<1.\] Then the exact probabilities of the three weak patterns are \[P(Y_1<Y_2)=\tau(1-\tau), \qquad P(Y_2<Y_1)=\tau(1-\tau),\] and \[P(Y_1=Y_2)=\tau^2+(1-\tau)^2.\] These probabilities depend on the marginal threshold \(\tau\), whereas the corresponding extended-rank values are \(1/2\), \(1/2\), and \(1\). The extended rank likelihood is therefore margin-free, but it is not an exact likelihood for the weak-rank/tie statistic.

Example S.2: two binary margins and two observations. Take \(n=2\) and \(J=2\), with binary margins \[Y_{ij}=1\{U_{ij}>\tau_j\}, \qquad j=1,2.\] Suppose the observed data are \[\boldsymbol{Y}_1=(0,0), \qquad \boldsymbol{Y}_2=(1,0).\] The first margin implies \(Y_{11}<Y_{21}\), while the second margin is tied, \(Y_{12}=Y_{22}\). Hence \[\mathcal{D}_n(\boldsymbol{Y}^{(n)}) = \{\boldsymbol{u}:u_{11}<u_{21}\},\] with no constraint on \(u_{12}\) or \(u_{22}\). Therefore \[L_n^R(\psi;\boldsymbol{Y}^{(n)}) = P_\psi(U_{11}<U_{21}) = \frac{1}{2},\] because \(U_{11}\) and \(U_{21}\) are independent uniforms from two independent subjects. By contrast, the exact probability of the observed binary data under fixed thresholds is \[\begin{align} &P_{\psi,F}\{\boldsymbol{Y}_1=(0,0),\boldsymbol{Y}_2=(1,0)\} \\ &\qquad = P_\psi(U_{11}\le\tau_1,U_{12}\le\tau_2) P_\psi(U_{21}>\tau_1,U_{22}\le\tau_2) \\ &\qquad = C_\psi(\tau_1,\tau_2) \{\tau_2-C_\psi(\tau_1,\tau_2)\}. \end{align}\] This exact probability depends on both the marginal thresholds and the copula. The extended rank likelihood deliberately avoids the unknown thresholds, but it no longer coincides with the exact likelihood of the observed data or of the observed weak-rank/tie pattern.

8.3.2 Coarsened-rank identifiability↩︎

For a fixed true margin vector \(F=(F_{Y_1},\ldots,F_{Y_J})\), define \[T_{F_j}(u)=F_{Y_j}\{F_{Y_j}^{-1}(u)\}, \qquad T_F(u_1,\ldots,u_J)=\{T_{F_1}(u_1),\ldots,T_{F_J}(u_J)\},\] and \[Q_{\psi,F}=\mathcal{L}_\psi\{T_F(\boldsymbol{U})\}, \qquad \boldsymbol{U}\sim C_\psi.\] For a measurable map \(T\) and probability measure \(\mu\), let \(T\#\mu\) denote the pushforward measure. Thus \(Q_{\psi,F}=T_F\# C_\psi\).

8.3.3 Proof of Theorem 3↩︎

9 Detailed Algorithms↩︎

9.1 Tempering to Reduce Multimodality↩︎

We first review how to apply tempering, which is mentioned in Section 4.2. For the latent binary updates, all log-odds \(\Delta_{ik}^{(d),t}\) are rescaled as \(\tilde{\Delta}_{ik}^{(d),t} \;\leftarrow\; \tau \, \Delta_{ik}^{(d),t}\), which shrinks probabilities toward \(1/2\) for \(\tau < 1\) and promotes exploration across configurations. For the Gaussian latent variables under the rank likelihood, tempering corresponds to inflating the conditional variances, \(Z_{ij} \mid (\text{--}) \sim \mathrm{TN}(\mu_{ij}, \gamma_j / \tau; L_{ij}, U_{ij}),\) which further smooths the latent landscape while preserving rank constraints.

In the maximization step, tempering rescales the contribution of the likelihood in \(\widehat Q(\boldsymbol{\theta} \mid \boldsymbol{\theta}^t)\) by \(\tau\), which is equivalent to applying the same scaling within the regression-based updates for \(\boldsymbol{B}^{(d)}\) and \(\boldsymbol{\gamma}\).

We employ a deterministic schedule \(\{\tau_t\}_{t \geq 0}\) that gradually increases from a moderate value (e.g., \(\tau_0 \approx 0.7\)) to \(\tau=1\), at which point the algorithm targets the approximation to the RL posterior provided in Algorithm 3. Early iterations thus emphasize exploration of the posterior landscape, while later iterations recover the untempered objective for accurate estimation. Empirically, this strategy substantially improves convergence and reduces sensitivity to initialization in deeper DDE models. 5.

9.2 Coordinate Ascent Monte Carlo EM for the DDE Copula↩︎

We provide detailed updates and sampling steps for Algorithm 3 in the main paper, as well as the rank likelihood SAEM algorithm described in the simulations of Section Here, we describe the algorithm for a \(D = 2\) layer DDE copula. The steps may be generalized for \(D > 2\). Let \[\boldsymbol{\theta} = \Bigl( \boldsymbol{\eta,\, B^{(1)},\, B^{(2)},\, \gamma}, \{c_{k}^{(1)}\}_{k=1}^{K_{\text{max}}^{(1)}},\, \{c_{k}^{(2)}\}_{k=1}^{K_{\text{max}}^{(2)}}\Bigr)\] denote the collection of model parameters. At iteration \(t\), given the current state \(\boldsymbol{\theta}^{(t)}\), the coordinate-ascent Monte Carlo EM algorithm alternates between a Monte Carlo E-step and a conditional maximization step. For each subject \(i=1,\dots,N\), let \[\boldsymbol{A}_i^{(2)} \in \{0,1\}^{K_{\text{max}}^{(1)}}, \qquad \boldsymbol{A}_i^{(1)} \in \{0,1\}^{K_{\text{max}}^{(2)}},\] denote the two layers of binary latent variables. We augment each layer with an intercept entry equal to one: \[\widetilde{\boldsymbol{A}}_i^{(1)} = (1, \boldsymbol{A}_i^{(1)}), \qquad \widetilde{\boldsymbol{A}}_i^{(2)} = (1, \boldsymbol{A}_i^{(2)}).\] The latent Gaussian copula variable for observed feature \(j\) is denoted \(Z_{ij}\), and the residual variance for feature \(j\) is \(\gamma_j\).

Monte Carlo E-step. Given \(\boldsymbol{\theta}^{(t)}\), perform the following updates. For all of our implementations, we draw \(C = 1\) samples in each step.

• Sample the binary latent variables \(\{\boldsymbol{A}_{i}^{(1)},\boldsymbol{A}_{i}^{(2)}\}_{i=1}^{N}\)

  • Top layer (\(d = D\)): \[\begin{align} \Delta_{ik}^{(D),t} &= \log \frac{\pi_k^t}{1-\pi_k^t} + \log p\!\left(\boldsymbol{A}_i^{(D-1),t-1} \mid A_{ik}^{(D)} = 1, \boldsymbol{A}_{i,-k}^{(D),t-1}, \boldsymbol{\theta}^t\right)\\ &- \log p\!\left(\boldsymbol{A}_i^{(D-1),t} \mid A_{ik}^{(D)} = 0, \boldsymbol{A}_{i,-k}^{(D),t-1}, \boldsymbol{\theta}^t\right). \end{align}\]

  • Intermediate layers (\(2 \le d \le D-1\)): \[\Delta_{ik}^{(d),t} = \log p\!\left(A_{ik}^{(d)} = 1 \mid \boldsymbol{A}_i^{(d+1),t-1}, \boldsymbol{\theta}^t\right) - \log p\!\left(A_{ik}^{(d)} = 0 \mid \boldsymbol{A}_i^{(d+1),t-1}, \boldsymbol{\theta}^t\right)\] \[\quad + \log p\!\left(\boldsymbol{A}_i^{(d-1),t-1} \mid A_{ik}^{(d)} = 1, \boldsymbol{A}_{i,-k}^{(d),t-1}, \boldsymbol{\theta}^t\right) - \log p\!\left(\boldsymbol{A}_i^{(d-1),t-1} \mid A_{ik}^{(d)} = 0, \boldsymbol{A}_{i,-k}^{(d),t-1}, \boldsymbol{\theta}^t\right).\]

  • Bottom layer (\(d = 1\)): \[\Delta_{ik}^{(1),t} = \log p\!\left(A_{ik}^{(1)} = 1 \mid \boldsymbol{A}_i^{(2),t-1}, \boldsymbol{\theta}^t\right) - \log p\!\left(A_{ik}^{(1)} = 0 \mid \boldsymbol{A}_i^{(2),t-1}, \boldsymbol{\theta}^t\right)\] \[\quad + \log p\!\left(\boldsymbol{Z}_i^{t} \mid A_{ik}^{(1)} = 1, \boldsymbol{A}_{i,-k}^{(1),t-1}, \boldsymbol{\theta}^t\right) - \log p\!\left(\boldsymbol{Z}_i^{t} \mid A_{ik}^{(1)} = 0, \boldsymbol{A}_{i,-k}^{(1),t-1}, \boldsymbol{\theta}^t\right).\]

• Sample the rank likelihood latent Gaussian variables \(\boldsymbol{Z}\) via Algorithm 2.

• Update the CUSP column-allocation variables for layer 2.

  For each column \(h=1,\dots,K^{(2)}_{\text{max}}\), let \(x_k^{(2),t}=B_{\cdot,h+1}^{(2),t}\) denote the \(h\)th non-intercept column of \(\boldsymbol{B}^{(2),t}\). Define \[\log p_{\mathrm{spike}}^{(2)}(x_k^{(2),t}) = -K^{(1)}_{\text{max}} \log(2\lambda_{1}^{(2)}) - \frac{\|x_k^{(2),t}\|_1}{\lambda_{1}^{(2)}},\] and \[\log p_{\mathrm{slab}}^{(2)}(x_k^{(2),t}) = \log 5 - K^{(1), t}_{\text{active}}\log 2 + \log\Gamma(1+K^{(1),t}_{\text{active}}) - \log\Gamma(1) - (1+K^{(1),t}_{\text{active}})\log\!\bigl(5+\|x_k^{(2),t}\|_1\bigr).\] This derived by marginalizing the Laplace slab density with rate \(\lambda^{(2)}_{0k}\), where \(\lambda^{(2)}_{0k}\sim\text{Gamma}(5,1)\). \(\Gamma(.)\) is the gamma function. Then for each allocation state \(k=1,\dots,K^{(2)}_{\text{max}}+1\), \[\log \Pr(c_k^{(2),t}=\ell \mid \boldsymbol{B}^{(2),t},\omega^{(2),t}) \propto \log \omega_k^{(2),t} + \begin{cases} \log p_{\mathrm{spike}}^{(2)}(x_k^{(2),t}), & \ell \le k,\\[0.3em] \log p_{\mathrm{slab}}^{(2)}(x_k^{(2),t}), & \ell > k. \end{cases}\] The implementation uses the MAP update \[c_k^{(2),t} = argmax_{\ell \in\{1,\dots,K^{(2)}_{\text{max}} + 1\}} \Pr(c_k^{(2),t}=\ell \mid \boldsymbol{B}^{(2),t},\omega^{(2),t}).\]

• Update the CUSP stick-breaking weights for layer 2. Let \[n_\ell^{(2),t} = \sum_{h=1}^{K_{\text{max}}^{(2)}}\mathbf{1}\{c_k^{(2),t}=\ell\}, \qquad n_{>\ell}^{(2),t} = \sum_{h=1}^{K_{\text{max}}^{(2)}}\mathbf{1}\{c_k^{(2),t}>\ell\}.\] Then for \(\ell=1,\dots,K^{(2)}_{\text{max}}\), \[v_\ell^{(2),t} \sim \mathrm{Beta}\!\bigl(1+n_\ell^{(2),t},\, \alpha^{(2)} + n_{>\ell}^{(2),t}\bigr), \qquad \alpha^{(2)} = K^{(2)}_{\text{max}},\] and set \(v_{K^{(2)}_{\text{max}} + 1}^{(2),t}=1\). The corresponding mixture weights are \[\omega_1^{(2),t}=v_1^{(2),t}, \qquad \omega_\ell^{(2),t} = v_\ell^{(2),t}\prod_{m=1}^{\ell-1}(1-v_m^{(2),t}), \quad \ell=2,\dots,K^{(2)}_{\text{max}}+1.\]

• Update the slab scales for layer 2. For each \(h=1,\dots,K^{(2)}_{\text{max}}\) such that \(c_k^{(2),t}>h\), sample \[\lambda_{1,k}^{(2),t} \sim \mathrm{Gamma}\!\left( 1+K^{(1),t}_{\text{active}},\, 5+\|\boldsymbol{B}_{\cdot,h+1}^{(2),t}\|_1 \right),\] where the second argument is the rate parameter. Then define the penalty matrix \[\lambda_{kh}^{(2),t} = \begin{cases} \lambda_{1}^{(2)}, & c_k^{(2),t}\le h,\\ \lambda_{1,k}^{(2),t}, & c_k^{(2),t}>h. \end{cases}\]

• Sample the CUSP column-allocation variables for layer 1. For each column \(k=1,\dots,K^{(1)}_{\text{max}}\), let \(x_k^{(1),t}=\boldsymbol{B}_{\cdot,h+1}^{(1),t}\). Define \[\log p_{\mathrm{spike}}^{(1)}(x_k^{(1),t}) = -J \log(2\lambda_{1}^{(1)}) - \frac{\|x_k^{(1),t}\|_1}{\lambda_{1}^{(1)}},\] and \[\log p_{\mathrm{slab}}^{(1)}(x_k^{(1),t}) = \log 5 - J\log 2 + \log\Gamma(1+J) - \log\Gamma(1) - (1+J)\log\!\bigl(5+\|x_k^{(1),t}\|_1\bigr).\] Then for \(\ell=1,\dots,K^{(1)}_{\text{max}}+1\), \[\log \Pr(c_k^{(1),t}=\ell \mid \boldsymbol{B}^{(1),t},\omega^{(1),t}) \propto \log \omega_\ell^{(1),t} + \begin{cases} \log p_{\mathrm{spike}}^{(1)}(x_k^{(1),t}), & \ell \le k,\\[0.3em] \log p_{\mathrm{slab}}^{(1)}(x_k^{(1),t}), & \ell > k. \end{cases}\] Again, the implementation uses the MAP choice \[c_k^{(1),t} = \arg\max_{\ell\in\{1,\dots,K^{(1)}_{\text{max}}+1\}} \Pr(c_k^{(1),t}=\ell \mid \boldsymbol{B}^{(1),t},\omega^{(1),t}).\]

• Update the CUSP stick-breaking weights and slab scales for layer 1. Let \[n_\ell^{(1),t} = \sum_{h=1}^{K_{\text{max}}^{(1)}}\mathbf{1}\{c_k^{(1),t}=\ell\}, \qquad n_{>\ell}^{(1),t} = \sum_{h=1}^{K_{\text{max}}^{(1)}}\mathbf{1}\{c_k^{(1),t}>\ell\}.\] Then for \(\ell=1,\dots,K^{(1)}_{\text{max}}\), \[v_\ell^{(1),t} \sim \mathrm{Beta}\!\bigl(1+n_\ell^{(1),t},\, \alpha^{(1)} + n_{>\ell}^{(1),t}\bigr), \qquad \alpha^{(1)} = K^{(1)}_{\text{max}},\] with \(v_{K^{(1)}_{\text{max}}+1}^{(1),t}=1\), and \[\omega_1^{(1),t}=v_1^{(1),t}, \qquad \omega_\ell^{(1),t} = v_\ell^{(1),t}\prod_{m=1}^{\ell-1}(1-v_m^{(1),t}), \quad \ell=2,\dots,K^{(1)}_{\text{max}}+1.\] For each active slab column \(h\) such that \(c_k^{(1),t}>h\), \[\lambda_{0,k}^{(1),t} \sim \mathrm{Gamma}\!\left( 1+J,\, 5+\|\boldsymbol{B}_{\cdot,k+1}^{(1),t}\|_1 \right),\] and \[\lambda_{jk}^{(1),t} = \begin{cases} \lambda_{1}^{(1)}, & c_k^{(1),t}\le k,\\ \lambda_{0,k}^{(1),t}, & c_k^{(1),t}>k. \end{cases}\]

Conditional M-step. Conditional on the Monte Carlo draws from the E-step, maximize the surrogate objective \(\widehat Q(\theta\mid \theta^{(t)})\) blockwise.

• Update the top-layer Bernoulli probabilities. The complete-data conditional log-likelihood for \(\pi=(\pi_1,\dots,\pi_{K_{\text{max}}}^{(2)})\) is \[\sum_{i=1}^N \sum_{l=1}^{K_{\text{max}}^{(2)}} \Bigl[ A_{il}^{(2)} \log \pi_l + (1-A_{il}^{(2)})\log(1-\pi_l) \Bigr].\] Replacing \(A_{il}^{(2)}\) by its Monte Carlo expectation \(q_{il}^{(t)}\) yields \[\pi_l^{(t+1)} = \arg\max_{\pi_l\in(0,1)} \sum_{i=1}^N \Bigl[ q_{il}^{(t)} \log \pi_l + (1-q_{il}^{(t)})\log(1-\pi_l) \Bigr] = \frac{1}{N}\sum_{i=1}^N q_{il}^{(t)}.\]

• Update the first-layer loading matrix \(B^{(1)}\). For each observed variable \(j=1,\dots,J\), define the Monte Carlo Gaussian objective for \(C\) simulations from the conditional posterior in the MC E-step as \[\ell_j^{(1)}(b) = \frac{1}{2\gamma_j^{(t)}} \sum_{i=1}^N \sum_{c=1}^C \left( Z_{ij}^{(t)} - \widetilde{\boldsymbol{A}}_i^{(1),(c),t} \boldsymbol{B}_{j.}^{(1),t^{\top}} \right)^2.\] The adaptive weighted \(\ell_1\) penalty is \[\mathcal{P}_j^{(1)}(b) = \sum_{k=1}^{K_{\text{max}}^{(1)}} |\boldsymbol{B}_{j k+1}| 1/\lambda_{jk}^{(1),t}.\] Thus the row update is \[\boldsymbol{B}_{j,\cdot}^{(1),t+1} = \arg\min_{b\in\mathbb{R}^{K^{(1)}_{\text{max}}+1}} \left\{ \ell_j^{(1)}(b) + \mathcal{P}_j^{(1)}(b) \right\}.\]

• Update the residual variances \(\gamma_j\). Given the updated row \(\boldsymbol{B}_{j,\cdot}^{(1),t+1}\), define \[\eta_{ij}^{(c),t+1} = \widetilde{A}_i^{(1),(c),t} \boldsymbol{B}_{j,\cdot}^{(1),t+1\top}.\] With \(\gamma_{j} \sim \mathrm{Gamma}(1,1)\), the conditional MAP update is \[\gamma_j^{(t+1)} = \frac{ 1 + \frac{1}{2} \sum_{i=1}^N \sum_{c=1}^C \bigl(Z_{ij}^{(t)}-\eta_{ij}^{(c),t+1}\bigr)^2 }{ 1 + \frac{1}{2} NC + 1 }.\]

• Update the second-layer loading matrix \(B^{(2)}\). For each hidden unit \(k=1,\dots,K^{(1)}_{\text{max}}\), define the Monte Carlo logistic objective \[\ell_k^{(2)}(b) = - \sum_{i=1}^N \sum_{c=1}^C \left[ A_{ik}^{(1),(c),t} \,\widetilde{A}_i^{(2),(c),t}B^{(2),t^{\top}}_{k.} - \log\!\left(1+\exp\bigl(\widetilde{A}_i^{(2),(c),t} \boldsymbol{B}_{k\cdot}^{(2),t\top}\bigr)\right) \right].\] The adaptive weighted \(\ell_1\) penalty is \[\mathcal{P}_k^{(2)}(b) = \sum_{l=1}^{K_{\text{max}}^{(2)}} |B^{(2)}_{kl+1}|1/\lambda_{kl}^{(2),t}.\] Hence \[\boldsymbol{B}_{k,\cdot}^{(2),t+1} = \arg\min_{b\in\mathbb{R}^{K^{(2)}_{\text{max}}+1}} \left\{ \ell_k^{(2)}(b) + \mathcal{P}_k^{(2)}(b) \right\},\] again solved numerically with fmincon.

• Threshold small coefficients. After the numerical updates, apply elementwise thresholding to the non-intercept coefficients: \[\boldsymbol{B}_{\cdot,2:(K^{(1)}_{\text{max}}+1)}^{(1),t+1} \leftarrow \operatorname{thres}\!\left( \boldsymbol{B}_{\cdot,2:(K^{(1)}_{\text{max}}+1)}^{(1),t+1},\, \xi \right),\] \[\boldsymbol{B}_{\cdot,2:(K^{(2)}_{\text{max}}+1)}^{(2),t+1} \leftarrow \operatorname{thres}\!\left( \boldsymbol{B}_{\cdot,2:(K^{(2)}_{\text{max}}+1)}^{(2),t+1},\, \xi \right).\] Here, the \(\operatorname{thres}(.)\) function maps any entry of the matrix less than \(\xi\) to zero, which speeds up convergence of the algorithm. We find this step necessary in order to minimize the noise injection from the Monte Carlo E-step. In practice, we tune \(\xi\) based on the predictive evaluation procedure outlined in Section 6. See Section 9.5 for suggested default values.

• Enforce hierarchical gating implied by the CSP activations in layer 1. Let \[\mathcal{I}_2^{(t)} = \{k : c_k^{(1),t}\le h\},\] Then, \(\boldsymbol{B}_{k,.}^{(2),t+1}=0 \quad \text{for } k\in\mathcal{I}_2^{(t)}.\)

• Update binary point estimates. After drawing \(C\) Monte Carlo replicates, the binary point estimates are updated by majority vote: \[A_{ik}^{(1),t+1} = \mathbf{1}\!\left( \frac{1}{C}\sum_{c=1}^C A_{ik}^{(1),(c),t} > \frac{1}{2} \right), \qquad A_{il}^{(2),t+1} = \mathbf{1}\!\left( \frac{1}{C}\sum_{c=1}^C A_{il}^{(2),(c),t} > \frac{1}{2} \right).\]

• Update the effective number of active nodes. The active widths of the two latent layers are \[K^{(1)^{*,t+1}_{\text{max}}} = \sum_{h=1}^{K_{\text{max}}^{(1)}}\mathbf{1}\{c_k^{(1),t}>h\}, \qquad K^{(2)^{*,t+1}}_{\text{max}} = \sum_{h=1}^{K_{\text{max}}^{(2)}}\mathbf{1}\{c_k^{(2),t}>h\}.\]

Convergence and temperature schedule. The algorithm monitors the relative Frobenius changes \[\frac{\|\boldsymbol{B}^{(1),t+1}-\boldsymbol{B}^{(1),t}\|_F}{\|\boldsymbol{B}^{(1),t}\|_F+\varepsilon}, \qquad \frac{\|\boldsymbol{B}^{(2),t+1}-\boldsymbol{B}^{(2),t}\|_F}{\|\boldsymbol{B}^{(2),t}\|_F+\varepsilon},\] and stops when their rolling-window variability is sufficiently small. Here \(\epsilon\) is a constant so that the quantity remains finite. After burn-in, the temperature is annealed upward according to \[\tau^{(t+1)} = \min\!\left\{1,\, \tau^{(t)} + 0.01\bigl((t-\mathrm{burn})-1\bigr) \right\}.\] Section 9.5 has recommendations default starting values.

9.3 Rank likelihood stochastic approximate EM↩︎

[9] introduce a stochastic approximate EM (SAEM) algorithm for estimation of parametric DDEs, which is adopted for DDE copula estimation in Section 5 (referred to as SAEM RL in the main text). We refer readers to their Algorithm 2 in the main text supplement S.3.2 for detail. To extend their algorithm to the maximum likelihood estimation under the rank likelihood DDE copula requires only one additional step. Namely, in addition to sampling binary latent variables \(\{\boldsymbol{A}^{(d)}\}_{d=1}^{D}\), we also sample rank-consistent \(\boldsymbol{Z}\) using Algorithm 2. Then, conditional on \(\boldsymbol{Z} \in \mathcal{R}(\boldsymbol{Y})\), the subsequent M-step, which maximizes a stochastically averaged objective, is identical to the M-step for a Gaussian DDE presented in [9]. The procedure is outlined in Algorithm 10.

In this case, the DDE copula parameters are \(\boldsymbol{\theta} = (\{\boldsymbol{B}^{(d)}\}, \boldsymbol{\gamma})\), since dimension is inferred through the sparsity pattern in the weights. We use the default hyperparameter selections for Normal DDEs recommended by [9], Supplement S.3.4, which include thresholding values limits for small entries in the weight matrices, as well as penalty parameters for the truncated LASSO penalty function [25] used in their optimization.

9.4 Addressing latent variable permutations for simulations↩︎

As mentioned in [9], resolving latent variable permutation is necessary to accurately compute errors between estimated weight matrices and data generating weight matrices. Given an estimated weight matrix \(\hat{\boldsymbol{B}}^{(d)}\), we resolve the permutation ambiguity by rearranging orderings of latent variables (and thus the rows and columns of each layer-specific weight matrix) and solving an assignment problem for each layer in a bottom-up fashion.

First, we construct a \(K^{(1)}_{\text{max}} \times K^{(1)}_{\text{max}}\) cost matrix, where each entry is the squared \(L^{2}\) norm between corresponding column vectors of \(\boldsymbol{B}^{(1)}\) and \(\hat{\boldsymbol{B}}^{(1)}\). Next, using the Hungarian algorithm of [36], we find the optimal permutation of the columns of \(\hat{\boldsymbol{B}}^{(1)}\) that minimizes assignment cost. In a recursive manner, we then permute the rows of \(\hat{\boldsymbol{ B}}^{(2)}\) based on the optimal permutation of the columns in \(\hat{\boldsymbol{B}}^{(1)}\), and solve the assignment problem Hungarian for the optimal permutation of the columns of \(\hat{\boldsymbol{B}}^{(2)}\). This procedure may be generalized for deeper models by solving the assignment problem for \(\hat{\boldsymbol{B}}^{(d)}\), permuting the rows of \(\hat{\boldsymbol{B}}^{(d+1)}\) based on the solution in the previous layer, and then again solving the assignment problem.

9.5 Hyperparameter selections↩︎

In our simulations and real data analysis, we find that there are four parameters that meaningfully effect DDE Copula estimation: the layer-specific rates of the spike-Laplace kernel in the CSP prior \(\{\lambda_{1}^{(d)}\}_{d=1}^{D}\), the temperature \(\tau\), and \(\xi\), the truncation limit for thresolding small coefficients.

• \(\boldsymbol{\{\lambda_{1}^{(d)}\}_{d=1}^{D}}\) For 2-layer DDE copulas, we generally specify \(\lambda_{1}^{(1)} <\lambda_{1}^{(2)}\), with the of each depending on the level of sparsity. For larger \(\lambda_{1}^{(d)}\), there is higher prior probability that a given column will come from the spike. Thus, we encourage more shrinkage in layer 2, which is affected by noise propagation from layer 1 latent variables. For \(J = 50\), we set \(\lambda_{1}^{(1)} = .02\) and \(\lambda_{1}^{(2)} = .04\), since there were fewer irrelevant columns in each weight matrix. For \(J = 100,150\), we set \(\lambda_{1}^{(1)} = .05\) and \(\lambda_{1}^{(2)} = .1\), as stronger shrinkage was required. For deeper DDE copula we recommend the ascending pattern for the penalty parameters for the first two layers, but suggest \(\lambda_{1}^{(d)} < .01, d>2\). We find that, while noise propagation through the network can create false positives in shallow layers, as the DDE becomes deeper, it is increasingly difficult to detect signal. In our deeper simulation (Section 10.3), we set \(\lambda_{1}^{(1)} = .01, \lambda_{1}^{(2)} = .03,\lambda_{1}^{(3)} = .005\)

• \(\boldsymbol{\tau}\) For the temperature, lower values will encourage more sparsity. The initial temperature should also balance the depth of the DDE. In our two-layer simulations we use \(\tau = .7\), which is subsequently raised in each iteration of the fitting algorithm. We note that we also apply tempering to the RL SAEM algorithm. For the deeper simulation in Section 10.3, the initial temperature of .7 prevented us from estimating any non-zero entries in \(\boldsymbol{B}^{(3)}\), so we used \(\tau = .9\) for those simulations.

• \(\boldsymbol{\xi}\) We follow the recommendations of [9] and use \(\xi = \max(.3, 3*N^{-.3})\) for all simulations and real data analyses

10 Additional Simulation Results↩︎

We complete the results presented in Section 5 by including the simulations for \(J = 50\) and \(J = 150\). We note largely consistent results to what is presented in the main text: RL CSP provides more accurate weight, sparsity structure, and dimension estimation, regardless of the dimension of \(\boldsymbol{Y}\). First we include the average estimate of \(\boldsymbol{B}^{(1)}\) across sample sizes and methods for \(J = 100\) in Figure 11, analogous to what we presented in Figure 8 in the main text.

We note that all rank-likelihood methods estimate \(\boldsymbol{B}^{(1)}\) accurately in increasing dimensions, but that the Poisson DDE performs substantially worse. Next we provide complete results for the remaining sample sizes.

10.1 J = 50↩︎

Analogous to the reuslts in the main text, in table Figure 6, we include the data generating weight matrices for the \(J = 50\) case. In Table 4, we provide dimension and sparsity structure estimation results, while in Figures 13-15 we provide weight matrix estimation statistics and visualizations.

10.2 J = 150↩︎

Analogous to the reuslts in the main text, in table Figure 16, we include the data generating weight matrices for the \(J = 50\) case. In Table 5, we provide dimension and sparsity structure estimation results, while in Figures 17-19we provide weight matrix estimation statistics and visualizations.

10.3 Deeper Simulations↩︎

We also ran simulations for deeper DDE Copulas (\(D = 3\)). We generated samples from a three layer DDE copula with \(J = 108\), \(K^{(1)} = 18\), \(K^{(2)} = 6\), and \(K^{(3)} = 2\). Data generating weight matrices are available in the right panel of Figure 21. We follow the pattern in Section 5 for the marginals \(F_{Y_{j}}\), alternating between irregular, discrete distributions. Like in Section 5, we generate 100 simulated data sets for each sample size and fit the DDE copula using Algorithm 3, first initializing the parameters based on the maximal allowable DDE network size of \(K^{(1)}_{\text{max}} = 36, K^{(2)}_{\text{max}} = 12, K^{(3)}_{\text{max}} = 4\). All models were fit under the same hyperparameter settings; \(\tau = .9\), \(\lambda_{1}^{(1)} = .01, \lambda_{1}^{(2)} = .03, \lambda_{1}^{(2)} = .005\).

MAP estimates for latent dimension, as well as the accuracy for the estimated sparsity structure, are available in Table 6. As the sample size increases, the latent dimension and sparsity structure is estimated with more accuracy. Similar to what was observed by [9], estimation accuracy is generally inferior for deeper layers, though convergence is apparent. This model is considerably more dense – marginally the induced mixture has \(2^{18+6+2} = 67108864\) components – and so larger sample sizes are required for accurate estimation.

We also include evaluations of the estimated weight matrices in Figures 21-20 via entry-wise MSE and the average MAP estimate of each layer-specific \(\boldsymbol{B}^{(d)}\), respectively. In Figure 20, we point out that the apparent lack of convergence for the deepest layer is due to the fact that, for simplicity, we fit all models under the same hyperparameter settings. In many cases, this resulted in the estimated dimension for the third dimension to either 1 or 4, which severely skews entry-wise MSE. In practice, one would conduct predictive evaluations as outlined in Section 6 to tune these parameter. In addition, at the two smallest sample sizes, graph recovery is uniformly poor, that the apparent better performance in these regimes is purely an artifact of the simulation design and sampling variability. After removing these outliers from simulations and focusing on sample sizes where the graph recovery begins to stabilize \((n \geq 2000)\), the transition toward consistent recovery is most clearly visible.

11 Comparisons to iVAE↩︎

We provide additional comparisons by analyzing the Big5 survey with identifiable variational auto-encoders (ivae; [13]). Motivated by our discovery of \(K^{(1)^{*}} = 9\) and \(K^{(2)^{*}} = 1\), we specify an iVAE with latent dimension \(K = 1\), while both encoder and decoder neural network architectures (two-layer perceptrons) have 9 hidden variables in each layer. The iVAE requires auxiliary information to construct “label" priors and we supply the model with the political ideology item used in Figure 9 in the main text. We utilize default settings provided by the authors for the algorithm.

We first compare the estimated encoder and decoder weights, each of which are \(J \times 9\), to the estimated \(\boldsymbol{B}^{(1)}\) under the DDE copula in Figure 22. The estimated \(\boldsymbol{B}^{(1)}\) is considerably more sparse than both the encoder and decoder weights. Next, we compute key items based on both the encoder and decoder weight matrices. We note that for the decoder weights, several columns were nearly entirely negative, and there were no items that surfaced based on the metric \(\max\{\min_{l \neq k} (\beta_{jk} - \beta_{jl}),\, 0\}\). The results are available in Table 9.

For the iVAE weights, \(A^{(1)}_{k}\) corresponds to the hidden variables in the respective perceptrons. To further facilitate comparisons, we include key items obtained from \(\boldsymbol{B}^{(1)}\) from the DDE copula, which are presented in Table 3 in the main text. Neither the encoder nor decoder weights capture structure in the data; key items are not meaningfully grouped across latent variables.

2pt

References↩︎

---

1. feldmanjr@wustl.edu↩︎

2. yuqi.gu@columbia.edu↩︎