Pelvic floor pressure distribution profile in urinary incontinence: a classification study with feature selection

Data acquisition

A fully instrumented non-deformable probe with capacitive transducers was used to acquire the vaginal spatiotemporal pressure distribution. The probe consisted of an Ertacetal cylinder (tensile modulus of elasticity = 2,800 MPa) covered by a 10 × 10 matrix of individually calibrated capacitive sensors (MLA-P1, Pliance System; novel; Munich, Germany). The cylinder was 23.2 mm in diameter and eight cm in length, and its sensing area was 70.7 × 70.7 mm (10 × 10 sensing elements of 7.07 × 7.07 mm, with 1.79 mm gaps between them) (Fig. 1). The capacitive sensors had a measurement range of 0.50–100.00 kPa and a resolution of 0.42 kPa. Reliability and testing capacities of the instrument were described by (Cacciari et al., 2017a; Cacciari et al., 2017b).

We divided the sensor matrix into various sections. The first was five planes with 20 sensors each. Each plane (36° apart from each other) represented a sum vector of pressures from two 10-sensor lines diametrically opposed along the cylinder. The second was 10 rings, with each one created by the 10-sensor perimeters surrounding the cylinder. The third division had three lines, each with a 10-sensor perimeter surrounding the cylinder: cranial (corresponding to the first three lines of sensors from the vaginal opening), medial (four mid-lines of sensors), and caudal (three last lines of sensors). The final division was left, posterior, right, and anterior sections. Figure 2 illustrates the instrument with the sensor matrix disposition and a scheme of the sensor subsets.

Data acquisition was performed on 24 adult continent women 35.3 ±  10.0 years 23.4 ±  4.2 kg/m²) and 24 diagnosed with stress UI 48.2 ±  8.1 years 27.5 ±  3.6 kg/m²; incontinence impact 66.7/100; severity measures 41.7/100 on the King’s Health Questionnaire). To be eligible, women had to be continent or urinary incontinent, not virgins, with no history of pregnancy within the past year, in premenopausal status with monthly menstrual cycles, with a body mass index (BMI) lower than 30 kg/m², and with no history of pelvic floor muscle (PFM) training or of any medical conditions that could interfere with PFM function. During clinical evaluation, participants were excluded if they presented above Stage II on the Pelvic Organ Prolapse Quantification (POPQ) scale (Haylen et al., 2016) and if they were not able to voluntarily contract their PFM. The King’s Health Questionnaire was used to classify the participants into the continent or incontinent group (Tamanini et al., 2003). Women reporting no symptoms of any type of UI were included in the control group. This study was approved by the Ethics Committee of the School of Medicine of the University of São Paulo (protocol n.023/14), and all participants provided written informed consent prior to participation.

The probe was always inserted with the same orientation and at a depth of seven cm from the hymenal caruncle according to references marks. After a 1-minute accommodation period, the participants were asked to accomplish four maneuvers (in the same order) with a 1-minute rest between them: (1) the maximum contraction maneuver, in which the participants had to lift and squeeze their pelvic floor as hard as possible for 3 s (Cacciari et al., 2017a); (2) the Valsalva maneuver, which consists of executing maximum intra-abdominal pressure effort leading to a downward movement of the pelvic floor for 5 s (Cacciari et al., 2017a); (3) the endurance maneuver, in which participants had to sustain pelvic floor contraction for 10 s while breathing normally (Cacciari et al., 2017b); and (4) the wave maneuver, in which women were instructed to contract their PFMs in a caudal-cranial direction for 2 s and then relax them in a cranial-caudal direction for 2 s (Cacciari et al., 2017b). All participants received standardized verbal support to encourage them to perform maximal PFM contractions throughout the maximal and endurance maneuvers. The start of the data acquisition was manually synchronized a few seconds before the verbal command. The sampling frequency during pressure data acquisition was set at 50 Hz (Cacciari et al., 2017a).

Feature extraction

Feature extraction was based in previous studies (Cacciari et al., 2017a; Cacciari et al., 2017b) where the extracted features were capable of distinguishing vaginal sub-regions, planes, rings and maneuvers, important aspects in pelvic floor assessment, and presented excellent inter- and intra-rater reliability and intra-trial repeatability.

Prior to any feature extraction, each sensor time series was filtered by a zero-lag, 8th order, low-pass, Butterworth filter with a cut-off frequency of 8 Hz. Then, the variable extraction process was performed over sensor sets, which were defined either by a maximum or by a sum operator, leading to two different time series: a peak pressure (Eq. 2.1) and a sum pressure time series (Eq. 2.2):

(2.1)${\displaystyle Y_{peak}^S\left[t\right]=\underset{s^{{}^{\prime }}}{max}S\left[t\right]}$ (2.2)${\displaystyle Y_{sum}^S\left[t\right]=\sum _{s^{{}^{\prime }}\in S\left[t\right]}S\left[t\right]}$ where S[t] is the pressure reading at time instant t of a sensor set S, and s′  is an element of this set.

For the maximum contraction and Valsalva maneuvers, the extracted variables were the maximum pressure (Eq. 2.3), maximum sum (Eq. 2.4), instant of maximum pressure (Eq. 2.5), instant of maximum sum (Eq. 2.6), and instant of activation (Eq. 2.7):

(2.3)${\displaystyle ma{x}_p=max{Y}_{peak}^S\left[t\right]}$ (2.4)${\displaystyle ma{x}_s=max{Y}_{sum}^S\left[t\right]}$ (2.5)${\displaystyle t_{ma{x}_p}=\underset{t}{\text{argmax}}{Y}_{peak}^S\left[t\right]}$ (2.6)${\displaystyle t_{ma{x}_s}=\underset{t}{\text{argmax}}{Y}_{sum}^S\left[t\right]}$ (2.7)${\displaystyle t_{activ}=t|Y_{peak}^S\left[t\right]>\delta }$ where δ is a threshold corresponding to twice the standard deviation of the peak pressure time series base value (Cacciari et al., 2017b).

All variables were computed over six distinct supersets of sensor groupings: S_long, S_lat, plane, ring, S_long ∩ S_lat, and whole matrix. The first superset, S_long, includes the posterior, anterior, right, and left groupings (Guaderrama et al., 2005), exactly as depicted in Fig. 2. S_lat includes the cranial, medial, and caudal groupings (Guaderrama et al., 2005); hence, S_long ∩ S_lat is the superset containing the results of the intersection between S _long and S _lat. The remaining two supersets include planes (Eq. 2.8) and rings (Eq. 2.9) of the sensors (Cacciari et al., 2017a):

(2.8)${\displaystyle Plan{e}_i=C_i\cup C_{i+5}}$ (2.9)${\displaystyle Rin{g}_j=R_j}$ where C_i is the group of sensors corresponding to the i-th column and i ϵ [1, 5] interval. R_j is the group of sensors corresponding to the j-th row and j ϵ [1, 10] interval.

Aside from the features extracted on the maximum contraction and Valsalva maneuvers, three other variables were extracted for the endurance maneuver. The integral of pressure (Eq. 2.10) and integral of sum (Eq. 2.11) were computed over the Plane and Ring supersets, while the plateau duration (Eq. 2.12) was computed only for the whole matrix of sensors (Cacciari et al., 2017b):

(2.10)${\displaystyle in{t}_p=\sum _t{Y}_{peak}^S\left[t\right]}$ (2.11)${\displaystyle in{t}_s=\sum _t{Y}_{sum}^S\left[t\right]}$ (2.12)${\displaystyle {\Delta }_{plateau}=max\left(\Delta t|Y_{peak}^S\left[t\right]\ge 0.9\cdot ma{x}_p,\forall t\in \Delta t\right)}$ where Δt is any time interval within the maneuver execution.

In the wave maneuver, the feature’s maximum pressure, instant of maximum pressure, integral of pressure, integral of sum, and instant of activation were extracted from all six aforementioned supersets. Moreover, the rate of contraction (Eq. 2.13) and rate of relaxation (Eq. 2.14) were computed only for the whole matrix (Cacciari et al., 2017b):

(2.13)${\displaystyle CR=\frac{max{Y}_{peak}^S\left[t^{{}^{\prime }}\right]-Y_{peak}^S\left[0\right]}{\underset{t^{{}^{\prime }}}{\text{argmax}}{Y}_{peak}^S\left[t^{{}^{\prime }}\right]}}$ (2.14)${\displaystyle RR=\frac{max{Y}_{peak}^S\left[1\right]-max{Y}_{peak}^S\left[t^{{}^{\prime }}\right]}{1-\underset{t}{\text{argmax}}{Y}_{peak}^S\left[t^{{}^{\prime }}\right]}}$ where t′  is the normalized time of the interval [0, 1].

Finally, sample covariances were computed from the sensor groupings pertaining to the S_long ∩ S_lat superset for the Valsalva, maximum contraction, and endurance maneuvers. The elements of the main and inferior diagonals of the sample covariance matrix were used as features.

Feature selection

The applied process of feature selection was composed of a ranking stage followed by a complete feature subset search (FSS), both independent of the classification process. Two FSS algorithms were tested: extended branch and bound algorithm (BB) and recursive feature elimination algorithm (RFE).

The extended BB FSS algorithm (Duc & Andrianasolo, 1998), with the exception of the functionality of determining variables to be disregarded in the search process, was implemented in the Python programming language using the Mahalanobis distance as the subset evaluation metric. This algorithm, despite performing a complete FSS without exhausting every feature combination, may imply a high computational cost since both the number of evaluated subsets and dimension of the covariance matrix to be inverted during the Mahalanobis distance calculation grow exponentially as a function of the number of features (Webb, 2002).

Given this computational burden, prior to the application of the BB algorithm, the N most relevant features (ranked according to the Pearson and RELIEF criteria) were kept for FSS. Incremental values of N were evaluated using the execution time as a stopping criterion, thus resulting in the value N = 25. Finally, in a similar procedure but with the mean accuracy as the stopping criteria, the range of dimensions of the subsets selected by the BB method was determined to be (Kao & Wei, 2011; Bø, 2004).

The Recursive Feature Elimination (RFE) algorithm iteratively fits a model, computes the model dependent ranking criterion and discards the last M features ranked by the criterion. This model dependency implies that the top ranked features are not necessarily the ones that are individually most important, thus working as a feature subset ranking method. The base model used to compute the criterion was the Random Forest, with the feature importance being its criterion. The model commonly used for this FSS method, the SVM, was not used since it is already part of the evaluated models, and it could create bias in its favor. Finally, the RFE parameters, as well as random forest parameters, both provided by the python library scikit-learn, were selected using an unreplicated regression analysis, with the target being the test accuracy, having rfe_step (M features to remove at a time), n_estimators, min_samples_split and max_features as 2 levels factors. The activeness of which factor was evaluated used Lenth method (Lenth, 2008), from which none was deemed active. This way, the sign of the coefficients was used to determine the values of each factor: rfe_step at 5, n_estimators at 10, min_samples_split at 2 and max_features at 0.5.

Classifiers configuration, selection, and evaluation

The Python package scikit-learn was used to implement the classifiers k-nearest neighbors (k-NN), LR, and SVM, which were evaluated on the classification of the intravaginal pressure data. The higher the numbers of hyperparameters to be explored, the higher the chances of overfitting due to configuration selection process, compromising the classification performance over independent samples (Cawley & Talbot, 2010). Thus, the number of hyperparameters with scanning range was fixed to 1 for each model.

For the k-NN classifier, the hyperparameter corresponding to the number of neighbors (exploring range) was set to the interval (Mumtaz et al., 2017; Bø, 2004). LR and SVM both have a hyperparameter that corresponds to the inverse of the regularization strength (C parameter), the scanning ranges of which were set to [0.0001, 0.001, 0.005, 0.01, 0.05, 0.1] and [0.0001, 0.001, 0.005, 0.01, 0.05], respectively. Then, with exception of the SVM configuration, which had its kernel type fixed to linear, no other default settings of the scikit-learn models were overridden.

Different SVM kernels were initially tested, without a conclusive result of which would be better. Thus, the kernel type of SVM was fixed to linear to reduce the computational cost in searching the better kernel during the hyperparameters selection, besides being less prone to overfitting the data.

Both the evaluation and selection of the configuration were based on leave-one-out (LOO) cross-validation (Längkvist, Karlsson & Loutfi, 2014; Shao, Meng & Wang, 2016). This procedure consisted of two nested LOO stages (Fig. 3), with the inner one used for ordering the models in terms of validation accuracy, whereas the outer one tracked the performance of the best of all three classifiers as well as their individual best configurations.

The whole procedure was repeated for each of the four data acquisition maneuvers, with and without the subset search (Fig. 3, Step 2). Moreover, a combined search scheme using all data acquisition maneuvers was performed; prior to Step 2, each activity underwent the ranking and branch and bound methods, with the output subset dimension fixed at 6. Then, the subsets of each maneuver were concatenated, thus constituting a 24-features subset, which was then fed to the remaining selection procedure.

Accuracy ((True positive + True negative)/all subjects), precision (True positive/(True positive + False positive)), and recall (True positive/(True positive + False negative)) for each model are presented. In this study, accuracy and recall should be carefully checked because there is a high cost associated with patients set as False Negative.

The non-parametric McNemar’s test (Stapor, 2017) was used to check for statistically significant differences (p ≤ 0.05) between the test accuracies of the three classifiers (Häfner et al., 2009) with FSS application. Additionally, the Shapiro–Wilk (for normal distribution testing) and Mann–Whitney U tests were applied to the distributions of the features that presented the highest selection frequencies, as selected by the combined search scheme.

Finally, a PCA was applied on the combined selection scheme with the aim of visualizing the separability between the two classes’ (continent and incontinent) for two situations, using all features of the four maneuvers and using only the features that presented the highest selection frequencies.

Individual maneuvers search

For the endurance maneuver features the highest accuracy was 68.8%, which was obtained by both the best configurations out of the three classifiers and the best configuration of the LR. The highest precision was 70.6%, which was obtained by the best configuration of the k-NN without FSS, and the highest recall was 83.3%, which was obtained by the best configuration out of the three classifiers with FSS (Table 1). LR reached the highest selection frequency for this maneuver, being selected 42 times out of 48 according to the validation accuracy estimated by the inner LOO (Fig. 3), regardless of the FSS application. In addition, with the exception of the k-NN classifier, the FSS application increased the accuracy, the precision, and the recall achieved by the best configuration of each classifier. On the other hand, with the exception of the SVM classifier, the RFE application decreased the accuracy, the precision, and the recall achieved by the best configuration of each classifier (Table 1).

For the Valsalva maneuver features, the highest accuracy achieved was 58.3%, which was obtained by both the best configuration out of the three classifiers without FSS and the best configuration of the LR. The highest precision was 77.8%, which was obtained by the best configuration of the LR with RFE, and the highest recall was 58.3%, which was obtained by the best configuration out of the three classifiers without FSS, LR without FSS, and LR with RFE (Table 1). LR also reached the highest selection frequency for this maneuver, being selected 43 times out of 48 without FSS, although k-NN was selected 40 times out of 48 with RFE. In addition, the application of FSS contributed to an overall accuracy, precision, and recall decrease, except for the k-NN classifier, which had a slight increase. On the other hand, with the exception of the best configuration out of the three classifiers, the RFE application produced an overall accuracy, precision, and recall increase (Table 1).

The highest accuracy achieved for the maximum contraction maneuver features was 60.4%, which was only reached without FSS. This value corresponds to the performance obtained by the best configurations of SVM and LR, with the latter being selected by the inner LOO 48 times out of 48. The highest precision was 60.0%, which was obtained by the best configuration of the SVM without FSS, and the highest recall was 66.7%, which was obtained by the best configuration out of the three classifiers without FSS, LR without FSS, and LR with RFE (Table 1). When the FSS and RFE procedures were employed, however, k-NN reached the highest selection frequency for this maneuver, being selected 40 times out of 48 and 38 times out of 48, respectively. However, both FSS and RFE also contributed to an accuracy, precision, and recall decrease for all classifiers (Table 1).

The highest accuracy achieved for the wave maneuver features was 79.2%, which was obtained only by the best configurations of the k-NN classifier with RFE application. This classifier, presented the highest selection frequency, being selected 43 times out of 48 with RFE, although LR was selected 41 times out of 48 without FSS and 37 times out of 48 with FSS (Table 1). The highest precision was 88.9%, which was obtained by the best configuration of the k-NN with RFE, and the highest recall was 79.2%, which was obtained by the best configuration of the LR with RFE (Table 1). In addition, the FSS and RFE application increased the accuracy, precision, and recall achieved by the best configurations of all classifiers (Table 1).

Regarding classifier performance with BB application, for the endurance maneuver, there was significant difference between LR and k-NN performance as well as between SVM and k-NN performance (Table 2). No significant differences were observed for the other maneuvers. The McNemar test results suggest statistical inferiority of the k-NN when applied to the endurance maneuver features. Regarding classifier performance with RFE application, for the Valsava maneuver, there was significant difference between SVM and LR performance as well as between k-NN and LR performance (Table 2). No significant differences were observed for the other maneuvers. The McNemar test results suggest statistical inferiority of the LR when applied to the Valsava maneuver features.

Combined maneuvers search

With the combined maneuvers search scheme, the highest achieved accuracy was 77.1%, which was obtained by the best configurations out of the three classifiers with FSS, the best configurations of the LR with FSS, the best configuration of the LR with RFE, and the best configuration of SVM with RFE (Table 3). The highest precision was 81.8%, which was obtained by the best configuration of the k-NN without FSS, and the highest recall was 83.3%, which was obtained by the best configuration of the LR with RFE, and the best configuration of SVM with RFE (Table 3). K-NN attained the highest selection frequency amongst the classifiers, being selected 41 times out of 48 with RFE. Overall, the FSS and RFE application increased the accuracy, precision, and recall achieved by the best configurations of all classifiers, with the greatest increase for the LR (Table 3). Moreover, no statistical differences were observed between the classifiers’ performances (Table 4).

From the six extracted features that presented the highest selection frequencies, i.e., the most frequently selected features, significant differences between the classes (continent vs. incontinent women) were observed (Mann–Whitney U test, p < 0.05), except for the instant of activation of the sensor grouping Ring₂ extracted from the Valsalva maneuver (Fig. 4). With these six features, the class separability gain, observed through the two-principal component plot, is quite visible if compared to the one attained with all features of the four activities (Figs. 5 and 6).

When applying the configuration selection procedure (Fig. 3) only on these six more frequently selected features with no further FSS, it was possible to achieve accuracy as high as 97.9% for the best configurations of the LR and SVM classifiers as well as 95.8% for the best configuration out of the three classifiers.

Finally, the test accuracies of the best models out of the three classifiers are presented in Table 5.