Survival Fully-Collapsed Copula Mixed Membership Blockmodel

1 Introduction

Survival models in networks describe the time-evolving dynamics of interaction systems by treating each edge as a timeto-event outcome rather than a purely topological link. Unlike traditional graph approaches that focus on who connects to whom, these probabilistic frameworks explicitly model waiting times and event timing, offering a natural way to explain response delays, interaction durations, and future event occurrence [Perry and Wolfe, 2013; Trivedi et al., 2019;

Kumar et al., 2019]. Recent advances have extended classical survival analysis to network data via multivariate point processes and hazard-based formulations, enabling continuoustime modeling of social and biological systems [Sit et al., 2021; Rastelli and Fop, 2020]. Within this paradigm, interaction timing is often shaped by latent roles of the participating actors, where different role combinations can produce distinct temporal patterns.

However, existing network survival frameworks still suffer from two key limitations. First, many approaches implicitly assume that the censoring mechanism is independent of the underlying event process once observed covariates are accounted for. This assumption is difficult to justify in real interaction logs and can introduce substantial bias when censoring is informative [Gharari et al., 2023; Kalbfleisch and Prentice, 2002]. Second, point-process-based models on networks commonly emphasize the initiation dynamics of senders and treat receivers primarily as context, which makes it hard to capture the shared latent structure that jointly governs pairwise behavior [Perry and Wolfe, 2013; Trivedi et al., 2019; Song et al., 2024]. These issues become especially pronounced at scale, where observation windows create pervasive censoring and interaction behavior is highly heterogeneous [Rastelli and Fop, 2020; Wang and Wong, 1987].

Recent work has begun to address the first limitation by explicitly relaxing the standard independent-censoring assumption. In particular, Gharari et al. [2023] propose a copulabased deep survival framework that couples the event-time and censoring-time mechanisms through a survival copula, allowing their dependence to be learned from right-censored observations, where the event time is only known to exceed the observation window. This approach builds on the classical role of copulas in separating marginal behavior from dependence structure [Nelsen, 2006]. However, this line of work is primarily developed for independent instances with covariates, and it does not explicitly model network topology or the latent roles that organize interaction patterns across node pairs [Perry and Wolfe, 2013; Sit et al., 2021].

A separate line of research focuses on the network setting and addresses the sender-centric bias common in pointprocess formulations. Song et al. [2024] propose the Survival Mixed Membership Blockmodel (SMMB), which combines mixed-membership blockmodeling with survival analysis to model edge-level response processes through role-dependent hazards [Song et al., 2024; Airoldi et al., 2008]. This provides a structured way to account for heterogeneity across role pairs and to accommodate censoring in time-to-event network data [Kalbfleisch and Prentice, 2002; Song et al., 2024]. SMMB follows the standard mixed-membership assumption that sender and receiver roles are drawn independently from their membership vectors [Airoldi et al., 2008; Song et al., 2024]. However, in many social and biological networks, structural patterns such as reciprocity, hierarchy, and assortativity imply that interacting roles are coupled rather than independent [Hoff, 2005; Karrer and Newman, 2011]. As a result, important role correlations can be missed when the prior cannot represent such dependence directly.

Copulas also provide a direct way to introduce dependence into mixed-membership modeling without changing the marginal membership interpretation [Nelsen, 2006]. Building on this idea, Fan et al. [2016] introduce the CopulaMMSB (CMMSB), which uses a copula layer to model the joint distribution of sender and receiver roles and thereby relaxes the conditional independence built into standard MMSB-style priors [Fan et al., 2016]. The modeling gain is clear, but it comes with a practical cost: once a copula is introduced, the Dirichlet-multinomial conjugacy that makes role-assignment updates simple no longer holds, so inference becomes non-conjugate and more expensive. As a result, common strategies rely on augmentation or auxiliaryvariable constructions to restore tractable updates [Damlen et al., 1999; Neal, 2000]. In CMMSB settings this typically leads to semi-collapsed samplers that only partially marginalize latent parameters, leaving strongly coupled latent variables in the chain; in practice, this coupling can yield high autocorrelation and slow mixing as network size and model complexity grow [Fan et al., 2016; Liu, 1994].

To bridge these gaps, we propose the first Survival FullyCollapsed Copula Mixed Membership Blockmodel, which jointly models time-to-event dynamics on networks and complex dependencies among latent roles. By combining a survival blockmodel with an improved Copula-MMSB role mechanism, it captures dependent sender-receiver role pairs while preserving the mixed-membership interpretation. Another key implication of this construction is that it resolves the computational bottleneck of copula-based mixedmembership blockmodels. Such models are typically nonconjugate, and existing inference procedures are often only semi-collapsed, which leads to slow mixing and poor convergence in large networks. Under our formulation, however, we can analytically integrate out both the node membership variables and the copula latent variables, yielding an exact closedform discrete distribution over role pairs. This makes it possible to use a fully collapsed Gibbs sampler that draws role assignments directly from the marginal distribution and avoids inner-loop Metropolis-Hastings updates or auxiliary-variable constructions. Empirically, on the Enron benchmark our approach improves predictive likelihood and yields clearer roledependence structure than survival blockmodel baselines, and in simulations it achieves a 100× gain in sampling efficiency over semi-collapsed alternatives.

2 Model

We consider a directed interaction network with N nodes, denoted by U = {1, . . . ,N}. For each ordered pair (i, j) ∈ E, we observe an observed time tij and an event indicator δij , where δij = 1 indicates that an event occurred and δij = 0 indicates right censoring. Each edge may also have an associated covariate vector xij . To capture both temporal dynamics and structural dependence between interacting nodes, we define a generative model that couples latent sender and receiver roles through a copula while driving event times through a role-dependent survival process. Key notations are summarized in Table 1.

Generative process. The model is specified hierarchically as follows. Figure 1 shows the resulting fully collapsed dependency structure used by the Gibbs sampler.

1. Node memberships. For each node i, draw a mixedmembership vector

where πik represents the propensity of node i to participate in role k.

1. Coupled role assignments. For each edge (i, j), draw a pair of latent uniform variables

from a bivariate copula with dependence parameter θ. These uniforms are mapped to discrete sender and receiver roles via the inverse CDFs of the membership vectors:

where Πi(k) = ∑km=1 πim. This construction ensures the correct marginals, sij ∼ Cat(πi) and rij ∼ Cat(πj ), while inducing dependence between them through Cθ.

1. Time-to-event outcome. Conditional on the role pair (sij, rij) = (k, ℓ), the edge (i, j) follows a role-specific survival model. Let λkℓ g denote the baseline hazard for role pair (k, ℓ) in time segment g, and let βkℓ be the regression coefficients for covariates xij. The observed time tij and event indicator δij are then generated according to this role-dependent hazard.

Collapsed Role Distribution. To enable efficient inference, we follow the principle of collapsed Gibbs sampling [Griffiths and Steyvers, 2004] and analytically integrate out the continuous latent variables (π, uij, vij). Let (uij, vij) ∼ Cθ denote the latent copula variables for edge (i, j). Conditioned on the membership vectors πi and πj , the probability that this edge is assigned to role pair (k, ℓ) is given by

Geometrically, this probability represents the volume under the copula density. The integration region is a rectangle defined by the membership boundaries. This leads to the copula rectangle mass:

The core of our contribution is the full collapse. By integrating out the Dirichlet-distributed membership vectors π, we obtain the marginal likelihood of the role pair (k, ℓ) in closed form:

This quantity p(ij) kℓ acts as a "collapsed prior" for the role pair, incorporating both the community structure and the dependency θ.

Survival Likelihood. Conditional on the latent role pair (sij, rij) = (k, ℓ), edge (i, j) follows a piecewise-constant survival model. Let the time axis be partitioned into g segments indexed by g = 1, . . . ,G, and let exposureij,g denote the time that edge (i, j) is at risk in segment g. The baseline hazard for role pair (k, ℓ) in segment g is λkℓ g , giving a cumulative hazard contribution

Covariates enter multiplicatively through

where xij is the covariate vector and βkℓ is the role-pairspecific coefficient. The survival likelihood is then

Marginalizing over role pairs yields the collapsed edge likelihood

which provides the discrete weights for sampling (sij , rij) in the fully collapsed Gibbs sampler.

3 Model Identifiability

Identifiability is a prerequisite for interpretability in latent variable models: distinct parameter settings should not induce the same observed data distribution. In our setting, this means that the role-specific survival mechanisms and the copula-induced role dependence are uniquely determined from the observed pairwise survival data. We adopt the general identifiability framework for finite mixtures and latent variable models in Allman et al. [2009].

Recall that the collapsed edge likelihood in (5) can be written as a finite mixture over role pairs:

where Lkℓ surv is the role-pair-specific survival kernel defined in Section 2, G is the number of time segments in the piecewiseconstant hazard specification and pkℓ(θ, α, η) denotes the collapsed role-pair mass induced by the copula parameter θ and the Dirichlet prior Dir(αη).

Our identification strategy proceeds in two steps: (i) identify the set of mixture components {Lkℓ surv } and their mixing weights {pkℓ} from the observed distribution; (ii) recover (θ, α, η) from the identified weights.

Regularity conditions. We impose the following sufficient conditions [Teicher, 1963; Allman et al., 2009; Kalbfleisch and Prentice, 2002]:

• A1. The number of roles K is known and finite. The Dirichlet hyperparameters satisfy α > 0 and η ∈ {v ∈ RK :

A2. The mapping

is locally one-to-one on the parameter space. Equivalently, the Jacobian of {pkℓ} with respect to (θ, α, η) has full column rank in a neighborhood of the true parameter.

• A3. The time grid contains at least two segments, and the baseline hazards {λkℓ g } are not identical across all role pairs. An identifiability constraint is imposed on {βkℓ} to fix the overall scale.
• A4. The family of kernels {Lkℓ surv (·)}k,ℓ is linearly independent on the support of (t, δ, x) under the design.
• A5. Both events and right censoring occur with positive probability.

Main result. Let Θ = (θ, α, η, λ, β) collect the model parameters, where λ = {λkℓ g : k, ℓ = 1, . . . ,K; g = 1, . . . ,G} and β = {βkℓ : k, ℓ = 1, . . . ,K}.

Theorem 1. Under assumptions A1-A5, if two parameter sets Θ and ˜Θ induce the same observed distribution in (6) for all (t, δ, x), then they are equivalent up to a permutation σ of the role labels {1, . . . ,K}. In particular:

1. the survival parameters match up to relabeling, βkℓ = ˜ βσ(k)σ(ℓ) and λkℓ g = ˜λσ(k)σ(ℓ) g for all g = 1, . . . , G;
2. the mixing weights match up to the same relabeling, pkℓ(θ, α, η) = pσ(k)σ(ℓ)(˜ θ, ˜α, ˜η);
3. by A2, the hyperparameters governing the collapsed weights are locally identified, so (θ, α, η) = (˜ θ, ˜α, ˜η).

Proof sketch. By A3-A5, the component kernels {Lkℓ surv } form a distinct and linearly independent family (A4), so standard identifiability results for finite mixtures imply that the observed distribution uniquely determines the set of components and their mixing weights {pkℓ} up to a permutation of labels [Teicher, 1963; Allman et al., 2009]. Given the identified weights, A2 ensures that the system of equations linking (θ, α, η) to {pkℓ} admits a unique (local) solution, completing the argument.

For the FGM copula used in our experiments, A2 reduces to a concrete non-degeneracy requirement: the Jacobian of {pkℓ(θ, α, η)} with respect to (θ, α, η) should have full column rank almost everywhere. In practice, we fix the label permutation by ordering roles according to an identifiable summary of their survival behavior, such as the average baseline hazard [Stephens, 2000].

4 Inference

We establish a coherent path from a variational lower bound to a fully collapsed Gibbs sampler. The formulation clarifies the origin of the bound, following standard variational arguments based on Jensen's inequality and KL divergences [Bishop, 2006; Wainwright and Jordan, 2008; Blei et al., 2017]. It also shows how the closed form arises and why it improves mixing, in the spirit of collapsed sampling for latent discrete structures [Liu, 1994; Van Dyk and Park, 2008].

After analytically integrating out the node memberships (πi,πj) and the copula uniforms (uij, vij), each edge (i, j) induces a finite mixture over role pairs with collapsed masses

Here p(ij) kℓ defines a valid bivariate probability mass function over role pairs, satisfying ∑Kk=1 ∑Kℓ=1 p(ij) kℓ = 1.

Combining the collapsed role mechanism with the survival layer yields the marginal likelihood contribution of edge (i, j),

Consequently, the role pair (sij, rij) can be updated from a discrete distribution with unnormalized weights

which eliminates the auxiliary-variable or inner-loop Metropolis-Hastings steps required by semi-collapsed copula blockmodels.

To compute the collapsed masses in (7), we consider polynomial copulas admitting an expansion

Under this representation, the copula rectangle volume factorizes term-wise,

Since πi ∼ Dir(αη), the cumulative membership Πi(k) = ∑m≤k πim follows a Beta distribution,

Standard Dirichlet-Beta identities give the moments

where (x)(a) denotes the rising Pochhammer symbol. Substituting these expressions into (7) yields the closed form

For the FGM copula [Nelsen, 2006], whose expansion is C(u, v) = uv - θu2 v - θuv2 + θu2 v2 , this reduces to

To define a valid copula, the function C(u, v) must satisfy the standard boundary conditions and be 2-increasing on [0, 1]2 [Joe, 2014]. These constraints ensure that the induced copula rectangle volumes VCθ (·) define a proper bivariate distribution over role indices.

Although arbitrary polynomials do not in general satisfy these conditions, several well-established constructions yield polynomial expansions of the form (9) while preserving copula validity. In particular, Bernstein copulas provide a basis in which the copula constraints reduce to linear inequalities on the coefficients {γab(θ)} [Sancetta and Satchell, 2004; Choudhuri et al., 2004]. More generally, copulas defined through valid densities or spline representations admit equivalent polynomial or series expansions after integration [Joe, 2014]. Crucially, regardless of the construction, any copula admitting a separable expansion of the form (9) yields collapsed role-pair masses through (10). Hence, enforcing copula validity does not interfere with the closed-form marginalization underlying our fully collapsed sampler.

We construct a fully collapsed Markov chain Monte Carlo algorithm that integrates out the node memberships (πi,πj) and copula uniforms (uij , vij ), and updates role pairs and survival parameters directly from their marginal posteriors. The number of roles K, the number of time segments G, and the segment cutpoints are fixed.

At each iteration, the sampler proceeds as follows.

1. Update role pairs (sij, rij). For each edge (i, j), remove its current contribution from the global sufficient statistics and compute the collapsed role-pair masses p(ij) kℓ using (10) or (11). The joint posterior update is then

This joint update avoids the auxiliary-variable or inner-loop Metropolis-Hastings steps required by semicollapsed copula blockmodels.

1. Update baseline hazards {λkℓ g }. With independent Gamma priors λkℓ g ∼ Gamma(κλ, νλ) for all (k, ℓ, g), where κλ > 0 is the shape and νλ > 0 is the rate,

where

Here g(t) denotes the time-segment index induced by fixed cutpoints 0 = τ0 < τ1 < · · · < τG, and exposureij,g is the time-at-risk contributed by edge (i, j) within segment g.

1. Update regression coefficients {βkℓ}. Each βkℓ is updated using a Metropolis-Hastings step targeting its conditional posterior.
2. Update concentration α. We propose α on the log scale and accept or reject using the collapsed Dirichlet-multinomial marginal likelihood.
3. Update base weights η. We employ a mixture of local Dirichlet proposals and global independence moves.
4. Update copula parameter θ. For low-dimensional copula families, θ is updated by grid-based sampling or treated as fixed.

The survival component follows a piecewise-constant hazard model, which admits conjugate Gamma updates and is standard in Bayesian event-history analysis [Sinha and Dey, 1997]. By analytically integrating out the node memberships and copula latents, the sampler operates on a discrete marginal over role pairs, dramatically reducing posterior dependence and enabling rapid traversal of the multimodal state space [Gelman et al., 1997].

4.1 Computational Complexity

The dominant computational cost of the proposed algorithm arises from the joint update of the role pair (sij, rij) for each edge. For a network with E observed edges and K latent roles, evaluating the fully collapsed conditional distribution requires computing K2 role-pair weights p(ij) kℓ Lkℓ surv . As a result, the per-sweep complexity of the collapsed sampler is O(EK2 ). In contrast, the SMMB baseline updates sender and receiver roles through conditionals that scale linearly in K, leading to an O(EK) per-sweep cost. However, such semi-collapsed or uncollapsed updates typically retain stronger posterior dependence between latent variables, which can translate into slower mixing and higher autocorrelation in large networks. The fully collapsed sampler makes a deliberate trade-off: it replaces separate conditional updates by a single joint update over (k, ℓ), thereby reducing posterior coupling between sender and receiver roles at the level of each edge. This trade-off is favorable in the regime most common for blockmodels, where K is small (e.g., K ≤ 10) while E may be large. In this setting, the extra factor of K has limited impact on runtime, whereas the gain in mixing can be substantial.

As shown in Section 5.1, the fully collapsed sampler yields markedly lower autocorrelation, which reduces the number of sweeps required to reach a target Monte Carlo precision. Consequently, despite the higher per-sweep complexity, the total wall-clock time to effective convergence can be lower than that of semi-collapsed baselines. Moreover, the algorithm scales linearly in the number of edges E, making it suitable for large sparse networks. This effect is reflected in the large ESS/s gains reported in Section 5.1.

5 Experiments

5.1 Simulation: MCMC Sampling Efficiency

To evaluate the practical benefits of the proposed fully collapsed inference scheme, we compare its sampling efficiency with two semi-collapsed baselines under the FGM copula specification. The goal is to isolate the effect of marginalizing out both the membership vectors (πi, πj) and the copula uniforms (U, V ) on the mixing behavior of the Markov chain.

We consider three samplers that differ only in the set of latent variables that are analytically integrated out. The first two correspond to the partially collapsed schemes studied in Fan et al. [2016], while the third is the fully collapsed method proposed in this paper:

1. CMMSBπ (node-membership marginalization). The copula uniforms (U, V ) are integrated out, but the node membership vectors πi are sampled explicitly. Each πi is updated by Metropolis-Hastings based on the product of likelihood contributions from all edges incident to node i, after which (sij, rij) is sampled from the conditional rectangle probabilities p(ij) kℓ (πi,πj).
2. CMMSBuv (copula-uniform marginalization). The membership vectors (πi,πj) are integrated out, while the copula uniforms (U, V ) are sampled. Each pair (uij, vij) is updated by Metropolis-Hastings given (sij, rij), after which (sij, rij) is resampled from the induced prior distribution with leave-one-out counts.
3. CMMSBFull (full marginalization). Both (πi,πj) and (U, V ) are integrated out, and each role pair (sij, rij) is sampled directly from the closed-form discrete distribution p(ij) kℓ derived in Section 4.

All samplers are run under identical settings. We simulate networks with K = 4 latent roles and N = 36 nodes, yielding E = 1260 directed edges. The Dirichlet concentration is set to α0 = 1, and the FGM copula parameter is fixed at θ = 0.5. Each chain is run for 5000 sweeps with 500 burnin iterations, and three independent chains are generated for each method.

To assess mixing, we monitor the match rate

together with the log prior. Such low-dimensional summaries of the latent configuration are commonly used as diagnostics for convergence and mixing in latent variable models [Brooks and Gelman, 1998]. Sampling efficiency is quantified using the Integrated Autocorrelation Time (IACT) [Geyer, 1992], the Effective Sample Size (ESS) [Kass and Raftery, 1995], and the Monte Carlo Standard Error (MCSE) computed via batch means [Flegal et al., 2008; Jones et al., 2006].

Figure 2 highlights a substantial difference in mixing behavior. For CMMSBFull, the autocorrelation of the match rate drops essentially to zero after the first lag, indicating that successive draws are nearly independent. By contrast, both semicollapsed samplers display pronounced autocorrelation over many lags, reflecting slow exploration of the posterior.

These differences are quantified in Table 2. The fully collapsed sampler attains the largest effective sample size (ESS = 4336) and an IACT close to one (1.04), indicating near-ideal mixing. This improvement directly translates into higher estimation precision, with CMMSBFull achieving the smallest Monte Carlo standard error.

In terms of computational throughput, CMMSBFull executes 400.2 sweeps per second, compared with 3.1 and 1.8 for the two semi-collapsed schemes. When combined with the gains in mixing, this yields more than two orders of magnitude improvement in effective samples per second (ESS/s).

Taken together, these results show that analytically integrating out both membership vectors and copula coordinates removes the auxiliary dependencies that dominate the posterior geometry in semi-collapsed copula blockmodels. In this experiment, the resulting fully collapsed sampler attains more than a 102 -fold increase in effective samples per second relative to the two semi-collapsed baselines, reflecting both substantially faster mixing and lower per-sweep overhead.

5.2 Real Data Analysis: Enron Corpus

We apply the proposed model to the Enron e-mail corpus, a standard benchmark for temporal and relational network analysis [Klimt and Yang, 2004]. Following the preprocessing protocol of Song et al. [2024], each directed edge (i, j) is associated with an observed event time tij, a right-censoring indicator δij, and a vector of covariates xij. Event times are linked to latent role pairs through the piecewise-constant hazard specification described in Section 2.

We fix the number of roles to K = 4 and the number of time segments to G = 5, using the same cutpoints as in Song et al. [2024]. For the copula family, we use the FGM copula with parameter θ = 0.30, chosen to reflect moderate positive dependence. The fully collapsed CMMSBFull sampler is compared against the standard MMSB baseline in terms of both sampling behavior and statistical fit.

Convergence and Mixing. We first evaluate robustness of posterior exploration using short MCMC runs (1000/500/10), where the notation denotes total iterations, burn-in, and thinning. Table 3 reports standard multi-chain diagnostics. The MMSB baseline exhibits severe lack of convergence, with a Gelman-Rubin statistic ˆR = 4.39 [Gelman and Rubin, 1992], indicating that chains remain trapped in distinct posterior modes. In contrast, the collapsed CMMSBFull sampler achieves ˆR = 1.26, reflecting substantially improved mixing. Although CMMSBFull has a higher per-iteration cost, its effective sampling rate is nearly three times larger (median ESS/s 0.069 versus 0.024), showing that the cost of collapsing is more than offset by the reduction in autocorrelation.

Model Fit and Predictive Performance. Using longer chains (2400/1200/10), we evaluate model fit and out-ofsample predictive performance. Table 4 reports the Deviance Information Criterion (DIC) [Spiegelhalter et al., 2002] and the posterior-averaged log-likelihood on a held-out test set. The fully collapsed copula model attains a substantially lower DIC than MMSB, with a difference exceeding 3000, indicating a markedly better tradeoff between goodness of fit and model complexity. The hold-out log-likelihood improves by +0.0221 nats per edge, which accumulates to approximately 280 nats over the full test set. This improvement reflects the fact that allowing dependence between sender and receiver roles captures systematic structure that the independence assumption of MMSB cannot represent.

Posterior Structure. To examine how the copula affects posterior uncertainty, we summarize the latent role distributions using entropy, top-1 posterior mass, and pairconditional dependence measures, reported in Table 5. The MMSB baseline yields higher posterior entropy and lower top-1 mass, indicating greater uncertainty in role assignments. By contrast, CMMSBFull produces more concentrated posteriors, with lower entropy and a higher probability assigned to the dominant role. The copula also reduces pairconditional mutual information, showing that role assignments rely less on idiosyncratic pair-specific interactions and more on the global dependence structure captured by θ. This leads to a more stable and interpretable representation of the network.

6 Discussion

We proposed a fully collapsed copula mixed-membership blockmodel that integrates dependence between sender and receiver roles into a dynamic survival network while retaining exact inference. By analytically marginalizing both node memberships and copula coordinates, the resulting Gibbs sampler operates on a purely discrete state space, eliminating the auxiliary-variable and inner-loop sampling that typically dominate copula-based blockmodels. This structural simplification empirically demonstrates several strengths. Across both simulations and the Enron benchmark, full marginalization yields orders-of-magnitude gains in effective sampling rate and substantially improves posterior stability and predictive fit. These results indicate that much of the apparent uncertainty and poor mixing in standard blockmodels is induced by the independence assumption rather than by intrinsic ambiguity in the data. The current formulation assumes a fixed node set and a parametric copula family. Extending the fully collapsed construction to settings with evolving node populations or more flexible dependence structures is an important direction for future work. More broadly, the framework illustrates how non-conjugate dependence in dynamic networks can be handled through exact marginalization rather than approximation, providing a principled route to scalable Bayesian inference in complex relational models.