A Machine Learning Model Based on Counterfactual Theory for Treatment Decision of Hepatocellular Carcinoma Patients

Introduction

Different treatment guidelines were proposed to provide appropriate treatment options for hepatocellular carcinoma (HCC) patients, such as China Liver Cancer (CNLC) staging¹ and Barcelona Clinic Liver Cancer (BCLC) staging.² However, they might be inadequate for clinical decision-making. First, there is disagreement among treatment guidelines over the recommended treatment regimens. In CNLC, hepatectomy is recommended for patients with Ia, Ib, and IIa stages, and transarterial chemoembolization (TACE) is recommended for patients with IIb, IIIa, and some IIIb stages.¹ In BCLC, patients with Ia and Ib stages are defined as 0 or A stage, IIa and IIb are defined as B stage, and IIIa and IIIb are defined as C stage.² Hepatectomy, TACE, and systemic treatments are recommended for different stages respectively.² Second, there exists heterogeneity among patients at the same stage, and survival times for patients even with the same treatment regimen are different.^3,4 Finally, in clinical practice, 50% of patients have deviated from their original recommendations based on the treatment guidelines.^4,5 Therefore, the treatment regimens for HCC patients are still unclear.

Several studies^6–10 have found that TACE and hepatectomy have similar survival rates for patients with early-stage HCC. The question of which treatment is superior remains unresolved. Previous studies comparing the efficacy of hepatectomy versus TACE in patients with single HCC have shown that patients in the hepatectomy group had higher survival rates at 1, 3, and 5 years than those in the TACE-treated group.¹¹. However, the study of Lee et al¹² presented contradictory results, with similar 1-, 3-, and 5-year survival rates for hepatectomy and TACE. In addition, TACE treatment was considered the standard treatment option for patients with stage BCLC-B in 2012, with the highest level of recommendation being 1A.¹³ When the guidelines were proposed, there was insufficient data to conclusively support the efficacy differences between TACE treatment and other treatments such as hepatectomy. In addition, hepatectomy has been shown to extend the survival of some patients with intermediate-stage HCC.^6,14–18 Hepatectomy is the treatment of choice for patients with type I/II portal vein tumor thrombus when the lesion is resectable.^19,20 In patients with advanced-stage HCC, TACE is relatively safe and has significantly better survival rates.^21–24 Thus, it is controversial which patients could benefit from hepatectomy or TACE.

To measure the causal effect of hepatectomy over TACE, the outcomes between the two treatments in the same individual would have to be compared.^25,26 So, each patient with liver cancer had two outcomes: the outcome after receiving the actual treatment (observed) and the counterfactual outcome if the patient had received a different treatment (unobserved).²⁶ However, counterfactual outcomes are difficult to observe and counterfactual cases may be missed even in randomized controlled trials.²⁶ This study attempted to estimate the counterfactual outcome of patients treated with transhepatic resection and TACE by machine learning techniques. Therefore, a machine learning model built upon the counterfactual theory can effectively compare the efficacy of hepatectomy with TACE based on retrospective data.

In this study, we first estimate the counterfactual outcomes of HCC patients treated with hepatectomy and TACE using retrospective data based on machine learning model, and construct a prognostic model to predict the overall survival of patients.

Material and Methods

Patients

This study was approved by the Ethics Review Committee of the Affiliated Hospital of North Sichuan Medical College (Number: 2022ER059-2), and the informed consent of the participants was waived for this retrospective study. All the patient identity information has been anonymized. The research was conducted in accordance with the principles of the Declaration of Helsinki.

The workflow of the study is shown in Figure 1. Patients with HCC whose first treatment option was hepatectomy or TACE at the Affiliated Hospital of North Sichuan Medical College from June 2016 to July 2021 were retrospectively collected in the cohort study. The inclusion and exclusion criteria are described in supplementary text 1.

Figure 1 The workflow of the study. (A) Collecting clinical features and extracting the radiomics features from CT images. (B) Predicting the patients undergoing the hepatectomy or TACE death counterfactual probability. (C) Building prognostic model to predicting the prognosis of HCC patients. (D) Comparing the performance of the models. (E) Analyzing the predictive ability of the model in different subgroups. Abbreviations: ROI, region of Interest; TACE, transarterial chemoembolization; HCC, hepatocellular carcinoma.

We created a standardized form for collecting the clinical variables (supplementary text 2). The initial treatment was decided based on CNLC guidelines and the patient’s willingness.

To ensure a causal effect of efficacy and treatment regimen, patients with similar baseline characteristics were screened by propensity score matching (PSM) and used as an independent testing cohort (supplementary text 3). The supplementary table 1 shows the correlation coefficients of baseline characteristics with the treatment regimen, survival time, and survival status. For the remaining unmatched patients, the inverse probability of treatment weighting (IPTW) was used to balance the confounders between treatment groups, and the weighted patients were used as the training cohort (supplementary text 3).

The primary endpoint was the overall survival (OS). The definition of OS and the follow-up program is described in supplementary Text 4. As of the last visit, 36 patients were lost to follow-up. Patients lost to follow-up were censored for survival.

Pre-Processing, Segmentation, and Extraction of Radiomics Features

All the enhanced CT examinations were completed within one month prior to treatment. The scanning parameters and acquisition protocols are shown in supplementary table 2 and supplementary text 5. Supplementary text 6 shows the image pre-processing information.

Two radiologists manually delineated the regions of interest (ROIs) along each tumor margin and peritumor regions using 3D Slicer (http://www.slicer.org). Details of annotation and its reproducibility assessment are described in supplementary text 6. Six ROIs were obtained for each subject, including tumoral ROI at the arterial phase (A), tumoral ROI at the portal phase (P), peritumoral regions within 3mm/5mm distances around the tumor at the arterial phase (A3/A5), and peritumoral regions within 3mm/5mm distances around the tumor at the portal phases (P3/P5).

Radiomics features were extracted from the pre-processed images based on these annotated ROIs by using the Pyradiomics package (https://pyradiomics.readthedocs.io/en/v3.0.1/), see supplementary text 7 for details. Z-score normalization was conducted for the radiomics features to reduce the deviance of dimension.

Treatment Efficacy Prediction

For patients in two treatment regimens (hepatectomy and TACE group), two survival prediction models (Death_hepatectomy and Death_TACE) were constructed respectively, so that patients’ survival status can be predicted from their clinical characteristics and imaging radiomics characteristics. Supplementary text 8 shows the details of the model construction. The Death_hepatectomy and Death_TACE models were then used to predict the survival (probability of death) and the counterfactual survival for each patient assuming their treatment regimens were switched. For patient l in the hepatectomy group, the probability of death predicted by the Death_hepatectomy model () was regarded as the proper or fitting prediction, whereas the probability of death predicted by applying the Death_TACE model () was considered as the counterfactual prediction. and were calculated as follows:

(1)

where T_i is the treatment regimen of group i (i=1 means hepatectomy group, i=0 represents TACE group), Y_i is the survival status of group i, _i represents predictive variable (features) for group i. Z_i is the difference of the probabilities of death between the two predictions. If Z_i is greater than a prescribed constant value C (C=0.1),²⁷ the patient can be considered to be more favorable for survival with the TACE, otherwise the favorable treatment remains the same.

Subsequently, we compared the prognosis of the TACE group using the same machine learning method, and the favorable treatment for each patient can be re-determined. Finally, all the patients were regrouped into hepatectomy- preferred and TACE-preferred groups. Then, we used stratified analysis to verify effectiveness of the new grouping, and the clinical characteristics of each group are more distinguishable and interpretable and are also more consistent with the guidelines.

Feature Selection of Prognostic Model

Clinical characteristics were selected by univariate and multivariate analysis of cox proportional hazards regression (COX). Then, prognostic models were constructed based on these selected clinical characteristics, which exhibited p-values below 0.05 in multivariate Cox analysis.

Correlation analysis and random survival forest (RSF) were used to select radiomics features from CT images. The former was used to eliminate variables with collinearity, and the latter was applied to rank the importance of the remaining non-collinearity features. Top 30 features were selected for each RSF operation, and the RSF feature selecting was performed 100 times, so that the top features appeared 80 times and above were finally maintained. The weighting of each feature was assigned based on its ranking of importance, and Radscores were calculated by linear multiplication of the features and weights (Radscore_A, Radscore_P, Radscore_A3, Radscore_P3, Radscore_A5, and Radscore_P5).

Prognostic Model Construction

We used the accelerate failure time model (AFT) to estimate the effect of covariates on patient survival time (supplementary text 9). The standard semiparametric AFT model correlates covariates to the logarithm of survival time through the following regression model:

(4)

where t₁ represents the survival time of patient l, β is the parameter vector, φ₁ is the covariate vector, and ε is a constant (8.455 in this study). Minimization of mean squared error can be used to estimate the regression parameter given a number of patient data.

Statistical Analysis

Open-source R (version 4.1.1) software was used for data analysis. Confounding factors between the two treatment groups were determined by correlation analysis. In the description of the characteristic distribution, the Kolmogorov–Smirnov test was used to test the normal distribution of the continuous variables. For normally distributed variables, t-tests were used and expressed as mean ± SD (standard deviation). Otherwise, the Mann–Whitney U-test was used and expressed as the median (inter-quartile range). For categorical variables, the chi-square test or Fisher test was used and expressed as counts (%). p < 0.05 was considered statistically significant. In addition, stratified analysis using clinical factors was performed to verify the predictive performance of machine learning models and the prognostic performance of the AFT model.

So that other researchers can independently replicate our findings, all experimental steps have been described in detail and completely. All the experimental analyses were performed with open-source software, including R software (http://www.r-project.org), 3D slicer (http://www.slicer.org), and Pyradiomics package (https://pyradiomics.readthedocs.io/en/v3.0.1/). And the data and code analyzed during the study are publicly available through the Treatment-efficacy-prediction GitHub repository (https://github.com/lalala111-hai/Treatment-efficacy-prediction).

Results

Patients

A total of 418 patients were enrolled in the study, of whom 267 (mean age, 58 years ± 11 [standard deviation]; 228 men) were treated with hepatectomy, and 151 (mean age, 59 years ± 13 [standard deviation]; 127 men) were treated with TACE. The flow chart of patient recruitment is shown in the supplementary figure 1. The median of the follow-up intervals was 864 days for the entire cohort (inter-quartile range [IQR]: 657–1004), 919 days for the training cohort (IQR: 643–1037), and 856 days for the testing cohort (IQR: 616–1253). The clinical characteristics of the patients are shown in the supplementary table 3.

One hundred and twenty patients were matched by using PSM to form an independent testing cohort (TACE group, 60 patients; hepatectomy group, 60 patients). The remaining 298 patients (TACE group, 91 patients; hepatectomy group, 207 patients) were used as the training cohort, and confounding variables were treated by IPTW. The PSM and IPTW results demonstrated there is no significant difference in clinical characteristics between the two treatment groups (P >0.05) (supplementary table 4, supplementary table 5).

Treatment Efficacy Prediction

Supplementary text 8 and supplementary table 6 shows the details of the model construction.Table 1 shows the performance of different models in the hepatectomy and TACE groups. The result shows that the RFS model (Death_hepatectomy) built by combining clinical characteristics, Radscore_A3 and Radscore_P has the best performance (training cohort: AUC=0.839 [0.781–0.896]; testing cohort: AUC=0.838 [0.727–0.948]). And the logistic regression model (Death_TACE) developed by combining clinical characteristics and Radscore_A5 was effective in predicting the probability of death (training cohort: AUC=0.790 [0.661–0.919]; testing cohort: AUC=0.803 [0.690–0.916]).

Table 1 Performances of the Different Models in the Hepatectomy and TACE Groups

When applying the Death_TACE model to predict the probability of death for the hepatectomy (counterfactual treatment) and the TACE group (fitting group), Figure 2 shows the results. Nineteen patients in the hepatectomy group were counterfactually predicted to be treated with TACE; 89 patients in the TACE group were suggested to be treated with hepatectomy (supplementary table 7). Patients were divided into the hepatectomy-preferred group and TACE-preferred group based on the results predicted by the machine learning model. In the hepatectomy-preferred group, hepatectomy had a prolonged survival time than TACE in both the training cohort (P < 0.001) and the testing cohort (P < 0.001) (Figure 3A–C). In the TACE-preferred group, patients treated with TACE had longer overall survival than hepatectomy in both the training cohort (P = 0.011) and the testing cohort (P < 0.001), indicating that TACE might be a better treatment option than hepatectomy (Figure 3B–D). In addition, in the hepatectomy group, hepatectomy-preferred patients had a prolonged survival time than TACE-preferred patients in both the training cohort (P < 0.001) and the testing cohort (P < 0.001) supplementary Figure 2 A–C). There were similar results in TACE treatment (supplementary figure 2 B–D).

Figure 2 Real and counterfactual treatment regimens. The bar plot shows the real treatment regimens, and the dotted lines indicates the counterfactual treatment regimens. Abbreviations: LR, Liver resection; Tace, transarterial chemoembolization.

Figure 3 Kaplan-Meier curves of prognostic differences between hepatectomy and TACE treatment. Comparisons of the hepatectomy and TACE groups in the Hepatectomy-preferred group (A), and the TACE-preferred group (B) in the training cohort. Comparisons of the hepatectomy and TACE groups in the Hepatectomy-preferred group (C), and TACE-preferred group (D) in the testing cohort. TACE, transarterial chemoembolization. Abbreviation: TACE, transarterial chemoembolization.

Feature Selection of Prognostic Model

Univariate and multivariate COX regression results showed (supplementary table 8) that treatment regimen (p=0.033, HR=1.440, 95% CI=1.030–2.010), ALBI grade (p=0.045, HR=1.420, 95% CI= 0.998–2.030), Tumor number (p= 0.018, HR=1.230, 95% CI= 1.040–1.460) and GGT (p=0.029, HR=1.230, 95% CI=1.020–1.480) were significantly associated with survival status in HCC patients. Covariance analysis and RFS were used for radiomics features selection, and the results of selected features for the six ROIs are shown in supplementary figure 3.

Prognostic Model Construction

Using the minimum Akaike information criteria (AIC) as the model selection criterion, the AFT model built by combining clinical characteristics, Radscore_A and Radscore_P5 had the best performance (AIC=73.391). Using the results of the machine learning model as a treatment regimen, the AIC value of the model was reduced to 69.753. The specific variables of the prognostic model are shown in Table 2 (all variables p<=0.05). In the validation cohort, the time-ROC curves (Figure 4) show that compared to the actual treatment regimen (Figure 4A and B), the counterfactual treatment regimen (Figure 4C and D) improved the 1–5 years AUC by 7.6% (0.728 vs 0.652), 7.7% (0.814 vs 0.737), 8.2% (0.789 vs 0.707), 4.3% (0.811 vs 0.768) and 4.3% (0.811 vs 0.768). The integrated Brier scores of the prognostic model are 0.167 in the training cohort and 0.168 in the testing cohort.

Table 2 Variables of the Prognostic Model

Figure 4 Time-ROC curves of the prognostic model. The time-ROC curves for the prognostic model of the actual treatment regimen (A), and the counterfactual treatment regimen (B) in the training cohort. The time-ROC curves for the prognostic model of the actual treatment regimen (C), and the counterfactual treatment regimen (D) in the testing cohort. Abbreviation: AUC, areas under the receiver operating characteristic curves.

The median risk score, 6.748, predicted by the prognostic model was used as a threshold to classify HCC patients into high-risk and low-risk groups. The cumulative mortality at 1- to 3-year for high-risk and low-risk patients were 57.02%, 67.03%, 67.03% and 16.78%, 30.87%, 34.9%, respectively. The Kaplan-Meier survival curve of the OS is shown in Figure 5.

Figure 5 Kaplan-Meier curves of high-risk and low-risk patients. Kaplan-Meier curves of overall survival in the training cohort (A) and the test cohort (B).

Stratified Analysis

Stratified analysis of patients by Cirrhosis (normal or abnormal), CNLC (I or II/ III), and Child-Pugh classification (grade A or B) showed that surgical treatment had significantly longer survival than TACE treatment in the hepatectomy-preferred group with p < 0.001 (supplementary figures 4–6). In contrast, for the TACE-preferred group, TACE treatment was associated with significantly longer survival (p < 0.05, supplementary figures 4–6). The prognostic model enables for risk stratification in different populations, with survival times significantly lower in the high-risk group than in the low-risk group (p≤0.001, supplementary figure 7).

Discussion

In this study, a predictive model for hepatectomy versus TACE treatments in HCC patients was developed based on counterfactual theory. The prognosis of patients with real treatment regimens versus counterfactual treatment regimens was compared. Our machine learning models estimated patients who might be beneficial from the preferred treatment regimen (hepatectomy and TACE), and thus recommend the treatment option for the patients that may be of greater benefit. In addition, a prognostic model for predicting the overall survival time of HCC patients was constructed, which showed high discrimination in both the training and testing cohorts.

Traditional decision support systems such as the CNLC staging and the BCLC staging classify patients into different grades and recommend corresponding treatment strategies. Although clinical guidelines are important for prognosis and treatment, the decision-making capacity may be inadequate.^2,15,28–30 Thus, an increasing number of studies have compared the efficacy of hepatectomy versus TACE treatment in patients with the same stage, but even randomized controlled trials would not immune from the causal inference problem.²⁵ Although previous radiomics or deep learning studies on HCC have built prognostic models, they have not taken into account the impact of treatment regimens.^31,32 In contrast, Cucchetti et al³³ considered the effect of treatment regimens based on the counterfactual principle, but only estimated the mean effect of the treatments. The study only evaluated the overall efficacy of different treatments, so it could not give individualized treatment recommendations according to the patient’s actual situation. The above studies may only provide limited help to clinicians in making treatment decisions. Therefore, this study inferred the prognosis of HCC patients taking different treatments based on the counterfactual principle. In addition, we provided a prognostic model that took into account treatment factors, which not only had high prognostic ability but also allowed for risk stratification. There is the first study designed to assess the efficacy and the prognosis of HCC patients in different treatments.

In this study, the probability of death was calculated for each patient after hepatectomy and TACE treatment. The results showed that 89 patients underwent TACE but the counterfactual predicted that hepatectomy was recommended, while only 19 patients who underwent surgery differed from the counterfactual results. As a whole, there may be a greater survival benefit with hepatectomy. This study confirmed previous findings, a large multicenter study comparing the benefits of hepatectomy at different stages of BCLC showed that hepatectomy was more beneficial to HCC patient survival in patients without liver dysfunction, regardless of the patient’s BCLC stage.²⁸ Similarly, Cucchetti et al³³ took a counterfactual approach to assess the difference in efficacy between hepatectomy and TACE, showing that patients treated with hepatectomy rather than TACE would have increased survival time from 18.3 months to 38.0 months and survival rates by more than 50%. In addition, the difference in efficacy between hepatectomy and radiofrequency ablation was considered in the study by Cucchetti et al,³³ but was not discussed in our study. In subsequent studies, we will focus on the other treatment options, such as ablation, immunotherapy and chemotherapy.

In this study, there were patients with early, intermediate and advanced HCC whose actual surgical procedures differed from the counterfactual recommended treatment, suggesting that there is a large heterogeneity among patients with the same clinical stage and that a reasonable treatment should be selected for each patient. Similar conclusions were reached in a study by Golfieri et al⁴ which showed that patients at the same clinical stage exhibited different tumor burden and liver function but were still recommended for the same treatment. The large variation in survival of patients at the same stage suggests that clinical staging guidelines do not fulfill the role of dividing patients into subgroups with significantly different prognosis.⁴ In our study, each patient was further analyzed, combining clinical factors and quantitative imaging features to provide a more personalized treatment opinion. Specifically, when confronted with patients who have difficulties with decision-making choices, the probability of death for each patient treated with hepatectomy and TACE was calculated and compared using the model from this study. If the probability of death is similar for those treated with hepatectomy and TACE, then TACE is recommended for pain control. Conversely, whichever treatment received had a lower probability of death was recommended. The predictive value of the probability of death provided in this study may provide a reference for physicians when making treatment decisions. In our study, a prognostic model to predict the OS of patients with HCC was also constructed. The prognostic model not only has high predictive ability but also could classify HCC patients into subgroups with significant differences in survival time. For patients in the high-risk group, a more rigorous follow-up strategy or consideration of additional adjuvant therapy should be adopted after surgery. For patients in the low-risk group, a more rational follow-up protocol could be adopted considering the risk-to-benefit ratio.

However, there are still some limitations of our study. First, the study only considered hepatectomy and TACE, and other treatment modalities should be included for comparison in future studies. Second, this retrospective study was observational and relied on statistical extrapolation to address the issue of causal inference, our results require further validation with other hospitals to check for the generalizability. Third, the sample size of the single-center study might be relatively small. Considering the complexity of the causal inference and treatment selection issues, the subsequent studies with the larger database to validate the results is needed. Finally, the initial treatment regimen was a carefully considered decision by the clinicians, and we used IPTW and PSM to control for baseline differences between patients to minimize the influence of the subjective judgment of clinicians and patients on the treatment regimen. But the possibility of selection bias cannot be completely eliminated.

Conclusion

In conclusion, we predicted the counterfactual outcomes in HCC patients treated by hepatectomy and TACE based on the counterfactual principle and predict the overall survival of HCC patients. This study might provide a new direction for individualized patient management and the selection of potentially appropriate treatments for patient.