## 2 Materials and Methods
### Statistical Analysis  
Before describing the analytical setup applied in this work, it is important to mention that the dataset did not contain any missing values, and therefore all data samples were used. Furthermore, even though we used six clinical rating scales when describing the dataset (see Section 2.1), only fourof these scales were used for the analysis, specifically: UPDRS III, UPDRS IV, RBDSQ, and FOG-Q. The reason is that previous studies have already shown that non-motor manifestations of PD are not linked with the phonatory aspects of HD, but rather with the impairments of prosody and articulation [44] that are commonly being quantified using a sentence reading task, free speech (monologue), etc. Since this study is focused on the phonatory aspects of HD, clinical rating scales describing only motor symptoms of PD were used. To reveal and assess the strength of a relationship between the computed acoustic features and patients' clinical data (UPDRS III, UPDRS IV, RBDSQ, and FOG-Q), [PERSON]'s correlation coefficient was computed (the statistical assumptions for [PERSON]'s correlation coefficient were satisfied as: (a) the acoustic features as well as the the clinical data are both variables that are measured on at least an ordinal scale, and (b) there is a monotonic relationship between the two variables). Since age, gender, and probably L-dopa, are manifested in a voice of PD patients [9], for the purpose of this work, we employed partial [PERSON]'s correlation controlling for the effect of the following confounding factors (also known as covariates): patients' age, gender [29; 45], and dopaminergic medication [32; 46]. The significance level of correlation was set to 0.05. More specifically, two correlation scenarios were considered: (a) correlation between the acoustic features at the baseline and the change in values of the selected clinical rating scales, and (b) correlation between the change in the acoustic features and the change in the values of the selected clinical rating scales. With this approach, we aimed at identifying those acoustic features that are significantly correlated with the specific motor and non-motor symptoms assessed by the selected clinical rating scales in both scenarios. Next, to evaluate the power of the acoustic features at the baseline to predict the change of the patients' clinical data in the horizon of two years, we used the acoustic features computed for the recordings acquired in session 1 (baseline examination) and built mathematical models predicting the change in the selected clinical rating scales (\(\Delta\)). For this purpose, we employed Gradient Boosted Trees (more specifically, the famous XGBoost algorithm [47]) in a supervised learning setup: 10-fold cross-validation with 20 repetitions [48]. The XGBoost algorithm belongs to the state-of-the-art in machine learning, which is supported by the fact that it has been recently used to win competitions on Kaggle. It works well even on small datasets (where it outperforms deep learning approaches), it is robust to outliers and it is able to model complex interdependencies. For these reasons, it has been used by many researchers in various biomedical fields, e.g., [49; 50; 51], etc. The performance of the models (precision of the predictions) was evaluated by MAE and estimation error rate (EER). These measures are defined as:  
\[\text{MAE} =\frac{1}{n}\sum_{i=1}^{n}|y_{i}-\hat{y_{i}}|,\] \[\text{EER} =\frac{1}{n\cdot r}\sum_{i=1}^{n}|y_{i}-\hat{y_{i}}|\cdot 100\,[ \%],\]  
where \(y_{i}\) stands for the true label of \(i\)-th observation, \(\hat{y_{i}}\) represents the predicted label of the \(i\)-th observation, \(n\) denotes the number of observations, and finally \(r\) stands for the range of values in the predicted clinical rating scale (not the range that can be theoretically reached, but the actual range of the values in the dataset).