The heterogeneity of the symptoms and the use of subjective clinical ratings make it challenging to diagnose PD with any degree of dependability at an early stage of the condition. Despite the high accuracy of unimodal AI methods based on speech, gait, or handwriting, they do not tend to be generalized across all populations and have no interpretability, which makes them less appealing to clinical users. In addition, the majority of current models are susceptible to noisy or incomplete data and do not offer much information regarding their decision-making process. This poses a critical gap towards the creation of a robust, adaptive, and explainable multimodal framework that incorporates various PD biomarkers, can withstand imperfect data, and can show interpretable information that can be used to make reliable clinical decisions.

As discussed in Figure 1, the initial step of the pipeline is to put the raw gait signals that are usually of different lengths and magnitudes across subjects into a uniform form. All recordings are broken up into temporal blocks of certain fixed lengths of 2 s (200 samples at 100 Hz) and still maintain the temporal continuity. In case one of the segments is shorter, it is zero-padded to preserve the dimension. A normalization of these windows is then made based on z-score scaling which rescales each feature by removing the mean and dividing by the standard deviation. This modification makes sure that the different sensors or trials with different absolute magnitudes do not overtake the training process. This normalization formula can be expressed as in Equation 1:Vnorm=q−meanSD(1)where Vnorm is the normalized value, q is the raw input value, mean is the average of the windowed interval, and SD is the standard deviation of the windowed interval.

This is applied to normalize the data so that its variance is 1 and its mean is 0, which improves the stability of training and makes the features of other sensors comparable.

The second step is to use temporal convolutions to get local fine gait dynamics. Convolution is nothing but a weighted average of time samples around it, which has gone through a non-linear activation (ReLU). This enables the network to automatically learn filters which identify stride changes, tremor bursts, or asymmetries. The temporal convolution at a time step is of a basic form which is expressed as in Equation 2:ht=a(∑i=0k−1⁡pi.el+i+(bias))(2)where pi are kernel weights, el+i are input samples within the convolutional window, bias is a bias term, and a is the activation function.

Instead of using local patterns, dilated convolutions are used, in which gaps are inserted between sampled inputs. This increases the receptive field (but does not raise the cost of computation) to allow the network to trace longer gait cycles. The convolution of the layer is dilated. l is formalized as in Equation 3:ht(l)=f(l)(hv−c(l−1),…,ht(l−1))(3)where c is the dilation factor, meaning that the filter skips every c samples from the previous layer.

By stacking layers with different dilations, the model can capture both short-term stride irregularities and long-term gait cycles, as defined in Equation 4.yt=ht(l)+p(id)(4)where p(id) is the residual mapping of the input.

This mechanism ensures that the network does not lose fine-scale gait information while learning higher-order abstractions.

The model generates rich temporal features for the maps after several convolutional layers. These maps are, however, different in length depending on the size of the window and convolution parameters. Global average pooling is used to compress this onto a fixed-length vector that is a representation of the gait segment. This is utilized to compute the average activation at each time step giving a compressed embedding. The operation is constituted as in Equation 5:f=1Q∑t=1Q⁡ht(5)where Q is the temporal length of the segment and f is the resulting embedding.

This latent vector captures dominant gait dynamics (e.g., stride frequency or imbalance) in a low-dimensional space, making it suitable for clustering and classification.

The model is trained, as an autoencoder, to make sure that important information is not lost in the learned embedding. The encoder is fed by inputs and transforms them into embeddings, and the decoder tries to rebuild the original input signals based on the embeddings. The reconstruction goal causes the embeddings to capture gait features of significance as opposed to noise. This is formalized as in Equation 6:m^k=gd(u)(6)where gd is the decoder function, u is the latent embedding, and m^k is the reconstructed signal.

Minimizing the difference between m^k and the true mk ensures that embeddings encode physiologically relevant gait dynamics.

After generating embeddings, they are clustered into groups to determine natural clusters of gait patterns. The redundancy is eliminated using dimensionality reduction (through PCA) after which the K-means clustering is carried out.

To estimate the best number of clusters, K, the silhouette score is calculated which weighs the intra-cluster cohesion and inter-cluster separation. The definition of this metric is shown in Equation 7:s=c−dmax(c,d)(7)where d is the mean distance of points of the same cluster and c is the mean distance between the current cluster and the other nearest cluster.

When the silhouette score is high, the clustering will be done in a better way, and it helps to find meaningful subgroups such as the case of healthy gait and the case of Parkinsonian gait.

A shallow neural network classifier is trained on the embeddings with pseudo-labels obtained through the process of clustering to transform the unsupervised cluster assignments into a predictive tool. The softmax function is the calculation of the probability of classes by the classifier, as shown in Equation 8:m^l=exp⁡(p⋅gl+bl)∑j=1l⁡exp⁡(p⋅gl+bl)(8)where g_c and b_c are classifier parameters, l is the total number of clusters (classes), and m^l is the probability that input p belongs to class c.

The model is optimized by reducing the non-entropic loss which penalizes variance between actual labels and probabilities. This is defined in Equation 9:L=−∑c=1C⁡ohlog(y^c)(9)where o_h is the actual text (one-hot encoded).

This guarantees that the classifier learns decision boundaries that match gait clusters and thus it can classify unseen gait segments with a lot of strength.

SHapley Additive exPlanations (SHAP) is used to provide interpretability. SHAP is based on the cooperative game theory, and the score of significance to each feature is calculated by measuring its marginal contribution to any feature subset. This is regularized in Equation 10:pi=∑S⊆F{l}∣A∣!(∣W∣−∣A∣−1)!∣W∣where p_i is the importance of feature l, W is the full feature set, A is a subset excluding feature l, and w(⋅) is the model prediction.

As a complementary interpretability technique, Integrated Gradients (IG) explains the prediction of a neural network using the input features of the network by summing gradients along a linear path between a baseline (e.g., zero input) and the actual input. IG contrasts with SHAP, which is a combinatorial-based method, and operates on a direct measure of the sensitivity of outputs to the change of input values. Equation 11 defines the attribution as follows:IGi(k)=(ki−ki′)⋅∫01∂f(k′+v(k−k′))∂kidv(11)where k is the actual input, k′ is the baseline input, and v is an integration parameter.

This intrinsic accumulation builds gradients along the path from input to baseline, pointing out which sensors and time points were most responsible for the ultimate decision. Within gait analysis, IG has been used to localize particular phases of a stride cycle (such as toe-off or heel-strike) that were most important for the classification of Parkinsonian gait.

In Parkinson’s disease (PD), some of the first impairments include inappropriate loudness, monotone tone, inaccurate articulation, or a hoarse voice. The incorporation of raw audio recordings of sustained vowel phonation into the proposed system is because of these characteristics to achieve the highest diagnostic capability. As discussed in Figure 1, these sound signals are uniformly sampled (e.g., 16 kHz) and represent a non-accented vocal prolongation (e.g., prolonged/a/sound) by each subject, so that any phonatory deviations can be isolated without lexical or linguistic variations.

The one-dimensional signal x(t) is transformed to a two-dimensional time–frequency representation through the short-time Fourier transform (STFT) which takes the Fourier transform of windowed slices of the signal. Mathematically, the STFT of a signal is computed as shown in Equation 12:o(g,v)=∑n=−∞∞⁡o(l)⋅w(l−g)⋅e−jwt(12)where o(g,v)is the STFT output, o(l)is the input signal, w(l−g)is the windowing function, and w is the angular frequency.

The result is a spectrogram in complex values showing the frequency content vs. time. To place a frequency scale more in accord with human perception of sound, the frequency axis on the spectrogram is converted using the Mel scale. The Mel-frequency m corresponding to a linear frequency f (in Hz) is calculated according to Equation 13:mel=2595⋅log10(1+l700)(13)where mel is the Mel frequency and l is the linear frequency in Hz.

Such perceptual scaling focuses on low-frequency bands of voice pathologies. The magnitude spectrogram is then squeezed in the logarithmic form as shown in Equation 14:Elog(h,k)=log(E(h,k)+b)(14)where Elog(h,k)represents the log-Mel spectrogram; E(h,k), magnitude Mel spectrogram; and b, small constant for stability.

The resulting log-Mel spectrogram is a suitable two-dimensional input to the EfficientNet-B0 architecture which is a convolutional feature extractor. Instead of derived handcrafted shallow convolutional kernels, EfficientNet-B0 uses an optimized sequence of convolutional blocks that uses depth-wise separable convolutions and squeeze and excitation modules. The feature maps obtained after the convolutional processing are pooled by using global average pooling to obtain a compact latent embedding, z. The resulting classification output is hence obtained as shown in Equation 15:K^=d(Oa+f)(15)where K^represents the predicted class probabilities; a, latent feature vector obtained after global average pooling; O, weight matrix of the classifier; f, bias term; and d, softmax activation function.

To make it more robust and insensitive to local spatial variations, the extracted feature maps from EfficientNet-B0 are aggregated with global average pooling. This function takes each feature map and compresses it into one representative value leading to a small latent embedding vector. z that captures the general discriminative information of the input spectrogram. The classification output is then optimized by using the cross-entropy loss function, which is written as in Equation 16:C=−∑k∈{PD,HC}⁡hklogh^k(16)where C represents the cross-entropy loss; h_k, ground-truth label for class c; and h_k, predicted probability for class k.

EfficientNet-B0 uses alternating layers of convolutional and pooling operations to capture information that is relevant to PD into voice features. In addition, to promote the interpretability of predictions, Grad-CAM is used. At a minimum, Grad-CAM calculates the importance weight αk of each feature map. The global average A^k of the gradients given by Equation 17:Ok=1Z∑i⁡∑j∂pl∂Aa,bk(17)where O_k represents the importance weight of feature map k; p^l, score for class l; Aa,bk, activation at location (a,b); and Z, number of pixels.

The class activation map is then computed using Equation 18:L=ReLU(∑k⁡olmg)(18)where L represents the heatmap; o_l, weights from Equation 12; m^g, feature map; and ReLU keeps only positive contributions. By means of this approach, clinicians will be able to understand predictions based on frequency divisions and time points that reflect the presence of any abnormalities in voice.

As discussed in Figure 1, motor disturbance in handwriting in the form of tremor, micrographia, and uneven stroke patterns is another identified symptom of PD. To capture such impairments, the proposed system will entail a handwriting analysis pipeline that will process drawn spirals that are submitted as images. The collection of these spirals is normally done on a digital tablet or scanned on paper where the participant is asked to trace over an Archimedean spiral. The resultant images capture the neuromuscular coordination, such as fine motor abilities, stroke fluency, and tremor strength. Raw images are preprocessed by conversion to grayscale, resizing (to 224,224 pixels), and contrast normalization queries. Before analysis, the images have been processed as follows: converted to grayscale, resized (to 224,224 pixels), and contrast has been normalized.

As an ancillary geometrical landmark, the curvature, kappa, of the spiral stroke can be derived to measure local tremor. In a two-dimensional curve given parametrically [x(t),y(t)], the curvature is computed as in Equation 19:h(k)=w′(k)g2(k)−g′(k)w2(k)(w′(k)2+g′(k)2)3/2(19)where h(k)represents the curvature; w′(k),g′(k), first derivatives; and w′′(k),g′′(k), second derivatives.

The now preprocessed spiral's image is then fed into a ResNet-50-based deep convolutional neural network (CNN) that is a stack of residual blocks. In every block, several layers of convolution that extract spatial features are used. The convolutional layer is a 2D convolution layer, and outputs (i,j) activation is calculated according to Equation 20:La,b=S(∑q,wWq,w×xi+q,j+w+e)(20)where La,brepresents the output at (a,b); xa+q,b+w, input patch; Wq,w, kernel weights; e, bias term; and S, activation function.

ResNet adds identity shortcuts: The output of a layer can be combined directly with the input of the next layer in the chain, as shown in Equation 21:r=F(g,{Ai})+o(21)where r represents the block output; o, input; and F(g,{Ai}), residual mapping.

This architecture alleviates the vanishing gradient issue. The non-linearity in each block is due to the ReLU activation function as indicated in Equation 22:A(d)=max(0,d)(22)where A(d)represents the ReLU activation and d represents the input.

The last convolutional image is then compressed using global average pooling to give a condensed feature vector. As an input, this vector is fed into the multimodal fusion and classification module. To interpret the decision process of the model, we use Grad-CAM on the final convolutional layer. Grad-CAM interprets the classification score by calculating the gradient of the class label y_c relative to feature maps A^k to identify the significance weight of each map, as mentioned in Equation 17. The resultant heatmap is a weighted summation of feature maps, as indicated in Equation 18. When used on spiral drawings, Grad-CAM can localize motor abnormalities found in stroke with sharp curvature changes or vibration.

The handwriting, gait, and speech datasets used in this study originate from independent sources and do not contain subject-level correspondence. Consequently, direct subject-wise fusion is not feasible. To address this, multimodal fusion is performed at the feature-distribution level. For each modality, deep feature representations are independently extracted and grouped according to diagnostic labels (PD or healthy control). Within each class, normalized statistical descriptors are computed to represent modality-specific feature distributions. These class-consistent representations are then concatenated to form trimodal feature vectors, which are used as inputs for XGBoost-based classification.

This fusion strategy ensures label consistency across modalities while enabling effective multimodal learning without requiring subject-level alignment.

Although unimodal analysis gives good insights into individual domains of symptoms in PD, they do not give the complete picture of the heterogeneity of Parkinson’s disease. To counter this, our system uses a multimodal fusion architecture to combine the voice, handwriting, and gait information. As discussed in Figure 1, a deep learning model that has been trained on feature vectors specific to each modality (voice, hand, gait) produces feature vectors v_voice, v_hand, and v_gait separately. All these are joined together in one single vector as presented in Equation 23:vfused=[vvoice∥vhand∥vgait](23)where vfusedrepresents fused feature vector; vvoice, voice features; vhand, handwriting features; vgait, gait features; ∥, concatenation.

This approach allows the model to be cross-modal learning. The fused vector is next forwarded to extreme gradient boosting (XGBoost).

XGBoost is an additive booster. At iteration t, the forecast is added to all trees, as in Equation 24:g(f)=∑u=1f⁡du(l)(24)where g(f)represents prediction at iteration f and du(l)represents weak learner (tree) at iteration u.

The full regularized objective function is given in Equations 25 and 26:T=∑i=1n⁡E(hi,h^i)+∑k=1t⁡R(Su)(25)Ω(f)=oM+12b∥g∥2(26)where T is the objective; E, loss; hi, true label; h^i, prediction; R(S), regularization term; o, complexity parameter; M, number of leaves; b, regularization strength; and g, leaf weights.

The preponderance of the 10 alternatives in the three models, which can be calculated by averaging class probabilities (soft voting) as in Equation 27, is as follows:PC=13∑m=13⁡P(m)(27)where PCis the final probability for class c and P(b)is the probability from model b.

The transparency is guaranteed by hiding SHapley Additive exPlanations (SHAP). To explain each prediction, it is necessary that SHAP values be assigned to each feature, as the attribution score using the Shapley values formula, mentioned in Equation 10.

The detailed training and evaluation process of the framework is further explained using the structured pseudocode, which explains data loading, encoder building, multimodal fusion, classifier training, validation, and explanation. Together, the architectural design (Figure 1), representation of the flow (Figure 2), and description of the pseudocode convey the methodology in detail for robust and interpretable prediction of Parkinson's disease in the three modalities of speech, gait, and handwriting.

The algorithm is based on two algorithms. (Unimodal) Algorithm 1 trains separate models of the three modalities: speech, gait, and handwriting; each modality is preprocessed (spectrograms, sliding windows, or images), encoded with EfficientNet-B0/AE + multilayer perceptron (MLP)/ResNet50, and optimized with cross-entropy loss. (Trimodal) Algorithms 2 fuses the embeddings of those three modalities (using PCA in the speech case and SMOTE in the balance case) and concatenates them into a common feature.

The interaction between the proposed system components is integrated into the total execution flow. As depicted in Figure 2, the framework sequentially goes through the acquisition of data from multiple modalities, extraction of features from the individual modalities by specialized deep learning pipelines, multimodal fusion, and classification to PD or healthy control. This end-to-end loop is augmented with explainability modules that give modality-specific and fusion-level insights that help in improving clinical explainability. Together with the workflow representation (Figure 2) and the description of the structured pseudocode necessary for the overall methodology, the architectural design (Figure 1) fully outlines our methodology for robust and explainable PD detection across speech, gait, and handwriting modalities.

To further ensure generalization and reduce dependency on a single train–test split, a fivefold stratified cross-validation strategy was employed for all unimodal and multimodal models. In each fold, class distributions were preserved to avoid bias. The reported performance metrics represent the mean and standard deviation across folds. The trimodal fusion model consistently achieved high accuracy with low variance across folds, indicating strong generalization and resistance to overfitting.

In this study, we employed early feature fusion followed by an XGBoost classifier for PD detection from multimodal handwriting, gait, and speech features. Early fusion was chosen to allow the model to learn joint representations across all modalities simultaneously, thereby capturing complementary interactions between heterogeneous features such as handwriting images, gait time-series, and speech acoustic signals. Our ablation studies indicated that early fusion outperformed late fusion or hybrid strategies in terms of classification accuracy and macro F1-score, highlighting its effectiveness in leveraging cross-modal correlations. The XGBoost classifier was selected due to its robustness and ability to handle tabular heterogeneous features efficiently, especially with limited sample sizes typical in biomedical datasets. XGBoost provides fast training, strong generalization, and inherent feature importance analysis, which complements the explainable AI (SHAP) approach employed in this work. We note that, in the current implementation, the fusion mechanism involves static concatenation of features; no dynamic weighting or gating is applied.

This study utilizes three publicly available benchmark datasets corresponding to handwriting, speech, and gait modalities for PD analysis. These datasets are modality-specific and independently collected, with no subject-level overlap across modalities, reflecting realistic clinical scenarios where complete multimodal data from the same subject may not always be available.

Handwriting dataset: Handwriting data were obtained from the Handwritten PD Spiral Dataset available on Kaggle. The dataset consists of 3,264 spiral drawing samples acquired from PD patients and healthy control subjects. Participants were instructed to trace Archimedean spiral patterns using digitized input, capturing fine motor impairments such as tremor, micrographia, and irregular stroke patterns. These spiral drawings serve as established biomarkers for assessing neuromuscular degradation associated with PD.

For all unimodal experiments, a subject-wise data splitting strategy was employed to prevent information leakage between training and evaluation sets. Specifically, 70% of the subjects were used for training, 15% for validation, and 15% for testing. Data were partitioned using a fivefold stratified cross-validation scheme, ensuring balanced representation of PD and healthy control subjects in each fold. Model performance metrics were reported as the average across all folds, providing a robust and unbiased evaluation. In the multimodal experiments, feature embeddings extracted independently from each modality were aligned at the feature level. Since datasets were modality-specific with no shared subjects, multimodal fusion was performed using concatenated latent representations rather than subject-level pairing. The same cross-validation protocol was consistently applied to unimodal and trimodal models to ensure fair comparison. A consolidated summary of dataset characteristics, including dataset size, feature types, acquisition conditions, and clinical relevance, is provided in Table 2.

The handwriting-based ResNet50 model confusion matrix in Figure 3 demonstrates an outstanding performance of the model in the classification of healthy individuals and patients with Parkinson. Among all of the healthy samples, 1,572 samples were recognized, and 60 samples were incorrectly recognized as Parkinson, that is, a very high specificity. Equally, 1,394 Parkinson samples were correctly identified whereas 238 were wrongly identified as healthy which shows high sensitivity but with a marginally higher misidentification than Healthy cases. Altogether, the model is capable of measuring biomarkers that are related to handwriting, including tremor-induced deviations and curvature, with a high degree of accuracy, making handwriting a powerful modality to detect Parkinson.

The accuracy curve shown in Figure 4 shows a gradual improvement in accuracy after 50 epochs, and both training and validation accuracy is 88%–89%. The fact that the two curves are very close suggests that the handwriting-based model is always able to work on the unknown validation data. These findings confirm that the handwriting signatures including spiral curve and irregularities caused by tremors are a powerful and consistent marker in detecting the existence of Parkinson’s disease.

The loss curve demonstrated in Figure 5 shows that both training and validation sets steadily decrease, and before reaching the 0.29 level, they leveled off after approximately 40 epochs. The intersection between training and validation loss emphasizes the fact that the model has good generalization without experiencing much variance or underfitting. This gradual transition is another confirmation that the selected CNN design can capture minor impairments of handwriting, whereas data augmentation must have contributed to a decreased risk of overfitting.

The areas where the model focuses the most when predicting, as shown in Figure 6, are emphasized by the spiral handwriting with the Grad-CAM overlay. The heatmap outlines unusual strokes and tremor-induced anomalies in following the spiral line that are unique features of the fine motor impairment caused by Parkinson. Such a visualization ensures that the model is also paying attention to clinically significant handwriting properties, which validates not only the reliability of the learnt features but also the readability of the detection structure.

The speech training loss plot presented in Figure 7 indicates a definite downward slope over 30 epochs, and it begins at approximately 0.73 and decreases to about 0.55. This gradual reduction indicates the capability of the model to incrementally train its parameters and discriminative speech attributes, including pitch variation, jitter, shimmering, and spectral patterns. The progressive enhancement indicates that the speech-based model is learning significant vocal biomarkers successfully; thus, it is an efficient single-modality component which can be used to detect Parkinson’s disease.

The gait data feature correlation heat map illustrated in Figure 8 shows that there are close relationships between various temporal and force-based measurements. Stride time, stance time, and swing time have high positive correlations with each other indicating that these parameters are naturally interdependent in the walking cycle. On the same note, mean force, impulse, and root mean square (RMS) values are found interacting in strong positive terms, meaning that they represent associated factors about gait dynamics and ground reaction forces. There are negative relationships between cadence and dominant frequency vs. stride time and stance time, which stands in agreement with the reality that the greater the cadence, the shorter the stride time is likely to be. On the whole, this correlation structure demonstrates effective biomechanical correlations between gait variables, which proves that the extracted variables represent mutually independent dimensions of motor performance that are useful in detecting Parkinson's Disease.

The loss curve as demonstrated in Figure 9 indicated by the autoencoder (AE) training records a significant reduction in mean squared error (MSE) in the early epochs, which declines between the mean squared error (MSE) of about 0.52 and 0.2 within the first five epochs. Following this rapid decrease, the loss decreases progressively and levels off at 0.07–0.08 in epoch 40. This gradual convergence actually signifies that the AE has successfully mastered compressed representations of the gait signals, which represent the key patterns of the gait signals at the lowest error. These low and constant reconstruction losses prove the appropriateness of the learned embeddings to downstream clustering and classification tasks in Parkinson’s disease detection.

The training and validation loss curves of the classifier shown in Figure 10 indicate that the first epoch refers to a sharp decrease in the training loss, which decreases above 0.20–0.03, and the validation loss follows almost the same pattern. The two losses then approach steadily to values near 0.01, and in all the further epochs, they are consistently at these values. Such a narrow gap indicates that the learned embeddings which the classifier is fine-tuning are not only effective in their optimization but also generalize well to unseen data. The high stability and extremely low loss rates support the claim that the representations obtained by the autoencoders are highly discriminative feature space that would result in good classification.

The curve of validation accuracy curve shown in Figure 11 demonstrates high results with the values of 99.7% and higher throughout all the epochs, and the maximum at approximately 99.7%. The little oscillations indicate the natural changes during training; however, they do not exceed a very slight margin, which proves stable generalization. In the meantime, the learning rate is fixed at about 0.001, which means that the model has converged with a marvelous performance without the need to adapt the learning rate. The combination of close accuracy of validation and stable optimization dynamics underscores the usefulness of the autoencoder-based embeddings in offering a deep separable feature space in which classification can be done.

The distribution of the cluster count illustrated in Figure 12 indicates two distinct clusters that are well separated with Cluster 0 having approximately 6,800 samples and Cluster 1 having approximately 8,700 samples. This means that the autoencoder-generated embeddings when clustered spontaneously form two major clusters which coincide with healthy controls and Parkinson patients. The comparatively well-balanced distribution implies that the learned feature space has significant differences in gait patterns between the two groups indicating the usefulness of the unsupervised clustering step in distinguishing motor signatures of Parkinsonism as compared to normal gait patterns.

The Integrated Gradients (IG) heatmap demonstrated in Figure 13 mentions the role played by various vertical ground reaction force (VGRF) channels and summed totals of the force signals throughout time during a gait trial that is identified as Parkinson's disease. The noticeable attributions are limited to particular left and right foot sensors (i.e., VGRF1left, VGRF2left, and VGRF1 irregularity), plus in particular time intervals, which implies that this model depends on slight biases and deviations in left and right foot pressure distribution to make its predictions. These patterns with emphasis are in line with clinically observed gait deficits in Parkinson patients including decreased force symmetry and non-regular step dynamics, which prove that the model reflects significant biomechanical characteristics associated with the condition.

The SHAP summary plot is used to show the relative significance of features in predicting the model using gait-based Parkinson classification as discussed in Figure 14. Characteristics including 115, 31, 68, and 3 have the highest SHAP values, meaning that they have the greatest contribution to the decision boundary. The color gradient indicates that high (pink/red) and low (blue) feature values have the potential to move the prediction one way or the other, implying that minor differences in gait signals are reflected by the model. The fact that the points around zero are concentrated in a smaller number in the case of less significant features also underlines the point that only a few extracted embeddings are significant in classification. Such an interpretation outcome also confirms that the model is not making predictions based on random noise but on particular and important gait-related characteristics which lends more credibility to its prediction of the ability to differentiate between the Parkinson patients and healthy controls.

As indicated in Figure 15 by the normalized confusion matrix, the trimodal fusion model has almost balanced and high performance on the two classes. In healthy subjects, 89% were rightly classified, and the rest were inaccurately classified as Parkinson and it was only 11%. In the same way, with the subjects of Parkinson, 90% of the subjects were correctly identified with only 10% misclassified as Healthy. This equal balance means that the model has a good sensitivity and specificity at the same time, meaning that it is robust in dealing with the two classes without high bias. The fact that the results were very similar in the categories indicates the validity of the trimodal framework in the clinical environment.

The heatmap of correlation shown in Figure 16 among modalities indicates a strong correspondence between the handwriting (image) and gait features, which have a value of +1 and imply that these modalities represent complementary structural and motor features of Parkinson's disease. Conversely, speech has a negative relationship with image (−1) and a less significant negative relationship with gait (−0.19) showing that speech is capturing distinct non-redundant biomarkers. This isolation has shown that multimodal fusion is worthwhile since the high-correlated modalities (gait + handwriting) are complemented with a different signal (speech) to increase the robustness and general diagnostic strength.

The trimodal early fusion approach as discussed in Table 3 had an accuracy of 92% with an equal level of precision and recall as when each of the three modalities was used separately. Although handwriting was a powerful feature (91%) and gait was well-accurate (90%) and poorly balanced, whereas speech was slower (74%), the trimodal model was more reliable and generalized, which is why it is the most efficient method of detecting Parkinson.

To assess the statistical reliability of the reported performance metrics, 95% confidence intervals (CIs) were computed for accuracy, precision, recall, F1-score, and AUC using a bootstrapping strategy with 1,000 resampling iterations on the test set. The trimodal fusion model achieved an accuracy of 92% (95% CI: 89.4%–94.1%), indicating stable and statistically reliable performance. The narrow confidence bounds across all evaluated metrics confirm that the observed performance improvements are not due to random variation and demonstrate the robustness of the proposed framework. Cross-validation results demonstrated consistent performance across all folds, with accuracy variation limited to ±1.2%, further validating the robustness of the proposed multimodal fusion strategy.

The uniform manifold approximation and projection (UMAP) of fused embeddings illustrated in Figure 17 clearly shows that there are clusters between healthy controls and Parkinson’s disease (PD) subjects. Although there is some overlap because of natural inter-patient variability, the clear separation points to the fact that the trimodal feature fusion allows the learner to learn discriminative boundaries. The diffusion of points through the embedding space indicates the way the multimodal integration achieves the variety of pathological patterns, tremors in handwriting, instability in walking, and speech difficulties, which would lead to a more significant and more accurate representation of Parkinson’s disease.

The ROC curve shown in Figure 18 indicates a great distinction capacity of the trimodal model among healthy and Parkinson subjects. The model has a very high discriminative power with an AUC of 0.95 indicating that it is capable of correctly classifying positive (Parkinson’s disease) and negative (Healthy) cases at various thresholds. This is because as the curve remains close to the upper-left corner, the sensitivity (true positive rate) and specificity (low false positive rate) are high.

The precision–recall curve illustrated in Figure 19 also demonstrates the strength of the model especially when the classes are imbalanced. The model has an average precision (AP) of 0.96, which means that it has a high precision and a high recall. It implies that the classifier is not only accurate on the vast majority of Parkinson cases, but also false alarms are kept to a minimum, which is essential when the clinical results would treat instances of missing a true case as a problem, as well as excessive false positives.

The SHAP summary plot demonstrated in Figure 20 indicates that the handwriting characteristics (hw1, hw2), as well as a number of speech characteristics (e.g., speech15, speech47, and speech9), have the most significant impact on the prediction displayed by the model, with a group of gait characteristics playing the supportive role. The diagnostic reliability of SHAP values of handwriting is through a large value of SHAP, representing tremor-induced anomalies of drawing tasks, and speech characteristics representing the instability of the voice and articulation variation. The outcomes of the fusion are that the model is not over-dependent on one modality; instead, it adopts complementary features between handwriting, speech, and gait, enabling interpretability and clinical confidence.

To evaluate the generalizability of the proposed framework, external validation experiments were conducted across handwriting, speech, and gait modalities using datasets collected under different acquisition conditions and participant demographics from the training data. For handwriting analysis, the trained ResNet50 model was evaluated on an independent publicly available Parkinson's handwriting dataset, demonstrating comparable classification trends with a marginal performance drop attributable to domain shift. Speech models were validated using external recordings obtained from heterogeneous acoustic environments, confirming robustness to recording variability. Gait embeddings learned by the autoencoder were evaluated on an independent cohort-based dataset, where clustering and classification behavior remained consistent with internal validation results. Overall, external validation results indicate that the proposed multimodal framework generalizes effectively beyond the original dataset, supporting its potential clinical applicability.

To quantify the contribution of each modality and the robustness of the trimodal fusion model, we conducted ablation experiments and statistical tests. In the ablation study, one modality was removed at a time from the trimodal fusion framework, and performance metrics (accuracy, precision, recall, F1-score) were evaluated on the test set. Table 4 summarizes the results. The removal of handwriting or gait significantly reduced accuracy (from 92% to 85% and 86%, respectively), while removing speech resulted in a smaller reduction (92% → 89%), indicating that handwriting and gait provide the most discriminative information, with speech acting as a complementary modality. Additionally, we assessed the statistical significance of the performance improvements of the trimodal model over unimodal baselines using the Wilcoxon signed-rank test. The results show that the differences in accuracy between the trimodal and each unimodal model are statistically significant (p < 0.01). An optional experiment was performed without XAI-guided feature selection. The model trained without interpretability-driven feature selection exhibited slightly lower performance (accuracy: 91% vs. 92%) and less clinically coherent decision patterns, confirming that XAI-driven feature interpretation enhances both model reliability and clinical explainability. These analyses collectively demonstrate that:

Each modality contributes meaningfully to the overall performance, with handwriting and gait being the most informative.

The proposed multimodal PD detection framework primarily contributes an integrated and explainable design, rather than introducing novel deep learning architectures. Ethical considerations are critical in deploying AI models in clinical contexts. Our datasets include subjects from diverse demographics; however, inherent biases in age, gender, and acquisition conditions may affect model generalization. To mitigate this, we applied normalization and cross-validation strategies and reported performance across all classes. Additionally, the framework ensures privacy-preserving practices by using de-identified datasets and following applicable data-sharing guidelines. Future work will focus on further bias auditing, fairness assessment, and clinician-informed validation to guarantee equitable and safe deployment of the model in real-world settings.

While the proposed multimodal framework demonstrates improved robustness compared to unimodal approaches within retrospective benchmark datasets, its adaptivity is limited to feature-level fusion and has not been validated under real-world dynamic clinical conditions.

Despite the robustness introduced through multimodal fusion, certain modality-specific generalizability challenges remain. Speech-based features may be influenced by language, accent, and recording conditions, potentially limiting cross-lingual applicability. Handwriting datasets are typically collected under controlled experimental settings, which may not fully reflect natural daily writing behavior. Similarly, gait signals are dependent on sensor type, placement, and acquisition protocols, which can vary across clinical and real-world settings. Although the fusion framework mitigates some of these issues by leveraging complementary modalities, complete generalization across diverse populations and environments cannot be guaranteed. Future studies should incorporate cross-lingual speech data, free-form handwriting samples, and protocol-independent gait recordings, along with domain adaptation and normalization strategies.

The proposed system integrates multiple deep feature extractors, a fusion and classification stage, and explainable AI modules, resulting in a relatively complex computational pipeline. While this design enables high accuracy and strong interpretability, it may pose challenges for direct deployment in resource-constrained clinical settings or real-time screening scenarios. The current implementation prioritizes methodological validation rather than deployment efficiency. However, the architecture is modular and can be optimized through encoder pruning, model compression, and modality-adaptive inference, allowing the system to operate with a subset of available modalities when necessary. Future work will explore lightweight model variants and edge–cloud hybrid deployment strategies to improve practicality without compromising diagnostic reliability.