Toward Conditional Distribution Calibration in Survival Prediction

Authors: shi-ang  qi
yakun  yu
russell  greiner
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 10:48:16

Abstract: Survival prediction often involves estimating the time-to-event distribution from censored datasets. Previous approaches have focused on enhancing discrimination and marginal calibration. In this paper, we highlight the significance of *conditional calibration* for real-world applications – especially its role in individual decision-making. We propose a method based on conformal prediction that uses the model’s predicted individual survival probability at that instance’s observed time. This method effectively improves the model’s marginal and conditional calibration, without compromising discrimination. We provide asymptotic theoretical guarantees for both marginal and conditional calibration and test it extensively across 15 diverse real-world datasets, demonstrating the method’s practical effectiveness and
versatility in various settings.

Main paper: 
 introduction: 
 Individual survival distribution (ISD), or time-to-event distribution, is a probability distribution that describes the times until the occurrence of a specific event of interest for an instance, based on information about that individual. Accurately estimating ISD is essential for effective decisionmaking and clinical resource allocation. However, a challenge in learning such survival prediction models is training on datasets that include censored instances, where we only know a lower bound of their time-to-event.Survival models typically focus on two important but distinct properties during optimization and evaluation: (i) discrimination measures how well a model's relative predictions between individuals align with the observed order [1, 2] , which is useful for pairwise decisions such as prioritizing treatments; (ii) calibration assesses how well the predicted survival probabilities match the actual distribution of observations [3, 4] , supporting both individual-level (e.g., determining high-risk treatments based on the probability) and group-level (e.g., allocating clinical resources) decisions. Some prior research has sought to improve calibration by integrating a calibration-specific loss during optimization [5, 6, 4] . However, these often produce models with poor discrimination [7, 8] , limiting their utility in scenarios where precise pairwise decisions are critical.Furthermore, previous studies have typically addressed calibration in a marginal sense -i.e., assessing whether probabilities align with the actual distribution across the entire population. However, for many applications, marginal calibration may be inadequate -we often require that predictions are correctly calibrated, conditional on any combination of features. This can be helpful for making more precise clinical decisions for individuals and groups. For example, when treating an overweight male, a doctor might decide on cardiovascular surgery using a model calibrated for both overweight and male. Note this might lead to a different decision that one based on a model that was calibrated for all patients. Similarly, a hospice institution may want to allocate nursing care based on a model that generates calibrated predictions for elderly individuals. This also aligns with the fairness perspective [9] , where clinical decision systems should guarantee equalized calibration performance across any protected groups.Contributions To overcome these challenges, we introduce the CSD-iPOT framework, a postprocessing approach built upon conformal prediction [10] [11] [12] 8 ] that uses the Individual survival Probability at Observed Time (iPOT) as conformity scores and generates conformalized survival distributions. The method has 3 important properties: (i) this conformity score naturally conforms to the definition in distribution calibration in survival analysis [3] ; (ii) it also captures the distribution variance of the ISD, therefore is adaptive to the features; and (iii) the method is computationally friendly for survival analysis models. Our key contributions are:• Motivating the use of conditional distribution calibration in survival analysis, and proposing a metric (Cal ws , defined in Section 4) to evaluate this property. • Developing the CSD-iPOT framework, to accommodate censorship. The method effectively solves some issues of previous conformal methods wrt inaccurate Kaplan-Meier estimation; • Theoretically proving that CSD-iPOT asymptotically guarantees marginal and conditional distribution calibration under some specified assumptions; • Conducting extensive experiments across 15 datasets, showing that CSD-iPOT improves both marginal and conditional distribution calibration without sacrificing discriminative ability; • Demonstrating that CSD-iPOT is computationally more efficient than prior conformal method on survival analysis.2 Problem statement and Related Work 
notation: 
 A survival dataset D = {(x i , t i , δ i )} n i=1 contains n tuples, each containing covariates x i ∈ R d , an observed time t i ∈ R + , and an event indicator δ i ∈ {0, 1}. For each subject, there are two potential times of interest: the event time e i and the censoring time c i . However, only the earlier of the two is observable. We assign t i ≜ min{e i , c i } and δ i ≜ 1[e i ≤ c i ], so δ i = 0 means the event has not happened by t i (right-censored) and δ i = 1 indicates the event occurred at t i (uncensored). Let I denote the set of indices in dataset D, then we can use i ∈ I to represent (x i , t i , δ i ) ∈ D.Our objective is to estimate the Individualized Survival Distribution (ISD), S(t | x i ) = P(e i > t | X = x i ), which represents the survival probabilities of the i-th subject for any time t ≥ 0.Calibration measures the alignment between the predictions against observations. Consider distribution calibration at the individual level: if an oracle knows the true ISD S(t | x i ), and draws realizations of e i | x i (call them e (1) i , ei , . . .), then the survival probability at observed time {S(e (m) i | x i )} m should be distributed across a standard uniform distribution U [0,1] (probability integral theorem [13] ). However, in practice, for each unique x i , there is only one realization of e i | x i , meaning we cannot check the calibration in this individual manner.To solve this, Haider et al. [3] proposed marginal calibration, which holds if the predicted survival probabilities at event times e i over the x i in the dataset, { Ŝ(e i | x i )} i∈I , matches U [0,1] . Definition 2.1. For uncensored dataset, a model has perfect marginal calibration iff ∀ [ρ 1 , ρ 2 ] ⊂ [0, 1],EQUATIONWe can "blur" each censored subject uniformly over the probability intervals after the survival probability at censored time Ŝ(c i | x i ) [3] (see the derivation in Appendix A):P ( Ŝ(e i | x i ) ∈ [ρ 1 , ρ 2 ] | δ i = 0) = ( Ŝ(t i |x i )-ρ 1 )1[ Ŝ(t i |x i )∈[ρ 1 ,ρ 2 ]]+(ρ 2 -ρ 1 )1[ Ŝ(t i |x i )≥ρ 2 ] Ŝ(t i |x i ). Figure 1 : Two notions of distribution calibration: marginal and conditional, illustrated using 3 bins separated at 1 3 and 2 3 . The curves in (a, d) represent the predicted ISDs. The colors of the stars distinguish the six subjects, with horizontal coordinates indicating the true event time (consistent across all panels) and vertical coordinates representing predicted survival probability at event time. Note the two groups (orange for x = 0 and blue for x = 1) correspond to the colors of the curves and histograms in (a, b, d, e). Note that all three P-P lines in the conditional case (f) coincide.Figure 1 (a) illustrates how marginal calibration is assessed using 6 uncensored subjects. Figure 1(b, c ) show the histograms and P-P plots, showing the predictions are marginally calibrated, over the predefined 3 bins, as we can see there are 6/3 = 2 instances in each of the 3 bins. However, if we divide the datasets into two groups (orange vs. blue -think men vs. women), we can see that this is not the case, as there is no orange instance in the [ 2 3 , 1] bin, and 2 orange instances in the [ 1 3 , 1] bin. In summary, individual calibration is ideal but impractical. Conversely, marginal calibration is more feasible but fails to assess calibration relative to certain subsets of the population by features. This discrepancy motivates us to explore a middle ground -conditional calibration. A conditionally calibrated prediction, which ensures that the predicted survival probabilities are uniformly distributed in each of these groups, as shown in Figure 1(d, e, f ), is more effective in real-world scenarios. Consider predicting employee attrition within a company: while a marginal calibration using a Kaplan-Meier (KM) [14] curve might reflect overall population trends, it fails to account for variations such as the tendency of lower-salaried employees to leave earlier. A model that is calibrated for both high and low salary levels would be more helpful for predicting the actual quitting times and facilitate planning. Similarly, when predicting the timing of death from cardiovascular diseases, models calibrated for older populations, who exhibit more predictable and less varied outcomes [15] , may not apply to younger individuals with higher outcome variability. Using age-inappropriate models could lead to inaccurate predictions, adversely affecting treatment plans.Methods based on the objective function [5, 6, 4] have been developed to enhance the marginal calibration of ISDs, involving the addition of a calibration loss to the model's original objective function (e.g., likelihood loss). However, while those methods are effective in improving the marginal calibration performance of the model, their model often significantly harms the discrimination performance [6, 4, 7] , a phenomenon known as the discrimination-calibration trade-off [7] .Post-processed methods [12, 8] have been proposed to solve this trade-off by disentangling calibration from discrimination in the optimization process. Candès et al. [12] uses the individual censoring probability as the weighting in addition to the regular Conformalized Quantile Regression (CQR) [11] method. However, their weighting method is only applicable to Type-I censoring settings where each subject must have a known censored time [16] -which is not applicable to most of the right-censoring datasets.Qi et al. [8] developed Conformalized Survival Distribution (CSD) by first discretizing the ISD curves into percentile times (via predefined percentile levels), and then applying CQR [11] for each percentile level (see the visual illustration of CSD in Figure 6 in Appendix B). Their method handles right-censoring using KM-sampling, which employs a conditional KM curve to simulate multiple event times for a censored subject, offering a calibrated approximation for the ISD based on the subject's censored time. However, their method struggles with some inherent problems of KM [14] e.g., KM can be inaccurate when the dataset contains a high proportion of censoring [17] . Furthermore, we also observed that the KM estimation often concludes at high probabilities (as seen in datasets like HFCR, FLCHAIN, and Employee in Figure 9 ). This poses a challenge in extrapolating beyond the last KM time point, which hinders the accuracy of KM-sampling, thereby constraining the efficacy of CSD (see our results in Figure 3 ).Our work is inspired by CSD [8] , and can be seen as a percentile-based refinement of their regressionbased approach. Specifically, our CSD-iPOT effectively addresses and resolves issues prevalent in the KM-sampling, significantly outperforming existing methods in terms of improving the marginal distribution calibration performance. Furthermore, to our best knowledge, this is the first approach that optimizes conditional calibration within the survival analysis that can deal with censorship.However, achieving conditional calibration (also known as conditional coverage in some literature [10] ) is challenging because it cannot be attained in a distribution-free manner for non-trivial predictions. In fact, guarantees of finite sample for conditional calibration are impossible to achieve even for standard regression datasets without censorship [10, 18, 19] . This limitation is an important topic in the statistical learning and conformal prediction literature [10, 18, 19] . Therefore, our paper does not attempt to provide finite sample guarantees. Instead, following the approach of many other researchers [11, [20] [21] [22] [23] [24] , we provide only asymptotic guarantees as the sample size approaches infinity.The key idea behind this asymptotic conditional guarantee is that the construction of post-processing predictions relies on the quality of the original predictions. Thus, we aim for conditional calibration only within the class of predictions that can be learned well -that is, consistent estimators.We acknowledge that this assumption may not hold in practice; however, (i) reliance on consistent estimators is a standard (albeit strong) assumption in the field of conformal prediction [21] [22] [23] , (ii) to the best of our knowledge, no previous results have proven conditional calibration under more relaxed conditions, and (iii) we provide empirical evidence of conditional calibration using extensive experiments (see Section 5) . 
methods: 
 This section describes our proposed method: Conformalized Survival Distribution using Individual survival Probability at Observed Time (CSD-iPOT), which is motivated by the definition of distribution calibration [3] and consists of three components: the estimation of continuous ISD prediction, the computation of suitable conformity scores (especially for censored subjects), and their conformal calibration.For simplicity, our method is motivated by the split conformal prediction [25, 11] . We start the process by splitting the instances of the training data into a proper training set D train and a conformal set D con . Then, we can use any survival algorithm or quantile regression algorithm (with the capability of handling censorship) to train a model M using D train that can make ISD predictions for D con -see Figure 2(a) .With little loss of generality, we assume that the ISD predicted by the model, ŜM (t | x i ), are right-continuous and have unbounded range, i.e., ŜM (t | x i ) > 0 for all t ≥ 0. For survival algorithms that can only generate piecewise constant survival probabilities (e.g., Cox-based methods [26, 27] , discrete-time methods [28, 29] , etc.), the continuous issue can be fixed by applying some interpolation algorithms (e.g., linear or spline). 6 provide a side-by-side visual comparison between CSD and our method.We start by sketching how CSD-iPOT deals with only uncensored subjects. Within the conformal set, for each subject i ∈ I con , we define a distributional conformity score, wrt the model M, termed the predicted Individual survival Probability at Observed Time (iPOT):EQUATIONHere, for uncensored subjects, the observed time corresponds to the event time, t i = e i . Recall from Section 2.2 that predictions from model M are marginally calibrated if the iPOT values follow U [0,1] -i.e., if we collect the distributional conformity scores for every subject in the conformal set Γ M = {γ i,M } i∈I con , the ρ-th percentile value in this set should be equal to exactly ρ. If so, no post processing adjustments are necessary.In general, of course, the estimated Individualized Survival Distributions (ISDs) Ŝ(t | x i ) may not perfectly align with the true distributions S(t | x i ) from the oracle. Therefore, for a testing subject with index n + 1, we can simply apply the following adjustment to its estimated ISD:EQUATIONHere, Percentile(ρ;Γ M ) calculates the ⌈ρ(|D con |+1)⌉ |D con |-th empirical percentile of Γ M . This adjustment aims to re-calibrate the estimated ISD based on the empirical distribution of the conformity scores.Visually, this adjustment involves three procedures:(i) It first identifies the empirical percentiles of the conformity scores -Percentile( 13 ; Γ M ) and Percentile( 23 ; Γ M ), illustrated by the two grey lines at 0.28 and 0.47 in Figure 2 (d), respectively -which uniformly divide the stars according to their vertical locations; (ii) It then determines the corresponding times on the predicted ISDs that match these empirical percentiles (the hollow circles, where each ISD crosses the horizontal line); (iii) Finally, the procedure shifts the empirical percentiles (grey lines) to the appropriate height of desired percentiles ( 1 3 and 2 3 ), along with all the circles. This operation is indicated by the vertical shifts of the hollow points, depicted with curved red arrows in Figure 2(d) .This adjustment results in the post-processed curves depicted in Figure 2(e) . It shifts the vertical position of the green star ☀ from the interval [ 1 3 , 2 3 ] to [ 2 3 , 1], and the pink star ☀ from [0, 1 3 ] to [ 1 3 , 2 3 ]. These shifts ensure that the calibration histograms and P-P plots achieve a uniform distribution (both marginally on the whole population and conditionally on each group) across the defined intervals.After generating the post-processed curves, we can apply a final step, which involves transforming the inverse ISD function back into the ISD function for the testing subject:EQUATIONThe simple visual example in Figure 2 shows only two percentiles created at 1 3 and 2 3 . In practical applications, the user provides a predefined set of percentiles, P, to adjust the ISDs. The choice of P can slightly affect the resulting survival distributions, each capable of achieving provable distribution calibration; see ablation study #2 in Appendix E.6 for how P affects the performance.It is challenging to incorporate censored instances into the analysis as we do not observe their true event times, e i , which means we cannot directly apply conformity score in (3) and the subsequent conformal steps. Instead, we only observe the censoring times, which serve as lower bounds of the event times.Given the monotonic decreasing property of the ISD curves, the iPOT value for a censored subject, i.e., ŜM (t i | x i ) = ŜM (c i | x i ), now serves as the upper bound of ŜM (e i | x i ). Therefore, given the prior knowledge that ŜM (e i | x i ) ∼ U [0,1] , the observation of the censoring time updates the possible range of this distribution. Given that ŜM (e i | x i ) must be less than or equal to ŜM (c i | x i ), the updated posterior distribution follows ŜM (ei | x i ) ∼ U [0, ŜM (ci|xi)] .Following the calibration calculation in [3] , where censored patients are evenly "blurred" across subsequent bins of Ŝ(c i | x i ), our approach uses the above posterior distribution to uniformly draw R potential conformity scores for a censored subject, for some constant R ∈ Z + . Specifically, for a censored subject, we calculate the conformity scores as:γ i,M = ŜM (c i | x i ) ⋅ u R , where u R = [ 0/R , 1/R , . . . , R/R ] .Here, u R is a pseudo-uniform vector to mimic the uniform sampling operation, significantly reducing computational overhead compared to actual uniform distribution sampling. For uncensored subjects, we also need to apply a similar sampling strategy to maintain a balanced censoring rate within the conformal set. Because the exact iPOT value is known and deterministic for uncensored subjects, sampling involves directly drawing from a degenerate distribution centered at ŜM (e i | x i ) -i.e., just drawing ŜM (e i | x i ) R times. The pseudo-code for implementing the CSD-iPOT process with censoring is outlined in Algorithm 1 in Appendix B.Note that the primary computational demand of this method stems from the optional interpolation and extrapolation of the piecewise constant ISD predictions. Calculating the conformity scores and estimating their percentiles incur negligible costs in terms of both time and space, once the rightcontinuous survival distributions are established. We provide computational analysis in Appendix E.5.Here we discuss the theoretical properties of CSD-iPOT. Unlike CSD [8] , which adjusts the ISD curves horizontally (changing the times, for a fixed percentile), our refined version scales the ISD curves vertically. This vertical adjustment leads to several advantageous properties. In particular, we highlight why our method is expected to yield superior performance in terms of marginal and conditional calibration compared to CSD [8] . Table 1 summarizes the properties of the two methods.Calibration CSD-iPOT differs from CSD in two major ways: CSD-iPOT (i) essentially samples the event time from ŜM (t | t > c i , x i ) for a censored subject, and (ii) subsequently converts these times into corresponding survival probability values on the curve. CSD [8] X X X ✓ X O(N ⋅ |P| ⋅ R) CSD-iPOT ✓ ✓ ✓ X ✓ O(N ⋅ R)The first difference contrasts with the CSD method, which samples from a conditional KM distribution, S KM ( t | t > c i ), assuming a homoskedastic survival distribution across subjects (where the conditional KM curves have the same shape and the random disturbance of e i is independent of the features x i ). However, CSD-iPOT differs by considering the heteroskedastic nature of survival distributions ŜM (t | t > c i , x i ). For instance, consider the symptom onset times following exposure to the COVID-19 virus. Older adults, who may exhibit more variable immune responses, could experience a broader range of onset times compared to younger adults, whose symptom onset times are generally more consistent [30] . By integrating this feature-dependent variability, CSD-iPOT captures the inherent heteroskedasticity of survival distributions and adjusts the survival estimates accordingly, which helps with conditional calibration.Furthermore, by transforming the times into the survival probability values on the predicted ISD curves (the second difference), we mitigate the trouble of inaccurate interpolation and extrapolation of the distribution. This approach is particularly useful when the conditional distribution terminates at a relatively high probability, where extrapolating beyond the observed range is problematic due to the lack of data for estimating the tail behavior. Different extrapolation methods, whether parametric or spline-based, can yield widely varying behaviors in the tails of the distribution, potentially leading to significant inaccuracies in survival estimates. However, by converting event times into survival percentiles, CSD-iPOT circumvents these issues. This method capitalizes on the probability integral transform [13] , which ensures that regardless of the specific tail behavior of a survival function, its inverse probability values will follow a uniform distribution.The next results state that the output of our method has asymptotic marginal calibration, with necessary assumptions (exchangeability, conditional independent censoring, and continuity). We also prove the asymptotic conditional calibrated guarantee for CSD-iPOT. The proofs of these two results are inspired by the standard conformal prediction literature [11, 22] , with adequate modifications to accommodate our method. We refer the reader to Appendix C.1 for the complete proof.Theorem 3.1 (Asymptotic marginal calibration). If the instances in D are exchangeable, and follow the conditional independent censoring assumption, then for a new instance n+ 1, ∀ ρ 1 < ρ 2 ∈ [0, 1], ρ 2 -ρ 1 ≤ P ( SM (t n+1 | x n+1 ) ∈ [ρ 1 , ρ 2 ]) ≤ ρ 2 -ρ 1 + 1 |D con | + 1 .Theorem 3.2 (Asymptotic conditional calibration). In addition to the assumptions in Theorem 3.1, if (i) the non-processed prediction ŜM (t | x i ) is a consistent survival estimator; (ii) its inverse function is differentiable; and (iii) the 1st derivation of the inverse function is bounded by a constant, then the CSD-iPOT process will achieve asymptotic conditional distribution calibration.Monotonicity Unlike CSD, CSD-iPOT does not face any non-monotonic issues for the postprocessed curves as long as the original ISD predictions are monotonic; see proof in Appendix C.2.Theorem 3.3. CSD-iPOT process preserves the monotonic decreasing property of the ISD.CSD, built on the Conformalized Quantile Regression (CQR) framework, struggles with the common issue of non-monotonic quantile curves (refer to Appendix D.2 in [8] and our Appendix C.2). While some methods, like the one proposed by Chernozhukov et al. [31] , address this issue by rearranging quantiles, they can be computationally intensive and risk (slightly) recalibrating and distorting discrimination in the rearranged curves. By inherently maintaining monotonicity, CSD-iPOT not only enhances computational efficiency but also avoids these risks.Discrimination Qi et al. [8] demonstrated that CSD theoretically guarantees the preservation of the original model's discrimination performance in terms of Harrell's concordance index (C-index) [1] . However, CSD-iPOT lacks this property; see Appendix C.3 for details.As CSD-iPOT vertically scales the ISD curves, it preserves the relative order of survival probabilities at any single time point. This preservation means that the discrimination power, measured by the area under the receiver operating characteristic (AUROC) at any time, remains intact (Theorem C.4). Furthermore, Antolini's time-dependent C-index (C td ) [32] , which represents a weighted average AUROC across all time points, is also guaranteed to be maintained by our method (Lemma C.5). As a comparison, CSD does not have such a guarantee for neither AUROC nor C td . 
evaluation metrics: 
 We measure discrimination using Harrell's C-index [1] , rather than Antolini's C td [32] , as Lemma C.5 already established that C td is not changed by CSD-iPOT. We aim to assess our performance using a measure that represents a relative weakness of our method.As to the calibration metrics, the marginal calibration score evaluated on the test set D test is calculated as [3, 8] :EQUATIONwhere 1) and ( 2); see (8) in Appendix A.P( Ŝ(e i | x i ) ∈ [0, ρ] , i ∈ I test ) is calculated by combining (Based on the marginal calibration formulation, a natural way for evaluating the conditional calibration could be: (i) heuristically define a finite feature space set {S 1 , S 2 , . . .} -e.g., S 1 is the set of divorced elder males, S 2 is females with 2 children, etc.; and (ii) calculate the worst calibration score on all the predefined sub-spaces. This is similar to fairness settings, researchers normally select age, sex, or race as the sensitive attributes to form the feature space. However, this metric does not scale to higher-dimensional settings because it is challenging to create the feature space set that contains all possible combinations of the features.Motivated by Romano et al. [33] , we proposed a worst-slab distribution calibration, Cal ws . We start by partition the testing set into a 25% exploring set D explore and a 75% exploiting set D exploit . The exploring set is then used to find the worst calibrated sub-region in the feature space R d :S v,a,b = { x i ∈ R d ∶ a ≤ v ⊺ x i ≤ b } and P ( x i ∈ S v,a,b , i ∈ I explore ) ≥ κ,where v, a, b = arg maxv∈R d ,a<b∈R 1 |P| ∑ ρ∈P (P ( Ŝ(e i | x i ) ∈ [0, ρ] , i ∈ I explore , x i ∈ S v,a,b )) -ρ) 2 .In practice, the parameters v, a, and b are chosen adversarially by sampling i.i.d. vectors v on the unit sphere in R d then finding the arg max using a grid search on the exploring set. κ is a predefined threshold to ensure that we only consider slabs that contain at least κ% of the instances (so that we do not encounter a pregnant-man situation). Given this slab, we can calculate the conditional calibration score on the evaluation set for this slab:EQUATIONBesides the above metrics, we also evaluate using other commonly used metrics: integrated Brier score (IBS) [34] , and mean absolute error with pseudo-observation (MAE-PO) [35] ; see Appendix D. 
experiments: 
 The implementation of CSD-iPOT method, worst-slab distribution calibration score, and the code to reproduce all experiments in this section are available at https://github.com/shi-ang/ MakeSurvivalCalibratedAgain.Datasets We use 15 datasets to test the effectiveness of our method. Table 3 in Appendix E.1 summarizes the dataset statistics, and Appendix E.1 also contains details of preprocessing steps, KM curves, and histograms of event/censor times. Compared with [8] , we added datasets with high censoring rates (>60%) and ones whose KM ends with high probabilities (>50%).Baselines We compared 7 survival algorithms: AFT [36] , GB [37] , DeepSurv [38] , N-MTLR [39] , DeepHit [29] , CoxTime [27] , and CQRNN [40] . We also include KM as a benchmark (empirical lower bound) for marginal calibration, which is known to achieve perfect marginal calibration [3, 8] . Appendix E.2 describes the implementation details and hyperparameter settings.Procedure We divided the data into a training set (90%) and a testing set (10%) using a stratified split to balance time t i and censor indicator δ i . We also reserved a balanced 10% validation subset from the training data for hyperparameter tuning and early stopping. This procedure was replicated across 10 random splits for each dataset.Due to space constraints, the main text presents partial results for datasets with high censoring rates and high KM ending probabilities (HFCR, FLCHAIN, Employee, MIMIC-IV). Notably, CQRNN did not converge on the MIMIC-IV dataset. Thus, we conducted a total of 104 method comparisons (15 datasets × 7 baselines -1). Table 2 offers a detailed performance summary of CSD-iPOT versus baselines and CSD method across these comparisons. Appendix E.4 presents the complete results.Discrimination The upper panels of Figure 3 indicate minimal differences in the C-index among the methods, with notable exceptions primarily involving DeepHit. Specifically, CSD-iPOT matched the baseline C-index in 75 instances and outperformed it in 7 out of 104 comparisons. This suggests that CSD-iPOT maintains discriminative performance in approximately 79% of the cases.The lower panels of Figure 3 show significant improvements in marginal calibration with CSD-iPOT. It often achieved near-optimal performance, as marked by the red dashedMethods Baseline CSD CSD-iPOT lines. Table 2 also shows that CSD-iPOT provided better marginal calibration than the baselines in 95 (and significantly in 50) out of 104 comparisons (91%).CSD-iPOT's marginal calibration was better than CSD most of the time (68/104, 65%). The cases where CSD performs better typically involve models like DeepHit or CQRNN. This shows that our approach often does not perform as well as CSD when the original model is heavily miscalibrated, which suggests a minor limitation of our method. Appendix C.4 discusses why our method is sub-optimal for these models.Conditional Calibration For small datasets (sample size < 1000), in some random split, we can find a worst-slab region S v,a,b on the exploring set with P(x i ∈ S v,a,b , i ∈ I explore ) ≥ 33% but still no subjects in this region in the exploiting set. This is probably because we only ensure that the times and censored indicators are balanced during the partition, however, the features can still be unbalanced. Therefore, we only evaluated conditional calibration on the 10 larger datasets, resulting in 69 comparisons. Among them, CSD-iPOT improved conditional calibration in 64 cases (93%) compared to baselines and in 51 cases (74%) compared to CSD.Case Study We provide 4 case studies in Figure 13 in Appendix E. 4 , where CSD leads to significant miscalibration within certain subgroups, and CSD-iPOT can effectively generate more conditional calibrated predictions in those groups. These examples show that CSD's miscalibration is always located at the low-probability regions, which corresponds to our statement (in Section 3.4) that the conditional KM sampling method that CSD used is problematic for the tail of the distribution.Other Metrics Results in Table 2 and Appendix E.4 show that CSD-iPOT also showed improvement in both IBS and MAE-PO, outperforming 63 and 54 out of 104 comparisons, respectively.Computational Analysis Appendix E.5 shows the comprehensive results and experimental setup. In summary, CSD-iPOT significantly reduces the space consumption and running time.We conducted two ablation studies to assess (i) the impact of the repetitions value (R) and (ii) the impact of predefined percentiles (P) on the method; see Appendix E.6. 
conclusions: 
 Discrimination and marginal calibration are two fundamental yet distinct elements in survival analysis. While marginal calibration is feasible, it overlooks accuracy across different groups distinguished by specific features. In this paper, we emphasize the importance of conditional calibration for practical applications and propose a principled metric for this purpose. By generating conditionally calibrated Individual Survival Distributions (ISDs), we can better communicate the uncertainty in survival analysis models, enhancing their reliability, fairness, and real-world applicability.We therefore define the Conformalized Survival Distribution using Individual Survival Probability at Observed Time (CSD-iPOT) -a post-processing framework that enhances both marginal and conditional calibration without compromising discrimination. It addresses common issues in prior methods, particularly under high censoring rates or when the Kaplan-Meier curve terminates at a high probability. Moreover, this post-processing adjusts the ISDs by adapting the heteroskedasticity of the distribution, leading to asymptotic conditional calibration. Our extensive empirical tests confirm that CSD-iPOT significantly improves both marginal and conditional performance without diminishing the models' discriminative power. Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes]Justification: Section 5 states how we split the data in a stratified manner. Further experimental details and hyperparameters information are included in the Appendix E.2 and E.3.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes]Justification: The paper uses violin plots to show the density distribution for the results from 10 random split of the experiments. We also use two-sided t-test to check whether the differences are significant (Table 2 ).Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. 
 
 Appendix: Distribution calibration (or simply "calibration") examines the calibration ability across the entire range of ISD predictions [3] . This appendix provides more details about this metric.A.1 Why the survival probability at event times should be uniform?The probability integral transform [13] states: if the conditional cumulative distribution function (CDF) F (t | x i ) is legitimate and continuous in t for each fixed value of x i , thenF (t | x i ) has a standard uniform distribution, U [0,1] . Since S(t | x i ) = 1 -F (t | x i ), then S(t | x i ) ∼ U [0,1]. Definition 2.1 defines the marginal calibration for the uncensored dataset. In this section, we expand the definition to any dataset with censored instances. Note that Haider et al. [3] proposed the following method, here we just reformulate their methodology to fit the language in this paper, for a better presentation purpose.Given an uncensored subject, the probability of its survival probability at event time in a probability interval of [ρ 1 , ρ 2 ] is deterministic, as:P ( Ŝ(e i | x i ) ∈ [ρ 1 , ρ 2 ] | δ i = 1 ) = 1 [ Ŝ(e i | x i ) ∈ [ρ 1 , ρ 2 ], δ i = 1 ] .For censored subjects, because we do not know the true event time, so there is no way we can know whether the predicted probability is within the interval or not. We can "blur" the subject uniformly to the probability intervals after the survival probability at censored time Ŝ(c i | x i ) [3] .P ( Ŝ(e i | x i ) ∈ [ρ 1 , ρ 2 ] | δ i = 0) = P ( Ŝ(e i | x i ) ∈ [ρ 1 , ρ 2 ], δ i = 0) P(δ i = 0) = P ( Ŝ(e i | x i ) ∈ [ρ 1 , ρ 2 ], Ŝ(e i | x i ) < Ŝ(c i | x i )) P( Ŝ(e i | x i ) < Ŝ(c i | x i )) = P ( Ŝ(e i | x i ) ∈ [ρ 1 , ρ 2 ], Ŝ(e i | x i ) < Ŝ(c i | x i ), Ŝ(c i | x i ) ≥ ρ 2 ) P( Ŝ(e i | x i ) < Ŝ(c i | x i )) + P ( Ŝ(e i | x i ) ∈ [ρ 1 , ρ 2 ], Ŝ(e i | x i ) < Ŝ(c i | x i ), Ŝ(c i | x i ) ∈ [ρ 1 , ρ 2 ]) P( Ŝ(e i | x i ) < Ŝ(c i | x i )) + P ( Ŝ(e i | x i ) ∈ [ρ 1 , ρ 2 ], Ŝ(e i | x i ) < Ŝ(c i | x i ), Ŝ(c i | x i ) ≤ ρ 1 ) P( Ŝ(e i | x i ) < Ŝ(c i | x i )) = P (ρ 1 ≤ Ŝ(e i | x i ) ≤ ρ 2 ) P ( Ŝ(c i | x i ) ≥ ρ 2 ) Ŝ(t i | x i ) + P (ρ 1 ≤ Ŝ(e i | x i ) ≤ Ŝ(c i | x i )) P ( Ŝ(c i | x i ) ∈ [ρ 1 , ρ 2 ]) Ŝ(t i | x i ) + P(∅) Ŝ(t i | x i ) = (ρ 2 -ρ 1 )1 [ Ŝ(t i | x i ) ≥ ρ 2 ] + ( Ŝ(t i | x i ) -ρ 1 ) 1 [ Ŝ(t i | x i ) ∈ [ρ 1 , ρ 2 ]] Ŝ(t i | x i ),where the decomposition of probability in the second to last equality is because of the conditional independent censoring assumption. Therefore, for the entire dataset, considering both uncensored and censored subjects, the probabilityEQUATIONTherefore, given the above derivation, we can provide a formal definition for marginal distributional calibration for any survival dataset with censorship. Definition A.1 (Marginal calibration). For a survival dataset, a model has perfect marginal calibration iff∀ [ρ 1 , ρ 2 ] ⊂ [0, 1], P ( Ŝ(e i | x i ) ∈ [ρ 1 , ρ 2 ] , i ∈ I) = ρ 2 -ρ 1 .where the probability P is calculated using (8) .Here we present more details for the algorithm: Conformalized Survival Distribution using Individual survival Probability at Observed Time (CSD-iPOT). The pseudo-code is presented in Algorithm 1. 6: for i ∈ I con do 7:γ i,M = ŜM (t i | x i ) ▷ calculate the iPOT value 8:if δ i == 0 then 9:Γ M = Γ M + {u ⋅ γ i,M | u ∈ u R } ▷ sample R times from U [0, ŜM (ti|xi)]10:else 11:Γ M = Γ M + {γ i,M , ⋯, γ i,M R times } ▷ repeat the iPOT value, R times 12:end if 13: end for 14: ∀ ρ ∈ P, S-1M (ρ | x n+1 ) = Ŝ-1 M ( Percentile(ρ; Γ M ) | x n+1 ) 15: SM (t | x n+1 ) = inf{ρ ∶ S-1 M (ρ | x n+1 ) ≤ t} ▷ transform the inverse ISD into a ISD curveOur method is motivated by the split conformal prediction [25] . The algorithm starts by partitioning the dataset in a training and a conformal set (line 1 in Algorithm 1). Previous methods [11, 12] recommend mutually exclusive partitioning, i.e., the training and conformal set must satisfy D train ∪ D con = D and D train ∩ D con = ∅. However, this can cause one problem: reducing the training set size, resulting in underfitting for the models (especially for deep-learning models) and sacrificing discrimination performance. Instead, we consider the two partitioning policies proposed by Qi et al. Another algorithm detail is the interpolation and extrapolation (line 4 in Algorithm 1) for the ISDs. This optional operation is required for non-parametric survival algorithms (including semiparametric [26] and discretize time/quantile models [28, 29] ). For interpolation, we can use either linear or piecewise cubic Hermit interpolating polynomial (PCHIP) [41] which both can maintain the monotonic property of the survival curves. For extrapolation, we extend ISDs using the starting point [0, 1] and ending point[t end , ŜM (t end | x i )].We also include a side-by-side visual comparison of our method to CSD [8] in Figure 6 . Both CSD and CSD-iPOT approaches use the same general framework of conformal prediction, but they differ in how to calculate the conformity score and adjust the positions of the predictions. Note that in Figure 6 , the horizontal positions of the circles remain unchanged for CSD while the vertical positions of the circles remain unchanged for CSD-iPOT. This unique approach for calculating the conformity score, and the downstream benefit for handling censored subjects (thanks to the conformity score design), together provide a theoretical guarantee for the calibration, as elaborated in the next section.This appendix offers more details on the theoretical analysis in Section 3.4.This section presents the complete proof for Theorem 3.1 and Theorem 3.2.For the completeness, we start by restating Theorem 3.1:Theorem C.1 (Asymptotic marginal calibration). If the instances in D are exchangeable, and follow the conditional independent censoring assumption, then for a new dataset nEQUATIONProof. This proof is inspired by the proof in Theorem D.1 and D.2 from Angelopoulos et al. [42] , the main difference is that our conformity scores are different, and we also have censoring subjects.Given the target percentile levels ρ 1 and ρ 2 , As we remember in the CSD-iPOT procedure, we essentially vertically shift the prediction at Percentile(ρ; Γ) to ρ (see ( 4)). Therefore, the condition in ( 9)SM (t n+1 | x n+1 ) ∈ [ρ 1 , ρ 2 ],can be transformed into this equivalent conditionŜM (t n+1 | x n+1 ) ∈ [Percentile(ρ 1 ; Γ M ), Percentile(ρ 2 ; Γ M )].As we recall, the conformity score set consists of the iPOT score for each subject in the conformal set. Let us consider the easy case of an uncensored dataset, thenΓ M = { ŜM (e i | x i )} i∈I conWithout loss of generality, we assume that the conformity scores in Γ M are sorted. This assumption is purely technical as it aims for simpler math expression in this proof, i.e., ŜM (e 1 | x 1 ) < ŜM (e 2 | x 2 ) < ⋯ < ŜM (e |D| con | x |D con | ). Therefore, by the exchangeability assumption (the order between subject n + 1 and subjects i ∈ I con do not matter), the iPOT score of subject n + 1 is equally likely to fall into any of the |D con | + 1 intervals between the [0, 1], separated by the |D con | conformity scores. Therefore,P ( ŜM (t n+1 | x n+1 ) ∈ [Percentile(ρ 1 ; Γ M ), Percentile(ρ 2 ; Γ M )]) = ⌈(ρ 2 -ρ 1 )(|D con | + 1)⌉ (|D con | + 1)This value is always higher than ρ 2ρ 1 , and if the conformity scores in Γ M do not have any tie 1 , we can also see that the above equation is less than ρ 2ρ 1 + 1 |D con |+1 . Now let's consider the censored case. Based on the probability integral transform, for a censored subject j, the probability that its iPOT value falls into the percentile interval, before knowing its censored time, is (from here we omit the subscript M for simple expression)P ( Ŝ(e j | x j ) ≤ ρ) = ρ, and P (ρ 1 ≤ Ŝ(e j | x j ) ≤ ρ 2 ) = ρ 2 -ρ 1Now given its censored time c j , we can calculate the conditional probabilityP (ρ 1 ≤ Ŝ(e j | x j ) ≤ ρ 2 | e j > c i ) = P (ρ 1 ≤ Ŝ(e j | x j ) ≤ ρ 2 | Ŝ(e j | x j ) ≤ Ŝ(c j | x j )) = P (ρ 1 ≤ Ŝ(e j | x j ) ≤ ρ 2 , Ŝ(e j | x j ) ≤ Ŝ(c j | x j )) P ( Ŝ(e j | x j ) ≤ Ŝ(c j | x j )) = P (ρ 1 ≤ Ŝ(e j | x j ) ≤ ρ 2 , Ŝ(e j | x j ) ≤ Ŝ(c j | x j ), Ŝ(c j | x j ) < ρ 1 ) Ŝ(c j | x j ) + P (ρ 1 ≤ Ŝ(e j | x j ) ≤ ρ 2 , Ŝ(e j | x j ) ≤ Ŝ(c j | x j ), ρ 1 ≤ Ŝ(c j | x j ) ≤ ρ 2 ) Ŝ(c j | x j ) + P (ρ 1 ≤ Ŝ(e j | x j ) ≤ ρ 2 , Ŝ(e j | x j ) ≤ Ŝ(c j | x j ), Ŝ(c j | x j ) > ρ 2 ) Ŝ(c j | x j ) = 0 Ŝ(c j | x j ) + P (ρ 1 ≤ Ŝ(e j | x j ) ≤ Ŝ(c j | x j ), ρ 1 ≤ Ŝ(c j | x j ) ≤ ρ 2 ) Ŝ(c j | x j ) + P (ρ 1 ≤ Ŝ(e j | x j ) ≤ ρ 2 , Ŝ(c j | x j ) > ρ 2 ) Ŝ(c j | x j ) = ( Ŝ(c j | x j ) -ρ 1 ) 1 [ρ 1 ≤ Ŝ(c j | x j ) ≤ ρ 2 ] + (ρ 2 -ρ 1 )1 [ Ŝ(c j | x j ) > ρ 2 ] Ŝ(c j | x j ) ,where the probability decomposition in the last equality is because the conditional independent censoring assumption, i.e., e j c j | x j . This above derivation means the Ŝ(e j | x j ) follows the uniform distribution U [0, Ŝ(ci|xj )] . Therefore, if we do one sampling for each censored subject using U [0, Ŝ(cj |xj )] , the above proof asymptotically converges to the upper and lower bounds for uncensored subjects, for any survival dataset with censorship.For the conditional calibration, we start by formally restating Theorem 3.2. Theorem C.2 (Asymptotic conditional calibration). With the conditional independent censoring assumption, if we have these additional three assumptions:(i) the non-processed prediction Ŝ(t | x i ) is also a consistent survival estimator with:EQUATION) (ii) the inverse ISD estimation Ŝ-1 (t | x i ) is differentiable, (iii) there exist some M such that inf ρ d Ŝ-1 (ρ|xn+1) d ρ ≥ M -1 ,then the CSD-iPOT process can asymptotically achieve the conditional distribution calibration:S(t | x n+1 ) = S(t | x n+1 ) + o p (1).Proof. In order to prove this theorem, it is enough to show that (i)S(t | x n+1 ) = Ŝ(t | x n+1 ) + o p (1), and (ii) Ŝ(t | x n+1 ) = S(t | x n+1 ) + o p (1).The second equality is obvious to see under the consistent survival estimator assumption (10) . Now let's focus on the first equality. We borrow the idea from Lemma 5.2 and Lemma 5.3 in Izbicki et al. [22] . It states under the consistent survival estimator assumption, we can haveŜ-1 (Percentile(ρ; Γ) | x n+1 ) = Ŝ-1 (ρ | x n+1 ) + o p (1).The only difference between the above claim and the original claim in Izbicki et al. [22] , is that they use the CDF while we use the survival function, i.e., the complement of the CDF.According to (4), the above claim can be translate to S-1(ρ | x n+1 ) = Ŝ-1 (ρ | x n+1 ) + o p (1).Then if | S-1 (ρ | x n+1 ) -Ŝ-1 (ρ | x n+1 )| = o p (1), if the Ŝ is differentiable, we can have | S(t | x n+1 ) -Ŝ(t | x n+1 )| = o p (1) (inf ρ d Ŝ-1 (ρ | x i ) d ρ ) -1.Therefore, if there exist some M such that infρ d Ŝ-1 (ρ|xn+1) d ρ ≥ M -1 .We can get the first equality proved.For the censored subjects, we can use the same reasoning and steps in Theorem 3.1 under the conditional independent censoring assumption. This will finish the proof.We are aware that linear interpolation and extrapolation will make the survival curves nondifferentiable. Therefore, we recommend using PCHIP interpolation [41] for CSD-iPOT. And extrapolation normally does not need to apply for the survival prediction algorithm because the iPOT value ŜM (t i | x i ) can be obtained within the curve range for those methods. However, for quantile-based algorithms, sometimes we need extrapolation. And we recognize this as a future direction to improve.CSD is a conformalized quantile regression based [11] method. It first discretized the curves into a quantile curve, and adjusted the quantile curve at every discretized level [8] . Both the estimated quantile curves (q(ρ)) and the adjustment terms (adj(ρ)) are monotonically increasing with respect to the quantile levels. However, the adjusted quantile curve -calculated as the original quantile curve minus the adjustment -is no longer monotonic. For example, if q(50%) = 5 and q(60%) = 6, with corresponding adjustments of adj(50%) = 2 and adj(60%) = 4, the post-CSD quantile curve will be 5 -2 at 50% and 6 -4 at 60%, demonstrating non-monotonicity. For a detailed description of their algorithm, readers are referred to Qi et al. [8] .However, CSD-iPOT has this nice property. Here we restate the Theorem C.4 Theorem C.3. CSD-iPOT process preserves the monotonic decreasing property of the ISD, s.t.,EQUATIONProof. The proof of this theorem is straightforward. The essence of the proof lies in the monotonic nature of all operations within CSD-iPOT.First of all, the percentile operation Percentile(ρ; Γ M ) is a monotonic function, i.e., for all ρ 1 < ρ 2 , Percentile(ρ 1 ; Γ M ) < Percentile(ρ 2 ; Γ M ).Second, because the non-post-processed ISD curves Ŝ(t | x i ) are monotonic. Therefore, the inverse survival function Ŝ-1 (ρ | x i ) is also monotonic. Therefore, after the adjustment step as detailed in ( 4), ( 4), for all ρ 1 < ρ 2 , it follows that S-1M (ρ 1 | x n+1 ) < S-1 M (ρ 2 | x n+1 ).S-1 M (ρ | x n+1 ) back to survival function SM (t | x n+1), the monotonicity is preserved. These steps collectively affirm the theorem's proof through the intrinsic monotonicity of the operations involved in CSD-iPOT.This section explores the discrimination performance of CSD-iPOT in the context of survival analysis. Discrimination performance, which is crucial for evaluating the effectiveness of survival models, is typically assessed using three key metrics:• Area under the receiver operating characteristic curve (AUROC)We will analyze CSD-iPOT's performance across these metrics and compare it to the performance of CSD. The comparative analysis aims to highlight any improvements and trade-offs introduced by the CSD-iPOT methodology.Harrell's C-index C-index is calculated as the proportion of all comparable subject pairs whose predicted and outcome orders are concordant, defined asEQUATIONwhere ηi denotes the model's predicted risk score of subject i, which can be defined as the negative of predicted mean/median survival time(E t [ Ŝ(t | x i )] or Ŝ-1 (0.5 | x i )).CSD has been demonstrated to preserve the discrimination performance of baseline survival models, as established in Theorem 3.1 byin Qi et al. [8] . In contrast, CSD-iPOT does not retain this property.To illustrate this, we present a counterexample that explains why Harrell's C-index may not be maintained when using CSD-iPOT. As shown in Figure 7 (a), two ISD curves cross at a certain percentile level 25% < ρ * < 50%. Initially, the order of median survival times (where the curves cross 50%) for these curves indicates that patient A precedes patient B. However, after applying the adjustment as defined in (4) -which involves vertically shifting the prediction from empirical percentile level Percentile(50%; Γ M ) to desired percentile level 50%. The order of the post-processed median times for the two curves (indicated by the hollow circles) is patient B ahead of patient A. That means this adjustment leads to a reversal in the order of risk scores, thereby compromising the C-index.AUROC The area under the receiver operating characteristic curve (AUROC) is a widely recognized metric for evaluating discrimination in binary predictions. Harrell's C-index can be viewed as a special case of the AUROC [2] , if we use the negative of survival probability at a specified time t * -Ŝ(t * | x i ) as the risk score.The primary distinction lies in the definition of comparable pairs. In Harrell's C-index, comparable pairs are those for which the event order is unequivocally determined. Conversely, for the AUROC evaluation at time t * , comparable pairs are defined as one subject experiencing an event before t * and another experiencing it after t * . This implies that for AUROC at t * , a pair of uncensored subjects both having event times before (or both after) t * , is not considered comparable, whereas for the C-index, such a pair is indeed considered comparable.The AUROC can be calculated using:EQUATIONFrom this equation, because the part of δ i ⋅ 1[t i ≥ t * ] ⋅ 1[t j < t * ] is independent of the prediction (it only relates to the dataset labels). As long as the post-processing does not change the order of the survival probabilities, we can maintain the same AUROC score.Theorem C.4. Applying the CSD-iPOT adjustment to the ISD prediction does not affect the relative order of the survival probabilities at any single time, therefore does not affect the AUROC score of the model. Formally, ∀ i, j ∈ I and ∀ t * ∈ R + , givenŜM (t * | x i ) < ŜM (t * | x j ),we must haveSM (t * | x i ) < SM (t * | x j ).Here SM is calculated using 5 in Section 3.2.Proof. The intuition is that if we scale the ISD curves vertically. Then the vertical order of the ISD curves at every time point should not be changed.Formally, we first represent SM (t * | x i ) by ρ * i , and represent SM (t * | x j ) by ρ * j . Then, by applying (5) and then (4), we can haveS-1 M (ρ * i | x i ) = Ŝ-1 M (Percentile(ρ * i ; Γ M ) | x i ) S-1 M (ρ * j | x j ) = Ŝ-1 M (Percentile(ρ * j ; Γ M ) | x j )where Percentile(ρ * i ; Γ M ) is the original predicted probability at t * , i.e., ŜM (t * | x i ) = Percentile(ρ * i ; Γ M ). Because the Percentile operation and inverse functions are monotonic (Theorem 3.3), therefore, this theorem holds.CSD adjusts the survival curves horizontally (e.g., along the time axis). Hence, while the horizontal order of median/mean survival times does not change -as proved in the Theorem 3.1 from Qi et al. [8] -the vertical order, represented by survival probabilities, might not be preserved by CSD.Let's use a counter-example to illustrate our point. Figure 7 (b) shows two ISD predictions Ŝ(t | x A ) and Ŝ(t | x B ) for subjects A and B. Suppose the two ISD curves both have the median survival time at t = 100, and the two curves only cross once. Without loss of generality, we assume Ŝ(t * | x A ) < Ŝ(t * | x B ) holds for all t * < 100 and Ŝ(t * | x A ) > Ŝ(t * | x B ) holds for all t * > 100. Now, suppose that CSD modified the median survival time from t = 100 to t = 150 for both of the predictions. Then the order between these two predictions at any time in the range of t * ∈ [100, 150] is changed fromŜ(t * | x A ) > Ŝ(t * | x B ) to Ŝ(t * | x A ) < Ŝ(t * | x B ).It is worth mentioning that in Figure 2 , a blue curve is partially at the top in (a), intersecting an orange curve around 1.7 days, while the orange curve is consistently at the top in (e). This might raise concerns that the pre-and post-adjustment curves do not maintain the same probability ordering at every time point, suggesting a potential violation. In fact, this discrepancy arises from the discretization step used in our process, which did not capture the curve crossing at 1.7 days due to the limited number of percentile levels (2 levels at 1 3 and 2 3 ) used for simplicity in this visualization. The post-discretization positioning of the orange curve above the blue curve in Figure 2 (e) does not imply that the post-processing step alters the relative ordering of subjects. Instead, it reflects the limitations of using only fewer percentile levels. Note that other crossings, such as those at approximately 1.5 and 2.0 days, are captured. In practice, we typically employ more percentile levels (e.g., 9, 19, 39, or 49 as in Ablation Study #2 -see Appendix E.6), which allows for a more precise capture of all curve crossings, thereby preserving the relative ordering.Antolini's C-index Time-dependent C-index, C td , is a modified version of Harrell's C-index [32] . Instead of estimating the discrimination over the point predictions, it estimates the discrimination over the entire curve.C td ( Ŝ) = ∑ i,j∈I test δ i ⋅ 1[t i < t j ] ⋅ 1[ Ŝ(t i | x i ) < Ŝ(t i | x j )] ∑ i,j∈I test δ i ⋅ 1[t i < t j ].Compared with (12) , the only two differences are: the risk score for the earlier subject n i is represented by the iPOT value Ŝ(t i | x i ), and the risk score for the later subject n j is represented by the predicted survival probability for j at t i , Ŝ(t i | x j ). C td can also be represented by the weighted average of AUROC over all time points (Equation 7 and proof in Appendix 1 from Antolini et al. [32] ). Suppose given a series of time grid {t 0 , . . . , t k , . . . , t K },EQUATIONwhereEQUATIONThe physical meaning of ω(t k ) measures the proportion of comparable pairs at t k over all possible pairs. Therefore, we can have the following important property of our CSD-iPOT. Lemma C.5. Applying the CSD-iPOT adjustment to the ISD prediction does not affect the timedependent C-index of the model.Proof. Because the ω(t k ) in 15 is independent of the prediction, and also given Theorem C.4, AUROC( Ŝ, t k ) = AUROC( S, t k ), we can easily have C td ( Ŝ) = C td ( S) from the formula in 14.As in the previous discussion, CSD does not preserve the vertical order represented by survival probabilities. Therefore, it is natural to see that CSD also does not preserve the C td performance.Compared to CSD, our CSD-iPOT method exhibits weaker performance on the DeepHit baselines. This section explores the three potential reasons for this disparity.First of all, DeepHit tends to struggle with calibration for datasets that have high KM ending probability [43] and has poor calibration compared to other baselines for most datasets [7, 8] . This is because the DeepHit formulation assumes that, by the end of the predefined t max , every individual must already have had the event. Hence, this formulation incorrectly estimates the true underlying survival distribution (often overestimates the risks) for individuals who might survive beyond t max .Furthermore, apart from the standard likelihood loss, DeepHit also contains a ranking loss term that changes the undifferentiated indicator function in the C-index calculation in (12) with an exponential decay function. This modification potentially enhances the model's discrimination power but compromises its calibration.Lastly, Figure 8 (a) shows an example prediction using DeepHit on the FLCHAIN (72.48% censoring rate with KM curve ends at 68.16%). The solids curves represent the ISD prediction from DeepHit for 4 randomly selected subjects in the test set. And the dashed green curve represents the KM curve for the entire test set. It is evident that DeepHit tends to overestimate the subjects' risk scores (or underestimate the survival probabilities), see Figure 8(d) . Specifically, at the last time point (t = 5215), KM predicts that most of the instances (68.16%) should survive beyond this time point. However, the ISD predictions from DeepHit show everyone must die by this last point ( ŜDeepHit (5215 | x i ) = 0 for all x i -see Figure 8(a) ). This clearly violates the unbounded range assumption proposed in Section 3.1, which assumes ŜM (t | x i ) > 0 for all t ≥ 0. This violation is the main reason why CSD-iPOT exhibits weaker performance on the DeepHit baseline.CSD can effectively solve this overestimate issue (Figure 8 (b)), as it shift the curves horizontally, i.e., no upper limit for right-hand side for shifting. CSD-iPOT, on the other hand, scale the curves vertically. In such a case, the scaling must be performed within the percentile [0, 1]. Furthermore, CSD-iPOT does not have any intervention for the starting and ending probability (ρ = 1 and ρ = 0) of the curves. So no matter how the post-process changes the percentile in the middle of the curves, the starting and ending points should not be changed, just like the curves in Figure 8 (c), whereas the earlier parts of the curve are similar as CSD's, the last parts gradually drop to 0.Consequently, while CSD-iPOT significantly improves upon DeepHit, as shown in Figure 8 (d), it still underperforms compared to CSD when dealing with models that are notably miscalibrated like DeepHit.We use Harrell's C-index [1] for evaluating discrimination performance. The formula is presented in (12) . Because we are dealing with a model that may not have proportional hazard assumption, therefore, as recommended [2] , we use the negative value of the predicted median survival time as the risk score, i.e., ηi = Ŝ-1 (0.5 | x i ).The calculation of marginal calibration for a censored dataset is presented in Appendix A and calculated using ( 8) and (6) .Conditional calibration, Cal ws , is estimated using (7) . Here we present more details on the implementation. First of all, the evaluation involves further partitioning the testing set into exploring and exploiting datasets. Note that this partition does not need to be stratified (wrt to time t i and event indicator δ i ). Furthermore, for the vectors v, we sampling M i.i.d. vectors on the unit sphere in R d . Generally, we want to select a high value for M to enable all possible exploration. For small or medium datasets, we use M = 1000. However, due to the computational complexity, for large datasets, we gradually decrease the value of M ∈ [100, 1000] to get an acceptable evaluating time.We set κ = 33% in (7) for finding the S v,a,b , that means we want to find a worst-slab that contains a least 33% of the subjects in the testing set.Integrated Brier score (IBS) measures the accuracy of the predicted probabilities over all times. IBS for survival prediction is typically defined as the integral of Brier scores (BS) over time points:IBS( Ŝ; t max ) = 1 t max ⋅ ∫ tmax 0 BS(t) dt, = 1 |I test | ∑ i∈I test 1 t max ⋅ ∫ tmax 0 ( δ i ⋅ 1[t i ≤ t] ⋅ S(t | x i ) 2 G(t i ) + 1[t i > t] ⋅ (1 -S(t | x i )) 2 G(t) ) dt,where G(t) is the non-censoring probability at time t. It is estimated with KM on the censoring distribution (flip the event indicator of data), and its reciprocal 1 G(t) is referred to as the inverse probability censoring weights (IPCW). t max is defined as the maximum event time of the combined training and validation datasets. Mean absolute error calculates the time-to-event precision, i.e., the average error of predicted times and true times. Here we use MAE-pseudo observation (MAE-PO) [35] for handling censorship in the calculation.MAE PO ({ ti } i∈I ) = 1 ∑ i∈I test ω i ∑ i∈I test ω i × |(1 -δ i ) ⋅ e PO (t i , I test ) + δ i ⋅ t i -ti | , where e PO (t i , I test ) = { N × E t [S KM(I test ) (t)] -(N -1) × E t [S KM(I test -i) (t)] if δ i = 0, t i otherwise,andω i = { 1 -δ i ⋅ S KM(I test ) (t i ) if δ i = 0, 1 otherwise.Here S KM(I test ) (t) represents the population level KM curve estimated on the entire testing set I test , and S KM(I test -i) (t) represent the KM curves estimated on all the test subjects but exclude subject i.C-index, Cal margin (also called D-cal), ISB, and MAE-PO are implemented in the SurvivalEVAL package [44] . For Cal ws , please see our Python code for implementation.We provide a brief overview of the datasets used in our experiments.In this study, we evaluate the effectiveness of CSD-iPOT across 15 datasets. Table 3 summarizes the data statistics. Compared to the datasets used in [8] , we have added HFCR, WHAS, PdM, Churn, FLCHAIN, Employee, and MIMIC-IV. Specifically, we use the original GBSG dataset, as opposed to the modified version by Katzman et al. [38] used in [8] , which has a higher censoring rate and more features. For the rest of the datasets, we employ the same preprocessing methods as Qi et al. [8] -see Appendix E of their paper for details about these datasets. Below, we describe these newly added datasets:Table 3 : Key statistics of the datasets. We categorize datasets into small, medium, and large, based on the number of instances, using thresholds of 1,000 and 10,000 instances. The bolded number represents datasets with a high percentage of censorship (≥ 60%) or its KM estimation ends at a high probability (≥ 50%). Numbers in parentheses indicate the number of features after one-hot encoding. Heart Failure Clinical Record dataset (HFCR) [45] contains medical records of 299 patients with heart failure, aiming to predict mortality from left ventricular systolic dysfunction. This dataset can be downloaded from UCI Machine Learning Repository [46] .Worcester Heart Attack Study dataset (WHAS) [49] contains 500 patients with acute myocardial infarction, focusing on the time to death post-hospital admission. The data was already post-processed and can be downloaded from the scikit-survival package [56] .Predictive Maintenance (PdM) contains information on 1000 equipment failures. The goal is to predict the time to equipment failure and therefore help alert the maintenance team to prevent that failure. It includes 5 features that describe the pressure, moisture, temperature, team information (the team who is running this equipment), and equipment manufacturer. We apply one-hot encoding on the team information and equipment manufacturer features. The dataset can be downloaded from the PySurvival package [52] .The customer churn prediction dataset (Churn) focuses on predicting customer attrition. We apply one-hot encoding on the US region feature and exclude subjects who are censored at time 0. The dataset can be downloaded from the PySurvival package [52] .Serum Free Light Chain dataset (FLCHAIN) is a stratified random sample containing half of the subjects from a study on the relationship between serum free light chain (FLC) and mortality [53] . This dataset is available in R's survival package [57] . Upon downloading, we apply a few preprocessing steps. First, we remove the three subjects with events at time zero. We impute missing values for the "creatinine" feature using the median of this feature. Additionally, we eliminate the chapter feature (a disease description for the cause of death by chapter headings of the ICD code) because this feature is only available for deceased (uncensored) subjects -hence, knowing this feature will be equivalent to leaking the event indicator label to the model.Employee dataset contains employee activity information that can used to predict when an employee will quit. The dataset can be downloaded from the PySurvival package [52] . It contains duplicate entries; after dropping these duplicates, the number of subjects in the dataset is reduced from 14,999 to 11,991. We also apply one-hot encoding to the department information.MIMIC-IV database [58] provides critical care data information for patients within the hospital. We focus on a cohort of all-cause mortality data curated by [35] , featuring patients who survived at least 24 hours post-ICU admission. The event of interest, death, is derived from hospital records (during hospital stay) or state records (after discharge). The features are laboratory measurements within the first 24 hours after ICU admission.German Breast Cancer Study Group (GBSG) [51] contains 686 patients with node-positive breast cancer, complete with prognostic variables. This dataset is available in R's survival package [57] .While the original GBSG offers a higher rate of censoring and more features, [8] utilized a modified version of GBSG from [38] , merged with uncensored portions of the Rotterdam dataset, resulting in fewer features, a lower censor rate, and a larger sample size.Figure 9 shows the Kaplan-Meier (KM) estimation (blue curves) for all 15 datasets, alongside the event and censored histograms (represented by green and orange bars, respectively) in a stacked manner, where the number of bins is determined by the Sturges formula: ⌈log(|D|) + 1⌉.In this section, we detail the implementation of the seven baseline models used in our experiments, consistent with Qi et al. [8] .• Accelerate Failure Time (AFT) [36] with Weibull distribution is a linear parametric model that uses a small l 2 penalty on parameters during optimization. It is implemented in lifelines packages [59] .• Gradient Boosting Cox model (GB) [37] is an ensemble method that employs 100 boosting stages with a partial likelihood loss [60] for optimization and 100% subsampling for fitting each base learner. The model is implemented in scikit-survival packages [56] .• Neural Multi-Task Logistic Regression (N-MTLR) [39] is a discrete-time model which is an NN-extension of the linear multi-task logistic regression model (MTLR) [28] . The number of discrete times is determined by the square root of the number of uncensored patients. We use quantiles to divide those uncensored instances evenly into each time interval, as suggested in [61, 3] . We utilize the N-MTLR code provided in Qi et al. [8] .• DeepSurv [38] is a NN-extension of the Cox proportional hazard model (CoxPH) [26] . To make ISD prediction, we use the Breslow method [62] to estimate the population-level baseline hazard function. We utilize the DeepSurv code provided in Qi et al. [8] . • CoxTime [27] is a non-proportional neural network extension of the CoxPH. The model is implemented in pycox packages [27] .• Censored Quantile Regression Neural Network (CQRNN) [40] is a quantile regression-based method. We add the bootstrap-rearranging post-processing [31] to correct non-monotonic predictions. We use the CQRNN code provided in Qi et al. [8] .Full hyperparameter details for NN-Based survival baselines In the experiments, all neural network-based methods (including N-MTLR, DeepSurv, DeepHit, CoxTime, and CQRNN) used the same architecture and optimization procedure.• Training maximum epoch: 10000• Early stop patients: 50 In this section, we present the comprehensive results from our primary experiment, which focuses on evaluating the performance of CSD-iPOT compared to the original non-post-processed baselines and CSD.Note that for the MIMIC-IV datasets, CQRNN fails to converge with any hyperparameter setting possibly due to the extremely skewed distribution. As illustrated in Figure 9 , 80% of event and censoring times happen within the first bin, and the distribution exhibits long tails extending to the 17th bin.Discrimination In Figure 10 , we demonstrate the discrimination performance of CSD-iPOT compared to benchmark methods, as measured by C-index. The panels are ordered by dataset size, from smallest to largest. In each panel, the performance of the non-post-processed baselines is shown with blue bars, CSD with orange bars, and CSD-iPOT with green bars. We can see from the figure that the three methods exhibit basically the same performance across all datasets. Indeed, as we can see from the summary in Table 2 , CSD-iPOT ties with the baselines in 75 out of 104 times. In the remaining 29 times they do not tie, none of them are significantly different from each other. For the 22 times CSD-iPOT underperforms, most of them are wrt to DeepHit or CQRNN baselines (e.g., DeepHit for HFCR, PBC, GBSG, etc.). And CSD-iPOT can even outperform 7 times (e.g., GB in HFCR). Marginal calibration In Figure 11 , we present the marginal calibration performance of CSD-iPOT versus the benchmark models. The arrangement and color schemes of the panels are consistent with those used previously. For marginal calibration evaluation, we add a "dummy" model -Kaplan-Meier (KM) curve (depicted by red dashed lines in each panel) -to serve as the empirical lower limit. We calculate each KM curve using the training set and apply it identically to all test samples. It is called a "dummy" because it lacks the ability to discriminate between individuals. However, it asymptotically achieves perfect marginal calibration (see Appendix B in [8] ).Our results in Figure 11 indicate a significant improvement in calibration performance with CSD-iPOT over both baselines and CSD. Overall, our method outperforms the baselines in 95 out of 104 times (Table 2 ). In the remaining 9 times where it does not outperform the baseline, only once is the difference statistically significant, and the marginal calibration score in this single case (GB for SEER-liver) is still close to the empirical lower bound (red dashed line).For datasets characterized by higher censoring rates or high KM ending probabilities (HFCR, PBC, PdM, FLCHAIN, Employee, MIMIC-IV), our method shows superior performance. For those datasets, CSD tends to produce non-calibrated predictions compared with baselines. This outcome likely stems from the inaccuracies of the KM-sampling method under conditions of high censor rates or ending probabilities, as discussed in Section 2. Nonetheless, our approach still improves the marginal calibration scores for the baselines under these circumstances.Conditional calibration Figure 12 showcases the conditional calibration performance of CSD-iPOT versus the benchmarks, evaluated using Cal ws . Compared with marginal calibration, we cannot use any method to establish the empirical upper bound. One might think that once we establish the worst slab, we can calculate the KM curve on this slab and then use it as the empirical upper bound for the conditional calibration. However, identifying a universal worst slab across all the models is impractical. That means different models exhibit varying worst-slabs. e.g., N-MTLR can have the worst calibration for overweighted males, while DeepSurv might perform relatively good calibration for this group but poorly for disabled cardiovascular patients.The results in Figure 11 indicate a significant improvement in conditional calibration performance using our method over both baselines and CSD. Overall, our method improves the conditional calibration performance of baselines 64 out of 69 times, with significant improvements 29 times out of 64 (Table 2 ).Our method outperforms CSD in 51 out of 69 cases. Most instances where our method underperforms are relative to the DeepHit and CQRNN baselines (the reasons are explained in Appendix C.4). To evaluate the practical benefits of CSD-iPOT over CSD, we present four case studies in Figure 13 . The figure showcases 4 concrete examples where CSD (orange) leads to significant miscalibration within certain subgroups (i.e., elderly patients, women, high-salary, and non-white-racial), but CSD-iPOT (green) can effectively generate more conditional calibrated predictions which are closer to the optimal line. Moreover, all four examples illustrate that the miscalibration of CSD consistently occurs in low-probability regions, corroborating our assertion in Section 3.4 that the conditional Kaplan-Meier sampling method employed by CSD is problematic if the tail of the distribution is unknown. IBS Figure 14 illustrates the IBS performance of CSD-iPOT versus the benchmarks. According to DeGroot and Fienberg [63] , BS can be decomposed into a calibration part and a discrimination part. This implies that IBS, an integrated version of BS, assesses aspects of calibration. Our results in Figure 14 and Table 2 show that our method improves the IBS score in most of the cases (63 wins, 18 ties, and 23 losses).MAE-PO Our results in Figure 15 and Table 2 show that our method improves the MAE-PO in general (54 wins, 33 ties, and 17 losses). Space complexity Although most of the computational cost arises from ISD interpolation and extrapolation, most of the memory cost of our method stems from storing the conformity scores into an array and performing the Percentile(⋅) operation.Let's reuse the symbol N = |D con | as the number of subjects in the conformal set, P = {ρ 1 , ρ 2 , . . .} is the predefined discretized percentiles, so that |P| be the number of predefined percentiles. And R is the repetition parameter for KM-sampling.CSD [8] process first discretized the ISD curves into percentile times (PCTs) using P. Then it calculates a conformity score for each individual at every percentile level ρ (i.e., each individual will contribute for |P| conformity scores). Lastly, the "KM-sampling" process involves repeating each individual by R times. Therefore, the memory complexity of CSD is O(N ⋅ |P| ⋅ R). In contrast, our CSD-iPOT method only calculates one iPOT score (as the conformity score) for each duplicated individual. After sampling, the total memory complexity of CSD-iPOT is O(N ⋅ R).Let consider an example of using SEER-stomach dataset in our main experiment in Appendix E. Other memory costs, e.g., storing features and ISD predictions, incur negligible memory costs. This is because the number of feature d, and the length of the ISDs are much smaller than the repeat parameters R.The primary sources of time complexity in CSD-iPOT are two-fold:• ISD interpolation and extrapolation (line 4 in Algorithm 1); • An optional monotonic step for the CQRNN model (see discussion below).Note that other time complexity, e.g., running the Percentile operation or adjusting the ISD curves using lines 14-15 in Algorithm 1, incur negligible time cost.Here we analyze two kinds of time complexity: training complexity and inference complexity. In Figure 16 , we see that CSD-iPOT can significantly reduce the additional time imposed by the CSD step on all survival analysis models (AFT, GB, N-MTLR, DeepSurv, DeepHit, and CoxTime). The only exceptional is the quantile-based method CQRNN, where for the 4 small datasets (HFCR, PBC, WHAS, and GBM) CSD-iPOT actually increase the training time. This is because, for those datasets, the quantile curves predicted by CQRNN are not monotonic. To address this, we attempted three methods before we directly apply CSD-iPOT, including ceiling, flooring, and bootstrap rearranging Chernozhukov et al. [31] -with the bootstrap method proving most effective, albeit at a significant computational cost.Similarly, each point in Figure 17 represents an average inference time of the method, where for the non-post-processed baselines it is purely the inference time of the ISD predictions from the survival model, while for CSD and CSD-iPOT, it is the inference time of the ISDs plus the post-processing time. We observe the inference time follows the same trend as the training time. The extra cost for the 4 small datasets (HFCR, PBC, WHAS, and GBM) is still due to the non-monotonic issues. In those cases, the predicted curves by CQRNN become monotonic after applying CSD. However, surprisingly, after applying CSD-iPOT, the curve remains non-monotonic. That is why the inference time for baselines and CSD-iPOT is higher than CSD on those four datasets.Ablation Study #1: impact of repetition parameter R As we proved in Theorem 3.1, for a censored subject j, if we sample its iPOT value using U 0, ŜM (ci|xi) , we can asymptotically achieve the exact marginal calibration. However, empirically, due to limited sample sizes, we find that only making one sampling for each censored subject will not achieve a good calibration performance. Instead, we propose the method of repetition sampling, i.e., sampling R times from U 0, ŜM (ci|xi) .This ablation study tries to find how this repetition parameter R affects the performance, in terms of both discrimination and calibration. We gradually increase R, from 3 to 1000, and assume the performance should converge at a certain level. This ablation study uses the following hyperparameters to enable a fair comparison:• Interpolation: PCHIP TL;DR The repetition parameter value has barely any impact on the C-index, and increasing R can benefit the marginal and conditional calibration, with convergence observed around R = 100.In Figure 18 , we see that the C-index for CSD-iPOT using 5 different repetition numbers has no visible differences, for almost all baselines and all datasets. The only exception is there are slight differences for DeepHit baseline on HFCR and Employee datasets, where higher R will slightly decrease the C-index performance insignificantly.In Figure 19 , we can clearly see the repetition parameter has a great impact on the marginal calibration. For most datasets, a higher repetition will significantly improve (decrease) the marginal calibration score. This trend is more clear for high censoring rate datasets (FLCHAIN, Employee, MIMIC-IV, etc), while for low censoring rate datasets (GBM, NACD, SUPPORT, which have censoring rates less than 40%), the trends still exist but the differences are not very significant. The marginal calibration performance usually converges when R ≥ 100 (the pink and green violins are almost the same for all datasets and all baselines). As a higher repetition parameter will increase the memory usage and computational complexity (Appendix E.5), we suggest using R = 100.In Figure 20 , the trends in conditional calibration are similar to those in marginal calibration. Most of the differences are more significant for high censoring rate datasets and less significant for low censoring rate datasets. And the conditional calibration performance also converges at R ≥ 100 (the pink and green violins are almost the same for all datasets and all baselines).Ablation Study #2: impact of predefined percentiles P Different choices of P may lead to slightly different survival distributions, all of which allow us to obtain provable distribution calibration, as discussed next. Theoretically, discretizing a continuous curve into a series of discrete points may result in some loss of information. However, this can be mitigated by using a sufficiently fine grid for percentile discretization. Therefore, we anticipate that if we select P that contains finer grid percentiles, the performance (both calibration and discrimination) will be better.Previous studies have commonly employed 10 equal-width probability intervals for calculating distribution calibration [3, 6, 8] , making an intuitive starting choice for 9 percentile levels at P = {10%, 20%, . . . , 90%}3 . We compare this standard setting with 19 percentile levels (P = {5%, 10%, . . . , 95%}), 39 percentile levels (P = {2.5%, 5%, . . . , 97.5%}), and 49 percentile levels (P = {2%, 4%, . . . , 98%}). All the calibration evaluations (for both marginal and conditional) are performed on 10 equal-width intervals to maintain comparability, following recommendations by [3] . This ablation study uses the following hyperparameters to enable a fair comparison: TL;DR The number of percentiles has no impact on the C-index, and has a slight impact on marginal and conditional calibration.In Figure 21 , we see that the C-index for the 4 percentile numbers has no visible differences, for all baselines except CQRNN. It is worth mentioning that the difference for CQRNN is due to its quantile regression nature, requiring more trainable parameters as percentiles increase. In Figure 22 , we can see the number of percentiles has some impact on the marginal calibration. For example, the results for Employee and MIMIC-IV show some perturbation as we increase the number of percentiles. However, trends are inconsistent across models, e.g., fewer percentile levels are preferable by AFT, GB, N-MTLR, DeepSurv for Employee, while more percentile levels are preferable by DeepHit and CoxTime. Also, most differences in Figure 22 are insignificant.In Figure 23 , the trends in conditional calibration are similar to those in marginal calibration. Most of the differences are insignificant and the trends are inconsistent.These results suggest the choice of predefined percentiles P has minimal impact on the performance. In practical applications, the reader can choose the best-performing P at their preference. However, it is worth noticing that a higher percentile will result in more computational cost. Therefore, we suggest that it is generally enough to choose P = {10%, 20%, . . . , 90%}. Justification: The CSD-iPOT method does not have theoretically guarantee for the preservation of Harrell's C-index. This is mentioned in Section 3.4 and Table 1 , and theoretically justified in Appendix C.3. Section 5 measures this perspective empirically.• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.• While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark). • The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: //nips.cc/public/guides/CodeSubmissionPolicy) for more details. • The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc. • The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are omitted from the script and why. • At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable). • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.38th Conference on Neural Information Processing Systems (NeurIPS 2024).This is a technical assumption in conformal prediction. This assumption is easy to solve in practice because users can always add a vanishing amount of random noise to the scores to avoid ties.As recommend by Qi et al.[35], the performance reaches the optimal when |P| = 9 or 19. |P| = 19 means that the predefined percentile levels are P = {5%, 10%, . . . , 95%}.Note that 0% and 100% are excluded because S-1 M (1 | xi) = 0 and S-1 M (0 | xi) = tmax < ∞ are fixed endpoints. 


Decision: The paper studies conditional calibration in the context of survival analysis, proposing a method based on conformal prediction that uses the model's predicted individual survival probabilities at observed times. This method effectively improves that model's marginal and conditional calibration, without compromising discrimination. The authors provide theoretical guarantees for both marginal and conditional calibration and test it's performance across 15 real-world datasets.

The majority of the maintained their positivity towards the paper during the discussion phase. One reviewer remained negative, on the principle that the paper contains a significant flaw. In short, the theoretical guarantee provided by the authors holds under an assumption for asymptotic consistency of the distribution estimator; this negates the need for calibration. Moreover, conditional calibration is a hard problem, as cited by many works that the reviewer referenced, and the authors did not include in their original submission.  Hence, the paper should be rejected on principle: the claim by the authors that they solve a hard problem is sufficient grounds for rejection, tantamount to claiming that the paper solves an NP hard problem, and subsequently proceeds to provide experimental results on its performance.

The authors responded in their rebuttal that their theoretical claim does not contradict prior work, that shows hardness under finite samples (as their result is asymptotic). They also acknowledged that, nevertheless, this assumption is indeed strong, even if quite common in conformal prediction literature. They also indicated that their contribution (and therefore their emphasis) was on the experimental results, that indeed establish excellent performance on 15 datasets. 

Remaining reviewers, and the AC, agreed with the authors on this front. They all agreed that the authors have adequately addressed these concerns in their rebuttal and that, even if theoretical results are weak, they are (a) correct, in that there are no errors in the statements of theorems and proofs, (b) rely on standard (albeit strong) assumptions in the field of conformal prediction, and (c) they are supported by very strong and exhaustive experiments, that make the paper worthy of reporting to the community. There was also an agreement with the authors that no results can prove conditional calibration under more relaxed assumptions.

It was also pointed out that a similar argument about weakness of theoretical guarantees could be made for marginal calibration; for example, Kaplan-Meier is a consistent estimator for the survival function but cannot be used for time-to-event prediction. In practice, one is interested in accurate time-to-event prediction models, where calibration is of secondary concern. For any given problem, numerous solutions can yield well-calibrated predictions that are not necessarily accurate and, thus, not useful. Hence, focusing on a strong theoretical result that guarantees conditional calibration might not be of practical use if the models derived from such an estimator do not yield accurate predictions. As the paper demonstrates, via extensive experimental results across 15 datasets and 7 baselines, that CSD-iPOT significantly improves the calibration (both marginal and conditional) of baseline time-to-event models with minimal loss in the concordance index (C-index). Therefore, even if the theoretical results may not be robust, the empirical results are compelling. 

The authors are strongly encouraged to address the concerns raised by all reviewers, particularly the reviewer that remained negative. Papers showing hardness should be cited and contrasted to. The claim regarding theoretical guarantees should be toned down and stated in the appropriate context. Additional comments made by remaining reviewers, that were sufficiently addressed by the authors, should also make it to the paper, even if in the supplement.

Reviews: review: paper_summary: paper_summary: While previous work focuses on calibration in a marginal sense, this paper proposes a post-processing approach that also imposes conditional calibration on all individual features. Therefore, the proposed method can guarantee equal calibration for different groups, ensuring fairness.

Contributions:

1. The paper proposes the post-processing approach to ensure conditional distribution calibration in survival analysis (accommodating censorship) and develops the corresponding metric.
2. The paper provides asymptotic guarantees for both marginal and conditional calibration and demonstrates the computation complexity.
3. The paper validates the method by conducting experiments on 15 datasets.
strength_weakness: strength_weakness: 1. Missing related work in conformal prediction: the paper claims the method is based on conformal prediction but it is not mentioned and illustrated clearly in the main content. Thus, this knowledge gap hampers the understanding.

2. The method seems to be sensitive to the percentiles. As shown in Figures 2 (d) and (e), the post-processed curves are piecewise functions highly related to the percentiles (1/3, 2/3). Thus, the resulting curves would be sensitive to the choice of percentiles. Therefore, the number of percentiles should be included in the complexity and asymptotic behavior analysis.  

3. The discrimination properties of the method are unclear. The authors first claim it lacks the discrimination performance preservation property for Harrel's C-index, which conflicts with the results in Table 1.

4. The empirical exploration and validation of the group fairness is missing.
questions: questions: As mentioned in the weakness:

Q1: How does the choice of percentiles affect the performance? Are there any empirical results that can illustrate the effect?

Q2: Does the approach keep the Harrel C-index?

Q3: There are plenty of survival analysis methods exploring the heterogeneity of groups. How does the proposed method compare to the SOTA heterogeneity methods? Will it improve the performance?
limitations: limitations: It remains unknown if this method would improve/hurt fairness in survival analysis. The study is missing.
paper_summary: paper_summary: The paper enhances the Conformalized Survival Distribution (CSD) post-processing framework to account conditional calibration. The proposed framework, CSD-iPOT, utilizes a conformal set to adjust survival curves vertically, aligning them with predetermined percentiles at test time. Unlike CSD, which relies on Kaplan-Meier curves, CSD-iPOT leverages Individualized Survival Distributions (ISD) for censored events when constructing the conformal set. The paper provides theoretical guarantees for marginal calibration, conditional calibration, and the monotonicity of ISD. Comprehensive experimental results from 15 datasets demonstrate that CSD-iPOT enhances both the marginal and conditional calibration of baseline models with minimal impact on the concordance index (C-index).
strength_weakness: strength_weakness: *The paper highlights that CSD-iPOT lacks theoretical guarantees for preserving Harrell's concordance index. This is a major limitation of this work; nevertheless, experimental results demonstrate that the impact is minimal*

*The description in lines 142-162 requires some improvements:*

- Eqn. 4: Needs to be adjusted to something like $\tilde{S}^{-1}(p | x_{n+1}) = T(\hat{S}^{-1}(\text{Percentile}(\cdot) | x_{n+1}))) $, where $T$ is the proposed vertical transformation. Then provide all the necessary numbered steps to finally obtain the calibrated ISD (Eqn. 5).

*For completeness, I encourage the author(s) to also benchmark competitive baseline models that directly model event times, e.g., Chapfuwa et al. 2018 and Miscouridou et al. 2018*

*Minor:*
- In Table 1, CSD-iPOT should be marked with an 'X' under Harrell’s concordance index category.
- Line 297: Typo -> "Compare to CSD"
questions: questions: - Given that CSD-iPOT lacks theoretical guarantees for preserving Harrell's concordance index, is there a hyperparameter that explicitly controls this trade-off?
- It seems that CSD-iPOT achieves calibration lower than the empirical lower limit set by Kaplan-Meier. Are these the instances that result in an impact on the C-index?
- Why does CSD-iPOT struggle with both calibration and the C-index for models such as DeepHit?
- Are the calibration metrics obtained at different percentiles than those used for the post-processed ISD? What happens if these two sets are different?
- Why does CSD-iPOT improve marginal and conditional calibration models better than CSD? Is the difference due to how CSD-iPOT handles censored events?
limitations: limitations: I encourage the authors to discuss the limitations of their work, including any violations of modeling assumptions.
paper_summary: paper_summary: The authors study the problem of how to create individual survival distribution (ISD) models which are well-calibrated, both in a marginal and conditional sense, without negatively affecting the discriminative performance.

They refine the "Conformalized Survival Distribution" (CSD) approach and propose "Conformalized Survival Distribution using Individual survival Probability at Observed Time" (CSD-iPOT). Both are post-processing methods which utilize conformal prediction, and can be applied on top of various survival analysis models to improve their calibration. CSD-iPOT is however designed to not only improve marginal calibration, but also conditional calibration.

The method is evaluated on 15 datasets, comparing baseline (without neither CSD nor CSD-iPOT), CSD and CSD-iPOT versions of 7 survival models. The results are quite promising overall.
strength_weakness: strength_weakness: - I found parts of Section 3 and 4 a bit difficult to follow, the proposed method/metric could perhaps be described in a more intuitive way. Figure 2 is neat, but would be helpful to have a similar visualization to illustrate the difference between CSD and CSD-iPOT.

- The technical novelty/innovation compared to CSD is perhaps somewhat limited, both methods employ the same general approach (post-processing methods which utilize conformal prediction).

- The experimental results could be more convincing. CSD-iPOT usually improves the calibration compared to CSD, but far from always, and the gains also seem to be relatively small quite often. Overall, I find it quite difficult to judge how much added benefit CSD-iPOT actually would have compared to CSD in practice. Are there concrete examples where CSD leads to significant mis-calibration within a certain patient subgroup, which is addressed by CSD-iPOT? The current results/metrics are a bit abstract / difficult to interpret.



Summary:
- Well-written and very solid paper overall that studies an important/interesting problem. The technical novelty compared to CSD, and how much better the proposed method actually would be in practice, is however a bit unclear. I am leaning towards accept.




***
***
***
***

**Update after the rebuttal:**

I have read the other reviews and all rebuttals.

2/3 other reviews are also positive, and I think the authors respond well to the third.

All my questions were addressed, Figure 1 and 2 in the provided pdf are neat.

I will increase my score to "7: Accept", this is a very well-written and solid paper that I think should be accepted.
questions: questions: - 189: "Unlike CSD [8], which adjusts the ISD curves horizontally (changing the times, for a fixed percentile), our refined version scales the ISD curves vertically", could you perhaps visualize how both CSD and CSD-iPOT modify the ISD curves in an example? I think that could help illustrate how these two methods differ.

- I thought that CSD-iPOT shouldn't affect the relative ordering of subjects, but when comparing the ISD curves before and after in Figure 2(d) and (e), the top orange curve is always above the top blue curve in (e), whereas it is partially below the blue curve in (d)?

- 236: "Qi et al. [8] demonstrated that CSD theoretically guarantees the preservation of the original model’s discrimination performance in terms of Harrell’s concordance index (C-index) [1]. However, CSD-iPOT lacks this property", but both methods have a check mark for "Discrimination guarantee Harrell’s" in Table 1?




Minor things:
- 77: "In summary, individual calibration is ideal but not impractical", impractical --> practical?

- 219: "with adequate modifications to the accommodate our method", remove "the"?

- 297: "Compared CSD" --> "Compared to CSD"?
limitations: limitations: Yes.
paper_summary: paper_summary: This paper proposes a new postprocessing method, CSD-iPOT, for survival analysis based on conformal prediction.
strength_weakness: strength_weakness: This paper puts emphasis on achieving conditional calibration (not marginal calibration), but the conditional calibration is known to be hard to achieve even for datasets without any censored data point.  See, e.g.,
+ Lei, J. and Wasserman, L. (2014). Distribution-free prediction bands for non-parametric regression. Journal
of the Royal Statistical Society: Series B (Statistical Methodology), 76(1):71–96.
+ Rina Foygel Barber, Emmanuel J Candes, Aaditya Ramdas, Ryan J Tibshirani, The limits of distribution-free conditional predictive inference, Information and Inference: A Journal of the IMA, Volume 10, Issue 2, June 2021, Pages 455–482.
+ Zhao et al., Individual Calibration with Randomized Forecasting, ICML 2020.

This paper does not cite these papers, and this paper does not discuss the hardness of conditional calibration.  Lack of extensive discussion on the hardness of conditional calibration is a serious problem of this paper.

Despite the hardness, this paper claims that the conditional calibration can be achieved in Theorem 3.2.  However, this theorem actually shows nothing: it shows only that, if we have an estimator that achieves the conditional calibration (i.e., Eq. (10) is satisfied), the output of the proposed method achieves the conditional calibration, too.  I think that if we have an estimator that achieves the conditional calibration, the proposed method is not required.

Furthermore, there are many other problems in this paper:
+ While $\Gamma_{M}$ should be a set of scalar values according to line 9 of Algorithm 1, $\Gamma_{M}$ is a set of pairs of scalar values according to line 11 of Algorithm 1.
+ According to lines 180-183, $R$ copies of comformity scores $\gamma_{i,M}$ are generated for each uncensored data point.  If so, the equation between lines 537 and 538 is incorrect.
+ According to lines 180-183, $R$ copies of comformity scores $\gamma_{i,M}$ are generated for each uncensored data point.  Neverthless, this paper assumes that the set $\Gamma_{M}$ does not have any tie in line 538.
+ The proof of Theorem C.1 are not fully described.   Even though the goal of the proof is to prove equation (9) on $\rho_{1}$ and $\rho_{2}$, the proof for censored data points does not argue $\rho_{1}$ and $\rho_{2}$.  An equation for censored data points analogous to the equation for uncensored data points (between lines 537-538) must be presented.
+ The proof of Theorem C.1 completely ignores $R$, even though Theorem C.1 does not hold if $R=1$.  (I think that the main idea of the proposed method, CSD-iPOT, is to use a sufficiently large $R$ to "blur" the censored data points.)
+ The assumptions used in the proof of Theorem C.2 are not clearly stated before the proof: many implicit assumptions are used in lines 563 and 564.   In particular, during the proof of (i), a statement (i.e., an implicit assumption) similar to (i) is used in line 563.
+ In the proof of Theorem C.2, the alleged proof of (i) on $x_{n+1}$ does not include any discussion on $x_{n+1}$.


------------------------
# Additional Comments

I submitted the following comments in the Author-Reviewer discussion period, but it accidentally did not visible to the authors (probably due to the complex system on comments visibility of OpenReview).  I noticed this fact during the Reviewer-AC discussion period, and the AC allowed me to post the comments here.

==comments begin==

Thank you for your comments. I will keep my score.

This paper has several critial problems.

+ The authors' comments did not give any evidence that the proposed algorithm achieves conditional calibration (claimed in lines 47-48). Since the hardness of the conditional calibration is an important topic in machine learning as already studied by many researchers, the authors must pay careful attention when they discuss on this topic.
+ Huge discrepancy between the implemented algorithm (with R=1000) in the experimental section and the provided proof (valid only for R=1), even though the key idea of the proposed algorithm CSD-iPOT is to "blur" a censored subject with a large 
.
Regarding the presentation, all the assumptions must be clearly presented before the proof. The assumption is stated at the end of the proof (in lines 546-548).

Minor thing:

The authors violated the rule on the 1-page pdf in the rebuttal phase: "Please use this PDF only for figures (including tables) and captions that describe the figure."

==comments end==
questions: questions: Nothing.
limitations: limitations: This paper does not discuss the hardness of conditional calibration.

Comments: title: Re: Official Comment by Reviewer kAoa
Official Comment by Reviewer kAoa
comment: We are grateful to the reviewer for taking the time to respond to our rebuttal.

We are sorry we may have misinterpreted the request, as we thought was asking for a comparison of models which only generate time predictions. If the requirement is to compare models that can predict both ISDs and survival times, we confirm that all 7 baselines in the paper (AFT, GB, N-MTLR, DeepSurv, DeepHit, CoxTime, and CQRNN) qualify. Notably, the censored quantile regression neural network model (CQRNN, Pearce et al. 2022), a variant of quantile regression, directly predicts the median survival time for each individual, aligning with what the reviewer has requested.
Thanks for addressing most of my concerns. However, this statement "If a baseline model can generate only a scalar value, it is not immediately clear how to convert this scalar value into a distribution," is not necessarily accurate. Individual survival distributions can be obtained from parametric time-to-event models (e.g., AFT) since $f(t|x) = h(t|x) S(t|x)$. For non-parametric approaches (e.g., DATE), time-to-event distributions are implicitly defined through sampling. 

After reviewing the rebuttal and reviewer's response, I find that overall, this is a solid paper, and I am still leaning more towards acceptance. I am keeping my score as it is.
Thank you for the responses.

I have read the other reviews and all the rebuttals. The responses have addressed my concerns.

I decided to keep my score as it is.
Thank you for the response.

I have read the other reviews and all rebuttals.

2/3 other reviews are also positive, and I think the authors respond well to the third.

All my questions were addressed, Figure 1 and 2 in the provided pdf are neat.

I will increase my score to "7: Accept", this is a very well-written and solid paper that I think should be accepted.