First, we state a method that does not take into account any error. To accurately determine the TMB-positivity thresholds from multifaceted efficacy analyses, previous work (27, 28) has integrated two types of clinical outcomes: binary tumor response (e.g., responder or non-responder) and continuous time-to-event (TTE) endpoints (e.g., survival or progression-free survival). This integration accounts for the within-subject dependency between these two types of endpoints. For simplicity, we refer to this method as the Bayesian Naïve Method (Bayes-NM).

Specifically, for patient i , R i denotes the tumor response status ( R i = 1 ,   0 for response and non-response, respectively), Z i denotes clinical covariates (e.g., age, gender, stage of cancer), and T M B i denotes the error-free biomarker. The tumor response R i depends on both the covariates Z i and T M B i . We model this using a logistic regression model, which relates the probability of response to these factors:

where θ represents a vector of unknown parameters.; α z   and α m denotes regression coefficients for the covariates Z i and T M B i ; b i is a random effect for patient i (accounting for unmeasured individual variations).

Next, for time-to-event analysis, the event time T i represents the observed time until the occurrence of an event (e.g., tumor recurrence, progression, or death), which is taken as the minimum of the actual event time U i and the censoring time C i . Define the event indicator as Δ i = I ( U i ≤ C i ) , where I ( · ) is the indicator function that equals 1 if the event occurs and 0 if censored. For the time-to-event (TTE) data, we use the Cox proportional hazards (Cox-PH) model, which focuses on classifying patients by their survival risks:

where h i ( t ) describes the instantaneous risk for patient i at time t; h 0 ( t ) is known as the baseline hazard, typically modeled using the Weibull distribution; β z and β m is effect coefficients for the covariates Z i and T M B i ; the shared term b i is a random effect term accounting for correlation between event times and responses. The random effect b i is assumed to follow a normal distribution N ( 0 ,   σ b 2 ) , which represented the intra-subject correlation between event times and individual response.

The observed dataset for Bayes-NM is denoted as D n = { R i , T i , Δ i , Z i , T M B i } i = 1 n . Multiscale endpoints (response and survival) can be jointly modeled by incorporating random effects and adjusting for the dependencies between the response probabilities and event times. Formally, the joint likelihood for the data is given by:

The joint log-likelihood function is then:

Inference about parameters θ is typically based on the maximization of this log-likelihood. Since the likelihood involves random effects, we use Markov Chain Monte Carlo (MCMC) sampling to estimate the posterior distribution of the parameters A.

To begin Bayesian inference, we must specify prior distributions for the unknown parameters. For the regression coefficients and variance components, we use non-informative priors as follows:

These priors reflect weak prior knowledge and allow the data to predominantly drive the posterior estimates.

Once the priors are defined, we use MCMC to sample from the posterior distribution of the parameters. The parameters are estimated as the average of the posterior sample:

where θ k ′ is the k-th sample from the posterior distribution and K is the number of MCMC iterations.

Equations 1, 2, and 4 comprise the core of Bayes-NM, assuming the absence of measurement errors, with further elaboration available in the Supplementary Material .

In practice, the measured TMB (i.e. T M B * ) is an approximation of the true TMB, due to inherent measurement errors, i.e., T M B * = T M B + e , where both T M B and e are unobserved potential variables and mutually independent. While external or internal validation sets can sometimes provide information on error characteristics, they may not always be available. This complexity in real-world data necessitates different strategies for handling errors. Below, we describe several methods for addressing these errors in the context of TMB measurement.

In addition to errors arising from TMB measurements, paired endpoints, such as clinical response, may also suffer from measurement inaccuracies. While errors in continuous TTE endpoints typically do not have a significant effect on the consistency of the analysis, response misclassification in discrete endpoints can substantially affect model accuracy and predictions. To address this, we propose the Bayesian Misclassification Correction Method (Bayes-MCM), which aims to account for misclassification errors in the observed response.

In this case, the observed R i   is no longer an accurate reflection of the true tumor state but may be a misclassified version of the actual response Y i . To address this, we define the following misclassification parameters:

the probability of correctly classifying a true response (tumor shrinkage) as a positive response.

the probability of correctly classifying a true non-response (no tumor shrinkage) as a negative response.

To incorporate these parameters, the standard mixed-effects logistic regression model for the true response endpoint:

The corresponding model for the observed response is:

where P ( Y i = 1 ) is derived from the logistic regression model.

Previous research (37) has demonstrated that without additional constraints or information, the misclassification parameters η and δ   are not practically identifiable. To overcome this limitation, we updated the Bayesian framework by incorporating informative priors and imposing constraints on these parameters. Specifically, the priors are modeled as Beta distributions:

where ϵ 1 , ϵ 2 , ϵ 3 , ϵ 4 are constants chosen based on prior knowledge or assumptions about misclassification rates. For example, η ∼ B e t a ( 9 ,   1 ) assumes that the probability of correctly identifying a response is high, but not certain, with a mean of 0.9 and a 95% credible interval of [0.664, 0.997]. Conversely, δ ∼ B e t a ( 72.5 ,   2.5 ) reflects a strong belief in the reliability of non-response classification, with a mean of 0.98 and a 95% credible interval of [0.917, 0.995]. These priors are constructed to reflect existing knowledge or reasonable assumptions about misclassification rates.

Furthermore, misclassification probabilities may vary across subpopulations, depending on factors such as tumor type, disease stage, or patient demographics. To address this, Bayes-MCM will incorporate hierarchical structures for η and δ, enabling the model to estimate separate parameters for different subgroups. This flexibility allows for tailored misclassification rates that reflect variations based on patient characteristics, such as prior treatments or response patterns.

The core of the Bayes-MCM includes Equations 2, 4, 7, and 8. Full mathematical details and derivations are provided in the Supplementary Material .

While misclassification errors can occur in individual components (like TMB and response endpoints), real-world errors often involve both types of measurements simultaneously. The Bayesian Pairwise Error Correction Method (Bayes-PECM) extends the Bayes-MCM by simultaneously modeling both TMB errors and response misclassification.

In this comprehensive model, Equations 2, 4–8 jointly account for both sources of error. The model is flexible enough to adapt to the uncertainties in both TMB measurements and response misclassification, with details in the Supplementary Material .

The Bayesian framework and pairwise error correction in TMBocelot enable a convergent and comprehensive estimate of ICI benefits for individual patients. By utilizing the joint likelihood to characterize a patient’s factual condition post-treatment and employing TMB as a biomarker, we used TMBcat to identify a stratification threshold based on the cut-off value with the minimum p-value (38, 39). Patients could be categorized into two categories for treatment prognosis comparison based on such a threshold.

We developed a robust framework for determining TMB-positive thresholds, ensuring broad applicability across diverse error scenarios. Depending on the specific context, the corresponding methods outlined above can be applied for parameter estimation, followed by the use of joint probabilities to identify precise TMB thresholds. The complete framework is given in pseudocode in Algorithm 1 and Figure 1A . To facilitate adoption and reproducibility, we provide implementation resources. The source code for TMBocelot, including the full implementation of the Bayesian framework, pairwise error correction, and threshold determination methods, is publicly accessible at https://github.com/YixuanWang1120/TMBocelot. This open-access repository serves as a practical guide for researchers and clinicians. By making the codebase available, we aim to bridge the gap between technical complexity and practical usability, ensuring TMBocelot can be readily integrated into clinical decision-making workflows.

To validate the applicability of TMBocelot, we assembled four cohorts of 438 different patients from publicly available studies, encompassing 57, 68, 73, and 240 patients with NSCLC (26, 40–42). Our primary efficacy endpoints were tumor response and progression-free survival (PFS) evaluated according to the RECIST criteria. Figure 1B illustrates the data collection process. The specific patient information is listed in Supplementary Table S1 .

The decision to focus exclusively on NSCLC cohorts was driven by the extensive availability of publicly accessible datasets for this cancer type and the well-established clinical relevance of TMB as a biomarker in NSCLC. Additionally, NSCLC represents a heterogeneous disease with varying response patterns to immune checkpoint inhibitors, making it a valuable model for evaluating the performance of TMBocelot.

Furthermore, we assessed the performance of TMBocelot in addressing various errors in oncology trials through a series of simulations. The simulations involved 200 individuals, each characterized by an individual-specific random effect b i across multiple endpoints. b i was drawn from a normal distribution with a mean of zero and variance σ b 2 . Tumor response states were categorized as response ( Y i = 1 ) and non-response ( Y i = 0 ). R i = Y i when misclassifications did not occur; otherwise, adjustments were made based on misclassification parameters η  and  δ . Actual response states were generated based on logistic probability, and patient event times were derived from a survival density following the Weibull distribution with a shape parameter of 1.0. Censoring time C was randomly generated from a uniform distribution U (0, 10).

We set α z = − 0.8 ,   α m = 0.4 ,   λ = 1.0 ,   β z = 1.0 ,   β m = − 0.4 and σ b = 0.5 . Z i was generated from the uniform distribution U (0, 1). Various error scenarios were designed, including error-free, TMB errors only, endpoints misclassification only, and both. We used the corresponding approaches described above. In experiments focusing on TMB errors, we assigned different distributions to the true TMB and measurement errors to demonstrate the methods’ efficiency. We separately set η , δ to (0.70, 0.95) and (0.75, 0.98) for experiments considering response misclassification. Detailed settings can be found in the Supplementary Table S2 . Figure 1C roughly depicts the generation of the simulation data.

In the simulations, we report fitted values, average bias, standard deviation (SD), and standard error (SE) for each parameter. SD measures variability in the estimates across 500 simulations, while SE represents their average error.

Table 1 summarizes parameter estimates, with extensive simulations detailed in Supplementary Table S2 . These results highlight the robust performance of TMBocelot across various error scenarios. While TMBocelot performs slightly worse than the true-data estimator under measurement errors or misclassification, it significantly outperforms the naive estimator, which ignores such errors. This strongly supports TMBocelot’s effectiveness in handling diverse error types.

Additionally, distinguishing between the corrected-score and Bayesian robust methods is crucial when considering only TMB errors. The results indicate: i) Regardless of the distribution of actual TMB, the corrected-score outperforms the Bayesian method when TMB errors follow a normal distribution with a standard deviation of 1.0. However, the Bayesian method performs better for errors with a larger deviation. The corrected-score may be more effective for minor errors, while the Bayesian robust method is preferable for more significant errors. ii) If TMB errors follow the extreme value distribution, the corrected-score’s performance deteriorates compared to the normal distribution, while the Bayesian method remains stable. This suggests that the corrected scoring method’s normality assumption may mis-specify the error distribution, especially in asymmetric distributions. Conversely, the Bayesian method demonstrates robustness concerning TMB error and truth value distribution. iii) The SE and SD for the Bayesian robust method exceeded those for the corrected-score method, attributed to a lack of a priori information. These limitations become more pronounced when accurately and robustly estimating the distribution of TMB actual values. The corrected-score is recommended for minor errors and a known variance; the Bayesian robust method is preferable for more significant errors or asymmetric distributions.

To evaluate the effect of errors on TMB thresholds and TMBocelot’s stability, we simulated treatment efficacy across different thresholds, considering 500 patients with actual TMBs positively correlated to favorable outcomes. The coefficients were set to   α z = − 1.8 ,   α m = 0.4 ,   λ = 1.0 ,   β z = 2.2 ,   β m = − 0.4 and σ b = 0.5 . TMB errors were generated by the normal distribution with a standard deviation of 1.0 and misclassification parameters η = 0.75 ,   δ = 0.98 . We derived different thresholds for classifying patients and comparing the treatment efficacy of TMB-based subgroups.

It is worth noting that the selection of T M B i and R i in the thresholding process and the efficacy comparison. When accounting for TMB errors or tumor response misclassification, we used the corresponding posterior estimates. For scenarios ignoring errors, the observed values were used. The comparison outcomes are depicted in Figure 2 .

Significant efficacy disparities between TMB-low and TMB-high groups reflect the correlation between higher TMB and increased antitumor immunogenicity. Specifically, results show: i) Thresholds considering response misclassification result in greater disparities in tumor response ( Figure 2D , p < 0.001) compared to the naive method ( Figure 2B , p = 0.118), with similar survival curve disparities ( Figures 2A, C ). ii) The threshold considering TMB errors results in more significant disparities in both the survival curve and response ( Figures 2E-H ) compared to the threshold based on the naive method ( Figures 2A, B ) though tumor response disparities ( Figures 2F, H , p = 0.011, 0.010) are slightly smaller than those considering only misclassification ( Figure 2D , p < 0.001). iii) Thresholds considering both TMB errors and response misclassification yield the most significant disparities in survival and tumor response ( Figures 2I, J ) compared to all other methods ( Figures 2A–H ). These findings indicate that TMB errors impact survival and tumor response, while response misclassification primarily affects tumor response. The framework effectively establishes robust TMB thresholds, supporting clinical decision-making even with pairwise errors.

To emphasize posterior TMB’s role in patient classification, we used accurate response labels. Logistic predictors based on observed, posterior, and true TMB generated ROC curves in Figure 3 , where Figures 3A – D derive from the same data set, and the differences stem from the randomness of the MCMC algorithm. Results showed posterior and true TMB yielded similar AUC values, both significantly higher than those for observed TMB. Variability between posterior and true TMB AUCs arises from MCMC algorithm randomness and posterior TMB’s tendency toward maximum likelihood, which cannot fully match true TMB. While posterior TMB sometimes outperforms true TMB through parameter combination, it may also perform slightly worse. Overall, posterior TMB surpasses observed TMB in classification performance, underscoring its importance in patient stratification.

Incorporating error considerations, the multi-endpoint joint analysis marks a significant advancement by addressing real-world measurement challenges and surpassing previous studies. To validate TMBocelot, we applied the Bayes-ECM to each experimental cohort, considering TMB errors as recommended by the framework. Subsequently, we determined the TMB threshold by dividing patients into two subgroups using TMBcat and comparing it with the medians. The comparison outcomes are presented in Figure 4 .

The threshold calculated by the proposed method surpasses the medians in all cohorts. This superiority leads to more significant efficacy disparities, notably in nsclc_cohort1_57 ( Figure 4A ) and nsclc_cohort2_68 ( Figure 4B ), nsclc_cohort3_73 ( Figure 4C ) and nsclc_cohort4_240 ( Figure 4D ). Furthermore, to validate the proposed Bayes-PECM, we introduced artificial perturbations in tumor response for nsclc_cohort2_68 and nsclc_cohort4_240. This included a 25% error rate for patients with PR/CR and a 2% error rate for patients with PD/SD. The Bayes-PECM was tested on these two cohorts; the outcomes are depicted in Figure 5 . Figure 5 clearly shows pronounced efficacy disparities for nsclc_68 and nsclc_240, grouped by the Bayes-PECM-based threshold. Notably, the disparities in efficacy for experimental cohorts with added tumor response misclassification are similar to those achieved by experimental cohorts without misclassification, validating the effectiveness of the Bayes-PECM proposed in the framework.

A joint model with pairwise error control facilitates a more comprehensive and robust TMB subgrouping. This approach reveals more significant discrepancies between the efficacies of the TMB-low and TMB-high groups, showcasing the potential of the proposed framework in enhancing the precision and reliability of TMB subgroup classifications.