Individualized treatment effect inference of head and neck cancer with multimodal data Open Access

Predicting the outcome of a specific treatment for an individual patient, known as individualized treatment effect (ITE) inference, is a cornerstone of precision medicine (Lei and Candès, 2021; Pearl, 2001; Prosperi et al., 2020). In complex diseases like head and neck cancer (HNC), where multiple treatment options such as radiotherapy, chemotherapy or their combinations exist, the ability to answer the question, “would this patient have lived longer had an alternative treatment been applied?” is of profound clinical importance. The increasing availability of large-scale medical data, including imaging and electronic health records (EHRs), offers an unprecedented opportunity for data-driven deep learning (DL) algorithms to learn the intricate causal relationships between patient status, treatment choices and subsequent outcomes. However, developing a DL framework for practical ITE prediction from observational data faces several significant challenges.

First, accurately predicting ITE requires a foundation built on patient status, which is characterized by a combination of multimodal data, including 3D medical images (e.g. CT) and tabular clinical variables (Welch et al., 2024). In addition to tabular EHR (Lacroix et al., 2001), a growing body of evidence suggests that imaging data holds great potential for prognosis. Effectively fusing information from heterogeneous yet complementary radiomic, geometric and demographic factors remains a technical hurdle.

Second, although ITE varies across treatments, many existing DL models for medical image-based survival prediction are direct models, as shown in Figure 1(a), that map preoperative patient status $x$ to an outcome $y$ (Kaur et al., 2022). These models either focus on a single treatment (Chaddad et al., 2016; Chang et al., 2016; Nie et al., 2019) or simply ignore the treatment $t$ (Baid et al., 2020; Bakas et al., 2018; Huang et al., 2021). Such models can predict outcomes under a single treatment regimen but cannot be used to compare different treatments for planning.

Third, observational studies, unlike randomized controlled trials, are subject to treatment selection bias (Shalit et al., 2017). For example, patients with certain characteristics (e.g. younger age, better overall health) might be systematically chosen for more aggressive treatments. A naive model, as shown in Figure 1(b), that learns a mapping from patient status $x$ and treatment $t$ to outcome $y$ without explicitly considering the causal role of treatment will confound the effect of the treatment with the baseline prognosis of the patient, failing to generalize to counterfactual scenarios (Liu et al., 2024; Shalit et al., 2017). To build a model that is generalizable to both factual (observed) and counterfactual (unobserved) patient-treatment cases, it is necessary to learn a treatment-invariant representation of the patient’s status.

Recent DL studies for ITE on tabular data (Bica et al., 2020; Curth and van der Schaar, 2021; Shalit et al., 2017; Shi et al., 2019; Xu and Yadlowsky, 2022; Yoon et al., 2018) have established that the core challenge is extracting unbiased, treatment-invariant status features $fθ(x) ⊥ t$ (i.e. disentangling $t$ from $fθ(x)$⁠), a challenge that is underexplored for multimodal status and multiple treatments. With multiple treatments, however, disentanglement via adversarial training – a widely used strategy in binary-treatment ITE prediction (Bica et al., 2020; Curth and van der Schaar, 2021; Shalit et al., 2017; Shi et al., 2019; Xu and Yadlowsky, 2022; Yoon et al., 2018) – cannot effectively disentangle the information between $fθ(x)$ and multiple treatments $t$ (Liu et al., 2021a, 2021b; Roy and Boddeti, 2019). Moreover, adversarial training is notoriously unstable, which can be even more challenging for multimodal data compared to simple tabular data.

To address these challenges, we propose an end-to-end mutual information-regularized deep causal learning framework for ITE estimation in HNC. Instead of training independent models for each treatment, our approach uses all treatment data in a unified manner and explicitly considers their causal relationships. Specifically, we develop a bi-stage adaptive instance normalization (Bi-AdaIN) to inject relevant factors into their corresponding layers. This approach is robust to missing values and efficiently conditions the model on the multiple treatment factors. Our method is general and can be applied to modern DL backbones, such as the 3D Swin transformer (Liu et al., 2021c) and 3D ResNet (Turnbull, 2022). In our bi-stage injection, we fuse clinical variables after each transformer block in the Swin transformer or convolutional layer in ResNet before feature disentanglement, while conditioning on the treatment in fully connected (FC) layers after feature disentanglement. Furthermore, we propose an advanced mutual information (MI)-based disentanglement method over conventional adversarial training to ensure that factual and counterfactual patient-treatment cases share the same feature distribution. Consequently, a single model can be applied unbiasedly to both scenarios.

In summary, our work makes the following novel contributions:

• We present, to our knowledge, the first DL-based study using multimodal medical imaging for ITE estimation.

• We propose an efficient multimodal fusion and treatment conditioning network using Bi-AdaIN, allowing for a synergistic integration of information and causal correlation modeling.

• We develop a robust deep disentanglement learning framework based on MI regularization to extract a treatment-invariant patient representation. This enables a single, unified model to predict outcomes for all treatments, accounting for both factual and unseen counterfactual scenarios.

We demonstrate the efficacy of our approach on the large-scale RADCURE clinical data set (Welch et al., 2024), showing that our model significantly reduces treatment effect bias compared to baseline approaches.

The goal of ITE estimation is to predict the potential outcomes for a patient under different treatments. Let $x$ denote the pretreatment status of a patient (e.g. CT images and clinical variables), $t∈{t1,t2,…,tK}$ be the treatment received from a set of K possible treatments and $y$ be the outcome (e.g. survival length). In an observational data set, we only observe the factual outcome $y$ for the treatment $t$ that was actually administered. The outcomes for all other treatments $t′ ≠ t$ are unobserved counterfactuals.

As illustrated in Figure 1(a), direct models that ignore treatment or simply use it as an input feature fail to address the confounding issue where patient status $x$ influences both the treatment choice $t$ and the outcome $y$⁠. Moreover, in retrospective data sets, simply taking ${x,t}→y$⁠, as illustrated in Figure 1(b), cannot guarantee performance. We only observe partial factuals ${xn,tn,yn}n=1N$⁠, without the counterfactuals associated with the alternative treatments ${xn,tn¬,yn¬}n=1N$⁠, where $tn¬ ≠ tn$⁠. Since the underlying treatment assignment mechanism is $p(t|x)$⁠, the inputs for factuals ${xn,tn}n=1N∼pF(x,t)=p(x)p(t|x)$ and counterfactuals ${xn,tn¬}n=1N∼pCF(x,t)=p(x)p(¬t|x)$ follow different empirical distributions (Johansson et al., 2016), i.e. there is a covariate shift $pF(x,t) ≠ pCF(x,t)$ (Liu et al., 2022). Thus, a model learned from factuals may not perform well on counterfactuals, failing to compare ITEs across treatments.

Our framework follows Figure 1(c), which explicitly models this causal structure. The core idea is that to accurately estimate the effect of $t$ on $y$⁠, we need to learn a representation of the patient status, $z=fϕ(x)$⁠, that is disentangled from the treatment assignment, i.e. $p(z|t)=p(z)$⁠. This ensures that the representation captures the patient’s underlying prognosis, independent of the treatment they received due to selection bias. This is an underexplored area for multimodal patient status and multi-treatment scenarios. After training, the model is expected to be unbiased with respect to treatments and all treatments can be traversed to select the one with the best ITE.

The overall architecture is shown in Figure 2. It has three main components: a multimodal feature encoder with Bi-AdaIN conditioning, a disentanglement module (either adversarial or MI-based) and a survival prediction head. The gradient from the $Ladv$ or $LMI$ loss is used to update the main encoder, while the survival prediction loss trains both the encoder and the prediction head.

Patient status $x$ consists of 3D CT image volumes, which are preprocessed to $ximg∈R256×256×128$ with 128 slices to cover the head and neck region. In addition, the tabular clinical vector $xtab$⁠. The HNC treatment usually involves radiotherapy (RT) as the standard procedure, with the addition of either chemotherapy (ChemoRT) or an epidermal growth factor receptor inhibitor (EGFRI). Therefore, we have three classes of treatment $t$ (RT, ChemoRT and RT+EGFRI), which are encoded as a three-dimensional one-hot vector. We aim to predict the unbiased treatment effect $y$ (e.g. survival days), which can be censored data representing the time from the CT scan to either death or the end of the study.

Without loss of generality, we adopt the widely used 3D Swin transformer (Liu et al., 2021c) and 3D ResNet (Turnbull, 2022) from the MONAI package (Cardoso et al., 2022) as the backbone image feature extractor, $Eimg$⁠. Specifically, we used the transformer blocks in the 3D Swin transformer or the convolutional layers of 3D ResNet18 as $Eimg$⁠, which includes four transformer stages or four residual layers (each with two basic ResNet blocks). The extracted patient status feature is $z∈R8×8×8×512$⁠. We enforce the independence of $z$ and $t$ with either adversarial training or MI regularization. Then, the patient status feature $z$ is conditioned on the treatment $t$ for the final ITE inference.

Thus, a critical problem is how to seamlessly inject $xtab$ and $t$ into the feature extraction process before and after disentanglement, respectively. A simple solution is to concatenate the clinical variables $xtab$ and $t$ into two FC layers and enforce disentanglement between them. We note that the conventional concatenation of a 1D vector and a 4D (3D spatial plus channel) convolutional feature map can be challenging and inefficient. However, combining $xtab$ only at the end of the network in a single FC layer can hardly allow the model to learn their complementary relationship (Liu et al., 2024) and it is hard to learn their correlation with a limited number of subsequent parameters. In addition, adding more FC layers to sequentially incorporate $xtab$⁠, disentanglement and $t$ can introduce many more parameters and be less efficient for learning their correlations.

Therefore, we resort to adaptive instance normalization (AdaIN) layers and propose a Bi-AdaIN to inject the relevant factors into corresponding layers. Following previous successful conditional modeling works (Karras et al., 2019; Liu et al., 2020), we adopt AdaIN (Huang and Belongie, 2017) to inject $xtab$ in each transformer or convolutional block, while conditioning on $t$ in each FC layer.

Specifically, we process $xtab$ or $t$ with a three-layer multi-layer perceptron (MLP) to produce a 16-dim vector $ltab∈R16$ or $lt∈R16$ in an intermediate latent space following (Liu et al., 2024). The learned affine transformations then specialize $ltab$ or $t$ to subject-wise scale and bias scalars $tabi′=(tabis,tabib)$ and treatment-wise scale and bias scalars $ti′=(tis,tib)$ in the i-th network layer to control the AdaIN operators. Note that the dimensionality of $tabis$⁠, $tabib$⁠, $tis$ or $tib$ must match the number of feature maps in that layer. Specifically, we have the following AdaIN operation in each layer:

where $xi$ is the extracted feature map after the i-th layer. As a result, $xi$ is individually normalized and then scaled and biased using the corresponding scalar components from $tabi′=(tabis,tabib)$ and $ti′=(tis,tib)$⁠. By doing so, a stronger inductive bias of $t$ is injected into the DL model.

The model f is detailed in Figure 2. Specifically, we used four AdaIN layers for $xtab$ in both the 3D Swin Transformer and 3D ResNet18 models. In addition, the treatment $t$ is induced in two FC layers.

A core principle for mitigating treatment selection bias and covariate shift in ITE is to enforce the independence of the representation $fθ(x)$ from the treatment $t$ (Hassanpour and Greiner, 2019; Johansson et al., 2016; Shalit et al., 2017). To enforce the independence of the representation $z=fϕ(x)$ and the treatment $t$⁠, we explore two methods.

Adversarial training. This approach introduces a discriminator network $Dadv$ that attempts to predict the treatment $t$ from the representation $z$⁠. The feature encoder $fϕ$ is concurrently trained to produce representations that “fool” this discriminator, making them indistinguishable with respect to treatment. This is formulated as a min-max game. The discriminator is trained to minimize the treatment prediction loss (e.g. cross-entropy), while the encoder is trained to maximize it:

In practice, this is often implemented using a gradient reversal layer (GRL) (Ganin et al., 2016). The GRL passes the input through unchanged during the forward pass but multiplies the gradient by a negative constant $(−λ)$ during the backward pass, allowing the entire system to be trained with a single optimization step.

Mutual information regularization. While conceptually simple, adversarial training can be unstable and less effective for multiple treatments (Liu et al., 2021a, 2021b; Roy and Boddeti, 2019). To mitigate treatment selection bias, we regularize the encoder $fϕ$ by explicitly minimizing the MI between the representation and the treatment variable, $I(z; t)$⁠. The MI is defined as:

where $Pzt$ is the joint probability distribution and $Pz⊗Pt$ is the product of the marginals. Minimizing MI forces the joint distribution to approximate the product of marginals, rendering $z$ and $t$ statistically independent.

Calculating MI directly for high-dimensional continuous variables is intractable. Therefore, we use the mutual information neural estimator (MINE) (Belghazi et al., 2018), which provides a scalable and unbiased estimate of the MI by maximizing a lower bound based on the Donsker–Varadhan representation of the KL-divergence. A separate network, the MI estimator $Tψ(z,t)$⁠, is trained to provide this lower bound:

The critic $Tψ$ is trained to output higher values for samples from the joint distribution (paired $(z,t)$ samples) and lower values for samples from the product of marginals (shuffled, unpaired $(z,t′)$ samples). The expectations are approximated using Monte Carlo sampling on mini-batches:

This estimated MI, $LMI$⁠, serves as a regularization loss for the main encoder $fϕ$⁠.

The final representation $z$ is fed into a prediction head $hθ$ to estimate a risk score. Since our outcome data is right-censored, a standard regression loss like mean squared error (MSE) is inappropriate. Instead, we use a survival loss that can handle censored observations. We use the negative log partial likelihood loss from the Cox proportional hazards model (Cox, 1972).

Given a predicted risk score $y^i=hθ(zi)$ for patient i, the loss is calculated over all patients i for whom the event was observed (⁠$δi=1$⁠). The loss encourages the risk score of patient i to be higher than that of any patient j who was still at risk at the time of patient i’s event (i.e. $yj≥yi$⁠):

This loss effectively ranks patients by their risk of an event without requiring the exact event times for censored individuals.

The complete model is trained end-to-end by optimizing a composite objective function:

where $Ldisentangle$ is either the adversarial loss $Ladv$ or the MI loss $LMI$ and $λ$ is a hyperparameter balancing prediction and disentanglement.

We used the publicly available RADCURE data set (Welch et al., 2024), an extensive collection of data from 3,346 HNC patients. This data set provides pretreatment computed tomography (CT) volumes, a comprehensive set of clinical variables and long-term survival outcomes. The patient cohort has a median age of 63 years, consists of 80% males and has a median follow-up time of 5 years.

Our study focuses on comparing three primary treatment regimens. The patient cohort is divided as follows: 1,413 patients received RT with concurrent chemotherapy (ChemoRT), 72 patients received RT with an epidermal growth factor receptor inhibitor (RT+EGFRI) and 1,861 patients received RT alone. For our model, the treatment variable $t$ is encoded as a three-dimensional one-hot vector.

A standardized preprocessing pipeline was implemented for all CT volumes using the MONAI framework (Cardoso et al., 2022). All CT volumes were resampled to an isotropic voxel spacing of $1.0×1.0×1.0$ mm³ using third-order spline interpolation. Voxel intensities were normalized by clipping values to a window of [−1000, 1000] hounsfield units (HU) to focus on relevant soft tissue and bone structures. Subsequently, we applied Z-score normalization by subtracting the mean and dividing by the standard deviation of the non-zero voxels within the patient’s body mask. A fixed-size volume of $256×256×128$ voxels was cropped from the center of each normalized volume to serve as the input for our DL models.

$xtab$ includes age, sex, eastern cooperative oncology group performance status (ECOG PS), smoking pack-years, smoking status, disease site, tumor subsite, tumor size category (T), nodal involvement (N), metastasis status (M), clinical stage, pathological type and human papillomavirus (HPV) status. For continuous or ordinal variables, scalar encoding is used. For nominal categorical variables, one-hot encoding is applied. The detailed processing is given in Table 1.

To improve model generalization, we applied on-the-fly data augmentation to the training set. The augmentation techniques, also implemented with MONAI (Cardoso et al., 2022), included geometric transformations such as random rotations within a range of [−15, 15] degrees, random scaling (zooming) with a factor between [0.9, 1.1] and random flipping along the sagittal plane. We also applied intensity transformations with random adjustments to brightness and contrast and used noise injection with the addition of random Gaussian noise.

The data set was split into training (70%), validation (10%) and testing (20%) sets. Notably, we select a testing set with all non-censored data for more reliable evaluation. We also keep the distribution of treatment types consistent across all three sets. All models were implemented in PyTorch. We used the AdamW optimizer with an initial learning rate of $1×10−4$ and a weight decay of $1×10−5$⁠. A cosine annealing learning rate scheduler was used. Due to GPU memory constraints, a batch size of 4 was used. The models were trained for up to 200 epochs on NVIDIA A100 GPUs. The hyperparameter $λ$ was determined through a grid search.

Prediction accuracy: We measured the accuracy of factual outcome prediction using MSE and mean absolute error (MAE).

We use years as the unit of our survival time. Notably, these two metrics can only be applied to observed (factual) samples.

Causal effect estimation: Since ground-truth counterfactuals are unobservable, we evaluated the model’s ability to estimate causal effects using the bias-adjusted treatment effect (BATE) (Shalit et al., 2017). For any two treatments A and B, BATE measures the difference between the average treatment effect estimated by the model on the entire population and the true average effect observed in the factual data:

where $y^$ are the model predictions. A lower BATE indicates that the model is less biased by the treatment assignment mechanism.

For example, consider comparing ChemoRT (A) and RT alone (B). The second term, $(Ex|t=A[y]−Ex|t=B[y])$⁠, is the observed difference in average survival between the 1,413 patients who actually received ChemoRT and the 1,861 who received RT alone. This value is biased because these two groups of patients may have been different from the start. The first term, $Ex[y^A−y^B]$⁠, is the model’s estimated treatment effect, averaged over the entire population. To calculate this, the model takes every single patient and predicts their outcome for both treatments. BATE measures the gap between the model’s unbiased population estimate and the real-world biased observation. A small BATE means the model has successfully adjusted for the initial differences between the treatment groups.

We present a comprehensive evaluation of our proposed methods, assessing the impact of different architectural choices on both predictive accuracy and causal inference robustness. We compare our final proposed model (Bi-AdaIN + MI) against several ablative and baseline configurations using two different network backbones: 3D Swin transformer and 3D ResNet18.

The primary results are summarized in Table 2 for the 3D Swin transformer backbone and Table 3 for the 3D ResNet18 backbone. As shown, our proposed method (“Bi-AdaIN + MI”) achieves the best causal inference performance (lowest BATE) while maintaining strong predictive accuracy. Comparing the “Baseline (Concatenation)” with “Baseline + Bi-AdaIN,” we observe that incorporating our Bi-AdaIN fusion mechanism consistently improves predictive accuracy (lower MAE/MSE). This suggests that Bi-AdaIN provides a more effective way to integrate tabular clinical data. When adding disentanglement (“+ Adversarial” or “+ MI”), we observe the expected tradeoff: a slight degradation in MAE/MSE for a significant improvement in BATE. Notably, the MI-based approach demonstrates a more significant reduction in BATE than the adversarial method and the 3D Swin transformer backbone consistently yields better results than the 3D ResNet18 backbone.

To further understand the interplay between modalities and how our model makes decisions, we conducted a permutation importance analysis. This technique measures a feature’s importance by evaluating the drop in model performance when that feature’s values are randomly shuffled. A larger drop in the BATE after shuffling a feature therefore indicates higher feature importance. The results for the tabular-only model and our proposed multimodal model (using the Swin transformer backbone) are presented in Figure 3.

For the model trained only on tabular data [Figure 3(a)], the feature importance largely aligns with clinical intuition. Clinical staging information (“Stage,” “N”, “T”) and the biological marker “HPV Status” are the most predictive features. We note that Stage, T and N are the cornerstones of the TNM staging system, which is the universal standard for classifying the anatomical extent of cancer and is a primary determinant of prognosis. HPV status is a critical molecular biomarker, as HPV-associated oropharyngeal cancers have a distinct biology and a significantly better prognosis than HPV-negative cancers. ECOG Performance Status is explained as a measure of a patient’s functional well-being and their ability to tolerate aggressive treatments.

However, a dramatic shift occurs in the multimodal model [Figure 3(b)]. The “Image Features,” extracted by the deep network, emerge as the single most important predictor. Consequently, the importance of the tabular features derived from imaging, such as “Stage,” “T,” and “N,” decreases substantially. This indicates that the model learns to extract richer, more nuanced prognostic information directly from the CT scans, making the manually annotated staging information partially redundant. Conversely, complementary, non-visual features like “HPV Status,” “ECOG PS,” and “Age” retain their high importance.

We also performed an ablation study on the 3D Swin Transformer backbone to compare three strategies for injecting tabular data via AdaIN:

1. All stages (proposed): AdaIN layers applied after all four transformer stages.

2. Late stages only: AdaIN layers applied only after the last two stages (stage 3 and stage 4).

3. Final stage only: an AdaIN layer applied only after the final stage (stage 4). (Table 4)

The findings confirm that while late-stage fusion is beneficial, a deeper integration of tabular data across all stages provides superior performance. Our rationale is that early fusion allows the network to learn a hierarchical co-adaptation of features. For instance, tabular variables like tumor stage (T) and nodal status (N) have direct visual correlates that are best captured in the early-to-mid layers of the network. By modulating these features from the beginning, Bi-AdaIN facilitates a more profound and interactive fusion, enabling the model to learn richer, more robust representations that are synergistically informed by both modalities throughout the entire feature extraction hierarchy. This ultimately leads to a more significant reduction in the causal effect bias (BATE).

The hyperparameter $λ$ controls the tradeoff between the survival prediction loss and the disentanglement regularization loss. To analyze its impact, we trained our proposed “Bi-AdaIN + MI” model with varying values of $λ$ from 0.01–10.0, as shown in Figure 4. As $λ$ increases, the BATE consistently decreases, indicating more effective bias removal. Conversely, the predictive accuracy tends to degrade. We selected a value of $λ=1.0$ for the Swin transformer and $λ=0.8$ for the ResNet18, as these values offered the best balance between low causal bias and high predictive performance on the validation set.

In this work, we developed and validated a deep causal learning framework to estimate individualized treatment effects for HNC patients using multimodal data. Our primary finding is that the synergistic combination of our proposed Bi-AdaIN fusion mechanism and a MI regularization strategy can significantly mitigate treatment selection bias, leading to more robust ITE predictions compared to conventional DL models and those using adversarial disentanglement.

This synergistic effect is vividly illustrated by our feature importance analysis (Figure 3). The analysis reveals that when the model gains access to the raw image data, the deep-learned “Image Features” become the dominant predictor. Consequently, the importance of tabular features that are themselves derived from imaging, such as “Stage,” “T” and “N,” is substantially diminished. This is not because they are clinically unimportant, but because their predictive information is now more richly and directly captured by the model from the raw pixels. This explicitly demonstrates the model’s ability to handle information redundancy.

Conversely, non-visual, complementary features like “HPV Status” and “ECOG PS” retain their high importance in the multimodal model. This confirms that our framework learns to intelligently weigh and integrate different sources of information, leveraging tabular data for systemic biological and performance status while relying on the imaging data for detailed morphological and anatomical assessment. This moves beyond simple correlations to model the complex, interdependent relationships required for accurate, personalized prognosis.

A key observation from our results is the tradeoff between predictive accuracy on factual data (MAE/MSE) and causal inference robustness (BATE). This outcome, while seemingly counterintuitive, is a strong indicator that our causal framework is operating correctly. It demonstrates that the baseline models were leveraging the spurious correlations inherent in the biased observational data to improve predictions. By forcing the model to learn treatment-invariant representations, we remove its reliance on these confounding factors, which, while slightly hurting its ability to predict a biased reality, vastly improves its capacity to generalize and make accurate counterfactual predictions – the ultimate goal of ITE estimation.

Our technical innovations directly address two fundamental challenges in this domain. First, the Bi-AdaIN module proved superior to simple concatenation for integrating clinical variables with 3D imaging data. By dynamically modulating feature maps at multiple stages of the network, Bi-AdaIN facilitates a more profound and interactive fusion, which translates to improved predictive performance as seen in our ablation studies. Second, our results confirm that MI regularization is a more stable and effective approach for disentanglement than adversarial training in this complex, multi-treatment, multimodal setting.

The clinical implications of this work are significant. A reliable ITE estimation tool could revolutionize treatment planning in oncology. For any given HNC patient, our model could provide clinicians with predicted outcomes for each available treatment option (e.g. “predicted survival with ChemoRT is 3 yrs, versus 2 yrs with RT alone”). This would empower clinicians to move beyond population-level guidelines and make more personalized, data-driven decisions. The superiority of the Swin transformer backbone also suggests that advanced architectures capable of capturing long-range spatial dependencies are crucial for extracting subtle prognostic biomarkers from medical images.

We also articulate that the model’s ability to learn image features as the dominant predictor demonstrates that it successfully extracts nuanced prognostic information (related to tumor morphology, infiltration and volume) directly from the CT scan, thereby capturing the information contained in T, N and stage in a more comprehensive, data-driven manner. We emphasize that the retained importance of complementary, non-visual features like HPV Status (a systemic viral marker) and ECOG PS (a measure of patient fitness) highlights the framework’s ability to intelligently integrate different, synergistic sources of information.

Despite these promising results, we acknowledge several limitations. First, our study is retrospective and relies on the RADCURE data set. While large, this data set lacks granularity on treatment protocols, such as specific radiation dosage or chemotherapy agents. The small sample size of the RT+EGFRI group also limits the statistical power of comparisons involving this treatment. Second, while our causal framework is principled, the feature extraction process within the deep neural network remains a “black box,” which can be a barrier to clinical trust. Finally, the model requires external validation on independent data sets to ensure its generalizability before any clinical application can be considered.

Future work will proceed along several key paths. To address model interpretability, we plan to integrate visualization techniques like SHapley Additive exPlanations (SHAP) score (Lundberg and Lee, 2017) and gradient-weighted class activation mapping (Grad-CAM) (He et al., 2019) to highlight the specific anatomical regions in the CT scans that are most influential in predicting the effects of different treatments. Furthermore, we aim to extend the framework to incorporate uncertainty quantification, providing confidence intervals alongside ITE predictions to better inform clinical decision-making. Finally, the ultimate goal is to validate this framework in prospective studies, assessing its real-world impact on treatment selection and patient outcomes.

In this work, we introduced a novel deep causal learning framework for estimating individualized treatment effects in HNC from multimodal observational data. By leveraging a proposed Bi-AdaIN method for effective data fusion and a robust MI regularization strategy for disentanglement, our model successfully overcomes the challenge of treatment selection bias. Our comprehensive experiments demonstrate that this approach yields significantly more accurate and reliable estimates of causal treatment effects than baseline and adversarial methods. This study represents a critical step toward developing trustworthy AI-driven decision support tools for personalized treatment planning in oncology, paving the way for more precise counterfactual reasoning in clinical practice.

Yawen Wei and Zhen Li, both authors, contributed equally.