Prior to any augmentation or preprocessing, the full dataset was randomly split into a 70–30 ratio for training and testing. No separate validation set was used. This ensured that augmented samples derived from the training set did not leak into the evaluation pipeline.

The HAM10000 dataset exhibits significant class imbalance, with some classes (e.g., DF) containing as few as 115 images and others (e.g., NV) containing over 6,000. To mitigate this, we applied data augmentation exclusively on the training set using the following techniques:

After augmentation, each class in the training set contained 6,000 images, resulting in a balanced training dataset of 42,000 images.

Each dermatoscopic image underwent the following preprocessing steps:

Metadata was preprocessed in the following stages:

To ensure proper alignment between images and metadata, we used the image_id column as the primary key to merge records. This guaranteed that each dermatoscopic image was correctly paired with its corresponding metadata during training and evaluation.

The complete preprocessing workflow is visualized in Figure 1.

A Multi-Layer Perceptron (MLP) is a class of feedforward artificial neural networks composed of fully connected layers and nonlinear activation functions. MLPs are particularly well-suited for processing structured tabular data such as patient metadata, which lacks spatial structure and does not benefit from convolutional operations (Rumelhart et al., 1986).

In this study, we utilize a lightweight yet effective multi-layer perceptron (MLP) to process structured clinical metadata, including variables such as age, sex, and anatomical site of the lesion. The input is a 1D vector representation of all available metadata fields, which is first mapped to a 128-dimensional latent space via a fully connected layer with ReLU activation. This is followed by a second fully connected layer that expands the representation to 256 dimensions. This compact architecture is designed to extract meaningful latent embeddings from clinical features while being computationally efficient.

Figure 3 delineates the architectural pathway of the Clinical MLP, emphasizing its role in encoding patient metadata into latent space.

Residual Networks (ResNets), introduced to mitigate the degradation problem in deep architectures, extend conventional convolutional networks by learning residual mappings instead of direct functions (He et al., 2016). Rather than learning an unreferenced function, ResNets learn a residual mapping, which simplifies optimization and enables very deep architectures. Unlike AlexNet, ResNet is conceptually derived from the simpler and deeper VGG networks (Simonyan and Zisserman, 2015), but it introduces residual or skip connections that alleviate vanishing gradient issues during training.

A typical residual block includes two convolutional layers. If x is the input to a block, and W₁, W₂ are convolution kernels with non-linear activation σ, then the residual output is given by:

This results in the following expression for the block's output:

DermiResNet extends this formulation by introducing a learnable weight α for the skip connection:

Here, α ∈ ℝ is a learnable scalar parameter that determines the relative importance of the shortcut connection, enabling the model to adaptively weigh the residual and identity paths during training (Xu et al., 2024).

The architecture begins with a primary convolutional module (conv1), after which the network progresses through four sequential stages. Each stage consists of a downsampling convolutional unit and a corresponding residual unit. Across these stages, the number of feature channels is gradually increased (64 → 128 → 256 → 512), while spatial dimensions are reduced through stride-2 convolutions. Within each residual unit, two convolutional layers are used, each followed by batch normalization and LeakyReLU activation. To mitigate overfitting, dropout is applied in the later residual units (res3 and res4).

The final classification head includes an adaptive average pooling layer, flattening, a two-layer fully connected network, and a softmax activation to output a 512-dimensional feature vector. The complete architecture is visualized in Figures 4, 5.

In this layer, we fuse the extracted features from the DermiResNet (image features) and the Clinical MLP (clinical metadata features). The fusion process combines these multimodal features into a unified representation, which is then passed to a simple classifier for final prediction. The classifier consists of fully connected layers with ReLU activation functions, reducing the fused feature space to 7 output classes corresponding to the skin lesion types.

We explore several fusion techniques to combine the features effectively, including:

Detailed descriptions of these fusion techniques, including their mathematical formulations and implementation, are provided in Section 3.5

Simple concatenation involves merging features from different modalities by stacking them along a specific dimension. Given two feature vectors x ∈ ℝ⁵¹² (from the image model) and y ∈ ℝ²⁵⁶ (from the clinical model), the concatenated feature vector z ∈ ℝ⁷⁶⁸ is:

Equation 4 shows the simple concatenation of two feature vectors.

Weighted concatenation improves on simple concatenation by applying modality-specific weights. Instead of blindly stacking features, we scale each vector by a learnable scalar weight to reflect its importance (Ngiam et al., 2011; Kiela et al., 2018). Let w₁, w₂ ∈ ℝ be scalar weights applied to the 512D image vector x and 256D clinical vector y, respectively. The fused vector z ∈ ℝ⁷⁶⁸ is:

The Hadamard product, or element-wise multiplication, combines two modalities by interacting their elements multiplicatively. Since this requires equal dimensions, we first project both x ∈ ℝ⁵¹² and y ∈ ℝ²⁵⁶ into a common latent space ℝ²⁵⁶. The fused representation z ∈ ℝ²⁵⁶ is given by:

Here, x′=Linear512→256(x), and ⊙ represents element-wise multiplication (Kim et al., 2017).

Tensor fusion captures all possible interactions between modalities using an outer product, resulting in a second-order tensor (Zadeh et al., 2017). For x ∈ ℝ⁵¹² and y ∈ ℝ²⁵⁶, the fused tensor Z ∈ ℝ^{512 × 256} is:

This allows pairwise modeling of every feature from one modality with every feature from the other, at the cost of increased dimensionality.

Bilinear fusion is a feature interaction mechanism that combines information from two different modalities by explicitly modeling the pairwise multiplicative interactions between their respective features. Unlike simple concatenation or element-wise operations, bilinear fusion generates a richer and more expressive representation by learning how every feature from one modality interacts with every feature from the other. Conceptually, bilinear fusion captures second-order statistics between modalities—unlike first-order techniques such as concatenation, which only represent raw values. This is particularly valuable in multimodal tasks where the interplay between modalities is non-trivial and nonlinear (Fukui et al., 2016).

Given an image feature vector x ∈ ℝ⁵¹² and a clinical metadata feature vector y ∈ ℝ²⁵⁶, the bilinear fusion mechanism applies a set of bilinear mappings to produce a fused feature vector z ∈ ℝ²⁵⁶. Each dimension z_i of the output vector is computed using a learnable bilinear interaction:

Here, each Wi∈ℝ512×256 is a slice of the 3D learnable weight tensor W∈ℝ256×512×256, such that:

This results in a final fused representation z ∈ ℝ²⁵⁶, which is then passed through a fully connected layer for classification.

Gated fusion introduces a dynamic weighting mechanism that assigns different attention to modalities based on input features (Arevalo et al., 2020). The gating vector g ∈ ℝ²⁵⁶ is derived from a sigmoid function applied to a learnable affine combination of inputs:

The final fusion is:

Here, x′=Linear512→256(x), and ⊙ denotes element-wise multiplication. This fusion strategy adaptively prioritizes modalities at an instance level.

Self-attention enables the model to attend to the most relevant parts of a single modality. Applied independently on each modality, it transforms the sequence of features X ∈ ℝ^n×d via attention weights:

Where:

Here, WQ,WK,WV∈ℝd×d are learnable weight matrices. In this context, Q (Query) represents the features for which we want to find contextual relevance, K (Key) encodes the features to be compared against, and V (Value) holds the actual information to be aggregated. Self-attention enables each position in the feature sequence to attend to all positions, allowing the model to learn intra-modality dependencies. It is used before fusion to enhance modality-specific representations (Vaswani et al., 2017).

Cross-attention aligns features between modalities by using one modality as query and the other as key-value pairs. For image features X ∈ ℝ^{1 × 512} and metadata features Y ∈ ℝ^{1 × 256}, the query is derived from image features and the key/value from metadata:

Where:

In cross-attention, Q_x represents the query derived from the primary modality (e.g., image features), which is seeking relevant complementary information. K_y and V_y are the key and value vectors derived from the secondary modality (e.g., clinical metadata), where the key determines alignment and the value contributes the corresponding context. This fusion allows modality X to selectively attend to modality Y, creating cross-modal representations that are aligned and context-aware (Lu et al., 2019).

The hardware and software specifications used for the experiments are summarized in Table 4.

The specific training configuration, which outlines hyperparameters and other details, is documented in Table 5. Each model was trained for an estimated duration of 2 h.

Grad-CAM is a powerful technique for visualizing the regions of an image that are most influential in a model's decision. It extends the Class Activation Mapping (CAM) approach by using gradient information to weight the importance of feature maps, making it applicable to a wider range of architectures, including those without global average pooling layers.

Mathematical formulation: Let A^k be the activation map of the k-th channel in the target convolutional layer, and let y^c be the score for class c. The weight αkc for the k-th channel is computed as the global average of the gradients of y^c with respect to A^k:

where Z is the number of pixels in the activation map. These weights capture the importance of each feature map for the target class.

The Grad-CAM heatmap L^c is then obtained by a weighted combination of the activation maps, followed by a ReLU function:

The ReLU function ensures that only positive influences are considered, as negative values are not relevant for the target class (Selvaraju et al., 2017).

Grad-CAM indicates areas of the image that the model considers most significant for its prediction. For instance, in images obtained with dermatoscopy, these areas may be lesion boundaries, texture, or color transitions. By projecting these areas, Grad-CAM opens a window to the model's “thinking process,” allowing physicians to ensure that the AI system is concentrating on clinically relevant features (Jaworek-Korjakowska et al., 2021).

To understand the decision-making process of the top-performing model, which utilizes a cross-attention fusion approach, Grad-CAM visualizations were employed. This technique facilitates the interpretation of the model's predictions by identifying key areas in dermatoscopic images that most strongly influence the results. Such post-hoc explainability is especially valuable in medical contexts, as it allows clinicians to confirm that the model's predictions are grounded in clinically significant features.

The cross-attention fusion mechanism enhances the interpretability of the model by dynamically aligning and integrating features from different modalities or scales. Unlike simpler fusion techniques, cross-attention computes attention weights between modalities, enabling the model to emphasize the most diagnostically relevant interactions. While Grad-CAM visualizations help interpret the image-based decision process by highlighting regions of importance in dermatoscopic inputs, they are limited to visual modalities. To gain a holistic understanding of the model's reasoning, particularly for clinical metadata, we later discuss Clinical Feature Relevance Scores, which complement Grad-CAM by providing insight into how non-image features influence predictions. This combined interpretability offers a more comprehensive explanation of the model's diagnostic behavior.

Figure 8 provides an illustrative example of a Grad-CAM heatmap overlaid on an input dermatoscopic image. The model predicts the lesion as BCC with a confidence score of 1, focusing primarily on the lesion's irregular borders and regions of heterogeneous pigmentation. These features are clinically significant as BCC are characterized by asymmetry, border irregularity, and color variation (Puckett et al., 2025).

In this specific example, the heatmap demonstrates that the model successfully identifies features associated with malignancy while ignoring irrelevant artifacts, such as hair strands and uniform skin areas. The model's ability to focus on clinically relevant features is a direct outcome of the cross-attention fusion mechanism, which dynamically aligns and weights features from different modalities. By doing so, the model achieves a higher level of alignment with dermatological practices, where lesion borders and internal variations are critical for diagnosis.

As discussed earlier, the model occasionally misclassifies BKL as MEL due to overlapping visual characteristics. One contributing factor is that both lesion types can exhibit irregular pigmentation, asymmetric structures, and varying border definitions, which are also key diagnostic markers for MEL. This confusion is particularly evident in cases where keratosis presents with darker pigmentation and irregular borders, mimicking features of malignant lesions.

Figure 9 illustrates such a misclassification, where a BKL lesion has been incorrectly predicted as MEL with a confidence score of 0.46. The left panel shows the original image, while the right panel presents the corresponding Grad-CAM heatmap, highlighting the regions that influenced the model's decision.

From the heatmap, it is evident that the model assigns high importance (red/yellow regions) to darker pigmented areas and irregular structures within the lesion. This suggests that the model's decision boundary between benign and malignant lesions is influenced primarily by pigmentation and border irregularity, which are not always exclusive to melanoma.

Beyond this example, Grad-CAM visualizations across a wide range of test cases reveal consistent patterns in the model's behavior. The model often prioritizes irregular lesion borders and regions of color variation, which are crucial for distinguishing malignant lesions from benign ones. In cases of BCC, the heatmaps highlight central regions with ulceration or shiny surfaces, further reinforcing the model's alignment with clinical indicators. This interpretability is invaluable for real-world deployment, as it allows clinicians to confirm that the model's focus areas correspond to meaningful diagnostic features.

The insights provided by Grad-CAM are directly correlated with the model's strong quantitative performance. For example, the regions identified by the heatmaps frequently align with features responsible for the model's high sensitivity and specificity, particularly in challenging cases of melanoma and basal cell carcinoma. This combination of performance metrics and visual interpretability underscores the potential of the cross-attention fusion-based model as a trustworthy diagnostic aid.

Additional Grad-CAM visualizations for Cross-attention are provided in Figure 10, showcasing further examples of model attention across different lesion types and cases.

To assess the instance-specific relevance of clinical features, we employed an attribution-based approach using Integrated Gradients (IG) (Sundararajan et al., 2017). IG is a path-based attribution method that quantifies feature importance by computing the integral of gradients along an interpolation path from a baseline input to the actual input.

Given an input clinical feature vector x ∈ ℝ^d and a baseline vector x′ ∈ ℝ^d, the integrated gradient for the i^th feature is computed as:

where f(·) represents the model's prediction score for the target class. The integral is approximated using a summation over discrete steps:

where S is the number of interpolation steps.

For each instance, we set the baseline x′ as a zero vector, representing the absence of clinical features. The integrated gradients were computed over S = 50 steps, and the absolute values of the attributions were taken as feature importance scores. These scores were then normalized to the range [0, 1] to facilitate interpretability.

The resultant feature relevance scores highlight the contribution of individual clinical attributes to the model's prediction. Features with higher attributions indicate stronger influence on the classification decision, thereby providing insights into the model's reliance on clinical metadata.

For a specific instance of NV, we analyzed the relative importance of clinical features using our best-performing Cross-Attention fusion model. As shown in Table 8, the top features with high importance include localization (scalp, foot, neck), and diagnostic methods (consensus, follow-up, confocal). These features likely play a significant role due to the distinct characteristics of lesions in these locations and the reliability of the diagnostic methods. Features with moderate importance, such as localization (face, genital) and diagnosis type (histopathology), contribute to predictions but are less critical. Interestingly, age and certain locations (e.g., ear, lower extremity) show low importance, suggesting they have minimal influence on the model's predictions for this instance. The balanced impact of sex (both male and female) indicates that while it influences outcomes, it is not among the strongest predictors. Overall, localization emerges as the most relevant feature, aligning with real-world clinical practice where lesion location is a key diagnostic factor.

Across various lesion types, we observed that localization consistently emerged as the most important clinical feature, underscoring its critical role in skin lesion diagnosis. This aligns with real-world clinical practice, where the anatomical location of a lesion is a key diagnostic factor due to its correlation with sun exposure, skin type, and lesion characteristics.

In cases of MEL, localization (e.g., back, face) was the top predictor, reflecting the higher prevalence of malignant lesions in sun-exposed areas. Additionally, age and gender showed moderate influence, consistent with their established roles as risk factors for melanoma. Older patients and males exhibited a higher likelihood of malignancy, further validating the model's alignment with epidemiological trends (Raimondi et al., 2020; Garbe et al., 2021).

For BCC, the model again prioritized localization (e.g., face, neck), as these areas are most susceptible to UV damage (Marzuka and Book, 2015). Diagnostic methods such as histopathology also played a significant role, as BCC is often confirmed through biopsy. Interestingly, age showed a stronger influence for BCC compared to other lesion types, likely due to the cumulative effect of UV exposure over time (Seretis et al., 2025).

In cases of AKIEC, localization (e.g., face, scalp) remained the most important feature, as AKIEC lesions are strongly associated with chronic sun exposure (Reinehr and Bakos, 2019). The model also highlighted the importance of diagnostic methods (e.g., follow-up, confocal microscopy), reflecting the need for repeated evaluations to monitor these precancerous lesions.

For DF, a benign lesion, localization (e.g., lower extremities) was again the dominant feature, while age and gender had minimal influence. This suggests that the visual appearance and location of DF lesions are more critical for diagnosis than demographic factors (Han et al., 2011).

Finally, in cases of VASC, localization (e.g., face, trunk) was the most influential feature, as these lesions often appear in specific anatomical regions (Mulligan et al., 2014). The model also relied heavily on dermoscopic examination, highlighting the importance of visual data for diagnosing vascular lesions.