An analysis of 3D knee kinematic data complexity in knee osteoarthritis and asymptomatic controls

Three-dimensional (3D) knee kinematic data, measuring flexion/extension, abduction/adduction, and internal/external rotation angle variations during locomotion, provide essential information to diagnose, classify, and treat musculoskeletal knee pathologies. However, and so across genders, the curse of dimensionality, intra-class high variability, and inter-class proximity make this data usually difficult to interpret, particularly in tasks such as knee pathology classification. The purpose of this study is to use data complexity analysis to get some insight into this difficulty. Using 3D knee kinematic measurements recorded from osteoarthritis and asymptomatic subjects, we evaluated both single feature complexity, where each feature is taken individually, and global feature complexity, where features are considered simultaneously. These evaluations afford a characterization of data complexity independent of the used classifier and, therefore, provide information as to the level of classification performance one can expect. Comparative results, using reference databases, reveal that knee kinematic data are highly complex, and thus foretell the difficulty of knee pathology classification.


Introduction
Three-dimensional (3D) knee kinematic data, measuring knee flexion/extension, abduction/ adduction, and internal/external rotation, are increasingly used in gait analysis towards quantifying the knee function [1], understanding pathological knee alterations [2], and assessing the progression of knee pathologies and their impact on the gait [3,5]. They also offer opportunities for diagnostics [6], classification [7], and therapy of knee musculoskeletal pathologies [8]. Kinematic data are generally given in the form of normalized curves of measurements recorded on a treadmill at regular intervals of time during a gait cycle. Several studies have used knee kinematic data to obtain crucial information about musculoskeletal pathologies [9], to distinguish between knee osteoarthritis and asymptomatic subjects [10], and to further and transverse planes during a full gait cycle (Fig 1), serve knee pathology classification algorithms [22] and can subsequently guide OA therapy plans [8].
The data used in this study was collected, for research projects [3,23,24], using a state-ofthe-art KneeKG acquisition system. The accuracy and reproducibility of the system have been assessed in several studies [25][26][27]. For each subject, the positional angles are recorded during two walking trials of 25 sec on a treadmill. For each participant patient, current methods of analysis have used the average of about a dozen curves in the set of measurement curves, so chosen as to maximize a correlation criterion between classes, namely the Inter Class Correlation Index of the observations [28]. In our study, a series of fifteen (15) full gait cycle walks on the treadmill are averaged to obtain a representative mean pattern per subject based on an adjusted coefficient of multiple determinations [29]. This is followed by interpolation and resampling from 1% to 100% of the gait cycle, therefore giving 100 measurement points for each participant. This yields a significant amount of high-dimensional data, leading to "the curse of dimensionality," in classification [30]. Also, knee pathologies notoriously cause high withinclass variability and in-between class similarity (example in Fig 1).
The data collection was approved by institutional ethics committees of the University of Montreal Hospital Research Center (Reference numbers CE 10.001-BSP and BD 07.001-BSP) and of the École de technologie supérieure (Reference numbers H20100301 and H20170901). All subjects provided written informed consent before the studies began. We used 3D knee kinematic data from two groups. The first group included 40 patients with knee OA diagnosed according to clinical and radiographic criteria defined by Altman et al. [31]. The patients were excluded if they had vestibular, neurological, or musculoskeletal disorders, fracture of the lower extremity, rheumatoid arthritis, or generalized osteoarthritis, limping gait or any condition that could affect a treadmill walking evaluation. The second group, which formed the control group, contains 40 asymptomatic participants (AS). Table 1 summarizes the participants' demographic characteristics and their walking speed. A statistical t-test was performed to compare the two groups characteristics using SPSS 18.0 (Statistical Package for Social Sciences) (SPSS Inc. Released, 2009. PASW Statistics for Windows). A threshold of significance set at p = 0.05 shows that there is no statistical difference between the two groups except for the height variable.

Features Extraction on 3D Knee kinematic data
A set of 62 biomechanical features are extracted from these 3D kinematic patterns for data complexity analysis. The features chosen are biomechanical parameters routinely assessed in clinical biomechanical studies of knee OA populations; for instance, maximum, minimum, varus and valgus thrust, angles at initial contact, mean values and range of motion (ROM) throughout gait cycles or specific gait sub-cycles (i.e., loading, stance, swing, etc.) [2,3,32]. The 3D knee kinematic data were normalized to remove dependence on magnitude and scale.

Knee kinematic data complexity
Data complexity is evaluated on the descriptive biomechanical features according to both single features and global features complexity metrics. The first metric evaluates the overlap of individual feature measurements whereas the latter evaluates the separability of classes and the geometry of the manifolds spanned by each class when the features are considered simultaneously.
Single feature complexity. Complexity of a single feature informs us on that particular feature discriminant power, i.e., its capacity to distinguish data samples from distinct classes. Such a characterization is useful for clinical assessments. Indeed, it allows the study of specific kinematic points of interest, having a clinical meaning, in order to correlate them with typical knee function behavior.
For each individual feature, the following three metrics examine the range and spread of the values of instances from distinct classes. Specifically, let D ¼ fx 1 ; x 2 ; :::::; x n g be a data set containing n knee kinematic data measurements, each belonging to one of two distinct classes c 1 and c 2 . Each element x i is characterized by a feature vector (f 1 , f 2 , . . .f p ), where p is the dimension of the feature space. In our case p = 62. The complexity metrics are: • Fisher discriminant ratio (F1): for a two-class data set, ratios F1 = {F1 1 , F1 2 , . . ., F1 p } for feature f i , (i = 1, . . ., p), are defined as: where μ i1 , μ i2 , σ i1 , σ i2 are the means and variances of the two classes, respectively, according to the i th feature, (i = 1, 2, . . ., p). The range of F1 i is [0, +1[. Small values indicate strong class overlap. High values indicate good class separation when using the i th feature.
• Volume of overlap region (F2): This measure estimates the amount of relative overlap of the bounding regions of two classes [33].
) the maximum and minimum values of each feature f i in the classes c 1 and c 2 . For two-class data sets, the ratio F2 i for each feature f i , (i = 1, 2, . . ., p) is computed as: The range of F2 is [0, 1]. Larger values indicate larger class overlap, and higher data complexity thereof. A shortcoming of this measure is that a single feature with non-overlapping regions causes F2 to be zero.
• Feature efficiency (F3): this measure is particularly relevant when dealing with high-dimensional data. It informs on how much each feature contributes to the separation of the classes, and the contribution is called efficiency. For each feature, the ambiguous (overlapping) regions are removed so that only non-overlapping regions remain for consideration. Proportion of women 72% (24) 16% (10) each feature f i , let Q i be the following ratio: where | | denotes the number of non-overlapping elements and n is the total number of elements in both classes. When there is no overlap between classes, F3 i achieves the maximum value of 1.
Each single feature complexity metric requires the estimation of high and low thresholds to assess whether the feature is complex or not. We estimated the various thresholds based on the probability distributions of the feature complexity measures. For each metric, i.e., the Fisher discriminant ratio F1 i , the volume of overlap region F2 i , and the feature efficiency F3 i , the best-fitting distribution is determined by the maximum likelihood estimation. The thresholds are subsequently estimated using the rule of thumb of 5th and 95th quantiles of the probability distribution.
Global feature complexity. The global feature complexity is based on the following metrics: fraction of points on the class boundary (N1), ratio of average intra/inter class nearest neighbor distance (N2), error rate of a nearest neighbor distance (N3), fraction of maximum covering spheres (T1), average number of points per dimension (T2).
• Fraction of points on the class boundary (N1): This measure estimates the length of the class boundary. The minimum spanning tree (MST) connecting all the training samples is generated and the fraction of points of a class connected to the opposite class is computed [34,35].
The range of N1 is [0, 1]. Large values indicate that the majority of points lay closely to the class boundary, which can make it difficult to locate the class boundary accurately.
• The ratio of average intra/inter class nearest neighbor distance (N2): Let x i denote an input sample, i = 1, 2, . . ., n, where n is the number of samples. Let IntraDist(x i ) designate the distance between x i and its nearest neighbor within the class, and InterDist(x i ) the distance to its nearest neighbor of the other class. The ratio of average intra/inter class nearest neighbor distance (N2) is the ratio of the sum of the intra-class distances to the sum of the inter-class distances for each input example: • Fraction of maximum covering spheres (T1): This measure describes the shapes of class manifolds using the notion of adherence subset [36]. An adherence subset is a sphere centered on an sample of the data set which is grown as much as possible before touching any sample of another class. The T1 is measured based on the number of spheres normalized by the total number of points.
• Average number of points per dimension (T2): The average number of points per dimension is the ratio of the number of elements in the data set to the number of features. Its value is a rough indicator of sparseness of the data set.
Statistical analysis. Following individual feature complexity evaluation, we performed a t-test statistical analysis to assess feature dissimilarity between the two groups (osteoarthritis patients and control participants). Since each feature single complexity variable informs on that particular feature discriminant power, i.e., its capacity to distinguish data samples from distinct classes, our hypothesis is that the features retained by complexity analysis as discriminant are statistically different in the two groups. A P-value of 0.05 was set as the statistical level of significance.

Results and discussion
We recall that data complexity assessment used 40 osteoarthritis patients (OA) and 40 control participants (AS). A set of 62 biomechanical features (p = 62) routinely assessed in clinical biomechanical studies of knee OA populations have been selected (as described in section Features Extraction on 3D Knee Kinematic Data). The analysis was done using Matlab R2017a software (Mathworks, Massachusetts, United State).

Single feature complexity
The single feature complexity measures have been evaluated using Fisher discriminant ratio (F1), the volume of overlap region (F2), and the individual feature efficiency (F3). The threshold values have been identified, for each complexity metric, using a probability density which best fit the complexity data. Fig 2 indicates that the Fisher discriminant ratio F1 and the feature efficiency F2 best follow a Weibull distribution, and the volume of overlap region F3 follows a normal distribution. Table 2 summarizes the thresholds value determined using the 10th and 90th percentiles of each complexity measure empirical density.
The single complexity measures of each feature and the corresponding thresholds are presented in Fig 3. Fig 3A represents the discriminant ratio F1 i for each feature f i (i = 1, 2, . . ., 62) computed according to Eq 1. Its threshold was set to 0.33 because we are interested in high values of F1 i ( Table 2, line 1). Except for the set of features {f 1 , f 2 , f 3 , f 6 , f 7 , f 11 , f 57 }, the low values of F1 i demonstrate the low discriminative power of the 3D kinematic features.
In Fig 3B, we represented the ratio of the width of the overlap interval F2 i of each feature. The threshold here was set to 0.47 because, for this metric, we are interested in low value of F2 i ( Table 2,

Statistical analysis
Following the single feature complexity analysis, we performed a statistical t−test analysis to examine the differences between the two groups (the osteoarthritis patients (OA) and the control participants (AS)) of these 13 retained features. Table 3 summarizes their clinical meaning and their corresponding p-value. The results of the statistical analysis confirms those obtained using the complexity analysis. Indeed, the retained 13 features are statistically significantly different, and can therefore be considered discriminant of OA patients and AS participants (p<0.001).
Both statistical analysis and single complexity analysis support previous studies on biomechanical data and their association with knee pathologies. Indeed, the features f 1 , f 2 , f 3 and f 11 which are related to the varus or valgus thrust during loading phase, have been identified as the most useful parameter, serving as a diagnostic biomarker [3,21]. Also the feature f 57 , which corresponds to the range of motion during loading phase, has been identified as a part of burden of disease biomarkers to discriminate between moderates OA grades and severe OA  An analysis of 3D knee kinematic data complexity grades [3]. Another study compared a set of biomechanical parameters in patients categorized as moderate to severe OA grades [21,37]. It was reported that both peak knee adduction moment and knee adduction angular impulse increased with knee radiographic grade.

Global feature complexity
The global feature complexity metrics have been measured using: the Fraction of points on the class boundary (N1), the Ratio of average intra/inter class nearest neighbor distance (N2), the An analysis of 3D knee kinematic data complexity Error rate of a nearest neighbor distance (N3), the Fraction of maximum covering spheres (T1), the Average number of points per dimension (T2). Table 4 summarizes these metrics. It shows high values of the ratio of average intra/inter class nearest neighbor distance (N2 = 0.95); this indicates that samples from the same class are dispersed. This result is confirmed by the fraction of points on the class boundary (N1 = 0.54) which shows that an important number of samples are located near the class boundaries. As a result the nearest neighbor error rate (N3 = 0.42) approaches that of random classification. Moreover, the analysis of the global complexity metrics shows the sparsity of the features. Indeed, the fraction of maximum covering spheres is equal to 1 (T1 = 1), which means that the largest adherence subset is one sample. Moreover, the average number of points per dimension (T2 = 1.29) indicates high sparseness of the data set. This sparsity is an important problematic in classification [30]. Table 1, which describes the general subject characteristics, shows that the proportion of women with respect to that of men is higher in the OA group (72% of women against 38% of men) while it is lower in the AS group (16% of women against 84% of men). For further analysis, to take into account gender, we measured the data global complexity for men and women separately to determine if it is gender-dependent. Table 5 summarizes the complexity results per gender. Overall, global complexity is high in both disease groups when analyzed separately per gender, thus confirming the conclusions drawn from Table 4. Moreover, the complexity is higher for women than men. This explains recent literature [38,39] reporting more effective automatic classification in men's groups than in women's.

Comparison to the complexity of other types of data
We compared the knee kinematic data complexity measures to those obtained using the data in two public databases for classification tests: the Wind database, which contains 178 samples from two classes described by 13 features, and the IRIS database, which contains 150 samples from two classes described by 4 features ( Table 6).
The ratio of average intra/inter class Nearest neighbor distance (N2) shows that the discriminant power is higher for Wine and IRIS than for the knee kinematic data. Indeed, N2 is 0.9 for the kinematic data, 0.026 for Wine, and 0.1 for IRIS. Also, the average number of points per dimension (T2) is very low for knee kinematic data. This indicates that 3D knee kinematic data is significantly more complex compared to that of IRIS and Wind data.

Conclusion
This study investigated 3D knee kinematic data complexity using both single feature and global feature complexity metrics. The first evaluates the overlap of individual feature measurements, whereas the latter measures the separability of classes and the geometry of the manifolds spanned by each class when the features are considered simultaneously. The results reveal that knee kinematic data are sparse and highly complex and, thus, foretell the difficulty of knee pathology classification. Indeed, from the initial set of 62 biomechanical features, only 13 features have been identified as not complex. Such information might be useful in future work (1) in the selection of an appropriate classification algorithm and for the interpretation of classification rates and, (2) in clinical studies to identify biomarkers for knee pathologies.

Author Contributions
Data curation: Jacques A. de Guise.