arXiv CSopen access

Proper Scoring Rules for Right-Censored Survival Data

Proper Scoring Rules for Right-Censored Survival Data

arXiv:2606.06393v1 [cs.LG] 4 Jun 2026

Jef Jonkers IDLab Department of Electronics and Information Systems Ghent University, Belgium [email protected]

Glenn Van Wallendael IDLab Department of Electronics and Information Systems Ghent University - imec, Belgium

Luc Duchateau Biometrics Research Group Department of Morphology, Imaging, Orthopedics, Rehabilitation and Nutrition Ghent University, Belgium

Sofie Van Hoecke IDLab Department of Electronics and Information Systems Ghent University - imec, Belgium

Abstract Proper scoring rules provide a rigorous theoretical basis for the training and evaluation of probabilistic forecasts. However, in the presence of right censoring, the event time is only partially observed, rendering conventional scoring rules inapplicable in their standard form. We propose a framework for proper scoring of right-censored survival outcomes based on a simple idea: first, map the predictive distribution through the censoring mechanism, then apply the underlying proper score on the induced observed-data law. This yields localized scores for fixed censoring times and marginalized scores when the censoring time is random or only partially observed. The resulting construction recovers familiar right-censored likelihood and IPCW-type criteria within a coherent framework, while also yielding right-censored versions of the CRPS, pinball loss, Brier score, and energy score. We show that the marginalized score is proper under conditional independent censoring and strictly proper on the identifiable region. The same principle also leads to censored engression, a sample-based learning objective for multivariate right-censored survival modeling. In experiments, our scores correctly rank the oracle forecast across several censoring regimes, whereas forecast-dependent plugin weighted scores can exhibit ranking reversals. Censored engression likewise substantially improves over naive training on censored outcomes.

1

Introduction

Predicting time-to-event outcomes is central in many applications, including healthcare, reliability, customer churn, and risk modeling. In these settings, the goal is often not merely to estimate the mean survival time or a hazard, but to forecast the full conditional event-time distribution FX given covariates X. Proper scoring rules [1] provide principled criteria for training and evaluating probabilistic forecasts, and they ensure that the true predictive distribution is optimal in expectation. The main difficulty in survival prediction is that event times are only partially observed. Under right censoring, one observes Y = min(T, C) and ∆ = 1{T ≤ C}, rather than the latent event time T itself. As a result, standard proper scoring rules for uncensored outcomes cannot be applied directly. Preprint.

This raises a basic question: how should one define a proper scoring rule when the event time is coarsened by right censoring? Our answer is to score forecasts after applying the same coarsening operation that generated the data. Given a predictive distribution F for T and a censoring time c, we evaluate the induced distribution of (min(T, c), 1{T ≤ c}) under T ∼ F , rather than evaluating F against the unobserved T . Thus, the score is determined by the observation mechanism: the censoring map is fixed independently of the forecast being evaluated. For fixed c, this gives a localized proper score in the sense of de Punder et al. [2]. For random censoring, the realized censoring time may be only partially observed; for example, if an event is observed, one only knows that C ≥ Y . We therefore average the localized score over the conditional law of C given (Y, ∆, X) under the censoring mechanism. The resulting marginalized score remains forecast-independent once the censoring law is fixed, or estimated separately from the forecast being evaluated, in contrast to plug-in constructions [3] where censoring weights or tail probabilities are recomputed from the forecast itself. Beyond evaluation, the same principle naturally leads to a training objective for flexible survival models. We extend engression [4], an implicit distributional regression method based on the energy score, to right-censored multivariate survival data. The resulting censored engression method trains a noise-conditioned generator by comparing censored generated samples with observed censored outcomes, using the same localized or marginalized right-censored energy score. This makes it possible to learn flexible joint event-time distributions under censoring without requiring a tractable likelihood. Contributions. First, we present a unified construction of proper scoring rules under right censoring by censoring the predictive distribution prior to scoring. Second, we show that the resulting marginalized score is proper under conditional independent censoring and strictly proper on the identifiable region. Third, we derive concrete right-censored versions of the logarithmic score, CRPS, pinball loss, Brier score, and energy score, clarifying their relationship to existing survival criteria (see Table 1). Fourth, we use the same principle to define censored engression for multivariate right-censored survival modeling. Finally, in experiments, we separate three empirical questions: in controlled settings with known data-generating processes, we show that our forecast-independent scores correctly rank the oracle forecast across censoring regimes, whereas forecast-dependent plug-in weighted scores can exhibit ranking reversals; in synthetic multivariate learning experiments, we demonstrate that censored engression outperforms naive training and exhibits competitive performance relative to likelihood-based Weibull, copula, and low-dimensional grid benchmarks trained using censored log scores. Furthermore, in an intensive care unit (ICU) acute kidney injury (AKI) case study, we show that censored engression can flexibly model joint time-to-event distributions for clinically relevant endpoints, including creatinine-based and urine-output-based AKI.

2

Related work

Our proper scoring rule framework under censoring connects and clarifies several strands of prior work. For right-censored survival evaluation, Graf et al. [5] introduced the inverse probability of censoring weighting (IPCW) Brier score and its integrated version, and Gerds and Schumacher [6] showed, implicitly, that, with a known or correctly specified censoring distribution, these criteria target the proper population objectives in regions where censoring positivity holds. In particular, administrative censoring is not problematic per se; the difficulty, as emphasized by Kvamme and Borgan [7], arises from a positivity issue when evaluation extends beyond the identifiable follow-up region. Likewise, the non-propriety results in Rindt et al. [8] should be understood as applying to IPCW implementations using marginal, rather than covariate-conditional, censoring weights; under correct conditional censoring weights, the pathology disappears. The same work also shows that the usual right-censored negative log-likelihood is proper, and that our construction recovers this score as a special case. Beyond likelihood and IPCW-type criteria, Avati et al. [9] proposed an adapted continuous ranked probability score (CRPS) for censored outcomes, but Rindt et al. [8] gave counterexamples showing that this score is not proper. More recently, Yanagisawa [3] derived proper extensions of discrete pinball, Brier, and ranked probability scores under censoring, but their weights depend on unknown functionals of the latent event-time law. When such weights are replaced by plug-in estimates from the forecast itself, the resulting score becomes forecast-dependent. Our work instead keeps the weighting forecast-independent by basing it solely on the censoring mechanism. 2

3

Proper scoring under censoring

We propose a general framework for proper scoring rules for right-censored survival outcomes that recovers, up to a forecast-independent positive factor, previously proposed proper scoring rules, such as the censored logarithmic score [8] and the IPCW Brier score [5, 6]. Our starting point is a strictly proper scoring rule on the latent outcome space. Following de Punder et al. [2], we view censoring as a transformation of the predictive distribution itself: rather than scoring the latent law directly, we first map it to the censored outcome space and then apply the corresponding proper score to the induced observed-data law. This preserves the original score divergence after censoring and yields a score that is locally proper for the latent law on the part of the outcome space that remains identifiable under the censoring mechanism. Let X ∈ X ⊆ Rd denote covariates, and let T ∈ T ⊆ Rk+ denote the latent event-time outcome, with k = 1 in the univariate case and k ≥ 2 in the multivariate setting. For each covariate value x, we seek a predictive distribution F = L(T | X = x). Under right censoring, however, we observe only Y = min(T, C) and ∆ = 1{T ≤ C}. Equip T with a σ-algebra T , let P be a class of probability measures on (T , T ), and let S : P × T → R̄ be a scoring rule with divergence DS (P ∥F ) = EP [S(F, T )] − EP [S(P, T )]. We call S strictly proper if DS (P ∥F ) ≥ 0 for all P, F ∈ P, with equality only when P = F . This is the uncensored benchmark to which we compare the censored scores constructed below. 3.1

Localized proper scoring rules under fixed censoring

We begin with a fixed censoring configuration. Let A ⊆ T denote the part of the outcome space that remains identifiable under this censoring mechanism, and let ∗ denote an abstract censored outcome. The corresponding censored outcome space is TA♭ = A ∪ {∗}. For a latent predictive law F ∈ P, define its censored version by FA♭ = FA + F̄A δ∗ ,

FA (B) := F (B ∩ A),

F̄A := F (Ac ).

Thus, the law is left unchanged on the identifiable region A, while all mass outside A is collapsed into the censored outcome ∗. The corresponding localized score is obtained by applying the proper score on the censored outcome space,  t, t ∈ A, ♭ SA (F, t) := S(FA♭ , t♭A ), where t♭A = ∗, t ∈ Ac . Intuitively, following the localization perspective of de Punder et al. [2], the forecast is first censored in the same way as the observation, and the original proper score is then evaluated on the resulting censored outcome space. By construction, the induced divergence is the original divergence applied to the censored laws, DS ♭ (P ∥F ) = DS (PA♭ ∥FA♭ ), so equality of population score implies equality A of the induced censored distributions, not necessarily of the full latent laws. In this sense, the score is strict only to the extent that the latent law remains identifiable under the fixed censoring mechanism. For fixed right censoring at time c, the identifiable region depends on the boundary convention. With the convention ∆ = 1{T ≤ c}, the identifiable region is Ac = [0, c], and the censored outcome map is  t, t ≤ c, ψc♭ (t) = ∗, t > c. The corresponding censored law is Fc♭ = F (· ∩ [0, c]) + F ((c, ∞))δ∗ . Thus, an event at the boundary c is distinguished from a censored observation beyond c. We write Sc♭ (F ; Y, ∆) for the localized score induced by this fixed censoring time. When T | X is continuous, boundary events have probability zero, so the convention at T = c does not affect expected scores. For distance-based scores such as the CRPS and energy score, however, 3

the abstract censored state ∗ cannot be inserted directly into Euclidean distances. By a slight abuse of notation, we keep the notation ψc♭ for this numerical encoding of the same coarsened observation: ψc♭ (t) = min(t, c),

Fc♭ = L(ψc♭ (T )).

This representation should be read as an encoding of the coarsened observation: the value c represents the censored region (c, ∞), not an event known to occur at c. When the event indicator is needed, we keep the convention ∆ = 1{T ≤ c}, so that (c, 0) denotes censoring beyond c whereas (c, 1) denotes an event observed at the boundary. In the multivariate setting, with T = (T1 , . . . , Tk ) and componentwise censoring thresholds c = (c1 , . . . , ck ), we use ψc♭ (t) = (min(t1 , c1 ), . . . , min(tk , ck )),

Fc♭ = L(ψc♭ (T )).

The score is then applied to the joint law of the censored vector, rather than to separate marginal laws, so dependence that remains visible after censoring is retained. For a shared censoring time, this corresponds to the special case c1 = · · · = ck = c. 3.2

Marginalized proper scoring rules

We now turn to random right censoring. Suppose we observe Y = min(T, C) and ∆ = 1{T ≤ C}, and let H(t | x) = P(C ≤ t | X = x) denote the conditional censoring distribution function. When C is random, the censoring time is not always fully determined by the observed data. Let Π(dc | Y, ∆, X) := L(C ∈ dc | Y, ∆, X) be the conditional law of the censoring time given the observed censored outcome, induced by the true censoring distribution H(· | X) and the right-censoring mechanism. Under conditional independent censoring, T ⊥ C | X, the observed-data score is obtained by averaging the fixed-c localized score over this conditional law: Z S̄(F ; Y, ∆, X) := Sc♭ (F ; Y, ∆) Π(dc | Y, ∆, X). In the univariate case, this reduces to Z  Sc♭ (F ; Y, 1) H(dc | C ≥ Y, X), ∆ = 1, S̄(F ; Y, ∆, X) = [Y,∞)  ♭ SY (F ; Y, 0), ∆ = 0. For a shared censoring time in the multivariate case, C is observed exactly whenever at least one coordinate is censored. If all coordinates are observed, one only knows that C ≥ maxj Yj , so Π is the law of C conditional on C ≥ maxj Yj . For fixed x, let S̄x (F ) := EP [S̄(F ; Y, ∆, X) | X = x] denote the conditional expected score. For a shared scalar censoring time, the identifiable region is     Ixshared = t ∈ Rk+ : P C ≥ max tj | X = x > 0 . j=1,...,k

For componentwise censoring times C = (C1 , . . . , Ck ), the analogous region is  Ixcomp = t ∈ Rk+ : P(Cj ≥ tj for all j = 1, . . . , k | X = x) > 0 . We write Ix for the relevant region under the censoring mechanism considered. In the univariate case, both definitions reduce to Ix = {t ≥ 0 : P(C ≥ t | X = x) > 0}. The next proposition states that the marginalized score is minimized by the true latent law, with uniqueness only on Ix . Proposition 1 (Marginalized propriety under right censoring). Assume conditional independent censoring, T ⊥ C | X, and let H(· | X) denote the true conditional distribution of the censoring time. Suppose that the conditional law of C given (Y, ∆, X) used in the marginalized score is the one induced by H(· | X) and the right-censoring mechanism. If, for H(· | X = x)-almost every fixed censoring time c, the localized score Sc♭ is proper for the induced censored law, then for every x, S̄x (F ) ≥ S̄x (P ) 4

Table 1: Examples of uncensored, fixed-censoring localized, and marginalized scores under right censoring. All scores are written as losses, so lower values are better. We write G(t | x) = P(C > t | X = x) and G(t− | x) = P(C ≥ t | X = x). Type

Score

Logarithmic score (log) Uncensored Slog (F ; t, x) = − log f (t | x) Fixed local

♭ Slog,c (F ; y, δ, x) = −δ log f (y | x) − (1 − δ) log(1 − F (c | x))

Marginalized

S̄log (F ; y, δ, x) = −δ log f (y | x) − (1 − δ) log(1 − F (y | x))

Continuous ranked probability score (CRPS) Rt R∞ Uncensored SCRPS (F ; t, x) = 0 F (s | x)2 ds + t (1 − F (s | x))2 ds R Rc y ♭ Fixed local SCRPS,c (F ; y, δ, x) = 0 F (s | x)2 ds + δ y (1 − F (s | x))2 ds Ry R ∞ G(s|x) 2 Marginalized S̄CRPS (F ; y, δ, x) = 0 F (s | x)2 ds + δ y G(y − |x) (1 − F (s | x)) ds Brier score at horizon τ (BSτ )
2 Uncensored SBS,τ (F ; t, x) = F (τ | x) − 1{t ≤ τ }
2

Fixed local

♭ SBS,τ,c (F ; t, x) = Fc♭ (τ | x) − 1{min(t, c) ≤ τ }

Marginalized

G(τ |X) 2 S̄BS,τ (F ; Y, ∆, X) = 1{Y > τ }F (τ | X)2 + ∆1{Y ≤ τ } G(Y − |X) (1 − F (τ | X))

Pinball / quantile score at level α (Qα )

Uncensored SQ,α (F ; t, x) = α − 1{t < qα (F | x)} t − qα (F | x)

♭ (F ; t, x) = α − 1{min(t, c) < qα (Fc♭ | x)} min(t, c) − qα (Fc♭ | x) Fixed local SQ,α,c S̄Q,α (F ; Y, ∆, X) = α(Y − qα (F | X))+ Marginalized

+∆

1−α 1{qα (F | X) > Y } G(Y − | X)

Z qα (F |X) G(t | X) dt Y

Energy score (ES) Uncensored SES,β (F ; t, x) = E[∥Z − t∥β | X = x] − 21 E[∥Z − Z ′ ∥β | X = x] Fixed local Marginalized

♭ SES,β,c (F ; y obs , x) = E[∥ψc (Z) − y obs ∥β | X = x] − 12 E[∥ψc (Z) − ψc (Z ′ )∥β | X = x] h i ♭ S̄ES,β (F ; Y obs , X) = E SES,β,C (F ; Y obs , X) | Y obs , X

for all forecasts F . If, in addition, the localized scores are strictly proper and the collection of induced censored laws for censoring times with positive probability identifies the latent law on Ix , then S̄x (F ) = S̄x (P )

⇐⇒

Fx = Px on Ix .

Hence, the marginalized score is proper as an observed-data score for forecasts of the latent law, and strictly proper only on the identifiable region. Outside Ix , two latent distributions can differ while inducing the same observed censored law. Hence, no scoring rule based only on (Y, ∆, X) can be strictly proper for those non-identifiable tail differences. In the univariate continuous case, if the base score is the logarithmic score and the forecast F has density f , then the marginalized score recovers the usual right-censored negative log-likelihood, see Section B.1 and Table 1. For the logarithmic score, marginalization simplifies because the fixed-c localized contribution is either observed exactly or constant over the compatible censoring values. Table 1 summarizes several scores obtained from the localized and marginalized constructions under right censoring. The formulas are written under the continuous no-boundary-ambiguity convention, so that the abstract censored outcome may be represented by the boundary value c through the boundary encoding ψc♭ (t) = min(t, c). Detailed derivations are given in Appendix B. 5

3.3

Special cases: administrative censoring and observed censoring times

Administrative censoring is the degenerate case in which the censoring time is fixed by design. Conditional on X = x, the censoring configuration is known, say c(x), and the observed-data score ♭ reduces to the fixed-censoring localized score S̄(F ; Y, ∆, X = x) = Sc(x) (F ; Y, ∆). Equivalently, the censoring law is Hx (dc) = δc(x) (dc). The score is strictly proper for the induced censored law, or for the latent law, only on the administratively identifiable region; it cannot distinguish tail differences beyond the follow-up horizon. A second useful case occurs when the realized censoring time C is observed for every individual, including those who experience the event before censoring. Then, no averaging over compatible ♭ censoring times is needed: for each observation, one evaluates S̄(F ; Y, ∆, C, X) = SC (F ; Y, ∆). Under conditional independent censoring, the resulting expected score is proper on the identifiable region. This setting is relevant when the censoring event, such as discharge or death in an ICU study (Section 5.3) is recorded even if the clinical event of interest occurs earlier. In such applications, the assumption that this censoring process is conditionally independent of the event time still remains.

4

Censored engression for multivariate survival

We now introduce censored engression, a sample-based method for learning flexible, possibly multivariate, event-time distributions from censored observations. The method combines two ideas. First, following engression [4], we represent the conditional event-time distribution via an implicit generator gθ (x, ε) that maps covariates and exogenous noise to samples from the predictive distribution. This makes the approach suitable for flexible joint survival modeling, including settings in which the dependence among multiple event times is difficult to parameterize explicitly. Second, following the localization-by-censoring principle of Section 3.1, we do not score the latent predictive distribution directly when event times are censored. Instead, we apply the censoring map to both generated samples and observations, and train the generator using a localized energy score, marginalized over compatible censoring configurations when necessary. Engression [4] represents the conditional law through a generator Yg (x) = g(x, ε),

ε ∼ Pε ,

(1)

where Pε is a simple reference distribution. The generator induces the conditional predictive law Pg (· | x) := L(g(x, ε)), so unlike classical regression, engression directly models the full conditional distribution via samples. For two independent noise variables ε, ε′ , the energy score population objective is   1 RES (g) = E ∥Y − g(X, ε)∥ − ∥g(X, ε) − g(X, ε′ )∥ , 2 where (X, Y ) ∼ P . Under correct specification, minimizing this objective recovers the true conditional law; see Proposition 1 of Shen and Meinshausen [4]. 4.1

Censored engression

We now extend engression to multivariate right-censored survival outcomes. Let T = (T1 , . . . , Tk ) ∈ Rk+ , k ≥ 1, denote the latent event-time vector, and let the generator in (1) induce the latent predictive law Fg,x := L(g(x, ε)). Under right censoring, the event times
are not fully observed. For a fixed censoring time c, let ψc♭ (t) = min(t1 , c), . . . , min(tk , c) denote the shared censoring map, and let Y = ψc♭ (T ) denote the corresponding censored outcome. Following the localization-by-censoring principle, we apply the energy score after censoring, that is, to the induced censored law L(ψc♭ (g(x, ε))) rather than to the latent law Fg,x itself. The resulting fixed-c censored energy loss is 1 Rlocal-cens-ES (g; Y, x, c) = Eε [∥ψc (g(x, ε)) − Y ∥] − Eε,ε′ [∥ψc (g(x, ε)) − ψc (g(x, ε′ ))∥] , (2) 2 iid

where ε, ε′ ∼ Pε . 6

When the censoring time is not fully observed, we average the fixed-c loss over the compatible conditional law of C given the observed data: Rcens-ES (g) := E[Rlocal-cens-ES (g; Y, X, C) | Y, ∆, X] .

(3)

Under conditional independent censoring, this is exactly the marginalized right-censored energy score objective introduced above. When there is no censoring, ψc is the identity map and the loss reduces to the ordinary engression energy score. With censoring, the model is instead trained to match the law of the censored outcome, and therefore recovers the latent event-time law only to the extent that it remains identifiable under the censoring mechanism. We refer to this procedure as censored engression. We also present this approach with a population optimality guarantee; see Proposition 2.

5

Experiments

5.1

Proper scoring rules for evaluating forecasts

We evaluate whether the proposed scores rank forecasts correctly under censoring. Because this oracle-ranking property concerns the unknown latent law of T | X, it cannot be verified from a real censored dataset. We therefore use controlled DGP-known simulations to validate scoring rules. We compare our localized and marginalized scores with the forecast-dependent weighted scores of Yanagisawa [3]. For the latter, we report both an oracle-weight version, where weights are computed from the true event-time law F0 , and a plug-in version, where weights are computed from the candidate forecast itself. The simulation design is given in Appendix C.2. Results and discussion. Table 2 reports the rank of the oracle forecast under each scoring rule or score family, for four censoring Table 2: Rank of the oracle foreregimes, with rank 1 indicating the best mean score. The corre- cast. Rank 1 means best mean score. Regimes A–C compare sponding raw mean score values are reported in Appendix C.3. F0 against F1 , . . . , F4 ; regime The latent gold-standard scores rank the oracle first in every regime, (D) as expected. The same is true for all of our localized and marginal- D compares F0 against four ized scores: across regimes A–D, our CRPS, pinball, Brier, and exploit forecasts F5 (ϵ). censored log scores all assign rank 1 to the oracle forecast. Thus, when evaluation is based on the censoring mechanism rather than on the forecast being scored, the oracle forecast is consistently preferred.

Scoring rule

A

B

C

D

CRPS, latent Pinball, latent Brier, latent NLL, latent

1 1 1 1

1 1 1 1

1 1 1 1

1 1 1 1

NLL [3], plug-in 2 1 2 2 The oracle-weight scores of Yanagisawa [3] also rank F0 first in Brier MC [3], plug-in 1 1 1 5 Brier [3], plug-in 1 1 2 5 all regimes, showing that the underlying scores behave as intended Pinball [3], plug-in 1 1 1 1 when the weights are fixed at their correct oracle values. The RPS [3], plug-in 2 1 2 5 failures arise in the practically relevant plug-in version, where the Censored NLL, ours 1 1 1 1 weights depend on the forecast being evaluated. The plug-in NLL CRPS, ours 1 1 1 1 Pinball, ours 1 1 1 1 ranks the oracle second in regimes A, C, and D; the plug-in RPS Brier, ours 1 1 1 1 ranks the oracle second in regimes A and C and last in regime D; and the plug-in Brier scores also fail in regimes C or D. The raw values in Appendix C.3 show that these ranking reversals are caused by exploit forecasts obtaining lower plug-in scores than the oracle, whereas our scores keep assigning lower scores to the oracle throughout.

A separate practical issue is estimation of the censoring distribution. The sensitivity results in Appendix C.4 show that when G(· | X) is replaced by pooled estimates, the oracle ranking is preserved in regime B and remains largely stable in regime C, with only a mild degradation for the discrete Brier score under a misspecified pooled Weibull censoring model. This reflects ordinary censoring-model misspecification rather than forecast-dependent circularity. 5.2

Synthetic multivariate censored-engression experiment

We next study censored engression as a learning method for multivariate right-censored event times. We consider synthetic survival data with Gaussian covariates, a shared censoring time, and a multivariate event-time vector T = (T1 , . . . , Tk ). The main experiments use a covariate-dependent 7

mixture log-normal latent DGP, which induces multimodal conditional event-time distributions and is intentionally misspecified for the Weibull and copula baselines. We vary the dimension k ∈ {2, 3, 5, 10} to assess scalability. Full details are given in Appendix C.5. We compare censored engression with oracle references, naive observed engression, and likelihoodbased baselines. DGP denotes the true data-generating conditional law and is therefore an oracle reference, not a fitted method. Latent engression is an infeasible benchmark trained directly on the latent event times T . The feasible methods are trained only on (X, Y, ∆): naive observed engression treats the censored vector Y as fully observed, censored engression uses the localized (2) or marginalized censored energy score (3), and the likelihood-based baselines use the censored logarithmic score (see Table 1). The likelihood baselines include independent Weibull, conditional MLP Weibull, Clayton copula Weibull, and conditional MLP Clayton copula Weibull models. For k = 2, we additionally include bivariate Gumbel copula and joint discrete-grid likelihood baselines. Results and discussion. We report two criteria. The primary criterion is the censored energy score on held-out censored data, which is available for all methods and corresponds to the observed-data prediction task. We also report the latent energy score against the unobserved true event-time vector T , which is available only because the experiment is synthetic. A selected comparison of censored engression against it naive counterpart is depicted in Table 3, complete results of all settings and benchmarks are presented in Tables 7–10. The experiments show a consistent pattern across censoring regimes and dimensions. Censored engression obtains the lowest censored energy score among feasible methods in all settings. The improvement over naive observed engression is substantial, especially under random and covariate-dependent censoring. For example, with k = 2 under uniform censoring, naive observed engression obtains a censored energy score of 0.5524, compared with 0.4904 for censored engression. With k = 10, the corresponding scores are 1.6402 and 1.3382. Thus, treating the censored vector as fully observed leads to an increasingly poor observed-data predictive distribution as the multivariate problem becomes more challenging. For k = 10, the censored engression approach performs even better than the latent engression approach on censored ES. The likelihood-based baselines are stronger competitors than naive observed training because they are fitted with the proper censored log score. Among them, the conditional MLP Clayton copula baseline is generally strongest, confirming that flexible covariate-dependent marginals and dependence modeling are important. Nevertheless, under the mixture log-normal DGP, censored engression consistently achieves lower censored energy score. This supports the main motivation for using a sample-based distributional model: the mixture DGP is multimodal and lies outside the Weibullcopula family, whereas censored engression can represent more flexible joint predictive distributions. The latent energy score gives a complementary diagnostic rather than the main evaluation target. It measures recovery of the full latent event-time vector, including parts of the distribution that are only partially identified from censored observations. As expected, the DGP and latent engression provide oracle references. Among feasible methods, likelihood-based copula models can sometimes obtain lower latent energy score than censored engression, especially for larger k, even when they are worse on the censored energy score. This is not contradictory: censored engression is trained to optimize the proper observed-data criterion, while latent recovery beyond the identifiable region is not fully determined by the censored data. The main empirical conclusion is therefore that censored engression gives the best observed-data predictive performance among feasible methods while remaining applicable as the event-time dimension increases. As a sanity check, Appendix C.5 repeats the k = 2 experiment under the original unimodal lognormal DGP. This setting is comparatively favorable to the likelihood-based Weibull and copula baselines. As expected, these baselines are highly competitive, and in some cases slightly better on latent energy score, while censored engression remains close on the observed censored energy score and continues to improve substantially over naive observed engression. 5.3

Illustrative clinical use case: early AKI forecasting in the ICU

We use early acute kidney injury (AKI) forecasting in the intensive care unit (ICU) as a realistic illustration of censored engression on multivariate time-to-event data. At each prediction time, the goal is to forecast the joint distribution of the remaining time to two clinically meaningful endpoints: 8

Table 3: Selected results for censored engression. Lower scores are better. Synthetic rows report both observed censored energy score and latent energy score under uniform censoring for the mixture log-normal DGP. DGP is the true data-generating law and latent engression is an infeasible reference trained on uncensored event times. AKI rows report localized energy scores; no latent oracle is available in the real-data use case. Bold indicates the best feasible observed-data result among the displayed methods. Setting

Score

Synthetic, k = 2, uniform

Censored ES Latent ES Synthetic, k = 10, uniform Censored ES Latent ES

AKI, joint AKI, creatinine AKI, urine-output

Loc. ES Loc. ES Loc. ES

DGP

Latent engression Naive obs. engression Censored engression

0.4840 1.2468 1.3131 9.4613

0.4907 1.2601 1.6082 10.0111

0.5524 1.4440 1.6402 16.5122

0.4904 1.2887 1.3382 14.2576

– – –

– – –

43.442 23.591 32.981

30.933 15.439 23.769

creatinine-defined AKI and urine-output-defined AKI. The predictors are routinely available variables, including demographics, baseline renal function, elapsed ICU time, and summary features of recent creatinine and urine-output measurements over the preceding 72 hours. The analysis uses deidentified ICU data from MIMIC-IV [10]; additional cohort, endpoint, and evaluation details are given in Appendix C.6. The observed target is right-censored by ICU discharge or death. This exit time is observed in the database, so the localized version of our scoring framework can be used directly: for each patient and prediction time, generated event times and observed outcomes are both mapped through the same censoring time before the score is evaluated. We report localized energy scores for the joint bivariate outcome and for each marginal endpoint, together with localized Brier scores for the composite endpoint “either AKI” at clinically relevant early-warning horizons. The comparison in Table 3 isolates the effect of the censoring-aware objective by comparing censored engression with a naive observed engression baseline that uses the same type of generator but treats the censored time Y = min(T, C) as if it were the true event time. Censored engression improves all selected localized energy scores in Table 3: the joint score decreases from 43.442 to 30.993, the creatinine-AKI marginal score from 23.591 to 15.439, and the urineoutput-AKI marginal score from 32.981 to 23.769. The improvement is therefore not driven by only one endpoint; accounting for censoring improves both the joint predictive distribution and each marginal event-time distribution. Fixed-horizon results show the same direction. For the composite endpoint, the localized Brier score decreases from 0.0257 to 0.0212 at 6 hours, from 0.0429 to 0.0398 at 12 hours, and from 0.0730 to 0.0714 at 24 hours. We also computed threshold-based enrichment metrics for short-term risk stratification, such as positive predictive value and sensitivity among the patients with the highest predicted risk of either AKI. These metrics are useful for interpretation but are not the primary evidence: unlike the localized scores, they depend on the chosen risk cutoff and are computed only among patients whose horizonspecific outcome is evaluable. With this caveat, they are directionally consistent with the proper-score results. For example, at 12 hours, the top-5% positive predictive value for the composite endpoint increases from 0.727 to 0.782, and at 24 hours it increases from 0.759 to 0.840; sensitivity at the same cutoff also increases from 0.316 to 0.340 and from 0.167 to 0.185, respectively. Appendix C.6 compares censored engression with standard survival baselines. Censored engression obtains the best localized Brier score for the composite endpoint across the reported horizons (Table 13) and the best joint and marginal localized energy scores among the methods that generate full event-time distributions (Table 14). These results reinforce the main conclusion of the use case: explicitly modeling the censoring mechanism improves observed-data predictive performance relative to both naive training on censored outcomes and several likelihood-based survival baselines.

6

Conclusion

This work argues that censored probabilistic evaluation should be built on forecast-independent censoring adjustment. IPCW scores target the latent event-time distribution only when the weights are 9

determined by the censoring mechanism rather than by the forecast being evaluated. This preserves the oracle ranking across a range of censoring regimes, while forecast-dependent plug-in scores change the evaluation target and can induce ranking reversals. The same methodological principle yields a training objective for censored generative survival modeling. Censored engression uses the censored energy score to extend engression to right-censored multivariate outcomes, improving over naive training on censored observations. Limitations. In practice, the censoring distribution must usually be estimated. This is a nuisanceestimation problem, not a circularity, as long as the censoring model is fixed independently of the forecast being scored. Our sensitivity results suggest that such estimation is less damaging than forecast-dependent weighting, though misspecification can still degrade performance. Code availability. tion.

The code accompanying this work will be made publicly available upon publica-

Acknowledgments and Disclosure of Funding Jef Jonkers is funded by the Research Foundation Flanders (FWO, Ref. 1S11525N).

References [1] Tilmann Gneiting and Adrian E Raftery. Strictly Proper Scoring Rules, Prediction, and Estimation. Journal of the American Statistical Association, 102(477):359–378, March 2007. ISSN 0162-1459. doi: 10.1198/016214506000001437. URL https://doi.org/10.1198/ 016214506000001437. [2] Ramon F. A. de Punder, Cees G. H. Diks, Roger J. A. Laeven, and Dick J. C. van Dijk. Localizing Strictly Proper Scoring Rules. Journal of the American Statistical Association, 0(0):1–13, January 2026. ISSN 0162-1459. doi: 10.1080/01621459.2025.2576189. URL https://doi.org/10.1080/01621459.2025.2576189. [3] Hiroki Yanagisawa. Proper Scoring Rules for Survival Analysis. In Proceedings of the 40th International Conference on Machine Learning, pages 39165–39182. PMLR, July 2023. URL https://proceedings.mlr.press/v202/yanagisawa23a.html. [4] Xinwei Shen and Nicolai Meinshausen. Engression: extrapolation through the lens of distributional regression. Journal of the Royal Statistical Society Series B: Statistical Methodology, 87(3):653–677, July 2025. ISSN 1369-7412. doi: 10.1093/jrsssb/qkae108. URL https://doi.org/10.1093/jrsssb/qkae108. [5] Erika Graf, Claudia Schmoor, Willi Sauerbrei, and Martin Schumacher. Assessment and comparison of prognostic classification schemes for survival data. Statistics in Medicine, 18(17-18):2529–2545, 1999. ISSN 1097-0258. doi: 10.1002/(SICI)1097-0258(19990915/30)18:17/18<2529::AID-SIM274>3.0.CO;2-5. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/%28SICI%291097-0258% 2819990915/30%2918%3A17/18%3C2529%3A%3AAID-SIM274%3E3.0.CO%3B2-5. [6] Thomas A. Gerds and Martin Schumacher. Consistent Estimation of the Expected Brier Score in General Survival Models with Right-Censored Event Times. Biometrical Journal, 48(6):1029–1040, 2006. ISSN 1521-4036. doi: 10.1002/bimj.200610301. URL https: //onlinelibrary.wiley.com/doi/abs/10.1002/bimj.200610301. [7] Håvard Kvamme and Ørnulf Borgan. The Brier Score under Administrative Censoring: Problems and a Solution. Journal of Machine Learning Research, 24(2):1–26, 2023. ISSN 1533-7928. URL http://jmlr.org/papers/v24/19-1030.html. [8] David Rindt, Robert Hu, David Steinsaltz, and Dino Sejdinovic. Survival regression with proper scoring rules and monotonic neural networks. In Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, pages 1190–1205. PMLR, May 2022. URL https://proceedings.mlr.press/v151/rindt22a.html. [9] Anand Avati, Tony Duan, Sharon Zhou, Kenneth Jung, Nigam H. Shah, and Andrew Y. Ng. Countdown Regression: Sharp and Calibrated Survival Predictions. In Proceedings of The 35th 10

Uncertainty in Artificial Intelligence Conference, pages 145–155. PMLR, August 2020. URL https://proceedings.mlr.press/v115/avati20a.html. [10] Alistair E. W. Johnson, Lucas Bulgarelli, Lu Shen, Alvin Gayles, Ayad Shammout, Steven Horng, Tom J. Pollard, Sicheng Hao, Benjamin Moody, Brian Gow, Li-wei H. Lehman, Leo A. Celi, and Roger G. Mark. MIMIC-IV, a freely accessible electronic health record dataset. Scientific Data, 10(1):1, January 2023. ISSN 2052-4463. doi: 10.1038/s41597-022-01899-x. URL https://www.nature.com/articles/s41597-022-01899-x. [11] Section 2: AKI Definition. Kidney International Supplements, 2(1):19–36, March 2012. ISSN 2157-1716, 2157-1724. doi: 10.1038/kisup.2011.32. URL https://www.kisupplements. org/article/S2157-1716(15)31031-5/fulltext.

A

Theory and proofs

A.1

Proof of Proposition 1

We use Sc♭ (F ; ψc♭ (t)) and Sc♭ (F ; Y, ∆) interchangeably for the abstract and right-censoring encodings of the same fixed-c censored observation. Proof. Fix x. Let Px denote the true conditional law of T | X = x and let Fx denote a candidate forecast. Let H(· | x) be the true conditional distribution of the censoring time C | X = x. For each fixed censoring time c, define the induced censored laws ♭ Pc,x := L(ψc♭ (T ) | X = x),

♭ Fc,x := L(ψc♭ (Z) | X = x),

Z ∼ Fx .

By definition,   S̄x (F ) − S̄x (P ) = EP S̄(F ; Y, ∆, X) − S̄(P ; Y, ∆, X) | X = x . The marginalized score averages the fixed-c localized score over the conditional law of C given the observed data, Π(dc | Y, ∆, X) := L(C ∈ dc | Y, ∆, X), where this conditional law is induced by the true censoring distribution H(· | X). Hence, by the tower property, h h i i ♭ ♭ S̄x (F ) − S̄x (P ) = EP EP SC (Fx ; Y, ∆) − SC (Px ; Y, ∆) | Y, ∆, X X = x h i ♭ ♭ = EP SC (Fx ; Y, ∆) − SC (Px ; Y, ∆) | X = x . Under conditional independent censoring, T ⊥ C | X, we may condition on the realized censoring time and integrate over its true conditional law: Z h i S̄x (F ) − S̄x (P ) = EPx Sc♭ (Fx ; ψc♭ (T )) − Sc♭ (Px ; ψc♭ (T )) H(dc | x). ♭ For fixed c, the distribution of ψc♭ (T ) under the truth is exactly Pc,x . Therefore the inner expectation is the divergence of the localized score between the induced censored laws: h i   ♭ ♭ EPx Sc♭ (Fx ; ψc♭ (T )) − Sc♭ (Px ; ψc♭ (T )) = DSc♭ Pc,x ∥Fc,x .

Hence

Z S̄x (F ) − S̄x (P ) =

  ♭ ♭ DSc♭ Pc,x ∥Fc,x H(dc | x).

By propriety of the fixed-censoring localized score, the integrand is nonnegative for H(· | x)-almost every c. Consequently, S̄x (F ) ≥ S̄x (P ). If, in addition, the localized scores are strictly proper and the collection of induced censored laws for censoring times with positive probability identifies the latent law on Ix , then equality can hold only if ♭ ♭ Pc,x = Fc,x

11

for H(· | x)-almost every relevant c. By the assumed identifiability condition, this is equivalent to Fx = Px

on Ix .

Therefore, S̄x (F ) = S̄x (P )

A.2

⇐⇒

Fx = Px on Ix .

Population optimality of censored engression

Proposition 2 (Population optimality of censored engression). Assume conditional independent censoring, T ⊥ C | X, and assume that Π(dc | Y, ∆, X) is the true forecast-independent conditional law of the censoring time given (Y, ∆, X). Assume that, for H(· | X = x)-almost every censoring time c, the fixed-c censored energy score is well-defined and strictly proper for the corresponding censored law. Suppose there exists g0 ∈ M such that g0 (x, ε) ∼ P (T | X = x)

for PX -almost every x.

Let ge ∈ arg ming∈M Rcens-ES (g). Then L(ψc (e g (x, ε)) | X = x) = L(ψc (T ) | X = x) for H(· | X = x)-almost every c and PX -almost every x. If these censored laws identify the latent law on an identifiable region Ix , then L(e g (x, ε) | X = x) = P (T | X = x)

on Ix

for PX -almost every x. Proof. Fix x and a censoring time c. By assumption, the fixed-c censored energy score is strictly proper for the censored law L(ψc (T ) | X = x). Hence, the corresponding fixed-c expected loss is minimized when L(ψc (g(x, ε)) | X = x) = L(ψc (T ) | X = x). Since Π(dc | Y, ∆, X) is the true forecast-independent conditional law of the censoring time given the observed data, the marginalized objective is the corresponding average of these fixed-c proper score divergences. Under correct specification, g0 induces the true latent law and therefore attains zero integrated divergence, so it minimizes the population objective. Any other population minimizer ge must also attain zero integrated divergence. Since the integrand is nonnegative for H(· | X = x)almost every c, it follows that L(ψc (e g (x, ε)) | X = x) = L(ψc (T ) | X = x) for H(· | X = x)-almost every c and PX -almost every x. The final claim then follows from the stated identifiability condition.

B

Specific score derivations

In each example, we proceed in the same way: we start from the uncensored score, define a fixedcensoring localized score, marginalize over the censoring time to obtain an observed-data score, and then relate this observed-data score back to the original latent target. B.1

Logarithmic score

For the logarithmic score, the marginalization step simplifies substantially. Under right censoring, the fixed-c localized score depends on the censoring time only through the censored contribution, which is either independent of c or observed exactly from the data. For fixed censoring time c, the induced censored law on the observed space has a density-mass representation  f (y | x), δ = 1, y ≤ c, qcF (y, δ | x) = S(c | x), δ = 0, y = c, 12

where f (· | x) is the conditional density of T and 1 − F (t | x) = P(T > t | X = x) is the corresponding survival function. Hence, the fixed-c localized logarithmic score is ♭ Slog,c (F ; y, δ, x) = − log qcF (y, δ | x) = −δ log f (y | x) − (1 − δ) log(1 − F (c | x)).

(4)

The marginalized score is h i ♭ S̄log (F ; y, δ, x) = E Slog,C (F ; y, δ, x) | Y = y, ∆ = δ, X = x . We distinguish the two possible values of δ. If δ = 1, then conditioning on Y = y, ∆ = 1, X = x implies that T = y and C ≥ y. By (4), ♭ Slog,C (F ; y, 1, x) = − log f (y | x),

which does not depend on C. Therefore, S̄log (F ; y, 1, x) = − log f (y | x). If δ = 0, then conditioning on Y = y, ∆ = 0, X = x implies that C = y. Again by (4), ♭ Slog,C (F ; y, 0, x) = − log(1 − F (C | x)) = − log(1 − F (y | x)),

and thus S̄log (F ; y, 0, x) = − log(1 − F (y | x)). Combining the two cases yields S̄log (F ; y, δ, x) = −δ log f (y | x) − (1 − δ) log(1 − F (y | x)). In discrete or mixed settings, a boundary convention is needed at Y = C. With the convention ∆ = 1{T ≤ C}, the censored contribution is naturally written in terms of P(T > Y | X), i.e. 1 − F (Y | X). B.2

CRPS

For the CRPS, the marginalization step is nontrivial. Under right censoring, the fixed-c localized score depends on the censoring time through the upper integration limit, so averaging over the conditional law of C is essential. Recall that for a distribution function F (· | x) and observation t, the CRPS is Z t Z ∞ 2 SCRPS (F ; t, x) = F (s | x) ds + (1 − F (s | x))2 ds. 0

t

For fixed censoring time c, the censored law is the law of min(T, c), whose distribution function is  F (s | x), s < c, ♭ Fc (s | x) = 1, s ≥ c. ♭ Hence, the fixed-c localized CRPS is SCRPS,c (F ; y, δ, x) = SCRPS (Fc♭ ; y, x), which can be written as Z y Z c ♭ SCRPS,c (F ; y, δ, x) =

F (s | x)2 ds + δ

0

(1 − F (s | x))2 ds.

(5)

y

If δ = 1, then y ≤ c is observed exactly and the second term runs only up to the censoring time c, because Fc♭ (s | x) = 1 for s ≥ c. If δ = 0, then necessarily y = c, and the second term vanishes. The marginalized CRPS is h i ♭ S̄CRPS (F ; y, δ, x) = E SCRPS,C (F ; y, δ, x) | Y = y, ∆ = δ, X = x . If δ = 0, then conditioning on Y = y, ∆ = 0, X = x implies that C = y, so Z y S̄CRPS (F ; y, 0, x) = F (s | x)2 ds. 0

13

If δ = 1, then conditioning on Y = y, ∆ = 1, X = x implies that C ≥ y. Using (5), "Z Z C

y

(1 − F (s | x)) ds C ≥ y, X = x .

F (s | x) ds + E

S̄CRPS (F ; y, 1, x) =

y

0

By Fubini’s theorem, "Z

#

C

Z ∞

2

(1 − F (s | x)) ds C ≥ y, X = x =

E

#

2

2

P(C > s | C ≥ y, X = x) (1 − F (s | x))2 ds.

y

y

Writing G(t− | x) := P(C ≥ t | X = x)„ this becomes Z ∞ G(s | x) (1 − F (s | x))2 ds. − | x) G(y y Therefore, Z y

Z ∞

2

F (s | x) ds + δ

S̄CRPS (F ; y, δ, x) = 0

y

G(s | x) (1 − F (s | x))2 ds. G(y − | x)

(6)

This is the usual right-censored CRPS under conditional independent censoring. B.3

Energy score

For the energy score, the fixed-censoring localized score is most naturally written in its expectation form. Unlike the logarithmic score, the marginalization step does not simplify away in general, because the fixed-censoring score depends on the censoring configuration through both expectation terms. In the univariate case with β = 1, the energy score reduces to the CRPS. Let β ∈ (0, 2), and recall that for a predictive law F (· | x) on RK + , the energy score is    1  SES,β (F ; t, x) = E ∥Z − t∥β | X = x − E ∥Z − Z ′ ∥β | X = x , 2 iid

where Z, Z ′ ∼ F (· | x). For fixed componentwise censoring threshold c = (c1 , . . . , cK ), define ψc (t) = (min(t1 , c1 ), . . . , min(tK , cK )). If T ∼ F (· | X = x), the induced censored law is Fc♭ (· | x) := L(ψc (T ) | X = x), and the observed censored outcome is y obs = ψc (y). Hence the fixed-c localized energy score is ♭ SES,β,c (F ; y obs , x) = SES,β (Fc♭ ; y obs , x)   1   = E ∥ψc (Z) − y obs ∥β | X = x − E ∥ψc (Z) − ψc (Z ′ )∥β | X = x , (7) 2 iid

where Z, Z ′ ∼ F (· | X = x). If the censoring threshold is random, the marginalized energy score is obtained by averaging (7) over the conditional law of C given the observed data: h i ♭ S̄ES,β (F ; Y, ∆, X) = E SES,β,C (F ; Y, ∆, X) | Y, ∆, X . In contrast to the logarithmic score and the CRPS, no further simplification is available in general, since the fixed-c localized score depends explicitly on c through both expectation terms. 14

B.4

Brier score

Let τ ≥ 0 denote a fixed evaluation horizon. The uncensored Brier score at horizon τ is
2 SBS,τ (F ; t, x) = F (τ | x) − 1{t ≤ τ } .

(8)

Note that we can do similar derivations for the binomial log-likelihood for τ ≥ 0. For fixed censoring time c, let Yc := min(T, c), Fc♭ := L(Yc | X = x). The natural localized Brier score at horizon τ is the ordinary Brier score applied to the censored distribution:
2 ♭ SBS,τ,c (F ; t, x) := Fc♭ (τ | x) − 1{min(t, c) ≤ τ } . (9) Since Fc♭ (τ | x) = P(min(T, c) ≤ τ | X = x) = the score can be written more explicitly as ( ♭ SBS,τ,c (F ; t, x) =



F (τ | x), τ < c, 1, τ ≥ c,

2 F (τ | x) − 1{t ≤ τ } ,

τ < c,

0,

τ ≥ c.

(10)

Thus, for each fixed c, the localized Brier score coincides with the ordinary Brier score on the identifiable region τ < c and becomes vacuous beyond it. In particular, because the ordinary Brier score is strictly proper for Bernoulli outcomes, (9) is strictly proper for the censored Bernoulli target 1{min(T, c) ≤ τ }, and therefore locally proper for the latent law on horizons τ < c. We now marginalize over the censoring time. Let Y = min(T, C), ∆ = 1{T ≤ C}, and define h i ♭ S̄BS,τ (F ; Y, ∆, X) := E SBS,τ,C (F ; Y, X) | Y, ∆, X . Using (10), we obtain the observed-data score by distinguishing three cases. If Y > τ , then necessarily C > τ and T > τ , so S̄BS,τ (F ; Y, ∆, X) = F (τ | X)2 . If Y ≤ τ and ∆ = 0, then C = Y ≤ τ , hence τ ≥ C and S̄BS,τ (F ; Y, 0, X) = 0. If Y ≤ τ and ∆ = 1, then T = Y ≤ τ and C ≥ Y . In this case,  (1 − F (τ | X))2 , C > τ, SBS,τ,C (F ; Y, X) = 0, Y ≤ C ≤ τ. Therefore,

S̄BS,τ (F ; Y, 1, X) = P(C > τ | C ≥ Y, X) (1 − F (τ | X))2 .

Writing

G(t− | x) := P(C ≥ t | X = x),

G(t | x) := P(C > t | X = x), this becomes S̄BS,τ (F ; Y, 1, X) =

G(τ | X) (1 − F (τ | X))2 . G(Y − | X)

Combining the three cases yields G(τ | X) (1 − F (τ | X))2 . G(Y − | X) This is the observed-data version of the localized proper Brier score. S̄BS,τ (F ; Y, ∆, X) = 1{Y > τ }F (τ | X)2 + ∆ 1{Y ≤ τ }

(11)

A convenient rescaling recovers the usual IPCW Brier score: dividing (11) by G(τ | X) yields S̄BS,τ (F ; Y, ∆, X) ∆1{Y ≤ τ }(1 − F (τ | X))2 1{Y > τ }F (τ | X)2 = + , G(τ | X) G(Y − | X) G(τ | X) which is exactly the standard IPCW Brier score [5, 6]. Since this rescales the score by a strictly positive factor that does not depend on F , it preserves the same minimizer and hence the same propriety characteristics on the identifiable region. 15

B.4.1

Integrated Brier score

The integrated version fits especially naturally into the localization framework. Define the fixed-c localized integrated Brier score by Z ∞ ♭ ♭ SIBS,c (F ; t, x) := SBS,τ,c (F ; t, x) dτ. (12) 0

Using (10), this reduces to ♭ SIBS,c (F ; t, x) =

Z c

2 F (τ | x) − 1{t ≤ τ } dτ.

(13)

0

Thus the fixed-censoring localized IBS is simply the ordinary integrated Brier score truncated to the identifiable region [0, c). Marginalizing over the censoring time gives h i ♭ S̄IBS (F ; Y, ∆, X) = E SIBS,C (F ; Y, X) | Y, ∆, X . By Fubini’s theorem and (11), Z ∞ S̄IBS (F ; Y, ∆, X) = S̄BS,τ (F ; Y, ∆, X) dτ 0

Z Y =

F (τ | X)2 dτ + ∆

0

Z ∞ Y

G(τ | X) (1 − F (τ | X))2 dτ. G(Y − | X)

(14)

This is the integrated observed-data score induced by the localized proper Brier construction. In particular, (14) has the same structural form as the marginalized CRPS: an integral of squared CDF error up to the observed time, plus a tail term weighted by the conditional survival of the censoring time. In the univariate continuous setting, this is not merely a formal similarity: the integrated Brier score and the CRPS coincide, so the marginalized IBS and the marginalized CRPS are the same score up to the usual boundary convention. B.5

Pinball loss

Let α ∈ (0, 1) denote a quantile level, and let qα (F | x) := inf{t ∈ R+ : F (t | x) ≥ α} be the corresponding conditional α-quantile implied by the forecast F (· | x). The uncensored quantile score is

SQ,α (F ; t, x) = α − 1{t < qα (F | x)} t − qα (F | x) . (15) Equivalently,

SQ,α (F ; t, x) = α t − qα (F | x) + + (1 − α) qα (F | x) − t + . This score is strictly consistent for the α-quantile functional. For fixed censoring time c, let Yc := min(T, c), Fc♭ := L(Yc | X = x). The natural localized quantile score is the ordinary quantile score applied to the censored law:

♭ SQ,α,c (F ; t, x) := α − 1{min(t, c) < qα (Fc♭ | x)} min(t, c) − qα (Fc♭ | x) . (16) Since

qα (Fc♭ | x) = min qα (F | x), c ,

this can be written as

♭ SQ,α,c (F ; t, x) = α − 1{min(t, c) < min(qα (F | x), c)} min(t, c) − min(qα (F | x), c) . (17) Thus the fixed-censoring localized score is simply the ordinary quantile score for the censored outcome min(T, c) and the censored quantile min(qα (F | x), c). 16

Let Y = min(T, C), ∆ = 1{T ≤ C}, and define h i ♭ S̄Q,α (F ; Y, ∆, X) := E SQ,α,C (F ; Y, X) | Y, ∆, X . Write qα (X) := qα (F | X),

G(t | x) := P(C > t | X = x),

G(t− | x) := P(C ≥ t | X = x).

If ∆ = 0, then C = Y , so qα (FC♭ | X) = min(qα (F | X), Y ), and therefore


S̄Q,α (F ; Y, 0, X) = α−1{Y < min(qα (F | X), Y )} Y −min(qα (F | X), Y ) = α Y −qα (F | X) + . If ∆ = 1, then T = Y and C ≥ Y . In this case, 

 S̄Q,α (F ; Y, 1, X) = E α − 1{Y < min(qα (F | X), C)} Y − min(qα (F | X), C) | C ≥ Y, X . If qα (F | X) ≤ Y , then min(qα (F | X), C) = qα (F | X) for every C ≥ Y , so
S̄Q,α (F ; Y, 1, X) = α Y − qα (F | X) . If qα (F | X) > Y , then


α−1{Y < min(qα (F | X), C)} Y −min(qα (F | X), C) = (1−α) min(qα (F | X), C)−Y , and hence 
 1−α E min(qα (F | X), C) − Y 1{C ≥ Y } | X S̄Q,α (F ; Y, 1, X) = − G(Y | X) Z qα (F |X) 1−α = G(t | X) dt. G(Y − | X) Y Combining the two cases yields Z qα (F |X)
1−α S̄Q,α (F ; Y, ∆, X) = α Y − qα (F | X) + + ∆ G(t | X) dt, (18) G(Y − | X) Y with the convention that the integral is zero whenever qα (X) ≤ Y .

C

Additional experimental details and results

C.1

Compute resources

All synthetic experiments were run on a machine with an Apple M1 Pro chip. The AKI use case was run on a Tesla V100-SXM3-32GB GPU with 32 GB of memory. We note that the computational requirements are primarily driven by training the neural-network models. They therefore depend mainly on the dataset size and on the size and type of the chosen architecture, rather than on the use of the proposed loss functions or proper scoring rules per se. C.2

Forecast-evaluation simulation setup

This subsection gives the exact simulation design for the forecast-evaluation experiment reported in Section 5.1. It specifies the latent event-time law, candidate forecasts, censoring regimes, discretization grid, and scoring rules used in the comparison. Latent event-time model and grid. We simulate covariates X = (X1 , X2 , X3 ) ∼ N (0, I3 ) and latent event times from T | X = x ∼ Weibull(k, λ(x)), k = 1.5, log λ(x) = 0.3 + 0.8x1 − 0.5x2 + 0.3x3 . All reported results are based on a test sample of size N = 1000. To match the discrete forecast-based scores, we represent forecasts on a grid 0 = ζ0 < ζ1 < · · · < ζB = zmax , B = 50, where zmax = 20.5471 is the empirical 0.995-quantile of a large pilot sample from the oracle distribution. We also evaluate continuous versions of our CRPS, pinball, Brier, and log scores before discretization when a continuous forecast representation is available; however, for conciseness, we omit them from Table 2, as they have the same ranking as their discrete counterparts. Note that their scores are lower, since the data-generating process is continuous and the discretization is an approximation. 17

Censoring regimes.

We consider four regimes.

Regime A: administrative censoring. The censoring time is deterministic, C ≡ cadmin , with cadmin = 0.9833 in the realized simulation, giving an event rate of 50.00%. Because the censoring time is known exactly, we evaluate this regime using localized scores, as in Section 3.3. Regime B: random independent censoring. The censoring time is sampled independently of T and X: C ∼ Uniform(0, Cmax ),

Cmax = 0.4zmax = 8.2188.

The event rate is 79.70%. Since the censoring law is marginal and forecast-independent, we evaluate this regime with marginalized scores using the true marginal censoring distribution. Regime C: conditionally independent censoring. The censoring time depends on covariates but remains conditionally independent of T : C | X = x ∼ Weibull(k, µC (x)),

k = 1.5,

log µC (x) = 0.2 − 0.3x1 + 0.4x3 .

The event rate is 48.60%. We evaluate this regime using marginalized scores with the true conditional censoring law. Regime D: targeted forecast-dependent weights stress test. Let j = ⌊B/2⌋ = 25, so that ζj−1 = 9.8626 and ζj = 10.2736. Define aj = ζj−1 + 0.25(ζj − ζj−1 ) = 9.9653,

bj = ζj−1 + 0.50(ζj − ζj−1 ) = 10.0681.

The oracle event distribution places probability 12 on (bj , ζj ] and probability 12 on (ζB−1 , ζB ], where ζB−1 = 20.1362 and ζB = 20.5471. Censoring is generated by a mixture: with probability 0.6, C ∼ Uniform(ζj−1 , aj ), so censoring occurs before the interior event mass; with probability 0.4, C is set to a late administrative value beyond zmax , so the event is observed. The resulting event rate is 39.50%. This regime is designed to allow forecast-dependent completion to be gamed by concentrating mass in the far tail. Candidate forecasts. In regimes A–C, we compare five forecasts. F0 is the oracle Weibull forecast. F1 is a shifted Weibull forecast with λ1 (x) = e0.25 λ(x). F2 is an overdispersed and rightshifted discrete forecast obtained by smoothing the oracle pmf with a symmetric kernel proportional to (1, 2, 3, 2, 1) and shifting mass two bins to the right. F3 is underdispersed and right-shifted, constructed by blending the oracle pmf with a narrow peaked template centered two bins to the right of the oracle median bin, using width 1.25 and mixing weight 0.70. F4 is upper-tail-heavy, with p4,i (x) ∝ p0,i (x) exp(i/B). (D)

Regime D uses a targeted stress-test construction. There, the oracle F0 places half its mass in an interior bin interval and half its mass in the last grid bin. It is compared with exploit forecasts F5 (ϵ), with ϵ ∈ {10−3 , 5 × 10−3 , 10−2 , 5 × 10−2 }, that place only ϵ mass in the interior event interval and nearly all remaining mass in the final bin, up to a small numerical probability floor. Scoring rules and evaluation protocol. As infeasible gold standards, we evaluate latent CRPS, pinball, Brier, and negative log-likelihood scores with respect to the unobserved true event times. On censored data, we evaluate the forecast-independent scores proposed in this paper: localized scores in regime A and marginalized scores in regimes B–D. For direct comparison with Yanagisawa [3], we thus evaluate their right-censored forecast-dependent weighted log score, multiclass Brier score, binary Brier score, pinball score, and ranked probability score (RPS). Each is evaluated either with oracle weights computed from F0 , or with plug-in weights computed from the forecast being scored. The binary Brier score is averaged over event-time percentiles, using the targeted grid threshold in regime D. C.3

Raw score values for forecast-evaluation simulations

The main text reports oracle ranks to summarize whether each scoring rule identifies the oracle forecast as best. Here we report the corresponding raw mean score values. For regimes A–C, Table 4 reports the oracle score together with the best competing forecast score in parentheses. For regime D, Table 5 reports the oracle score and the best exploit forecast score separately. Lower values are better in all tables. Bold entries indicate ranking reversals, where a non-oracle forecast obtains a lower score than the oracle. 18

Table 4: Raw mean scores for regimes A–C. Each entry gives the oracle forecast score, with the best competing forecast score in parentheses. Lower is better. Bold entries indicate ranking reversals, where the best competing forecast obtains a lower score than the oracle. Scoring rule

A

B

C

CRPS, latent Pinball, latent Brier, latent NLL, latent

0.7766 (0.7905) 0.3030 (0.3057) 0.1016 (0.1021) 1.0745 (1.1389)

0.7766 (0.7905) 0.3030 (0.3057) 0.1016 (0.1021) 1.0745 (1.1389)

0.7766 (0.7905) 0.3030 (0.3057) 0.1016 (0.1021) 1.0745 (1.1389)

NLL [3], plug-in Brier MC [3], plug-in Brier [3], plug-in Pinball [3], plug-in RPS [3], plug-in

0.9841 (0.9800) 1.5163 (1.5171) 0.9086 (0.9056) 0.7721 (0.7735) 0.7689 (0.7694) 0.7676 (0.7690) 0.1041 (0.1045) 0.1034 (0.1040) 0.1020 (0.1019) 0.2947 (0.3400) 0.3286 (0.3477) 0.3935 (0.4197) 1.6839 (1.4671) 1.6877 (1.7565) 1.6683 (1.4206)

NLL [3], oracle Brier MC [3], oracle Brier [3], oracle Pinball [3], oracle RPS [3], oracle

0.9841 (0.9856) 0.7721 (0.7726) 0.1041 (0.1045) 0.2947 (0.3400) 1.6839 (1.7096)

1.5163 (1.5193) 0.7689 (0.7695) 0.1034 (0.1039) 0.3286 (0.3448) 1.6877 (1.7098)

0.9086 (0.9099) 0.7676 (0.7680) 0.1020 (0.1023) 0.3935 (0.4178) 1.6683 (1.6938)

Censored NLL CRPS, ours Pinball, ours Brier, ours

0.3811 (0.4036) 0.1129 (0.1133) 0.0416 (0.0417) 0.0457 (0.0460)

0.7350 (0.7702) 0.3729 (0.3753) 0.1418 (0.1429) 0.0873 (0.0880)

0.3445 (0.3688) 0.1268 (0.1278) 0.0473 (0.0476) 0.0494 (0.0497)

C.4

Table 5: Raw mean scores for regime D. Entries report the oracle score and the best exploit forecast score. Lower is better. Bold entries indicate ranking reversals, where an exploit forecast obtains a lower score than the oracle. Scoring rule

Oracle score

Best exploit score

CRPS, latent Pinball, latent Brier, latent NLL, latent

2.6189 0.9320 0.0250 −0.5448

4.8200 1.8684 0.0455 0.2944

NLL [3], plug-in Brier MC [3], plug-in Brier [3], plug-in Pinball [3], plug-in RPS [3], plug-in

0.6931 0.5000 0.0250 1.1289 6.2500

0.6716 0.3683 0.0184 1.4673 4.6039

NLL [3], oracle Brier MC [3], oracle Brier [3], oracle Pinball [3], oracle RPS [3], oracle

0.6931 0.5000 0.0250 1.1289 6.2500

1.4235 0.8438 0.0422 1.5317 10.5473

Censored NLL CRPS, ours Pinball, ours Brier, ours

−0.2050 1.0195 0.3741 0.0099

0.0832 1.6975 0.6843 0.0167

Sensitivity to estimating the censoring distribution

We use the same synthetic setup as in Section 5.1, with the exact simulation details given in Appendix C.2, but now replace the true censoring law in regimes B and C by estimated censoring distributions. The estimators are fit by treating censoring as the event. We consider a Kaplan–Meier estimator and a pooled Weibull model. In regime B, these pooled estimators are well matched to the data-generating process because censoring is independent of both T and X. In regime C, they are intentionally misspecified because they ignore the covariate dependence of C | X. Table 6 reports the resulting raw score values for our marginalized CRPS, pinball, and Brier scores, using both discrete and continuous implementations where available. The rank summary is as follows. In regime B, both the Kaplan–Meier and pooled Weibull censoring estimators preserve the oracle ordering for all three discrete scores. In regime C, the pooled Kaplan–Meier estimator also ranks the oracle first for all three discrete scores. The pooled Weibull estimator preserves the oracle ranking for discrete CRPS and pinball, but the discrete Brier score ranks the near-oracle tail forecast slightly ahead of the oracle. This reflects ordinary censoring-model misspecification rather than forecast-dependent circularity: once the censoring estimator is fixed, all candidate forecasts are still evaluated by the same forecast-independent score. C.5

Synthetic censored-engression setup

This subsection gives the data-generating processes and implementation details for the synthetic multivariate censored-engression experiments in Section 5.2. The experiments use simulated data so that the latent event times and the true data-generating conditional law are available for oracle evaluation. Covariates.

For each observation i, we draw Xi = (Xi1 , Xi2 , Xi3 , Xi4 )⊤ ∼ N (0, I4 ).

The latent event-time vector has dimension k, where k = 2 for the bivariate experiments and k ∈ {3, 5, 10} for the scaling experiments. Base coefficient functions. For both latent DGPs, we use the same covariate-dependent location and scale functions. Let b = (−0.6, −0.2, 0.2, 0.6)⊤ . For event coordinate j = 1, . . . , k, define µj (x) = 0.35 + x⊤ βj ,

σj (x) = 0.35 + 0.15 sigmoid(x⊤ γj ), 19

Table 6: Sensitivity of our scores to estimated censoring distributions. Entries are mean scores; lower is better. Regime B uses censoring independent of both T and X. Regime C uses covariate-dependent censoring, but the censoring estimator is pooled and therefore ignores the covariate dependence. Regime / estimator Forecast

CRPS disc. CRPS cont. Pinball disc. Pinball cont. Brier disc. Brier cont.

Regime B, KM censoring estimator F0 oracle 0.4326 F1 shifted 0.4863 F2 over 0.7960 F3 under 0.6281 F4 tail 0.4405

0.3748 0.4002 0.6497 0.5095 0.3766

0.1637 0.1827 0.2943 0.2346 0.1666

0.1420 0.1524 0.2427 0.1931 0.1432

0.1023 0.1093 0.2072 0.1470 0.1031

0.0875 0.0926 0.1819 0.1269 0.0881

Regime B, Weibull censoring estimator F0 oracle 0.4326 F1 shifted 0.4874 F2 over 0.7926 F3 under 0.6248 F4 tail 0.4416

0.3752 0.4037 0.6456 0.5061 0.3786

0.1640 0.1837 0.2935 0.2333 0.1674

0.1424 0.1535 0.2416 0.1919 0.1441

0.1028 0.1096 0.2070 0.1472 0.1035

0.0872 0.0922 0.1808 0.1261 0.0878

Regime C, KM censoring estimator, pooled F0 oracle 0.2402 F1 shifted 0.2618 F2 over 0.4316 F3 under 0.3301 F4 tail 0.2432

0.1295 0.1365 0.2708 0.1845 0.1303

0.0891 0.0965 0.1578 0.1198 0.0901

0.0482 0.0510 0.0998 0.0679 0.0486

0.1023 0.1084 0.2099 0.1466 0.1028

0.0496 0.0516 0.1073 0.0716 0.0499

Regime C, Weibull censoring estimator, pooled F0 oracle 0.2399 F1 shifted 0.2604 F2 over 0.4352 F3 under 0.3317 F4 tail 0.2424

0.1281 0.1334 0.2701 0.1838 0.1288

0.0888 0.0956 0.1587 0.1200 0.0896

0.0474 0.0495 0.0991 0.0674 0.0476

0.1279 0.1285 0.2300 0.1772 0.1273

0.0496 0.0517 0.1080 0.0720 0.0499

where
⊤ βj = (0.35 + 0.08(j − 1)) cos(b1 + 0.7(j − 1)), . . . , cos(b4 + 0.7(j − 1)) , and

⊤ γj = 0.10 sin(b1 − 0.7(j − 1)), . . . , sin(b4 − 0.7(j − 1)) .

Unimodal log-normal DGP. The original, easier DGP is a correlated log-normal model. Let εi ∼ N (0, Σρ ), where Σρ is the k × k equicorrelation matrix with diagonal entries 1 and off-diagonal entries ρ = 0.45. Conditional on Xi = x, the latent log-times are log Tij = µj (x) + σj (x)εij ,

j = 1, . . . , k.

We use this DGP only as a sanity check in the appendix, with k = 2. Mixture log-normal DGP. The main experiments use a covariate-dependent two-regime mixture log-normal DGP. Let
π(x) = sigmoid 0.8x1 − 0.6x2 + 0.4 sin(x3 ) + 0.25(x24 − 1) , and draw a latent regime indicator Zi | Xi = x ∼ Bernoulli(π(x)). The two regimes correspond to an early and a late event-time mode. Let ℓj be event-coordinatespecific loadings, evenly spaced between 0.85 and 1.15: j−1 (1.15 − 0.85), k−1 with ℓ1 = 1 when k = 1. The regime shift is ℓj = 0.85 +

aj = 1.15 ℓj . 20

j = 1, . . . , k,

(0)

(1)

Let εi ∼ N (0, Σ0 ) and εi ∼ N (0, Σ1 ), independently, where Σ0 and Σ1 are k ×k equicorrelation matrices with off-diagonal correlations 0.70 and 0.20, respectively. The early regime has stronger dependence and smaller marginal scale, while the late regime has weaker dependence and larger marginal scale: (0) (0) log Tij = µj (x) − aj + 0.85 σj (x)εij , (1)

(1)

log Tij = µj (x) + aj + 1.15 σj (x)εij . The observed latent event time is then (0) (1) log Tij = (1 − Zi ) log Tij + Zi log Tij . This DGP induces a nonlinear, covariate-dependent, multimodal conditional distribution and is intentionally outside the Weibull-copula family used by the parametric likelihood baselines. Censoring mechanisms. A single scalar censoring time Ci is shared by all event coordinates. The observed outcomes are Yij = min(Tij , Ci ), ∆ij = 1{Tij ≤ Ci }, j = 1, . . . , k. Thus, censoring is perfectly dependent across coordinates through the shared censoring variable Ci . We consider three censoring mechanisms: • Administrative censoring: C ≡ cadmin = 3.0. • Random independent censoring: C ∼ Uniform(0, cmax ), cmax = 5.0. • Conditionally independent censoring: C | X = x ∼ Uniform(0, cmax (x)), with cmax (x) = 5.0 (0.45 + 0.70 sigmoid(0.45x1 − 0.35x2 + 0.20x3 )) . The administrative regime corresponds to fixed localized censoring. The two random censoring regimes require the marginalized censored-energy objective because, for uncensored observations, the realized censoring time is only known to exceed the largest observed event time. Methods and baselines. We compare censored engression with oracle references, naive observed engression, and likelihood-based survival baselines. The DGP row evaluates samples from the true conditional law and is not a fitted method. Latent engression is an infeasible benchmark trained directly on latent event times T . Naive observed engression trains the same generator architecture on Y as if it were fully observed. Censored engression trains the generator with the localized or marginalized censored energy score. The likelihood-based baselines are fitted only on (X, Y, ∆) using the censored logarithmic score. We include independent Weibull, conditional MLP independent Weibull, Clayton copula Weibull, and conditional MLP Clayton copula Weibull baselines. For k = 2, we additionally include Gumbel copula Weibull and joint discrete-grid likelihood baselines. The latter is included only in low dimension because its output size scales as O(B k ) for B grid bins. Training setup. All engression models use the same noise-conditioned MLP generator with hidden width 128, noise dimension 128, depth 2, and a softplus output transform. Training, validation, and test sizes are 4000, 1000, and 1000, respectively. We report means and standard errors over five independent repetitions. Engression models are trained for at most 200 epochs with early stopping. For censored engression, we use 16 latent samples per minibatch update. For evaluation, we use 1024 latent samples. In the random censoring regimes, the marginalized score is approximated with 8 censoring draws during training and 512 censoring draws during evaluation. The likelihood baselines are trained by minimizing their observed-data censored negative loglikelihood. Validation negative log-likelihood is used for early stopping. For the conditional MLP Weibull and conditional MLP Clayton baselines, we use hidden width 128, depth 2, batch size 512, and AdamW optimization. 21

Table 7: Synthetic multivariate survival experiment under a covariate-dependent mixture log-normal latent DGP with k = 2 event times and one shared censoring time. Entries are mean ± standard error over five independent repetitions. Lower is better. DGP denotes the true data-generating conditional law and is therefore an oracle reference, not a fitted method. Latent engression is the infeasible benchmark trained directly on latent event times. Among feasible methods, the best result in each censoring regime is shown in bold.

C.6

Censoring

Method

Censored ES

Latent ES

Administrative Administrative

DGP Latent engression

0.6332 ± 0.0061 0.6424 ± 0.0055

1.2699 ± 0.0183 1.2862 ± 0.0182

Administrative Administrative Administrative Administrative Administrative Administrative Administrative Administrative

Naive observed engression Censored engression Independent Weibull MLP independent Weibull Clayton–Weibull MLP Clayton–Weibull Gumbel–Weibull Joint grid likelihood

0.6444 ± 0.0053 0.6401 ± 0.0053 0.7850 ± 0.0033 0.7029 ± 0.0044 0.7586 ± 0.0027 0.6745 ± 0.0055 0.7567 ± 0.0025 0.7143 ± 0.0069

1.3941 ± 0.0180 1.3237 ± 0.0193 1.4765 ± 0.0165 1.3772 ± 0.0175 1.4515 ± 0.0175 1.3388 ± 0.0190 1.4509 ± 0.0172 1.4761 ± 0.0156

Uniform Uniform

DGP Latent engression

0.4840 ± 0.0067 0.4907 ± 0.0061

1.2468 ± 0.0168 1.2601 ± 0.0169

Uniform Uniform Uniform Uniform Uniform Uniform Uniform Uniform

Naive observed engression Censored engression Independent Weibull MLP independent Weibull Clayton–Weibull MLP Clayton–Weibull Gumbel–Weibull Joint grid likelihood

0.5524 ± 0.0125 0.4904 ± 0.0073 0.5975 ± 0.0071 0.5333 ± 0.0054 0.5781 ± 0.0070 0.5152 ± 0.0070 0.5773 ± 0.0071 0.5245 ± 0.0091

1.4440 ± 0.0155 1.2887 ± 0.0207 1.4703 ± 0.0198 1.3532 ± 0.0208 1.4426 ± 0.0200 1.3233 ± 0.0207 1.4409 ± 0.0194 1.3586 ± 0.0183

Conditional uniform Conditional uniform

DGP Latent engression

0.3811 ± 0.0053 0.3868 ± 0.0049

1.2468 ± 0.0168 1.2601 ± 0.0169

Conditional uniform Conditional uniform Conditional uniform Conditional uniform Conditional uniform Conditional uniform Conditional uniform Conditional uniform

Naive observed engression Censored engression Independent Weibull MLP independent Weibull Clayton–Weibull MLP Clayton–Weibull Gumbel–Weibull Joint grid likelihood

0.4553 ± 0.0104 0.3871 ± 0.0055 0.4774 ± 0.0056 0.4178 ± 0.0048 0.4633 ± 0.0053 0.4051 ± 0.0057 0.4620 ± 0.0053 0.4124 ± 0.0066

1.5009 ± 0.0199 1.3088 ± 0.0236 1.4879 ± 0.0214 1.3698 ± 0.0217 1.4640 ± 0.0206 1.3399 ± 0.0198 1.4628 ± 0.0209 1.3582 ± 0.0219

Additional AKI use-case setup and results

This subsection describes the ICU AKI use case used in Section 5.3. The goal is methodological: we use a realistic, irregularly sampled ICU prediction problem to test whether a censoring-aware distributional objective improves observed-data predictive distributions. We therefore report processed cohort sizes and endpoint rates, but do not include a demographics table. Cohort construction and prediction rows. We use MIMIC-IV ICU stays and construct one prediction row per eligible ICU stay and timestamp. Prediction rows are kept only when they occur within the first 14 days after ICU admission and have positive remaining follow-up time before ICU discharge. The two event times are the remaining times from the prediction timestamp to KDIGO [11] stage-2 AKI according to the creatinine criterion and according to the urine-output criterion, see Table 12. The shared censoring time C is the remaining time to ICU discharge. For endpoint j ∈ {1, 2}, Yj = min(Tj , C), ∆j = 1{Tj ≤ C}, where Tj is the latent time to the corresponding AKI definition. Times are transformed during training as log(1 + t) , t 7→ log(1 + q0.95 ) where q0.95 is the 95th percentile of the training discharge-time distribution. All reported localized energy scores and fixed-horizon metrics are computed after transforming sampled times back to hours. The split is chronological at the ICU-stay level. Stays from anchor-year groups 2008–2016, plus one-half of the 2017–2019 stays, are used for training; the remaining 2017–2019 stays are used for validation; 2020–2022 stays are held out for testing. This yields 4,831,810 training rows, 575,500 22

Table 8: Synthetic multivariate survival experiment under a covariate-dependent mixture log-normal latent DGP with k = 3 event times and one shared censoring time. Entries are mean ± standard error over five independent repetitions. Lower is better. DGP denotes the true data-generating conditional law and is therefore an oracle reference, not a fitted method. Latent engression is the infeasible benchmark trained directly on latent event times. Among feasible methods, the best result in each censoring regime is shown in bold. Censoring

Method

Censored ES

Latent ES

Administrative Administrative

DGP Latent engression

0.9189 ± 0.0055 0.9457 ± 0.0044

3.2238 ± 0.0389 3.2923 ± 0.0416

Administrative Administrative Administrative Administrative Administrative Administrative

Naive observed engression Censored engression Independent Weibull MLP independent Weibull Clayton–Weibull MLP Clayton–Weibull

0.9486 ± 0.0050 0.9434 ± 0.0039 1.1043 ± 0.0067 1.0809 ± 0.0056 0.9905 ± 0.0045 0.9817 ± 0.0019

4.5247 ± 0.0667 3.7475 ± 0.0427 3.8241 ± 0.0559 4.0990 ± 0.0955 3.6364 ± 0.0520 3.5207 ± 0.0361

Uniform Uniform

DGP Latent engression

0.7451 ± 0.0073 0.7651 ± 0.0059

3.1883 ± 0.0268 3.2389 ± 0.0296

Uniform Uniform Uniform Uniform Uniform Uniform

Naive observed engression Censored engression Independent Weibull MLP independent Weibull Clayton–Weibull MLP Clayton–Weibull

0.9003 ± 0.0089 0.7584 ± 0.0071 0.8922 ± 0.0086 0.8724 ± 0.0085 0.8059 ± 0.0085 0.8014 ± 0.0078

4.4224 ± 0.0810 3.5878 ± 0.0608 3.6567 ± 0.0421 3.8675 ± 0.1367 3.4878 ± 0.0338 3.4073 ± 0.0359

Conditional uniform Conditional uniform

DGP Latent engression

0.5965 ± 0.0073 0.6137 ± 0.0060

3.1883 ± 0.0268 3.2389 ± 0.0296

Conditional uniform Conditional uniform Conditional uniform Conditional uniform Conditional uniform Conditional uniform

Naive observed engression Censored engression Independent Weibull MLP independent Weibull Clayton–Weibull MLP Clayton–Weibull

0.7489 ± 0.0113 0.6089 ± 0.0061 0.7155 ± 0.0087 0.7006 ± 0.0086 0.6463 ± 0.0084 0.6416 ± 0.0080

4.5448 ± 0.0852 3.6355 ± 0.0626 3.6283 ± 0.0399 3.7896 ± 0.0984 3.4849 ± 0.0396 3.4174 ± 0.0448

validation rows, and 835,661 test rows. The training row-level event rates are 18.3% for creatininedefined AKI and 39.7% for urine-output-defined AKI; however, note that this is an artificially high value because stays in which the event has already occurred are also counted as events. We do so because urine-defined AKI may already have occurred; however, creatinine-defined AKI may still occur, and we still want to model it. Additionally, one can exit the event at the next time point. During training, the loss function is set to zero for events that have already occurred, since this can be inferred from the observation at that moment. The median remaining discharge time is 54.7 hours. Covariates and preprocessing. Each prediction row has 22 standardized covariates. These include gender, age at admission, weight, chronic renal failure status, baseline serum creatinine, elapsed ICU time, and 72-hour summary features for serum creatinine and urine output: count, mean, minimum, maximum, latest value, and slope. We also add missingness indicators and slope-availability indicators for the two time-varying measurements. Imputation and scaling are fit on the training split only: missing counts are set to zero, missing measurement summaries are filled with training medians, missing slopes are set to zero, and all model features are standardized using training means and standard deviations. Training objectives and baselines. The censored engression model is a noise-conditioned MLP generator with noise dimension 500, hidden width 500, depth 2, and a softplus output transform. It is trained with AdamW, learning rate 2 · 10−4 , weight decay 10−5 , batch size 1024, gradient clipping at 1, and early stopping on the validation localized censored energy score. Each training update uses 16 generated samples. A stay-aware minibatch sampler prevents two rows from the same ICU stay from appearing in the same minibatch. The naive engression baseline uses the same generator class but replaces the localized censored objective by an ordinary energy-score objective that treats the observed censored vector Y as if it were the true event-time vector. We compare against five likelihood-based baselines trained on the same processed tensors: independent Weibull, conditional MLP independent Weibull, Clayton copula Weibull, conditional MLP Clayton copula Weibull, and a joint discrete-grid likelihood model with 32 bins. The conditional Weibull and conditional Clayton baselines use hidden width 128 and depth 2; the joint grid baseline is included only because this AKI task is bivariate. 23

Table 9: Synthetic multivariate survival experiment under a covariate-dependent mixture log-normal latent DGP with k = 5 event times and one shared censoring time. Entries are mean ± standard error over five independent repetitions. Lower is better. DGP denotes the true data-generating conditional law and is therefore an oracle reference, not a fitted method. Latent engression is the infeasible benchmark trained directly on latent event times. Among feasible methods, the best result in each censoring regime is shown in bold. Censoring

Method

Censored ES

Latent ES

Administrative Administrative

DGP Latent engression

1.2108 ± 0.0129 1.3204 ± 0.0272

5.3617 ± 0.0205 5.5121 ± 0.0363

Administrative Administrative Administrative Administrative Administrative Administrative

Naive observed engression Censored engression Independent Weibull MLP independent Weibull Clayton–Weibull MLP Clayton–Weibull

1.2445 ± 0.0129 1.2258 ± 0.0121 1.4678 ± 0.0119 1.4394 ± 0.0112 1.3084 ± 0.0122 1.3140 ± 0.0167

7.9679 ± 0.0529 7.1165 ± 0.0267 6.1825 ± 0.0446 6.6426 ± 0.1223 5.9206 ± 0.0317 5.8265 ± 0.0393

Uniform Uniform

DGP Latent engression

0.9536 ± 0.0100 1.0477 ± 0.0204

5.1394 ± 0.0078 5.3048 ± 0.0275

Uniform Uniform Uniform Uniform Uniform Uniform

Naive observed engression Censored engression Independent Weibull MLP independent Weibull Clayton–Weibull MLP Clayton–Weibull

1.1433 ± 0.0070 0.9692 ± 0.0110 1.1529 ± 0.0132 1.1326 ± 0.0132 1.0356 ± 0.0117 1.0441 ± 0.0120

7.6696 ± 0.0984 6.6095 ± 0.0641 5.9272 ± 0.0301 6.1641 ± 0.0837 5.6937 ± 0.0365 5.5976 ± 0.0432

Conditional uniform Conditional uniform

DGP Latent engression

0.7638 ± 0.0105 0.8536 ± 0.0192

5.1394 ± 0.0078 5.3048 ± 0.0275

Conditional uniform Conditional uniform Conditional uniform Conditional uniform Conditional uniform Conditional uniform

Naive observed engression Censored engression Independent Weibull MLP independent Weibull Clayton–Weibull MLP Clayton–Weibull

0.9440 ± 0.0064 0.7754 ± 0.0117 0.9230 ± 0.0137 0.9022 ± 0.0137 0.8308 ± 0.0124 0.8392 ± 0.0148

7.8476 ± 0.0988 6.7030 ± 0.0690 5.9580 ± 0.0362 6.2551 ± 0.0646 5.7525 ± 0.0428 5.6475 ± 0.0527

Evaluation. Since the discharge censoring time is observed for every prediction row, we evaluate using localized scores. For the bivariate energy score, generated samples are first localized by ♭ ψC (z1 , z2 ) = (min(z1 , C), min(z2 , C)),

and are then compared with the observed censored vector Y . We report the joint localized energy score for rows with at least one endpoint still at risk, and the marginal localized energy scores after excluding rows in which the corresponding endpoint is already present at prediction time. The test set contributes 802,657 rows to the joint score, 725,596 rows to the creatinine marginal score, and 699,383 rows to the urine-output marginal score. We also compute localized Brier scores at fixed horizons for the composite endpoint “either AKI”. For a horizon τ , the censored predicted risk is set to one when C ≤ τ , and otherwise uses the latent sampled risk, e.g.  1, C ≤ τ, ♭ FC,τ (x) = PrF (T1 ≤ τ or T2 ≤ τ | X = x), C > τ. The observed censored binary outcome is 1{Y1 ≤ τ or Y2 ≤ τ }. Top-5% and top-10% positive predictive value and sensitivity are reported only as secondary ranking summaries on the horizonspecific evaluable subset. They are not valid scores and should not be interpreted as estimates of deployment utility. Fixed-horizon results. Table 13 reports the composite endpoint at 6, 12, 24, and 48 hours. Censored engression obtains the best localized Brier score at all four horizons. The gap is smallest at 6 hours and becomes more pronounced from 12 to 48 hours, where censoring by ICU discharge is more consequential. The secondary enrichment metrics are directionally consistent with the score results: censored engression usually selects a higher-risk top tail than the naive generator and the likelihood baselines. These summaries are useful for clinical interpretation, but the localized Brier score is the primary fixed-horizon criterion. Distributional results. Table 14 compares methods that can generate full bivariate event-time samples. Censored engression has the lowest joint localized energy score and the lowest marginal 24

Table 10: Synthetic multivariate survival experiment under a covariate-dependent mixture log-normal latent DGP with k = 10 event times and one shared censoring time. Entries are mean ± standard error over five independent repetitions. Lower is better. DGP denotes the true data-generating conditional law and is therefore an oracle reference, not a fitted method. Latent engression is the infeasible benchmark trained directly on latent event times. Among feasible methods, the best result in each censoring regime is shown in bold. Censoring

Method

Censored ES

Latent ES

Administrative Administrative

DGP Latent engression

1.6530 ± 0.0092 2.0170 ± 0.0403

9.6954 ± 0.0843 10.2168 ± 0.0699

Administrative Administrative Administrative Administrative Administrative Administrative

Naive observed engression Censored engression Independent Weibull MLP independent Weibull Clayton–Weibull MLP Clayton–Weibull

1.7181 ± 0.0098 17.3182 ± 0.2313 1.6884 ± 0.0107 15.7127 ± 0.1847 2.0100 ± 0.0088 14.9379 ± 0.6241 1.9585 ± 0.0118 13.3928 ± 0.2939 1.8086 ± 0.0074 14.3108 ± 0.5964 1.8793 ± 0.0186 13.1867 ± 0.2063

Uniform Uniform

DGP Latent engression

1.3131 ± 0.0055 1.6082 ± 0.0190

Uniform Uniform Uniform Uniform Uniform Uniform

Naive observed engression Censored engression Independent Weibull MLP independent Weibull Clayton–Weibull MLP Clayton–Weibull

1.6402 ± 0.0228 16.5122 ± 0.3144 1.3382 ± 0.0052 14.2576 ± 0.3110 1.5901 ± 0.0086 13.9407 ± 0.4351 1.5491 ± 0.0066 12.4447 ± 0.1425 1.4397 ± 0.0073 13.3828 ± 0.4909 1.4783 ± 0.0200 12.7285 ± 0.5428

Conditional uniform Conditional uniform

DGP Latent engression

1.0562 ± 0.0067 1.3390 ± 0.0175

Conditional uniform Conditional uniform Conditional uniform Conditional uniform Conditional uniform Conditional uniform

Naive observed engression Censored engression Independent Weibull MLP independent Weibull Clayton–Weibull MLP Clayton–Weibull

1.3623 ± 0.0180 16.8903 ± 0.3222 1.0815 ± 0.0055 14.6636 ± 0.2549 1.2751 ± 0.0089 13.7644 ± 0.4093 1.2406 ± 0.0083 12.7779 ± 0.1484 1.1574 ± 0.0079 13.2792 ± 0.4230 1.1837 ± 0.0091 12.5259 ± 0.4063

9.4613 ± 0.0914 10.0111 ± 0.1553

9.4613 ± 0.0914 10.0111 ± 0.1553

localized energy scores (CRPS essentially) for both AKI definitions. The joint discrete-grid likelihood baseline is the closest competitor, which suggests that flexible joint modeling is important for this task; however, it remains slightly worse than censored engression and is tied to the low-dimensional bivariate setting. Conditional Weibull and conditional Clayton baselines substantially improve upon their non-conditional counterparts, confirming the importance of covariate-dependent event-time distributions, but they do not match the sample-based censored engression objective. Overall, the AKI use case supports the same methodological conclusion as the synthetic experiments: treating censored observations as fully observed event times can distort the learned distribution, whereas localizing the score at the observed censoring time yields better predictive distributions based on the observed data. Interpretational caveats. The row-level event rates above are prediction-row rates, not patient-level incidence estimates, because the same ICU stay can contribute multiple prediction times. In addition, ICU discharge is an observed exit process rather than a randomized censoring mechanism. The reported localized scores, therefore, evaluate the predictive distribution under the stated observed-data censoring convention; they should not be interpreted as causal estimates of latent AKI incidence after discharge. The use case is intended to demonstrate the proposed scoring and training methodology, not to define a final clinical deployment model.

25

Table 11: Synthetic multivariate survival experiment under a regular log-normal latent DGP with k = 2 event times and one shared censoring time. Entries are mean ± standard error over five independent repetitions. Lower is better. DGP denotes the true data-generating conditional law and is therefore an oracle reference, not a fitted method. Latent engression is the infeasible benchmark trained directly on latent event times. Among feasible methods, the best result in each censoring regime is shown in bold. Censoring

Method

Censored ES

Latent ES

Administrative Administrative

DGP Latent engression

0.4288 ± 0.0040 0.4305 ± 0.0048

0.7176 ± 0.0069 0.7196 ± 0.0063

Administrative Administrative Administrative Administrative Administrative Administrative Administrative Administrative

Naive observed engression Censored engression Independent Weibull MLP independent Weibull Clayton–Weibull MLP Clayton–Weibull Gumbel–Weibull Joint grid likelihood

0.4416 ± 0.0035 0.4308 ± 0.0041 0.4365 ± 0.0039 0.4380 ± 0.0049 0.4329 ± 0.0038 0.4360 ± 0.0051 0.4340 ± 0.0040 0.4506 ± 0.0037

0.8642 ± 0.0135 0.7892 ± 0.0123 0.7268 ± 0.0071 0.7403 ± 0.0081 0.7225 ± 0.0069 0.7299 ± 0.0091 0.7236 ± 0.0071 0.9149 ± 0.0167

Uniform Uniform

DGP Latent engression

0.3221 ± 0.0056 0.3238 ± 0.0060

0.7175 ± 0.0117 0.7259 ± 0.0119

Uniform Uniform Uniform Uniform Uniform Uniform Uniform Uniform

Naive observed engression Censored engression Independent Weibull MLP independent Weibull Clayton–Weibull MLP Clayton–Weibull Gumbel–Weibull Joint grid likelihood

0.3975 ± 0.0088 0.3257 ± 0.0055 0.3275 ± 0.0054 0.3290 ± 0.0063 0.3251 ± 0.0053 0.3280 ± 0.0064 0.3257 ± 0.0055 0.3407 ± 0.0057

0.9709 ± 0.0211 0.7886 ± 0.0249 0.7259 ± 0.0127 0.7421 ± 0.0138 0.7229 ± 0.0128 0.7395 ± 0.0126 0.7225 ± 0.0128 0.8191 ± 0.0204

Conditional uniform Conditional uniform

DGP Latent engression

0.2672 ± 0.0055 0.2687 ± 0.0057

0.7175 ± 0.0117 0.7259 ± 0.0119

Conditional uniform Conditional uniform Conditional uniform Conditional uniform Conditional uniform Conditional uniform Conditional uniform Conditional uniform

Naive observed engression Censored engression Independent Weibull MLP independent Weibull Clayton–Weibull MLP Clayton–Weibull Gumbel–Weibull Joint grid likelihood

0.3557 ± 0.0098 0.2694 ± 0.0057 0.2717 ± 0.0051 0.2729 ± 0.0061 0.2699 ± 0.0053 0.2729 ± 0.0064 0.2703 ± 0.0053 0.2798 ± 0.0057

1.0458 ± 0.0231 0.7927 ± 0.0198 0.7266 ± 0.0128 0.7504 ± 0.0147 0.7235 ± 0.0125 0.7478 ± 0.0132 0.7233 ± 0.0129 0.8203 ± 0.0199

Table 12: KDIGO Acute Kidney Injury (AKI) stage definitions for serum creatinine (sCr) and urine output (UO) used in this work. Stage

sCr criteria

UO criteria

Stage 1

sCr ≥ 1.5 − 1.9× baseline, or UO ≤ 0.5 mL/kg/h for 6–12 h increase in sCr ≥ 0.3 mg/dL

Stage 2

sCr ≥ 2.0 − 2.9× baseline

UO ≤ 0.5 mL/kg/h for ≥ 12 h

Stage 3

sCr ≥ 3.0× baseline, or sCr ≥ 4.0 mg/dL

UO ≤ 0.3 mL/kg/h for ≥ 24 h

26

Table 13: Secondary utility metrics for the composite endpoint (either AKI). Localized Brier is a score, so lower is better; PPV and sensitivity are utility metrics, so higher is better. These metrics are descriptive only. Bold entries indicate the best method within each horizon and metric. Horizon

Method

Loc. Brier ↓

PPV@5% ↑

Sens.@5% ↑

PPV@10% ↑

Sens.@10% ↑

6h 6h 6h 6h 6h 6h 6h

Naive Censored Joint grid Indep. Weibull Clayton Weibull Cond. indep. Weibull Cond. Clayton Weibull

0.023 0.022 0.028 0.223 0.206 0.041 0.047

0.485 0.501 0.465 0.177 0.179 0.403 0.366

0.528 0.545 0.505 0.193 0.194 0.439 0.398

0.338 0.348 0.341 0.157 0.159 0.313 0.290

0.735 0.757 0.742 0.342 0.345 0.680 0.632

12h 12h 12h 12h 12h 12h 12h

Naive Censored Joint grid Indep. Weibull Clayton Weibull Cond. indep. Weibull Cond. Clayton Weibull

0.045 0.040 0.043 0.216 0.201 0.051 0.054

0.727 0.782 0.728 0.367 0.369 0.655 0.633

0.316 0.340 0.317 0.160 0.160 0.285 0.275

0.615 0.664 0.622 0.327 0.330 0.560 0.547

0.535 0.578 0.542 0.284 0.287 0.488 0.476

24h 24h 24h 24h 24h 24h 24h

Naive Censored Joint grid Indep. Weibull Clayton Weibull Cond. indep. Weibull Cond. Clayton Weibull

0.081 0.071 0.074 0.190 0.179 0.077 0.078

0.759 0.840 0.805 0.514 0.515 0.762 0.757

0.167 0.185 0.177 0.113 0.114 0.168 0.167

0.696 0.751 0.724 0.470 0.473 0.690 0.689

0.307 0.331 0.319 0.207 0.208 0.304 0.304

48h 48h 48h 48h 48h 48h 48h

Naive Censored Joint grid Indep. Weibull Clayton Weibull Cond. indep. Weibull Cond. Clayton Weibull

0.095 0.080 0.081 0.138 0.132 0.083 0.083

0.769 0.929 0.911 0.676 0.678 0.878 0.886

0.094 0.113 0.111 0.083 0.083 0.107 0.108

0.763 0.882 0.862 0.648 0.649 0.831 0.840

0.186 0.215 0.211 0.158 0.158 0.203 0.205

Table 14: Additional AKI use-case benchmark results using localized energy scores. Lower is better. These scores evaluate overall fit of the predicted event-time distributions. Method Naive observed engression Censored engression Independent Weibull Conditional MLP independent Weibull Clayton copula + Weibull Conditional MLP Clayton copula + Weibull Joint discrete-grid likelihood

Joint localized ES

Creatinine localized ES

Urine localized ES

43.442 30.933 45.948 35.728 45.749 34.006 31.675

23.591 15.439 25.915 19.534 25.915 18.863 16.210

32.981 23.769 34.154 26.061 34.154 24.606 24.057

27