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 \(F_X\) 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 \(\Delta=\mathbb{1}\{T\le C\}\), rather than the latent event time \(T\) itself. As a result, standard proper scoring rules for uncensored outcomes cannot be applied directly. 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),\mathbb{1}\{T\le c\})\) under \(T\sim 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 [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\ge Y\). We therefore average the localized score over the conditional law of \(C\) given \((Y,\Delta,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.

1.0.0.1 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 [tab:right-censoring-scores]). 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, [5] introduced the inverse probability of censoring weighting (IPCW) Brier score and its integrated version, and [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 [7], arises from a positivity issue when evaluation extends beyond the identifiable follow-up region. Likewise, the non-propriety results in [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, [9] proposed an adapted continuous ranked probability score (CRPS) for censored outcomes, but [8] gave counterexamples showing that this score is not proper. More recently, [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.

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 [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\in\mathcal{X}\subseteq\mathbb{R}^d\) denote covariates, and let \(T\in\mathcal{T}\subseteq\mathbb{R}_+^k\) denote the latent event-time outcome, with \(k=1\) in the univariate case and \(k\ge2\) in the multivariate setting. For each covariate value \(x\), we seek a predictive distribution \(F=\mathcal{L}(T\mid X=x)\). Under right censoring, however, we observe only \(Y=\min(T,C)\) and \(\Delta=\mathbb{1}\{T\le C\}\).

Equip \(\mathcal{T}\) with a \(\sigma\)-algebra \(\mathscr T\), let \(\mathcal{P}\) be a class of probability measures on \((\mathcal{T},\mathscr T)\), and let \(S:\mathcal{P}\times\mathcal{T}\to\bar{\mathbb{R}}\) be a scoring rule with divergence \[D_S(P\|F)=\mathbb{E}_P[S(F,T)]-\mathbb{E}_P[S(P,T)].\] We call \(S\) strictly proper if \(D_S(P\|F)\ge0\) for all \(P,F\in\mathcal{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\subseteq \mathcal{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 \(\mathcal{T}_A^{\flat}=A\cup\{*\}\). For a latent predictive law \(F\in\mathcal{P}\), define its censored version by \[F_A^{\flat} = F_A+\bar F_A\,\delta_*, \qquad F_A(B):=F(B\cap A), \qquad \bar F_A:=F(A^c).\] 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, \[S_A^{\flat}(F,t) := S(F_A^{\flat},t_A^{\flat}), \qquad \text{where} \quad t_A^{\flat} = \begin{cases} t, & t\in A,\\ *, & t\in A^c. \end{cases}\]

Intuitively, following the localization perspective of [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, \(D_{S_A^{\flat}}(P\|F) = D_S(P_A^{\flat}\|F_A^{\flat})\), so equality of population score implies equality 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 \(\Delta=\mathbb{1}\{T\le c\}\), the identifiable region is \(A_c=[0,c]\), and the censored outcome map is \[\psi^\flat_c(t) = \begin{cases} t, & t\le c,\\ *, & t>c. \end{cases}\] The corresponding censored law is \(F_c^{\flat} = F(\cdot\cap[0,c]) + F((c,\infty))\delta_*\). Thus, an event at the boundary \(c\) is distinguished from a censored observation beyond \(c\). We write \(S_c^\flat(F;Y,\Delta)\) for the localized score induced by this fixed censoring time.

When \(T\mid 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, the abstract censored state \(*\) cannot be inserted directly into Euclidean distances. By a slight abuse of notation, we keep the notation \(\psi_c^\flat\) for this numerical encoding of the same coarsened observation: \[\psi^\flat_c(t)=\min(t,c), \qquad F_c^\flat=\mathcal{L}(\psi^\flat_c(T)).\] This representation should be read as an encoding of the coarsened observation: the value \(c\) represents the censored region \((c,\infty)\), not an event known to occur at \(c\). When the event indicator is needed, we keep the convention \(\Delta=\mathbb{1}\{T\le 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=(T_1,\dots,T_k)\) and componentwise censoring thresholds \(c=(c_1,\dots,c_k)\), we use \[\psi^\flat_c(t) = (\min(t_1,c_1),\dots,\min(t_k,c_k)), \qquad F_c^\flat=\mathcal{L}(\psi^\flat_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 \(c_1=\cdots=c_k=c\).

3.2 Marginalized proper scoring rules↩︎

We now turn to random right censoring. Suppose we observe \(Y=\min(T,C)\) and \(\Delta=\mathbb{1}\{T\le C\}\), and let \(H(t\mid x)=\mathbb{P}(C\le t\mid 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 \[\Pi(dc\mid Y,\Delta,X):=\mathcal{L}(C\in dc\mid Y,\Delta,X)\] be the conditional law of the censoring time given the observed censored outcome, induced by the true censoring distribution \(H(\cdot\mid X)\) and the right-censoring mechanism. Under conditional independent censoring, \(T\perp C\mid X\), the observed-data score is obtained by averaging the fixed-\(c\) localized score over this conditional law: \[\bar S(F;Y,\Delta,X) := \int S_c^\flat(F;Y,\Delta)\,\Pi(dc\mid Y,\Delta,X).\] In the univariate case, this reduces to \[\bar S(F;Y,\Delta,X) = \begin{cases} \displaystyle \int_{[Y,\infty)} S_c^\flat(F;Y,1)\, H(dc\mid C\ge Y,X), & \Delta=1,\\[2mm] S_Y^\flat(F;Y,0), & \Delta=0. \end{cases}\] 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\ge \max_j Y_j\), so \(\Pi\) is the law of \(C\) conditional on \(C\ge \max_j Y_j\).

For fixed \(x\), let \(\bar{\mathcal{S}}_x(F) := \mathbb{E}_P[\bar S(F;Y,\Delta,X)\mid X=x]\) denote the conditional expected score. For a shared scalar censoring time, the identifiable region is \[\mathcal{I}_x^{\mathrm{shared}} = \left\{ t\in\mathbb{R}_+^k: \mathbb{P}\!\left(C\ge \max_{j=1,\ldots,k} t_j \mid X=x\right)>0 \right\}.\] For componentwise censoring times \(C=(C_1,\ldots,C_k)\), the analogous region is \[\mathcal{I}_x^{\mathrm{comp}} = \left\{ t\in\mathbb{R}_+^k: \mathbb{P}\!\left(C_j\ge t_j \;\text{for all } j=1,\ldots,k \mid X=x\right)>0 \right\}.\] We write \(\mathcal{I}_x\) for the relevant region under the censoring mechanism considered. In the univariate case, both definitions reduce to \(\mathcal{I}_x=\{t\ge0:\mathbb{P}(C\ge t\mid X=x)>0\}\). The next proposition states that the marginalized score is minimized by the true latent law, with uniqueness only on \(\mathcal{I}_x\).

Outside \(\mathcal{I}_x\), two latent distributions can differ while inducing the same observed censored law. Hence, no scoring rule based only on \((Y,\Delta,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 8.1 and Table [tab:right-censoring-scores]. For the logarithmic score, marginalization simplifies because the fixed-\(c\) localized contribution is either observed exactly or constant over the compatible censoring values.

4pt

max width=

Table [tab:right-censoring-scores] 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 \(\psi_c^\flat(t)=\min(t,c)\). Detailed derivations are given in Appendix 8.

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 \(\bar S(F;Y,\Delta,X=x)=S_{c(x)}^\flat(F;Y,\Delta)\).

Equivalently, the censoring law is \(H_x(dc)=\delta_{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 \(\bar S(F;Y,\Delta,C,X)=S_C^\flat(F;Y,\Delta)\).

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_\theta(x,\varepsilon)\) 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 \[\label{eq:noise-model} Y_g(x)=g(x,\varepsilon), \qquad \varepsilon\sim P_\varepsilon,\tag{1}\] where \(P_\varepsilon\) is a simple reference distribution. The generator induces the conditional predictive law \(P_g(\cdot\mid x):=\mathcal{L}(g(x,\varepsilon))\), so unlike classical regression, engression directly models the full conditional distribution via samples. For two independent noise variables \(\varepsilon,\varepsilon'\), the energy score population objective is \[\mathcal{R}_{\mathrm{ES}}(g) = \mathbb{E}\left[ \|Y-g(X,\varepsilon)\| - \frac{1}{2}\|g(X,\varepsilon)-g(X,\varepsilon')\| \right],\] where \((X,Y)\sim P\). Under correct specification, minimizing this objective recovers the true conditional law; see Proposition 1 of [4].

4.1 Censored engression↩︎

We now extend engression to multivariate right-censored survival outcomes. Let \(T=(T_1,\ldots,T_k)\in\mathbb{R}_+^k, k \ge 1\), denote the latent event-time vector, and let the generator in 1 induce the latent predictive law \(F_{g,x}:=\mathcal{L}(g(x,\varepsilon)).\)

Under right censoring, the event times are not fully observed. For a fixed censoring time \(c\), let \(\psi^\flat_c(t)=\bigl(\min(t_1,c),\ldots,\min(t_k,c)\bigr)\) denote the shared censoring map, and let \(Y=\psi^\flat_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 \(\mathcal{L}(\psi^\flat_c(g(x,\varepsilon)))\) rather than to the latent law \(F_{g,x}\) itself. The resulting fixed-\(c\) censored energy loss is \[\label{eq:local-cens-engression} \mathcal{R}_{\mathrm{local\text{-}cens\text{-}ES}}(g;Y,x,c) = \mathbb{E}_{\varepsilon} \left[ \|\psi_c(g(x,\varepsilon))-Y\| \right] - \frac{1}{2} \mathbb{E}_{\varepsilon,\varepsilon'} \left[ \|\psi_c(g(x,\varepsilon)) - \psi_c(g(x,\varepsilon'))\| \right],\tag{2}\] where \(\varepsilon,\varepsilon'\stackrel{\mathrm{iid}}{\sim}P_\varepsilon\).

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: \[\label{eq:cens-engression} \mathcal{R}_{\mathrm{cens\text{-}ES}}(g) := \mathbb{E}\!\left[ \mathcal{R}_{\mathrm{local\text{-}cens\text{-}ES}}(g;Y,X,C) \mid Y,\Delta,X \right].\tag{3}\] Under conditional independent censoring, this is exactly the marginalized right-censored energy score objective introduced above. When there is no censoring, \(\psi_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\mid 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 [3]. For the latter, we report both an oracle-weight version, where weights are computed from the true event-time law \(F_0\), and a plug-in version, where weights are computed from the candidate forecast itself. The simulation design is given in Appendix 9.2.

5.1.0.1 Results and discussion.

r0.32 4.0pt

Table 1 reports the rank of the oracle forecast under each scoring rule or score family, for four censoring regimes, with rank \(1\) indicating the best mean score. The corresponding raw mean score values are reported in Appendix 9.3.

The latent gold-standard scores rank the oracle first in every regime, as expected. The same is true for all of our localized and marginalized scores: across regimes A–D, our CRPS, pinball, Brier, and 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.

The oracle-weight scores of [3] also rank \(F_0\) first in all regimes, showing that the underlying scores behave as intended when the weights are fixed at their correct oracle values. The failures arise in the practically relevant plug-in version, where the weights depend on the forecast being evaluated. The plug-in NLL ranks the oracle second in regimes A, C, and D; the plug-in RPS 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 9.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 9.4 show that when \(G(\cdot\mid 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=(T_1,\ldots,T_k)\). The main experiments use a covariate-dependent 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\in\{2,3,5,10\}\) to assess scalability. Full details are given in Appendix 9.5.

We compare censored engression with oracle references, naive observed engression, and likelihood-based 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,\Delta)\): 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 [tab:right-censoring-scores]). 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.

5.2.0.1 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 [tab:censored-engression-selected-results], complete results of all settings and benchmarks are presented in Tables 2–5. 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 Weibull-copula 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 9.5 repeats the \(k=2\) experiment under the original unimodal log-normal 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: 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 de-identified ICU data from MIMIC-IV [10]; additional cohort, endpoint, and evaluation details are given in Appendix 9.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 [tab:censored-engression-selected-results] 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.

4.0pt

max width=

Censored engression improves all selected localized energy scores in Table [tab:censored-engression-selected-results]: the joint score decreases from \(43.442\) to \(30.993\), the creatinine-AKI marginal score from \(23.591\) to \(15.439\), and the urine-output-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 horizon-specific 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 9.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 [tab:aki-utility-compact]) and the best joint and marginal localized energy scores among the methods that generate full event-time distributions (Table 7). 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 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.

6.0.0.1 Limitations.

In practice, the censoring distribution must usually be estimated. This is a nuisance-estimation 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.

6.0.0.2 Code availability.

The code accompanying this work will be made publicly available upon publication.

7 Theory and proofs↩︎

7.1 Proof of Proposition 1↩︎

We use \(S_c^\flat(F;\psi_c^\flat(t))\) and \(S_c^\flat(F;Y,\Delta)\) interchangeably for the abstract and right-censoring encodings of the same fixed-\(c\) censored observation.

7.2 Population optimality of censored engression↩︎

8 Specific score derivations↩︎

In each example, we proceed in the same way: we start from the uncensored score, define a fixed-censoring 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.

8.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 \[q_c^F(y,\delta\mid x) = \begin{cases} f(y\mid x), & \delta=1,\; y\le c,\\ S(c\mid x), & \delta=0,\; y=c, \end{cases}\] where \(f(\cdot\mid x)\) is the conditional density of \(T\) and \(1-F(t\mid x)=\mathbb{P}(T>t\mid X=x)\) is the corresponding survival function. Hence, the fixed-\(c\) localized logarithmic score is \[S^{\flat}_{\log,c}(F;y,\delta,x) = -\log q_c^F(y,\delta\mid x) = -\delta \log f(y\mid x)-(1-\delta)\log (1 - F(c\mid x)). \label{eq:fixed-c-log-score}\tag{4}\]

The marginalized score is \[\bar S_{\log}(F;y,\delta,x) = \mathbb{E}\!\left[ S^{\flat}_{\log,C}(F;y,\delta,x) \mid Y=y,\Delta=\delta,X=x \right].\] We distinguish the two possible values of \(\delta\).

If \(\delta=1\), then conditioning on \(Y=y,\Delta=1,X=x\) implies that \(T=y\) and \(C\ge y\). By 4 , \[S^{\flat}_{\log,C}(F;y,1,x)=-\log f(y\mid x),\] which does not depend on \(C\). Therefore, \(\bar S_{\log}(F;y,1,x)=-\log f(y\mid x)\).

If \(\delta=0\), then conditioning on \(Y=y,\Delta=0,X=x\) implies that \(C=y\). Again by 4 , \[S^{\flat}_{\log,C}(F;y,0,x)=-\log (1-F(C\mid x))=-\log (1 - F(y\mid x)),\] and thus \(\bar S_{\log}(F;y,0,x)=-\log (1 - F(y\mid x))\).

Combining the two cases yields \[\bar S_{\log}(F;y,\delta,x) = -\delta\log f(y\mid x)-(1-\delta)\log (1 - F(y\mid x)).\]

In discrete or mixed settings, a boundary convention is needed at \(Y=C\). With the convention \(\Delta=\mathbb{1}\{T\le C\}\), the censored contribution is naturally written in terms of \(\mathbb{P}(T>Y\mid X)\), i.e.\(1 - F(Y\mid X)\).

8.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(\cdot\mid x)\) and observation \(t\), the CRPS is \[S_\mathrm{CRPS}(F;t,x) = \int_0^t F(s\mid x)^2\,ds + \int_t^\infty (1-F(s\mid x))^2\,ds.\]

For fixed censoring time \(c\), the censored law is the law of \(\min(T,c)\), whose distribution function is \[F_c^{\flat}(s\mid x) = \begin{cases} F(s\mid x), & s<c,\\ 1, & s\ge c. \end{cases}\] Hence, the fixed-\(c\) localized CRPS is \(S_{\mathrm{CRPS}, c}^{\flat}(F;y,\delta,x) = S_\mathrm{CRPS}(F_c^{\flat};y,x)\), which can be written as \[S_{\mathrm{CRPS}, c}^{\flat}(F;y,\delta,x) = \int_0^y F(s\mid x)^2\,ds + \delta\int_y^c (1-F(s\mid x))^2\,ds. \label{eq:fixed-c-crps}\tag{5}\] If \(\delta=1\), then \(y\le c\) is observed exactly and the second term runs only up to the censoring time \(c\), because \(F_c^{\flat}(s\mid x)=1\) for \(s\ge c\). If \(\delta=0\), then necessarily \(y=c\), and the second term vanishes.

The marginalized CRPS is \[\bar{S}_{\mathrm{CRPS}}(F;y,\delta,x) = \mathbb{E}\!\left[ S_{\mathrm{CRPS}, C}^{\flat}(F;y,\delta,x)\mid Y=y,\Delta=\delta,X=x \right].\]

If \(\delta=0\), then conditioning on \(Y=y,\Delta=0,X=x\) implies that \(C=y\), so \[\bar{S}_{\mathrm{CRPS}}(F;y,0,x) = \int_0^y F(s\mid x)^2\,ds.\]

If \(\delta=1\), then conditioning on \(Y=y,\Delta=1,X=x\) implies that \(C\ge y\). Using 5 , \[\bar{S}_{\mathrm{CRPS}}(F;y,1,x) = \int_0^y F(s\mid x)^2\,ds + \mathbb{E}\!\left[ \int_y^C (1-F(s\mid x))^2\,ds \;\middle|\; C\ge y,X=x \right].\]

By Fubini’s theorem, \[\mathbb{E}\!\left[ \int_y^C (1-F(s\mid x))^2\,ds \;\middle|\; C\ge y,X=x \right] = \int_y^\infty \mathbb{P}(C> s\mid C\ge y,X=x)\, (1-F(s\mid x))^2\,ds.\]

Writing \(G(t^-\mid x):=\mathbb{P}(C\ge t \mid X=x)\),, this becomes \[\int_y^\infty \frac{G(s\mid x)}{ G(y^-\mid x)}(1-F(s\mid x))^2\,ds.\]

Therefore, \[\bar{S}_{\mathrm{CRPS}}(F;y,\delta,x) = \int_0^y F(s\mid x)^2\,ds + \delta\int_y^\infty \frac{ G(s\mid x)}{ G(y^-\mid x)} (1-F(s\mid x))^2\,ds. \label{eq:marginalized-crps-right-censoring}\tag{6}\] This is the usual right-censored CRPS under conditional independent censoring.

8.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 \(\beta=1\), the energy score reduces to the CRPS.

Let \(\beta\in(0,2)\), and recall that for a predictive law \(F(\cdot\mid x)\) on \(\mathbb{R}_+^K\), the energy score is \[S_{\mathrm{ES}, \beta}(F;t,x) = \mathbb{E}\!\left[\|Z-t\|^\beta\mid X=x\right] - \frac{1}{2}\mathbb{E}\!\left[\|Z-Z'\|^\beta\mid X=x\right],\] where \(Z,Z'\stackrel{\text{iid}}{\sim}F(\cdot\mid x)\).

For fixed componentwise censoring threshold \(c=(c_1,\dots,c_K)\), define \[\psi_c(t) = (\min(t_1,c_1),\dots,\min(t_K,c_K)).\] If \(T\sim F(\cdot\mid X=x)\), the induced censored law is \[F_c^{\flat}(\cdot\mid x):=\mathcal{L}(\psi_c(T)\mid X=x),\] and the observed censored outcome is \(y^{\mathrm{obs}}=\psi_c(y)\). Hence the fixed-\(c\) localized energy score is \[\begin{align} S_{\mathrm{ES},\beta,c}^{\flat}(F;y^{\mathrm{obs}},x) &= S_{\mathrm{ES},\beta}(F_c^{\flat};y^{\mathrm{obs}},x) \notag\\ &= \mathbb{E}\!\left[\|\psi_c(Z)-y^{\mathrm{obs}}\|^\beta\mid X=x\right] - \frac{1}{2}\mathbb{E}\!\left[\|\psi_c(Z)-\psi_c(Z')\|^\beta\mid X=x\right], \label{eq:fixed-c-energy-score} \end{align}\tag{7}\] where \(Z,Z'\stackrel{\text{iid}}{\sim}F(\cdot\mid 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: \[\bar{S}_{\mathrm{ES}, \beta}(F;Y,\Delta,X) = \mathbb{E}\!\left[ S_{\mathrm{ES}, \beta, C}^{\flat}(F;Y,\Delta,X) \mid Y,\Delta,X \right].\] 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.

8.4 Brier score↩︎

Let \(\tau\ge 0\) denote a fixed evaluation horizon. The uncensored Brier score at horizon \(\tau\) is \[S_{\mathrm{BS},\tau}(F;t,x) = \bigl(F(\tau\mid x)-\mathbb{1}\{t\le \tau\}\bigr)^2. \label{eq:uncensored-brier}\tag{8}\]

Note that we can do similar derivations for the binomial log-likelihood for \(\tau\ge 0\).

For fixed censoring time \(c\), let \(Y_c:=\min(T,c)\), \(F_c^{\flat}:=\mathcal{L}(Y_c\mid X=x)\). The natural localized Brier score at horizon \(\tau\) is the ordinary Brier score applied to the censored distribution: \[S_{\mathrm{BS},\tau,c}^{\flat}(F;t,x) := \bigl(F_c^{\flat}(\tau\mid x)-\mathbb{1}\{\min(t,c)\le \tau\}\bigr)^2. \label{eq:fixed-c-brier-direct}\tag{9}\] Since \[F_c^{\flat}(\tau\mid x) = \mathbb{P}(\min(T,c)\le \tau\mid X=x) = \begin{cases} F(\tau\mid x), & \tau<c,\\ 1, & \tau\ge c, \end{cases}\] the score can be written more explicitly as \[S_{\mathrm{BS},\tau,c}^{\flat}(F;t,x) = \begin{cases} \bigl(F(\tau\mid x)-\mathbb{1}\{t\le \tau\}\bigr)^2, & \tau<c,\\[1mm] 0, & \tau\ge c. \end{cases} \label{eq:fixed-c-brier-piecewise}\tag{10}\] Thus, for each fixed \(c\), the localized Brier score coincides with the ordinary Brier score on the identifiable region \(\tau<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 \(\mathbb{1}\{\min(T,c)\le \tau\}\), and therefore locally proper for the latent law on horizons \(\tau<c\).

We now marginalize over the censoring time. Let \(Y=\min(T,C)\), \(\Delta=\mathbb{1}\{T\le C\}\), and define \[\bar S_{\mathrm{BS},\tau}(F;Y,\Delta,X) := \mathbb{E}\!\left[ S_{\mathrm{BS},\tau,C}^{\flat}(F;Y,X) \mid Y,\Delta,X \right].\]

Using 10 , we obtain the observed-data score by distinguishing three cases. If \(Y>\tau\), then necessarily \(C>\tau\) and \(T>\tau\), so \[\bar S_{\mathrm{BS},\tau}(F;Y,\Delta,X)=F(\tau\mid X)^2.\]

If \(Y\le \tau\) and \(\Delta=0\), then \(C=Y\le \tau\), hence \(\tau\ge C\) and \[\bar S_{\mathrm{BS},\tau}(F;Y,0,X)=0.\] If \(Y\le \tau\) and \(\Delta=1\), then \(T=Y\le \tau\) and \(C\ge Y\). In this case, \[S_{\mathrm{BS},\tau,C}(F;Y,X) = \begin{cases} (1-F(\tau\mid X))^2, & C>\tau,\\ 0, & Y\le C\le \tau. \end{cases}\]

Therefore, \[\bar S_{\mathrm{BS},\tau}(F;Y,1,X) = \mathbb{P}(C>\tau\mid C\ge Y,X)\,(1-F(\tau\mid X))^2.\] Writing \[G(t\mid x):=\mathbb{P}(C>t\mid X=x), \qquad G(t^-\mid x):=\mathbb{P}(C\ge t\mid X=x),\] this becomes \[\bar S_{\mathrm{BS},\tau}(F;Y,1,X) = \frac{G(\tau\mid X)}{G(Y^-\mid X)}(1-F(\tau\mid X))^2.\] Combining the three cases yields \[\bar S_{\mathrm{BS},\tau}(F;Y,\Delta,X) = \mathbb{1}\{Y>\tau\}F(\tau\mid X)^2 + \Delta\,\mathbb{1}\{Y\le \tau\} \frac{G(\tau\mid X)}{G(Y^-\mid X)} (1-F(\tau\mid X))^2. \label{eq:marginalized-brier-right}\tag{11}\] This is the observed-data version of the localized proper Brier score.

A convenient rescaling recovers the usual IPCW Brier score: dividing 11 by \(G(\tau \mid X)\) yields \[\frac{\bar S_{\mathrm{BS},\tau}(F;Y,\Delta,X)}{G(\tau \mid X)} = \frac{\Delta \mathbb{1} \{ Y \le \tau\} (1 - F(\tau \mid X))^2}{G(Y^- \mid X)} + \frac{\mathbb{1} \{ Y > \tau \} F(\tau \mid X)^2}{G(\tau \mid 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.

8.4.1 Integrated Brier score↩︎

The integrated version fits especially naturally into the localization framework. Define the fixed-\(c\) localized integrated Brier score by \[S_{\mathrm{IBS},c}^{\flat}(F;t,x) := \int_0^\infty S_{\mathrm{BS},\tau,c}^{\flat}(F;t,x)\,d\tau. \label{eq:fixed-c-ibs-def}\tag{12}\] Using 10 , this reduces to \[S_{\mathrm{IBS},c}^{\flat}(F;t,x) = \int_0^c \bigl(F(\tau\mid x)-\mathbb{1}\{t\le \tau\}\bigr)^2 \,d\tau. \label{eq:fixed-c-ibs}\tag{13}\] 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 \[\bar S_{\mathrm{IBS}}(F;Y,\Delta,X) = \mathbb{E}\!\left[ S_{\mathrm{IBS},C}^{\flat}(F;Y,X)\mid Y,\Delta,X \right].\]

By Fubini’s theorem and 11 , \[\begin{align} \bar S_{\mathrm{IBS}}(F;Y,\Delta,X) &= \int_0^\infty \bar S_{\mathrm{BS},\tau}(F;Y,\Delta,X)\,d\tau \notag\\ &= \int_0^Y F(\tau\mid X)^2\,d\tau + \Delta\int_Y^\infty \frac{G(\tau\mid X)}{G(Y^-\mid X)} (1-F(\tau\mid X))^2\,d\tau. \label{eq:marginalized-ibs-right} \end{align}\tag{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.

8.5 Pinball loss↩︎

Let \(\alpha\in(0,1)\) denote a quantile level, and let \[q_\alpha(F\mid x):=\inf\{t\in\mathbb{R}_+:F(t\mid x)\ge \alpha\}\] be the corresponding conditional \(\alpha\)-quantile implied by the forecast \(F(\cdot\mid x)\). The uncensored quantile score is \[S_{\mathrm{Q},\alpha}(F;t,x) = \bigl(\alpha-\mathbb{1}\{t<q_\alpha(F\mid x)\}\bigr)\bigl(t-q_\alpha(F\mid x)\bigr). \label{eq:uncensored-quantile-score}\tag{15}\] Equivalently, \[S_{\mathrm{Q},\alpha}(F;t,x) = \alpha\bigl(t-q_\alpha(F\mid x)\bigr)_+ + (1-\alpha)\bigl(q_\alpha(F\mid x)-t\bigr)_+.\] This score is strictly consistent for the \(\alpha\)-quantile functional.

For fixed censoring time \(c\), let \(Y_c:=\min(T,c)\), \(F_c^{\flat}:=\mathcal{L}(Y_c\mid X=x)\). The natural localized quantile score is the ordinary quantile score applied to the censored law: \[S_{\mathrm{Q},\alpha,c}^{\flat}(F;t,x) := \bigl(\alpha-\mathbb{1}\{\min(t,c)<q_\alpha(F_c^{\flat}\mid x)\}\bigr) \bigl(\min(t,c)-q_\alpha(F_c^{\flat}\mid x)\bigr). \label{eq:fixed-c-quantile-direct}\tag{16}\] Since \[q_\alpha(F_c^{\flat}\mid x)=\min\!\bigl(q_\alpha(F\mid x),c\bigr),\] this can be written as \[S_{\mathrm{Q},\alpha,c}^{\flat}(F;t,x) = \bigl(\alpha-\mathbb{1}\{\min(t,c)<\min(q_\alpha(F\mid x),c)\}\bigr) \bigl(\min(t,c)-\min(q_\alpha(F\mid x),c)\bigr). \label{eq:fixed-c-quantile}\tag{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_\alpha(F\mid x),c)\).

Let \(Y=\min(T,C)\), \(\Delta=\mathbb{1}\{T\le C\}\), and define \[\bar S_{\mathrm{Q},\alpha}(F;Y,\Delta,X) := \mathbb{E}\!\left[ S_{\mathrm{Q},\alpha,C}^{\flat}(F;Y,X)\mid Y,\Delta,X \right].\]

Write \[q_\alpha(X):=q_\alpha(F\mid X), \qquad G(t\mid x):=\mathbb{P}(C>t\mid X=x), \qquad G(t^-\mid x):=\mathbb{P}(C\ge t\mid X=x).\]

If \(\Delta=0\), then \(C=Y\), so \(q_\alpha(F_C^{\flat}\mid X)=\min(q_\alpha(F\mid X),Y),\) and therefore \[\bar S_{\mathrm{Q},\alpha}(F;Y,0,X) = \bigl(\alpha-\mathbb{1}\{Y<\min(q_\alpha(F\mid X),Y)\}\bigr) \bigl(Y-\min(q_\alpha(F\mid X),Y)\bigr) = \alpha\bigl(Y-q_\alpha(F\mid X)\bigr)_+.\]

If \(\Delta=1\), then \(T=Y\) and \(C\ge Y\). In this case, \[\bar S_{\mathrm{Q},\alpha}(F;Y,1,X) = \mathbb{E}\!\left[ \bigl(\alpha-\mathbb{1}\{Y<\min(q_\alpha(F\mid X),C)\}\bigr) \bigl(Y-\min(q_\alpha(F\mid X),C)\bigr) \mid C\ge Y,X \right].\]

If \(q_\alpha(F\mid X)\le Y\), then \(\min(q_\alpha(F\mid X),C)=q_\alpha(F\mid X)\) for every \(C\ge Y\), so \[\bar S_{\mathrm{Q},\alpha}(F;Y,1,X) = \alpha\bigl(Y-q_\alpha(F\mid X)\bigr).\] If \(q_\alpha(F\mid X)>Y\), then \[\bigl(\alpha-\mathbb{1}\{Y<\min(q_\alpha(F\mid X),C)\}\bigr) \bigl(Y-\min(q_\alpha(F\mid X),C)\bigr) = (1-\alpha)\bigl(\min(q_\alpha(F\mid X),C)-Y\bigr),\] and hence \[\begin{align} \bar S_{\mathrm{Q},\alpha}(F;Y,1,X) &= \frac{1-\alpha}{G(Y^-\mid X)} \mathbb{E}\!\left[ \bigl(\min(q_\alpha(F\mid X),C)-Y\bigr)\mathbb{1}\{C\ge Y\} \mid X \right] \\ &= \frac{1-\alpha}{G(Y^-\mid X)} \int_Y^{q_\alpha(F\mid X)} G(t\mid X)\,dt. \end{align}\] Combining the two cases yields \[\bar S_{\mathrm{Q},\alpha}(F;Y,\Delta,X) = \alpha\bigl(Y-q_\alpha(F\mid X)\bigr)_+ + \Delta\,\frac{1-\alpha}{G(Y^-\mid X)} \int_Y^{q_\alpha(F\mid X)} G(t\mid X)\,dt, \label{eq:marginalized-quantile-right}\tag{18}\] with the convention that the integral is zero whenever \(q_\alpha(X)\le Y\).

9 Additional experimental details and results↩︎

9.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.

9.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.

9.2.0.1 Latent event-time model and grid.

We simulate covariates \(X=(X_1,X_2,X_3)\sim N(0,I_3)\) and latent event times from \[T\mid X=x \sim \mathrm{Weibull}(k,\lambda(x)), \qquad k=1.5, \qquad \log \lambda(x)=0.3+0.8x_1-0.5x_2+0.3x_3 .\] 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=\zeta_0<\zeta_1<\cdots<\zeta_B=z_{\max}, \qquad B=50,\] where \(z_{\max}=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 1, 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.

9.2.0.2 Censoring regimes.

We consider four regimes.

Regime A: administrative censoring. The censoring time is deterministic, \(C\equiv c_{\mathrm{admin}}\), with \(c_{\mathrm{admin}}=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\sim \mathrm{Uniform}(0,C_{\max}), \qquad C_{\max}=0.4z_{\max}=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\mid X=x \sim \mathrm{Weibull}(k,\mu_C(x)), \qquad k=1.5, \qquad \log \mu_C(x)=0.2-0.3x_1+0.4x_3 .\] 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=\lfloor B/2\rfloor=25\), so that \(\zeta_{j-1}=9.8626\) and \(\zeta_j=10.2736\). Define \[a_j=\zeta_{j-1}+0.25(\zeta_j-\zeta_{j-1})=9.9653, \qquad b_j=\zeta_{j-1}+0.50(\zeta_j-\zeta_{j-1})=10.0681 .\] The oracle event distribution places probability \(\frac{1}{2}\) on \((b_j,\zeta_j]\) and probability \(\frac{1}{2}\) on \((\zeta_{B-1},\zeta_B]\), where \(\zeta_{B-1}=20.1362\) and \(\zeta_B=20.5471\). Censoring is generated by a mixture: with probability \(0.6\), \(C\sim\mathrm{Uniform}(\zeta_{j-1},a_j)\), so censoring occurs before the interior event mass; with probability \(0.4\), \(C\) is set to a late administrative value beyond \(z_{\max}\), 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.

9.2.0.3 Candidate forecasts.

In regimes A–C, we compare five forecasts. \(F_0\) is the oracle Weibull forecast. \(F_1\) is a shifted Weibull forecast with \(\lambda_1(x)=e^{0.25}\lambda(x)\). \(F_2\) is an overdispersed and right-shifted 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. \(F_3\) 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\). \(F_4\) is upper-tail-heavy, with \(p_{4,i}(x)\propto p_{0,i}(x)\exp(i/B)\).

Regime D uses a targeted stress-test construction. There, the oracle \(F_0^{(D)}\) places half its mass in an interior bin interval and half its mass in the last grid bin. It is compared with exploit forecasts \(F_5(\epsilon)\), with \(\epsilon\in\{10^{-3},5\times 10^{-3},10^{-2},5\times 10^{-2}\}\), that place only \(\epsilon\) mass in the interior event interval and nearly all remaining mass in the final bin, up to a small numerical probability floor.

9.2.0.4 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 [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 \(F_0\), 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.

9.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 [tab:oracle-score-values-abc] reports the oracle score together with the best competing forecast score in parentheses. For regime D, Table [tab:oracle-score-values-d] 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.

9.4 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 9.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\mid X\).

Table [tab:estimated-censoring-scores] 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.

3.5pt

max width=

9.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.

9.5.0.1 Covariates.

For each observation \(i\), we draw \[X_i=(X_{i1},X_{i2},X_{i3},X_{i4})^\top \sim N(0,I_4).\] The latent event-time vector has dimension \(k\), where \(k=2\) for the bivariate experiments and \(k\in\{3,5,10\}\) for the scaling experiments.

9.5.0.2 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)^\top .\] For event coordinate \(j=1,\ldots,k\), define \[\mu_j(x)=0.35 + x^\top \beta_j, \qquad \sigma_j(x)=0.35 + 0.15\,\mathrm{sigmoid}(x^\top \gamma_j),\] where \[\beta_j = (0.35 + 0.08(j-1)) \bigl( \cos(b_1+0.7(j-1)),\ldots,\cos(b_4+0.7(j-1)) \bigr)^\top ,\] and \[\gamma_j = 0.10 \bigl( \sin(b_1-0.7(j-1)),\ldots,\sin(b_4-0.7(j-1)) \bigr)^\top .\]

9.5.0.3 Unimodal log-normal DGP.

The original, easier DGP is a correlated log-normal model. Let \(\varepsilon_i\sim N(0,\Sigma_\rho)\), where \(\Sigma_\rho\) is the \(k\times k\) equicorrelation matrix with diagonal entries \(1\) and off-diagonal entries \(\rho=0.45\). Conditional on \(X_i=x\), the latent log-times are \[\log T_{ij} = \mu_j(x)+\sigma_j(x)\varepsilon_{ij}, \qquad j=1,\ldots,k.\] We use this DGP only as a sanity check in the appendix, with \(k=2\).

9.5.0.4 Mixture log-normal DGP.

The main experiments use a covariate-dependent two-regime mixture log-normal DGP. Let \[\pi(x) = \mathrm{sigmoid} \left( 0.8x_1 - 0.6x_2 + 0.4\sin(x_3) + 0.25(x_4^2-1) \right),\] and draw a latent regime indicator \[Z_i\mid X_i=x \sim \mathrm{Bernoulli}(\pi(x)).\] The two regimes correspond to an early and a late event-time mode. Let \(\ell_j\) be event-coordinate-specific loadings, evenly spaced between \(0.85\) and \(1.15\): \[\ell_j = 0.85 + \frac{j-1}{k-1}(1.15-0.85), \qquad j=1,\ldots,k,\] with \(\ell_1=1\) when \(k=1\). The regime shift is \[a_j = 1.15\,\ell_j.\]

Let \(\varepsilon_i^{(0)}\sim N(0,\Sigma_{0})\) and \(\varepsilon_i^{(1)}\sim N(0,\Sigma_{1})\), independently, where \(\Sigma_0\) and \(\Sigma_1\) are \(k\times 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: \[\log T_{ij}^{(0)} = \mu_j(x) - a_j + 0.85\,\sigma_j(x)\varepsilon_{ij}^{(0)},\] \[\log T_{ij}^{(1)} = \mu_j(x) + a_j + 1.15\,\sigma_j(x)\varepsilon_{ij}^{(1)}.\] The observed latent event time is then \[\log T_{ij} = (1-Z_i)\log T_{ij}^{(0)} + Z_i\log T_{ij}^{(1)}.\] This DGP induces a nonlinear, covariate-dependent, multimodal conditional distribution and is intentionally outside the Weibull-copula family used by the parametric likelihood baselines.

9.5.0.5 Censoring mechanisms.

A single scalar censoring time \(C_i\) is shared by all event coordinates. The observed outcomes are \[Y_{ij}=\min(T_{ij},C_i), \qquad \Delta_{ij}=\mathbb{1}\{T_{ij}\le C_i\}, \qquad j=1,\ldots,k.\] Thus, censoring is perfectly dependent across coordinates through the shared censoring variable \(C_i\).

We consider three censoring mechanisms:

• Administrative censoring: \[C \equiv c_{\mathrm{admin}} = 3.0.\]

• Random independent censoring: \[C \sim \mathrm{Uniform}(0,c_{\max}), \qquad c_{\max}=5.0.\]

• Conditionally independent censoring: \[C\mid X=x \sim \mathrm{Uniform}(0,c_{\max}(x)),\] with \[c_{\max}(x) = 5.0 \left( 0.45 + 0.70\,\mathrm{sigmoid}(0.45x_1-0.35x_2+0.20x_3) \right).\]

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.

9.5.0.6 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,\Delta)\) 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.

9.5.0.7 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 log-likelihood. 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.

3.5pt

3.5pt

3.5pt

3.5pt

3.5pt

9.6 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.

[ht]
\centering
\caption{KDIGO Acute Kidney Injury (AKI) stage definitions for serum creatinine (sCr) and urine output (UO) used in this work.}
\label{tab:kdigo-stages}
\begin{tabular}{p{1.1cm} p{4cm} p{4cm}}
\toprule
Stage & sCr criteria & UO criteria \\
\midrule
\textbf{Stage 1} 
& \begin{tabular}[t]{@{}l@{}}
    sCr $\geq$ $1.5-1.9 \times$ baseline, or\\
    increase in sCr $\geq 0.3$ mg/dL
  \end{tabular}
& \begin{tabular}[t]{@{}l@{}}
    UO $\le 0.5$ mL/kg/h for 6--12 h
  \end{tabular}
\\[2mm]
\midrule
\textbf{Stage 2} 
& \begin{tabular}[t]{@{}l@{}}
    sCr $\geq$ $2.0-2.9 \times$ baseline
  \end{tabular}
& \begin{tabular}[t]{@{}l@{}}
    UO $\le 0.5$ mL/kg/h for $\geq 12$ h
  \end{tabular}
\\[2mm]
\midrule
\textbf{Stage 3} 
& \begin{tabular}[t]{@{}l@{}}
    sCr $\geq 3.0 \times$ baseline, or \\
    sCr $\geq 4.0$ mg/dL
  \end{tabular}
& \begin{tabular}[t]{@{}l@{}}
    UO $\le 0.3$ mL/kg/h for $\geq 24$ h
  \end{tabular}
\\
\bottomrule
\end{tabular}

9.6.0.1 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 [tab:kdigo-stages]. The shared censoring time \(C\) is the remaining time to ICU discharge. For endpoint \(j\in\{1,2\}\), \[Y_j=\min(T_j,C), \qquad \Delta_j=\mathbb{1}\{T_j\le C\},\] where \(T_j\) is the latent time to the corresponding AKI definition. Times are transformed during training as \[t \mapsto \frac{\log(1+t)}{\log(1+q_{0.95})},\] where \(q_{0.95}\) is the \(95\)th 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\) validation rows, and \(835{,}661\) test rows. The training row-level event rates are \(18.3\%\) for creatinine-defined 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.

9.6.0.2 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.

9.6.0.3 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\cdot 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.

9.6.0.4 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 \[\psi^\flat_C(z_1,z_2)=(\min(z_1,C),\min(z_2,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 \(\tau\), the censored predicted risk is set to one when \(C\le \tau\), and otherwise uses the latent sampled risk, e.g. \[F^{\flat}_{C,\tau}(x)= \begin{cases} 1, & C\le \tau,\\ \Pr_F(T_1\le \tau \text{ or } T_2\le \tau\mid X=x), & C>\tau. \end{cases}\] The observed censored binary outcome is \(\mathbb{1}\{Y_1\le\tau \text{ or } Y_2\le\tau\}\). Top-\(5\%\) and top-\(10\%\) positive predictive value and sensitivity are reported only as secondary ranking summaries on the horizon-specific evaluable subset. They are not valid scores and should not be interpreted as estimates of deployment utility.

3.5pt

max width=

9.6.0.5 Fixed-horizon results.

Table [tab:aki-utility-compact] 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.

4pt

9.6.0.6 Distributional results.

Table 7 compares methods that can generate full bivariate event-time samples. Censored engression has the lowest joint localized energy score and the lowest marginal 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.

9.6.0.7 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.

References↩︎