Analysis of Various Fractional Order Derivatives Approaches in Assessment of Graphomotor Difficulties

Graphomotor disabilities (GD) are present in up to 30% of school-aged children and are associated with several symptoms in the field of kinematics. Although the basic kinematic features such as velocity, acceleration, and jerk were proved to effectively quantify these symptoms, a recent body of research identified that the theory of fractional calculus can be used to even improve the objective GD assessment. The goal of this study is to extend the current knowledge in this field and explore the abilities of several fractional order derivatives (FD) approximations to estimate the severity of GD in the children population. We enrolled 85 children attending the 3rd and 4th grade of primary school, who performed a combined loop task on a digitizing tablet. Their performance was rated by psychologists and the online handwriting signals were parametrised by kinematic features utilising three FD approximations: Grünwald-Letnikov’s, Riemann–Liouville’s, and Caputo’s. In this study, we showed the differences across the employed FD approaches for the same kinematic handwriting features and their potential in GD analysis. The results suggest that the Riemann-Liouville’s approximation in the field of quantitative GD analysis outperforms the other ones. Using this approach, we were able to estimate the overall score with a low error of 0.65 points, while the scale range is 4. In fact, the psychologists tend to make the error even higher.


I. INTRODUCTION
Fractional calculus (FC) is a name of the theory of integrals and derivatives of an arbitrary order [1]. The concept of fractional operators has been introduced almost simultaneously with the development of the classical differential, integral or other well-known calculus [2]. It attracted the interest of many famous mathematicians, including Euler, Liouville, Laplace, Riemann, Grünwald, and Letnikov. The principles of FC have been used in modeling of many physical and chemical processes, as well as in modern engineering and The associate editor coordinating the review of this manuscript and approving it for publication was Donato Impedovo . science in general [3]- [5]. It has been advantageously used during the modeling of different diseases such as the human immunodeficiency virus (HIV) [6] or malaria [7]. Recently, the FC has been significantly examined in computer vision, particularly in image restoration, super-resolution, image segmentation or motion estimation [8]. In line with this trend, in our recent research, we developed new parametrisation techniques of online handwriting (a handwritten signal with temporal information) based on the application of the fractional order derivatives (FD) [9]- [13].
It has been estimated that approximately 10-30 % of children experience graphomotor difficulties [20] such as graphomotor production deficits, motor feedback difficulties (e. g. the pen's tip location tracking problems), motormemory dysfunctions, etc. Considering that children spend 31-60 % of their school-time performing handwriting [21], the early identification of graphomotor disabilities (GD) is crucial in the prevention of serious pedagogical and psychological consequences [22]. Otherwise, a child's every-day life can be greatly affected starting with a lack of motivation to write, a decrease in self-esteem in combination with poor emotional well-being continuing to bad attitude and behaviour, communication and social interaction problems [23]. In some cases, it may result in being diagnosed with a serious neurodevelopmental disorder such as developmental dysgraphia (DD) [24], [25]. To identify and evaluate GD in school-aged children, several well-established questionnaires or tests based on a visual inspection of the handwritten product have been developed [26]. Though, their utilization on a day-to-day basis is still limited due to the fact that the administration and coding are very timeconsuming. Furthermore, the perceptual abilities, experience, and subjective judgment of an examiner are limited as well.
To overcome the limitations of the perceptual GD analysis, researchers have been focusing on computerized quantitative analysis of online handwriting. Pen and paper have been replaced by digitizing tablets used to record a variety of signals describing the evolution of the handwritten product in time. It allowed quantification of kinematic (velocity, acceleration or jerk) as well as dynamic (pen pressure, tilt or azimuth) components of the handwritten signal. For instance, Pagliarini et al. [27] (2017) presented the potential of quantitative analysis to identify the development of handwriting difficulties (HD) at a very early age. Mekyska et al. [14] (2017) built a classifier (random forests; 54 children) identifying the presence of DD with 96 % sensitivity and specificity. Rosenblum and Dror [15] (2017) achieved 90 % sensitivity and specificity in DD diagnosis (support vector machines; 99 children) using various kinematic and dynamic features. Asselborn et al. [16] (2018) reported 96 % sensitivity and 99 % specificity (random forests; 268 children) using 53 handwriting features quantifying different dimensions of handwriting. Next, Mekyska et al. [17] (2019) proposed a model (based on XGBoost, 76 children) and achieved 50 % sensitivity and 90 % specificity in identification of GD presence using 7 basic graphomotor elements quantified by conventional temporal, spatial, kinematic, and dynamic parameters. In 2020, Galaz et al. [13] published a work dealing with advanced analysis (utilising modulation spectra, fractional order derivatives, and tunable-Q wavelet transform) of graphomotor elements in 53 children attending 3rd and 4th grade of elementary schools. Employing random forests they reached 83 % sensitivity and 81 % specificity. In the same year, Asselborn et al. [18] proposed new data-driven based approaches for an assessment of handwriting difficulties, that were divided into 4 dimensions: kinematic, pressure, tilt and static. This novel approach enables a detailed analysis in children having a very similar overall score of dysgraphia, but differing in specific difficulties. Finally, Garot et al. [19] (2020) enrolled 280 children who were performing the Concise Evaluation Scale for Children's Handwriting (BHK) while recorded by digitising tablets. Employing a cluster analysis, the authors were able to automatically discriminate among 3 groups of children associated with dysgraphia: 1) children with mild dysgraphia usually not identified in schools, 2) children with severe dysgraphia manifested in kinematics and pressure, and 3) children with severe dysgraphia manifested mainly in tilt. The overview of the mentioned current works and their achievemnts can be found in Table 1.
Considering the success of utilizing the FD (Grünwald-Letnikov approach) in Parkinson's disease dysgraphia analysis in our previous works [9]- [12], and in the assessment of GD in school-aged children [13], [28], this study, as a next logical step, has the following aims: • to extend our previous research by the employment of several FD-approaches instead of one (Grünwald-Letnikov approach), • to explore the differences of several FD approaches in the assessment of GD in the children population, • to compare the power of the FD-based handwriting features computed by several FD approaches to estimate the severity of GD.

II. DATASET & METHODOLOGY A. DATASET
For this study, we enrolled 85 children (31 girls and 54 boys) attending 3rd and 4th grade at several primary schools in the Czech Republic. The demographic data of the participants can be found in Table 2 and the resulting grade distribution in Table 3. Children were asked to perform drawings, writings, and several cognitive tests based on a protocol consisting of 31 tasks designed in cooperation with psychologists and special educational counselors. Every graphomotor task of the protocol has been evaluated by a wellexperienced psychologist and rated on the scale from 0 to 4 where: 0 -no graphomotor difficulties; 1 -mild graphomotor difficulties; 2 -graphomotor difficulties; 3 -dysgraphia; 4 -severe dysgraphia. Finally, an overall score has been assigned to each child based on a complex analysis of all the 31 tasks in the protocol (i.e. including the cognitive tests). Although the protocol contains 7 graphomotor tasks such as Archimedean spiral, loops, sawtooth, or rainbow, in this study, we focused on one graphomotor task (combined loops), which has been proved to discriminate well between children with/without graphomotor difficulties [17].The distribution of scores (the overall and the subscore for the combined loops task) is presented in Fig. 1. Correlation between the scores and the demographic data is visualized in Fig. 2. Parents of all children participating in this study signed an informed consent form approved by the Ethical Committee of the Masaryk University. Throughout the entire duration of this study, we strictly followed the Ethical Principles of Psychologists and Code of Conduct released by the American Psychological Association (https://www.apa.org/ethics/code/).

B. DATA ACQUISITION
At first, a template of the combined loop task was shown to a child and then he/she was asked to replicate it on an A4 paper that was laid down and fixed to a digitizing tablet. The drawing was acquired by the Wacom Intuos Pro L (PHT-80) digitizer with the sampling frequency of 150 Hz, and the Wacom Inking pen, which provides a feeling of writing by a regular pen and offers immediate visual feedback.  ). For more information, see our previous works [12], [14], [28]. An example of the selected combined loop task performed by a child with/without GD can be seen in Fig. 3.

C. FRACTIONAL ORDER DERIVATIVES
The essential of this study is the investigation of the several (non-equivalent) FD approximations as a new advanced approach of drawing/handwriting parameterisation. We developed this method to substitute the conventional differential derivatives in the feature extraction process (see our previous works [9]- [12], [28]) in order to improve the quantitative analysis of the GD. In the scope of this study, we utilized three FD approximations: Grünwald-Letnikov (GL), Riemann-Liouville (RL), and Caputo (C), implemented by Valério Duarte in Matlab [29]- [31].

1) GRÜNWALD-LETNIKOV
The FD definition by Grünwald-Letnikov is one of the first and basic approaches [2]. A direct definition of the derivation of the function y(t) by the order α -D α y(t) [1] is based on the finite differences of an equidistant grid in [0, τ ], assuming that the function y(t) satisfies certain smoothness conditions in every finite interval (0, t), t ≤ T , where T denotes the period. Choosing the grid with and using the notation of finite differences where The Grünwald-Letnikov definition from 1867 is defined as where GL D α y(t) denotes the Grünwald-Letnikov derivatives of order α of the function y(t), and h represents the sampling lattice.

2) RIEMANN-LIOUVILLE
Another classical form of the FD has been given by Riemann-Liouville. The left-inverse interpretation of D α y(t) by Riemann-Liouville [1], [3] from 1869 is defined as where RL D α y(t) denotes the Riemann-Liouville derivatives of order α of the function y(t), is the gamma function and n − 1 < α ≤ n, n ∈ N, t > 0.

3) CAPUTO
Nowadays, the most significant contributions to the field of FC are the results achieved by M. Caputo [32]. In contrast to the previous ones, the improvement hereabouts lie in the unnecessity to define the initial FD condition [1], [3]. The Caputo's definition from 1967 is where C D α y(t) denotes the Caputo derivatives of order α of the function y(t), is the gamma function and n − 1 < α ≤ n, n ∈ N, t > 0.

D. HANDWRITING FEATURES
Altogether, we extracted 3 sets of handwriting features, one feature set per one employed FD approach. Basic kinematic VOLUME 8, 2020 features from the input handwritten signal were extracted as well, namely velocity, acceleration, jerk and their horizontal and vertical variants. Due to rare omissions of 3-4 samples by the digitizing tablet during the acquisition, we performed the in-signal outliers removal (outliers were considered as elements more than three scaled median absolute deviations from the median). If not pre-processed, the differentiation of this gap would leave significant peaks in the output handwriting feature as illustrated in Figure 4. All handwriting features were computed for α in the range of 0.1-1.0 (with 0.1 step), where α = 1.0 is equal to the full derivation. Finally, the statistical properties of all extracted handwriting features were described by the mean and the relative standard deviation (relstd). To sum up, each feature set consists of 180 computed kinematic features.

E. STATISTICAL ANALYSIS
At first, we performed the normality test of the handwriting features using the Shapiro-Wilk test [33]. In the case of non-normally distributed features, we utilised the Box-Cox transformation [34]. Next, to assess the strength of the relationship between the feature values and the scores (the overall score and the sub-score), Spearman's and Pearson's correlation coefficients were computed (we considered the level of significance 0.05). The p-values were adjusted using the False Discovery Rate (FDR) method to address the issue of multiple comparisons.
During the statistical analysis, we controlled for the effect of several confounding factors (covariates), namely age, grade, and sex.
Finally, to evaluate the power of the handwriting features to support the estimation of scores assessing the GD, we performed a multivariate analysis. For this purpose, we employed the state-of-the-art algorithm XGBoost [35] (10-fold cross-validation with 20 repetitions). The XGBoost algorithm was selected, because of its ability to achieve good performance on a small data set. Moreover, it is able to compete with the deep learning methods that are still not being used in the case of the small dataset as they require much larger data to be trained on [36]. Hyper-parameter space optimization was performed by a random search strategy with following parameter values:   Table 4).Framed sub-areas in each cross-correlation matrix visually isolates the handwriting features computed by the same FD approach. FIGURE 6. Distribution of the FD order α for the features mostly correlating with the overall score, see Table 4 (GL -Grünwald-Letnikov; C -Caputo; RL -Riemann-Liouville).

III. RESULTS
The results of the correlation analysis can be seen in Table 4. The table shows the top 5 features per FD approximation according to p-values of the Spearman's correlation related to the overall score (upper part) and the sub-score (bottom part). The strongest correlation (after the FDR adjustment) with the overall score was identified in features extracted by the Caputo's FD. However, in the case of the sub-score, the Riemann-Liouville's FD arises as the most significant.
The correlation matrices (using the Spearman's correlation) are visualized in Fig. 5. Each matrix includes the top 5 features per FD approximation (i. e. 15 features in one matrix) identified in Table 4. The distribution of the FD order α of 20 best features regarding the Spearman's correlation per FD approximation is visualised in Fig. 6 for the overall score and in Fig. 7 for the sub-score.
Finally, the results of the multivariate analysis can be found in Table 5. In the case of the overall score estimation, the best results were achieved by the Riemann-Liouville FD.  Table 4 (GL -Grünwald-Letnikov; C -Caputo; RL -Riemann-Liouville).
In the case of the sub-score estimation, the lowest error was achieved when combining features of all the approximations. Hyper-parameters of the best XGBoost models can be found in Table 6.

IV. DISCUSSION
The main goal of this study is to explore the differences across various FD approximations utilized in the analysis of the GD. A comparison of an identical feature (i. e. velocity for α = 0.2) extracted from the handwritten product associated with the GD (the same sample as in the bottom part of Fig. 3) is shown in Fig. 8. It illustrates the differences across the involved FD approximations. The velocity function extracted by the Caputo's FD dominates by significant peaks in the positions, where a child interrupts the performance for a moment and then continues writing. These interruptions are also visible in the function computed by the Riemann-Liouville approach, though in the form of a constant line followed by elevated oscillations instead of peaks. On the VOLUME 8, 2020  other hand, the function based on the Grünwald-Letnikov approach seems to be a constant line, nevertheless after a scale normalization (min-max normalization), see Fig. 9, it is clear that the function has the oscillatory nature as well.
The differences across FD approaches are underlined by the comparison in Fig. 10, where the dependency of the relative standard deviation of the velocity on the FD order α is visualized. Feature values computed by the Grünwald-Letnikov approach are generally higher in comparison with the Caputo and Riemann-Liouville ones, which are more similar. On the other hand, the envelope of the velocity profile based on the Grünwald-Letnikov approach is more similar to the Riemann-Liouville one. Moreover, all functions meet at the point where α = 0.9 and continue simultaneously to the full derivation (α = 1.0), which is expected, because the full derivation has to be the same for all approaches.
Experts in the field of psychology need to understand and clearly interpret the results of the graphomotor analysis, i. e. to link them with specific symptoms or physiological processes. This is very challenging especially in the case of advanced signal parameterisation, which is also our case. Therefore, to bring credibility for a non-technical reader, we provide an illustration in Fig. 11. In this figure, we compare the vertical projection of the movement (y axis) and the vertical velocity (Grünwald-Letnikov approach, α = 0.8) in a child without graphomotor difficulties (same as in the upper part of Fig. 3). The function extracted by FD for α = 0.8 is difficult to be understood, but the relationship to the velocity is obvious.
Regarding the results of the correlation analysis (association with the overall score), the most significant features (after the FDR adjustment) are extracted by the Caputo's FD, where the top 5 have the p-value < 0.05. Most significant handwriting features are related to the variability of the jerk, which refers to the disturbances in the fluent handwriting performance of the child with GD. The values of the correlation coefficients are negative, which means that the handwriting performance of the subject is worse with the lower variability of the jerk. This may be confusing, because just the opposite effect may be expected. Nevertheless, this is specific for the combined loop task. A child without GD is less focused on the writing (the movement is more automatic), therefore the changes between loops are more dynamic, which results in higher jerk variability. Vice versa, a child with GD is more focused on his/her performance, therefore, the handwriting is associated with lower acceleration and jerk. In the case of Grünwald-Letnikov based features, 4 out of the 5 most significant ones are jerk  related too, what supports the results obtained by the Caputo's approach. In the view of the Riemann-Liouville FD, the most significant features are mostly acceleration and jerk related, this likewise supports the association with the smooth handwriting disabilities.
Considering the correlation with the sub-score, the most significant features (after the FDR adjustment) are extracted by the Riemann-Liouville FD, while 4 out of 5 features are acceleration-based. This again refers to the disruptions in continuous handwriting of a child with GD (i. e. less automatic and dynamic movements). In the case of the Grünwald-Letnikov approach, the variation of the velocity is observed to be the most significant, however, none of the features is significant after the p-value adjustment (similarly to the Caputo's approach). Due to the omission of the full derivations in best correlation results, the FD-based features outperform the conventional handwriting features in the scope of the sub-score correlation analysis for the connected loops task. In addition, this is in line with our previous results. [11], [12].
Regarding the cross-correlation of the top-ranked features strongly associated with the overall score (see the left matrix in Fig. 5), we did not observe any strong correlations among the features based on the Caputo's approach. In the case of the Riemann-Liouville's approximation, we identified a significant correlation between the mean of the vertical acceleration and the relstd of the horizontal acceleration, in both features α = 1, which means full derivation. Similarly, in the Grünwald-Letnikov's approach, we identified a strong association between the relstd of the horizontal acceleration, and the mean vertical jerk and the mean vertical acceleration. The last two mentioned features are in fact very close to each other, because the acceleration with α = 1 is very similar to the jerk with α = 0.2 We assume that the above-mentioned association is linked with the fact that the vertical movement, contrary to the horizontal one, requires coordinated movement and finer flexions/extensions of more joints (interphalangeal and metacarpophalangeal) and therefore, it is more complex than ulnar abductions of the wrist [37], [38]. Since the vertical movement is complex, it is strongly affected by psychological and muscular fatigue [39], which could be manifested in lower vertical acceleration in children with GD. Nevertheless, low relstd in the horizontal direction could mean monotonous and less dynamic movement too.
In the case of the cross-correlation matrix linked with the sub-score, we can observe significant correlations only in features that express the same information, e. g. the mean of the horizontal acceleration, but differ only in α, e. g. the difference is 0.1. Since this difference is very low, it is obvious that these features significantly correlate. Except for this, the features do not correlate much among themselves which means that they are not redundant, but still relevant (see Table 4).
Based on the distribution of α in the 20 top-ranked features, we can observe that those based on the Caputo's approach are mostly concentrated around 0.2 and 0.5 for the overall score and almost evenly distributed in the case of subscore correlation analysis. The Grünwald-Letnikov FD-based features associated with the overall score have α concentrated around 0.2 and 0.9. Those associated with the subscore are mainly around 0.2, 0.5 and 0.7. Finally, in the case of the Riemann-Liouville's approach, we can observe a higher concentration in the range [0.4; 0.6] for the overall score, and in the range [0.7; 0.9] for the sub-score. Since the distribution of the α varies per FD approximation and rat-  ing scale, we hypothesise that further and finer optimization of this parameter would bring even better quantification of the GD.
Concerning the multivariate analysis (Table 5), where we estimated the overall score, the best results were achieved by the Riemann-Liouville FD-based features. The resulting MAE was 0.65, and RMSE = 0.79. When estimating the sub-score, all approaches had a very similar MAE, nevertheless, the lowest RMSE (0.79) was reached by the Riemann-Liouville's approach too. A combination of all the approaches slightly decreased the error. These results suggest that the Riemann-Liouville's approximation in the field of quantitative GD analysis outperforms the other ones. In addition, using this approach we were able to estimate the scores with MAE = 0.65 and MAE = 0.66, respectively. If we take into account that the range of the first scale is 4, and of the second one 3, the error can be considered as very low. In fact, when assessing GD in children, psychologists tend to make the error even higher, e. g. two experts can frequently differ by 1 point (compare it to 0.65 or 0.66).

V. CONCLUSION
To the best of our knowledge, this is a unique study that performs an investigation of the various FD approaches in the computerized assessment of the GD in school-aged children. Therefore, it should be considered as being rather exploratory and pilot in nature. We can conclude that the employment of various FD approximations brings major differences in kinematic handwriting features. In the scope of the correlation analysis associated with the overall score, the Caputo's FD approach exceeds the rest of the analysed FD approximations. However, in the scope of the sub-score, the Riemann-Liouville gained the most significant features. Moreover, the results of the multivariate analysis suggest that the Riemann-Liouville's approximation in the field of the quantitative GD analysis outperforms the other ones (MAE = 0.65 for overall score and MAE = 0.66 for sub-score).
This study has several limitations and possible parts, that could be further improved. First of all, the dataset is relatively small in terms of the statistical validity of the results. To generalize the results, the larger dataset have to be acquired and more handwriting tasks should be included in the analysis. Next, a more granular FD α order search (step of 0.01 or even less) in order to find the optimal α range should be performed. Moreover, other feature types, such as temporal, spatial, and dynamic, should be included in future comparisons. The future study should be detailly focused on the comparison of the FD-based features with the conventionally used handwriting features. The different handwriting tasks have to be investigated separately for the best performing FD-based features. Besides, when comparing the several feature sets performance (regression, etc.) an ANOVA test should be performed in the future to analyze the differences between them. Finally, various machine learning models should be trained and compared in the future studies to get more information about the classification performance of the proposed features and to obtain the most robust models for GD identification.
ZOLTAN GALAZ was born in Galanta, Slovakia. He received the Ph.D. degree in telecommunications from the Brno University of Technology, Czech Republic, in 2018.
He is currently the Team Leader of the Brain Diseases Analysis Laboratory, Department of Telecommunications, Brno University of Technology. He has coauthored more than ten publications indexed by the Web of Science. In cooperation with neurologists/psychologists from different countries, he develops diagnostic/monitoring systems focused on disorders, such as Parkinson's disease and developmental dysgraphia. He deals with the research of advanced digital signal processing and machine learning algorithms applied in the field of computerized non-invasive neurodegenerative and neurodevelopmental disorders diagnosis, analysis, and monitoring. He has authored more than 50 articles indexed in the ISI Journal citation report, more than 100 conference papers, around ten books, and responsible of ten national and European research projects. He is a Spanish liaison of EURASIP. She has been with the University of Hradec Kralove (part-time) since 2005. She has been cooperating on a long-term basis with the University of South Bohemia, Ceske Budejovice, on the research of the relationship between personality and language (the CPACT project). She is currently managing the DAQUAo project focused on the study of psycholinguistic phenomena of therapeutic diaries. Within the Dysgraphy project, she participates in data collection and partial analysis, focusing on pre-school children.

VOJTECH ZVONCAK was born in
JIRINA BEDNAROVA was born in Jindrichuv Hradec, Czech Republic. She received the degree in psychology from the Faculty of Arts, Masaryk University, Brno, Czech Republic.
Since 1995, she has been a Special Educational Counselor. She deals with the diagnosis and intervention in children with inequalities of cognitive and graphomotor development, and specific learning disabilities. She also deals with the education of adults working with such children. She has also authored several diagnostic materials (standardized and dynamic diagnosis) and a number of publications aimed at targeted support of children's development.
ZDENEK SMEKAL (Senior Member, IEEE) was born in Brno, Czech Republic. He received the M.Sc. and Ph.D. degrees in 1973 and 1978, respectively. He is currently a full-time Professor with the Brno University of Technology. He has been an Associate Professor since 1986 in the field of electronics and communication engineering and a Professor since 2000 in the field of theoretical electrical engineering.
He is currently the Chairman of the board for doctoral study programme of teleinformatics and the a Leading Researcher with the Signal Processing Laboratory, Department of Telecommunications, University of Technology. He has coauthored more than 90 publications indexed by the Web of Science and authored five monographs on digital signal processing.
He is a member of Scientific Council with the Brno University of Technology. He was involved in European COST research projects, including OC