Plantar-pressure data were collected using a bilateral in-shoe wearable system, with 1 insole per foot. Each insole contained 8 discrete pressure sensors, yielding 16 channels sampled at 20 Hz. The detailed 16-channel sensor-to-anatomy mapping and region definitions are provided in Multimedia Appendix 2.

Figure 1 illustrates the leakage-controlled preprocessing and evaluation pipeline used in this study. The full workflow consisted of following 3 steps: participant-wise data splitting, construction of parallel signal representations, and sliding-window construction. All windows from the same participant were assigned to the same subset throughout, preventing repeated-participant leakage across model development and final testing and reducing optimistic bias in performance estimation [22,23]. We first performed a participant-level stratified split, yielding an independent test set of 20 participants (11 non-ROA and 9 ROA) and a development set of 72 participants. All subsequent model-development procedures were confined to the development set.

For each participant, a fixed-length steady-walking segment of 500 samples was retained from the original plantar-pressure recording. From this same segment, 2 signal representations were constructed for subsequent analyses—the raw plantar-pressure signal in its original scale, and a participant-wise channel-wise standardized representation. Let μc and σc denote the mean and SD of channel c computed from that participant’s steady-walking signal. We first standardize each channel within a participant:

where ϵ=10-6 ensured numerical stability. The same representation-construction procedure was applied to both the development and independent test sets. After representation construction, the retained steady-walking segment was segmented into overlapping windows to capture short-term dynamics. We set the window length to L=80 samples (4 s at 20 Hz) and the stride to s=40 samples (2 s), yielding 50% overlap. This duration was chosen to span approximately 3‐4 complete gait cycles at typical walking cadences (1.0‐1.2 Hz), balancing within-window stability for PSI estimation against temporal resolution. This procedure produced Nw=11 windows per participant, so that all participants contributed the same number of window-level samples. Each window was denoted as X(w)∈RL×C for w=1,…,Nw.

We extracted 2 complementary window-level feature sets from the windowed plantar-pressure signals (ie, conventional biomechanical features and PSI-based coordination features; Figure 2). In the PSI branch, we aimed to capture distributed coordination patterns across plantar regions that may not be adequately characterized by conventional biomechanical descriptors alone. The PSI-related analyses were based on the participant-wise standardized representation. Estimating PSI directly from the window X~(w) can be sensitive to sensor noise, gait-phase misalignment, and interparticipant amplitude variation, particularly in small cohorts. To improve robustness while preserving a lightweight formulation, we introduced a mixer-style representation transform inspired by the token-channel mixing strategy of multilayer perceptron (MLP)-mixer and its time-series adaptations [24,25].

Within each development split, this transform was pretrained using only non-ROA windows from the corresponding training subset in a self-supervised one-step-ahead prediction task, consistent with predictive representation learning for sequential data [26,27]. The resulting module alternates temporal and channel mixing to produce a stable intermediate representation before PSI construction. Specifically, the transformed signal was defined as:

where M1 and M2 implement temporal mixing and channel mixing, respectively, using lightweight MLPs shared across channels and time steps to control complexity. PSI matrices were then constructed for each window over a multiscale lag set T = {4, 8, 12}, chosen to span short-to-long coordination offsets under the current sampling rate. The sensitivity of this design choice was further examined in the design-variation analyses. Each PSI matrix was defined from lag-aligned interregional correlation, motivated by the lead-lag structure commonly observed in plantar-pressure signals [28,29]. For window w, let z_i^(w)(t) denote channel i at time t in Z^(w). For any τ ∈ T, we formed a lead-lag aligned signal pair:

and defined the PSI matrix C(w)(τ)∈RC×C as

where corr(⋅,⋅) is Pearson correlation over length L-τ.

From each PSI matrix C(w)(τ), we extracted a compact set of PSI-derived summary features to characterize global coupling structure at each lag scale. These features capture complementary aspects of PSI structure, including edge density, distributional summaries (mean, SD, maximum, minimum), lagged asymmetry, and dominant spectral organization. The full definitions, mathematical formulations, and naming rules for PSI are provided in Multimedia Appendix 3. Thus, PSI features from all 3 lag scales were included in the candidate representation space before downstream feature selection.

Meanwhile, the biomechanics branch computed conventional plantar-pressure features directly from the raw pressure window X(w), providing interpretable anchors such as overall statistics, center-of-pressure (CoP), load ratios, and asymmetry indices. The complete feature list (35 features), channel mappings, region definitions, and formulas are provided in the Multimedia Appendix 2. This branch outputs window-level biomechanics feature vector fBio(w)∈R35. For each window, the PSI and biomechanics features were concatenated to form the window-level feature vector f(w)=[fPSI(w);fBio(w)].

To obtain 1 risk score per participant, we aggregated window-level features into a participant-level representation using fixed summary statistics. For a scalar window feature sequence {uk(w)}w=1Nw, we computed 5 statistics (mean, SD, maximum, minimum, and median) and stacked them as

Applying this aggregation to PSI features fPSI(w)∈R21 and biomechanics features fBio(w)∈R35 yielded F∈R105 and G∈R175, respectively, and the final participant-level feature vector was x=[F;G]∈R280.

Given the limited cohort size relative to the 280-dimensional representation, the resulting feature space was prone to redundancy and overfitting. We therefore used Fisher score ranking as a simple filter-based selection step to retain the most discriminative participant-level features. For feature i, the score was defined as

Here, μi,c and σi,c denote the class-specific mean and standard deviation for class c∈{0,1}, μi denotes the pooled mean, and ϵ=10-6 is a small constant for numerical stability. The Fisher score quantifies between-class separability relative to within-class variance. The top-50 setting was fixed before locked independent test-set evaluation according to the development protocol. The top-k sensitivity analysis reported below was conducted to examine whether this feature count lay within a stable performance range, rather than to select the final feature count using the locked test set.

The downstream classifier was a random forest (RF), chosen for its robustness in small-sample settings and its ability to capture nonlinear feature interactions [30]. To avoid information leakage in threshold selection, the operating threshold was determined within the development set by participant-level 5-fold cross-validation. In each fold, all model-development steps, including Fisher score ranking, top-50 retention, and RF fitting, were performed using only the fold-specific training partition. Each validation participant contributed 1 aggregated participant-level feature vector and 1 corresponding RF probability as an out-of-fold estimate, ensuring that the probability used for threshold selection was generated without using that participant in representation learning, feature selection, or classifier fitting. The final threshold was selected from these pooled out-of-fold predictions to maximize the score.

After this internal development step, Fisher score ranking was reapplied once to the full development set to construct the final top-50 panel, and the RF classifier was refit on all 72 development participants using the locked feature panel and locked operating threshold before 1-time evaluation on the independent test set. Thus, cross-validation was used for internal model development rather than for reporting averaged final performance, whereas the independent test set was reserved for a single locked evaluation. Further implementation details are provided in Multimedia Appendix 4.

We additionally benchmarked the proposed fusion model against 3 representative end-to-end temporal baselines to examine how it compares with alternative modeling paradigms under the same fixed-length input, development or test split, and participant-level evaluation protocol. InceptionTime was included as a strong convolutional baseline for time-series classification [31], Medformer as a Transformer baseline whose multigranularity patching design is well-suited to medical time-series classification [32], and TodyNet as a graph-based baseline that directly models multivariate temporal dynamics through dynamic graph construction and temporal graph learning [33]. For fairness, each baseline was trained with its own tuned model-specific hyperparameter configuration, while the shared data handling and evaluation setup are summarized in Multimedia Appendix 4.

All evaluations were performed at the participant level, with 1 ROA risk score and 1 binary decision assigned to each participant. We used the F₁-score as the primary metric:

In addition to F1, we reported AUROC, sensitivity, specificity, and precision to provide a more comprehensive characterization of discrimination performance. Given the limited test-set size (n=20), point estimates alone could overstate differences between models. We therefore computed 2000-iteration bootstrap 95% CIs for all metrics, and supplemented the bootstrap analysis with paired DeLong tests to evaluate whether observed AUROC differences were statistically distinguishable from chance. To assess whether the primary model performance could be explained primarily by available demographic and anthropometric variables, we conducted a supplementary model-level covariate sensitivity analysis. First, we trained a demographics-only RF model using age, height, weight, and BMI under the same participant-wise evaluation protocol. Second, we added these variables to the Fusion feature set before Fisher-score selection and RF classification. The complete covariate sensitivity analysis protocol is detailed in Multimedia Appendix 6.

To better understand which participant-level descriptors contributed consistently under resampling, we examined feature-selection stability across the 5 development folds. Within each fold, Fisher score ranking was applied to the training partition using the same feature-selection procedure as in the main pipeline. For each candidate feature, we summarized the number of folds in which it was retained and its mean within-fold Fisher score. For interpretation, retained features were further grouped by descriptor family, lag scale, and anatomical sensor region according to the insole topology.

To characterize robustness, we evaluated the fusion, bio-only, and PSI-only models under the following 3 perturbation categories: temporal scaling, random channel dropout, and additive Gaussian noise. Temporal scaling was tested at 0.8×, 1.0×, and 1.2× relative to the native sampling rate. Random channel dropout was evaluated at masking rates (0.0, 0.1, 0.2, and 0.3), and additive Gaussian noise was evaluated at σ=0.00, 0.01, 0.03, and 0.05 applied to the standardized signal. Design sensitivity was further examined by varying the sparsity threshold in PSI matrix construction, the lag configuration (single-scale vs multiscale), and the window length using 3 s, 4 s, and 5 s windows with 50% overlap. Phase-resolved PSI patterns were visualized by aligning trajectories to the normalized stance phase using 5% bins and pooling them at the group level within each functional module and lag setting. These perturbation and design-sensitivity analyses were used to assess the stability of the selected configuration and were not used to reselect the feature count, window length, lag configuration, sparsity threshold, classifier, or operating threshold.

The study was approved by the Institutional Review Board (IRB) of Zhujiang Hospital, Southern Medical University, Guangzhou (IRB 2019-KY-016‐02). All participants were enrolled from the PROC study and provided written informed consent before data collection; the original consent permitted secondary analysis of deidentified data without additional consent. Participants received financial compensation of approximately US $30 per participant. No identifiable images of individual participants appear in this manuscript or its supplementary materials.

We first compared PSI-only, bio-only, and fusion models under the leakage-controlled development or test protocol. As summarized in Table 1, PSI-only showed limited standalone classification performance, with a low F₁-score, AUROC, accuracy, and precision despite relatively high recall. Bio-only provided a strong baseline with substantially improved discrimination and a balanced sensitivity-specificity profile. Fusion showed numerically higher F₁-scores (0.842 vs 0.778), accuracy (0.850 vs 0.800), precision (0.800 vs 0.778), and recall (0.889 vs 0.778) than bio-only, while AUROC and specificity were identical between the 2 models.

Bootstrap 95% CIs for all metrics are shown in Figure 3. The intervals were wide and overlapping across branches; for instance, the F₁-score spanned 0.500‐0.948 for bio-only and 0.615‐1.000 for fusion. Paired DeLong tests showed that both fusion and bio-only outperformed PSI-only at the AUROC level (P=.02 and P=.04, respectively), whereas fusion did not differ from bio-only in AUROC (P>.99). Together, these results suggest that PSI alone is insufficient as a standalone classifier, while fusion showed a directionally favorable threshold-dependent profile relative to bio-only. Given the wide and overlapping CIs, this pattern should be interpreted cautiously as preliminary evidence that PSI may complement biomechanical descriptors rather than as definitive superiority over bio-only. The complete bootstrap confidence intervals and paired DeLong comparisons are provided in Multimedia Appendix 7.

The covariate sensitivity analysis is summarized in Multimedia Appendix 6. The demographics-only RF model performed substantially worse than the fusion model, indicating that age, height, weight, and BMI alone could not reproduce the primary plantar-pressure–based classification performance. Adding these demographic and anthropometric variables to the fusion feature set did not materially improve the F₁-score, AUROC, or specificity. Within the augmented feature space, age and BMI were retained by Fisher score ranking, whereas height and weight were not retained. These findings suggest that the fusion model was not simply a demographic proxy, although residual confounding cannot be excluded.

To examine whether the locked top-50 feature-retention setting lay within a stable performance range, we conducted a top-k scan across 15 values; results are illustrated in Figure 4. Fusion performance at k=50 (F₁-score=0.842) lay within a high-performance plateau, with similarly strong F₁-score and accuracy at neighboring k values. At smaller k values, performance was limited by an insufficient number of retained features, whereas larger panels around k>=90‐100 showed degraded classification performance, with F₁-score dropping to 0.667 at k=100. These results suggest that k=50 was located within a relatively stable operating range rather than representing an isolated favorable point. However, this analysis was used to assess the robustness of the selected configuration and was not used for final feature-count tuning.

To evaluate whether the descriptors selected during the 5-fold development procedure were reproducible under sample variation, we analyzed feature-selection stability across the 5 development folds. In each fold, Fisher score ranking was applied to the fold-specific training partition using the same selection procedure as in the main pipeline, and the top 50 descriptors were retained. Each descriptor was then summarized by its cross-fold selection rate and mean within-fold Fisher score, and descriptors retained in at least 3 of 5 folds were regarded as stable. This cross-fold view was used for interpretation because, in the present small-sample setting, the single top-50 panel obtained from the full development set may be sensitive to sample composition, whereas repeatedly retained descriptors are more likely to reflect reproducible discriminative structure [34-36].

Figure 5 summarizes the cross-fold overlap structure of selected descriptors across the 5 development folds. A substantial shared overlap, together with several additional higher-overlap intersections, indicates that the retained descriptor set was not purely fold-specific but contained a reproducible core. Importantly, the higher-overlap intersections were dominated by biomechanical descriptors, supporting conventional biomechanics as the primary anchor of the stable feature set. PSI-derived descriptors appeared less frequently and in smaller numbers, but they were still recurrent across multiple higher-overlap intersections, indicating a complementary yet reproducible coordination-level contribution. Lower-overlap intersections were also present, suggesting some sensitivity of the selected panel to sample composition, although these intersections were generally smaller and more heterogeneous. Importantly, these stable descriptors were used for interpretation rather than as a fixed final feature panel. For final model construction, Fisher score ranking was reapplied once to the full development set, and the locked top-50 panel was selected from all candidate descriptors using all 72 development participants.

To further illustrate the main discriminative signals within this reproducible subset, Table 2 highlights representative descriptors that were both relatively stable across folds and highly ranked by Fisher score.

The representative examples make the reproducible structure more concrete and highlight two complementary, clinically meaningful patterns. First, the biomechanical examples are concentrated in interpretable domains, such as asymmetry, CoP, and gait timing, reinforcing conventional biomechanics as the primary anchor of the stable signal. This pattern is consistent with previous evidence that knee osteoarthritis is associated with abnormal plantar-pressure distribution, altered load redistribution, and CoP-related gait changes [6,16-18]. Second, PSI-derived descriptors were also retained among the stable high-ranking features, indicating that coordination-level information was recurrently represented alongside conventional biomechanical descriptors. Taken together, these representative descriptors indicate that ROA-related discrimination in the present cohort was anchored in clinically interpretable abnormalities of plantar loading, load redistribution, and whole-foot progression.

To further interpret the coordination-level information captured by PSI, we examined whether ROA-related PSI differences showed coherent organization across lag scales, stance phase, and functional plantar modules. Because PSI was designed to characterize interregional temporal coupling rather than local pressure magnitude, the analysis was structured as a progressive visualization from sensor-pair spatial organization to phase-resolved stance dynamics, and finally to module-level difference trajectories. This design allowed us to assess whether PSI captured a distributed coordination pattern associated with ROA rather than isolated or nonspecific correlation differences.

We first visualized sensor-pair delta-correlation matrices, defined as ROA minus non-ROA, at representative lag settings. As illustrated in Figure 6, the delta-correlation matrices revealed spatially distributed group differences across multiple sensor pairs, rather than a focal abnormality confined to a single plantar region. The pattern also varied across lag settings, indicating that ROA-related PSI alterations depended on the temporal offset used to characterize interregional coupling. At shorter lags, the observed differences may reflect altered near-synchronous load sharing among plantar regions, whereas longer lags may capture delayed coordination related to forward load transfer and bilateral stance progression. These findings suggest that PSI differences were not merely a consequence of local pressure magnitude changes, but reflected lag-dependent reorganization of plantar coordination.

The spatially distributed and lag-dependent delta-correlation patterns in Figure 6 suggested that ROA-related PSI differences were linked to the temporal organization of plantar coordination during walking. However, the heatmap representation does not indicate when during stance these coordination differences emerged. We therefore further aligned PSI trajectories to the normalized stance phase (0%‐100%) to examine the temporal expression of the observed lag-dependent coordination changes. As illustrated in Figure 7, group-level PSI trajectories differed between ROA and non-ROA across all examined lag settings. Across broad portions of stance, the ROA curves remained lower than the non-ROA curves, suggesting reduced or less coherent interregional plantar coordination in the ROA group. From a gait-phase perspective, this reduction was not limited to initial contact or early loading, but extended across mid-stance and late stance, where controlled CoP progression, forward load transfer, and preparation for push-off become increasingly important. The magnitude and temporal expression of this separation varied with lag scale; shorter lags emphasized more immediate interregional coupling, whereas longer lags showed greater fluctuation near late stance and push-off. These lag-dependent patterns suggest that the proposed PSI may capture ROA-related compensatory changes in plantar load redistribution, particularly during forward load transfer and propulsion, when coordinated forefoot and heel-to-midfoot interactions become biomechanically important.

Figure 8 further shows phase-resolved PSI difference trajectories at a representative lag (τ=0.4s), where the vertical axis represents ΔPSI=PSIROA-PSInon-ROA. Across most of stance, the major module-level trajectories remained below zero, indicating that PSI was generally lower in ROA than in non-ROA for these interregion couplings. The persistence of negative differences across broad stance intervals suggests that ROA-related coordination changes were not restricted to a single gait event, but were expressed throughout stance-phase load transfer. This group difference persisted across broad portions of stance and was especially evident in forefoot-related and heel-to-mid couplings, while larger fluctuations near late stance and push-off suggested greater between-participant heterogeneity when forward load transfer and propulsion became dominant [6,16,17].

Taken together, Figures 6-8 indicate a coherent progression of PSI evidence: ROA-related differences were spatially distributed across sensor pairs, temporally expressed across the normalized stance phase, and functionally organized across interregion plantar modules. Importantly, these phase-resolved analyses were intended to characterize the organization of PSI differences, rather than to imply that all stance bins, lag scales, or modules contributed equally to the final classifier. The final retained descriptors were still determined by the participant-level feature-selection procedure.

We further examined whether PSI-related information remained useful under perturbation and design variation. Across temporal scaling, channel dropout, and Gaussian noise, fusion generally maintained the highest or among the highest participant-level F₁-scores among the 3 models (Table 3). Performance was strongest at the native and mildly stretched temporal scales (1.0× and 1.2×, both F₁-score=0.842), with only a modest reduction under compression (F₁-score=0.824 at 0.8×). Under channel dropout and Gaussian noise, fusion also maintained a generally favorable profile relative to both PSI-only and bio-only models. These perturbation analyses suggest preliminary robustness of the combined representation under moderate signal variation, although the small locked test set limits strong comparative claims.

As shown in Table 4, across design variations, fusion showed F₁-scores values within a relatively narrow range over sparsity thresholds (F₁-score=0.824‐0.842) and lag settings (F₁-score=0.818‐0.842), while showing greater sensitivity to window length (F₁-score=0.761‐0.842). The 4-second window (80 samples at 20 Hz, spanning approximately 3‐4 complete gait cycles at typical walking cadences of 1.0‐1.2 Hz) showed the numerically highest F₁-score among the examined window lengths. Shorter 3-second windows, covering only 1‐2 gait cycles, may have yielded less reliable correlation estimates for PSI matrix construction because the reduced sample size per window limited the statistical stability of the lag-aligned Pearson correlations. Conversely, 5-second windows, although providing more samples per correlation estimate, may have smoothed phase-specific coordination changes across stance subphases, reducing the ability to resolve temporally localized PSI differences between ROA and non-ROA groups. The 4-second window therefore appeared to provide a suitable balance between within-window statistical stability for correlation-based PSI estimation and temporal resolution across stance subphases.

Within the lag analysis, the multiscale configuration (T={4, 8, 12}, corresponding to 0.20, 0.40, and 0.60 s at 20 Hz) showed the numerically highest performance among the examined lag settings (F₁-score=0.842), compared with each single-scale alternative. A single lag may not simultaneously capture rapid intrafoot pressure propagation between adjacent sensors and slower cross-region coordination reorganization spanning multiple functional modules. The multiscale design therefore provided a broader candidate coordination space—spanning short-to-intermediate coupling timescales relevant to mid-to-late stance gait dynamics—from which downstream Fisher score selection could retain informative lag-specific descriptors across all 3 scales. This exploratory pattern is consistent with the idea that the multiscale design may capture complementary coupling timescales, but it requires confirmation in larger independent cohorts. Overall, these exploratory perturbation and design-sensitivity analyses suggest that the selected configuration was not uniquely dependent on a single perturbation or design setting. These findings should be interpreted as preliminary robustness checks rather than as independent optimization or validation of the model design.

Table 5 shows that fusion achieved the strongest overall participant-level performance among the evaluated models. Among the baselines, InceptionTime [31] was the strongest competitor and showed a relatively balanced classification profile, but remained below fusion in accuracy, F₁-score, AUROC, precision, and recall, while matching the best specificity. Medformer [32] achieved high recall comparable with fusion, but with lower AUROC, lower precision, and markedly reduced specificity, indicating a more false-positive−prone decision profile. TodyNet [33] showed perfect recall but near-zero specificity, yielding the weakest overall discrimination and suggesting a strong bias toward positive predictions under the present data scale.

Overall, these results suggest that, under the current data scale and evaluation setting, fusion showed a directionally favorable threshold-dependent profile relative to the evaluated end-to-end baselines. Relative to InceptionTime, the strongest baseline, fusion showed numerically higher accuracy, F₁-score, AUROC, precision, and recall while matching specificity, but these differences should be interpreted cautiously given the small locked test set. Medformer preserved high recall but at the cost of lower specificity and lower overall discrimination, whereas TodyNet showed a marked false-positive tendency with poor ranking quality. Together, these comparisons suggest that the proposed representation-based framework remained competitive with representative end-to-end baselines in this proof-of-concept evaluation.

This study addressed a key representation gap in plantar-pressure–based ROA identification by introducing a coordination-aware framework that integrates local biomechanical loading with multiscale interregional temporal information. The combined framework achieved competitive participant-level ROA identification, but it did not show consistent superiority over bio-only across performance metrics. The biomechanics branch made the strongest individual contribution, whereas PSI showed limited standalone discrimination. However, PSI captured complementary coordination patterns that were not represented by conventional biomechanical descriptors, and a reproducible PSI component was retained within the combined representation. Robustness and design-sensitivity analyses provided preliminary evidence that the findings were not confined to a single evaluated configuration. Accordingly, the principal innovation of this study is representational and methodological rather than cohort- or device-based; it explicitly combines interpretable local loading descriptors with multiscale interregional temporal coordination under a participant-level, leakage-controlled evaluation protocol.

The branch-wise findings help explain the primary source of discrimination. The strong contribution of the biomechanics branch is consistent with previous evidence that ROA-related plantar-pressure abnormalities occur in anatomically meaningful regions, particularly the hallux and heel, and are also reflected in CoP stability and variability [17,18]. This interpretation is further supported by recent studies associating ROA-related functional impairment with altered plantar-pressure behavior, gastrocnemius muscle-tendon-unit stiffness, and 3D gait biomechanics [7,8]. These patterns are biomechanically plausible because structural degeneration, pain-related gait adaptation, and altered neuromuscular control at the knee may be associated with compensatory load redistribution and changes in stance progression [6-8,15-18]. Experimental modulation of tibialis anterior activity has also been shown to alter plantar pressure and gastrocnemius activity during gait, supporting the sensitivity of plantar loading to neuromuscular modulation [37]. The findings of this study therefore reinforce the importance of anatomically grounded loading descriptors for participant-level ROA identification.

The phase-resolved and fold-wise analyses help clarify what PSI may contribute beyond these local loading descriptors. The observed PSI differences were distributed across multiple functional modules and lag settings rather than being confined to isolated local extrema. In addition, a core group of PSI descriptors was repeatedly retained across the development folds, supporting the reproducibility of these selected patterns within the present dataset in accordance with established feature-stability principles [34,35]. Although PSI is not a direct measure of joint kinematics or electromyographic muscle activity, its distributed temporal patterns are compatible with studies reporting altered kinematic synergy, muscle synergy, and muscle-activation patterns in individuals with KOA during walking [38-40]. PSI may therefore provide an interpretable coordination-level description of bilateral load transfer and whole-foot progression. However, feature stability and phase-resolved interpretability do not by themselves establish incremental predictive value beyond the biomechanics branch. Given the limited standalone discrimination of PSI and the absence of consistent superiority of fusion over bio-only, the findings of this study should not be interpreted as definitive evidence of additional predictive or clinical value.

Comparisons with the end-to-end baselines place these findings in a broader modeling context. The evaluated models represented complementary learning paradigms but exhibited different operating profiles: InceptionTime was the strongest and most balanced competitor, Medformer favored sensitivity over specificity, and TodyNet showed lower specificity and a higher false-positive tendency in the present locked test set. The value of the proposed framework therefore lies less in definitive performance superiority over every alternative than in combining competitive participant-level identification with an explicit and interpretable organization of plantar-pressure information before classification. The robustness and design-sensitivity analyses further suggested that the main findings were not entirely dependent on a single analytical configuration, although these analyses remain internal and exploratory. Together with the structured comparison in Multimedia Appendix 5, these results support the positioning of the framework as a representational and methodological contribution rather than a cohort- or device-based advance.

To provide a screening-context interpretation of the locked-test operating characteristics, we translated sensitivity and specificity into positive predictive value and negative predictive value across assumed community ROA prevalence levels in Multimedia Appendix 8 [41]. This supplementary analysis illustrates how disease prevalence could affect the potential rule-in and rule-out value of the framework. It should not, however, be interpreted as evidence of clinical deployment readiness because the underlying operating characteristics were estimated from a small internal test set rather than from a representative prospective screening population.

Several limitations should be acknowledged. The sample size was modest, and all participants were drawn from a single previously reported cohort, limiting the precision and external generalizability of the performance estimates. Although the locked test set was separated from model development, it remained part of the same source cohort and therefore represents an internal held-out proof-of-concept evaluation rather than external validation. Contemporary prediction-model guidance emphasizes that independent validation should be adequately powered to estimate performance with sufficient precision [14,42]. The findings of this study should consequently be interpreted as preliminary and require confirmation in larger independent cohorts.

Important clinical and demographic factors, including age, BMI, and affected-side heterogeneity, may influence plantar-pressure patterns [15]. The supplementary covariate sensitivity analysis indicated that the available demographic variables alone did not reproduce the performance of the fusion model and that their inclusion did not materially improve its performance (Multimedia Appendix 6). Nevertheless, this analysis cannot exclude residual confounding or establish the independent and interactive effects of age and BMI in the present modest sample. These questions will require larger, more carefully stratified cohorts with richer clinical covariate information. Furthermore, to our knowledge, no publicly available dataset combined confirmed ROA labels with suitable in-shoe plantar-pressure recordings for direct external validation. Available public resources lacked either ROA annotations or compatible pressure measurements, precluding cross-dataset benchmarking [43,44]. Multicenter efforts to establish open datasets with imaging-grade clinical labels and standardized plantar-pressure recordings would substantially strengthen future evaluation.

This study also focused on binary ROA identification rather than disease-severity stratification or longitudinal progression. PSI depends on sensor layout, signal quality, and preprocessing choices, which may affect the validity and comparability of plantar-pressure measurements [10]; its transferability across devices and clinical settings therefore remains uncertain. Recent Soft Science studies demonstrate advances in biocompatible piezoelectric, gel-based, and electrospun flexible pressure sensors for walking, human-motion, gait, and continuous health monitoring [45-47]. These developments strengthen the technological basis for wearable plantar sensing, but they also underscore the need to verify signal stability, durability, environmental resilience, and cross-platform calibration before clinical translation [45-47]. In addition, the proposed coordination-level interpretation was not directly validated using concurrent joint-kinematic, electromyographic, or other neuromuscular measurements. Future studies should validate the framework prospectively in larger multicenter cohorts, evaluate its robustness across sensor platforms and patient subgroups, examine its association with direct biomechanical and neuromuscular measurements, and extend it toward severity modeling and longitudinal prediction.

This study introduces an interpretable participant-level representation that integrates local biomechanical loading with multiscale interregional temporal coordination, extending plantar-pressure approaches centered primarily on summary descriptors or less interpretable end-to-end models. Complementary evaluations provided convergent evidence of competitive participant-level ROA identification and preliminary robustness across the evaluated settings, while showing that local biomechanical descriptors constituted the main discriminative basis and that PSI offered a reproducible coordination-level component whose incremental predictive value remains to be established. By broadening plantar-pressure analysis from local loading to distributed gait coordination, the framework provides a methodological basis for more interpretable functional characterization of ROA. If prospectively validated in external multicenter cohorts, it could warrant further evaluation as an adjunct to imaging for functional assessment and longitudinal monitoring. Prospective multicenter validation is required before clinical translation.