Are Multilevel functional models the next step in sports biomechanics and wearable technology? A case study of Knee Biomechanics patterns in typical training sessions of recreational runners

This paper illustrates how multilevel functional models can detect and characterize biomechanical changes along different sport training sessions. Our analysis focuses on the relevant cases to identify differences in knee biomechanics in recreational runners during low and high-intensity exercise sessions with the same energy expenditure by recording $20$ steps. To do so, we review the existing literature of multilevel models, and then, we propose a new hypothesis test to look at the changes between different levels of the multilevel model as low and high-intensity training sessions. We also evaluate the reliability of measures recorded in three-dimension knee angles from the functional intra-class correlation coefficient (ICC) obtained from the decomposition performed with the multilevel funcional model taking into account $20$ measures recorded in each test. The results show that there are no statistically significant differences between the two modes of exercise. However, we have to be careful with the conclusions since, as we have shown, human gait-patterns are very individual and heterogeneous between groups of athletes, and other alternatives to the p-value may be more appropriate to detect statistical differences in biomechanical changes in this context.


General overview and motivation
Advances in biosensors and digital medicine are improving disease monitoring and detection. A promising field for implementing these novel strategies is sports training and biomechanics. These tools are a crucial element for optimizing athlete training and reducing the incidence of sports injuries. Multiple repeated measurements are collected from each individual over different sessions, weeks, or the entire season in these domains. So, it is essential to evaluate the changes produced along with relevant outcomes at different resolution levels scales among individuals. In addition, much of the information recorded is of a functional nature, such as the cycle of gait movement. Functional gait analysis enables a more accurate assessment of the effect of fatigue and the detection of potential injury risk factors. This paper illustrates how multilevel functional models can detect and characterize biomechanical changes along different training sessions. Besides, the multilevel models can provide a vectorial representation of different athlete training activities and feed supervised predictive models in various modeling tasks. Our analysis focuses on the relevant cases to identify differences in knee biomechanics in recreational runners during low and high-intensity exercise sessions with the same energy expenditure by recording 20 steps. To do so, we review the existing literature of multilevel models and then, we propose a new hypothesis test to look at the changes between different levels of the multilevel model as low and high-intensity training sessions. We also evaluate the reliability of measures recorded in three-dimension knee angles from the functional intra-class correlation coefficient (ICC) obtained from the decomposition performed with the multilevel funcional model taking into account 20 measures recorded in each test. The results show that there are no statistically significant differences between the two modes of exercise. However, we have to be careful with the conclusions since, as we have shown, human gaitpatterns are very individual and heterogeneous between groups of athletes, and other alternatives to the p-value may be more appropriate to detect statistical differences in biomechanical changes in this context.

Introduction
In recent years, there has been a big increase in the availability of powerful biosensors. These are now capable of monitoring an individual's energy expenditure with great accuracy and measuring various physiological and biomechanical variables in real-time. This provides the opportunity to have a unique assessment of an individual's physical capability [Lencioni et al., 2019] and performance and thus, be able to schedule optimal interventions over time [Kosorok andLaber, 2019, Buford et al., 2013]. One field that can benefit from the intensive use of these technologies is biomechanics [Ibrahim, 2021, Uhlrich et al., 2020. In both sports and general populations, abnormal movement patterns are synonymous with muscular and motor problems, risk of injury, or even the appearance of severe neurological diseases such as Parkinson's [Morris et al., 2001]. Therefore, the detection and characterization of these abnormal movement patterns in biological activities such as walking and running, are essential in areas beyond professional and sports medicine, such as clinical medicine [Chia et al., 2020].
Nowadays, with the growing boom of wearables, these technologies are be-ing democratized, and their use is increasingly more common among the general population, such as amateur runners. In this setting, the remote control of athlete training and even monitoring their daily routine out of sport activity is feasible. Although we are in the early stages of this technological revolution, the first research papers are appearing, which through high-resolution data gathered with biosensors, can begin to answer unknown and complex questions about the relationship between training load [Cardinale and Varley, 2017], daily biomechanical patterns [Karas et al., 2019], and injury prediction [Bittencourt et al., 2016, Malone et al., 2017. Furthermore, they may even enable us to build predictive models that support decision-making and help optimize the performance [Matabuena and Rodríguez-López, 2019, Hemingway et al., 2020, Piatrikova et al., 2021. For example, several contemporary works provide new epidemiological knowledge using biomechanical data of human locomotion [Karas et al., 2019], [Warmenhoven et al., 2020].
Other papers have tried to predict sports injuries [Rossi et al., 2018] or other motor or neurological diseases prematurely [Belić et al., 2019], or even the impact of therapy together with their prognosis in the recovery phase after surgery [Karas et al., 2020].
The rising proliferation of running as a sporting activity carries a substantial risk for recreational runners who often perform high-intensity training as interval sessions without a formal training schedule. For recreational runners there has been an increased prevalence of running related injures the most common of which is the knee [Van Gent et al., 2007, Messier et al., 2008.
To date, several works have studied the aetiology of running related knee injuries in recreational runners, some even using 3-dimensional analysis [Messier et al., 2008]. However, to the best of our knowledge, no studies have compared biomechanical changes during high intensity interval training (HIIT) compared to lower intensity continuous running. Moreover, some essential questions remain unanswered, for example the reliability of biomechanical measures at the knee in two or more HIIT training sessions.
Traditionally, gait analysis has been performed at fixed points within the gait cycle. A more detailed and meaningful analysis can be attained by using a complete stride cycle with functional data analysis (FDA) [Febrero Bande and Oviedo de la Fuente, 2012] techniques. These analyses can provide a greater insight into understanding the mechanics of locomotion especially as a runner fatigues. Under different fatigue conditions it would be possible to identify with greater clarity, changes that take place within each part of the gait cycle.
The predominant data-analysis practice in biomechanics is to summarize the curve recorded for each stride using several statistical metrics and applied standard multivariate techniques, although there is a loss of information with this approach. There is a rising popularity of using functional data in biomechanics, with the purpose too bridge this gap between complex statistical modeling with functional data and standard analysis between practitioners. An interesting paper [Warmenhoven et al., 2020] explains a global target audience using principal functional component analysis to more accurately analyze biomechanics data. However, this methodology does not take into account that we can obtain multiple strides per individual in each training session.
Here, the general procedure is to normalize the curve obtained for step and body segment to the [0, 1] interval and take the mean of the different curves recorded and create an average functional curve for analysis. Nevertheless, this procedure can be suboptimal because the constructed mean representation ig-nores the individual variability between the distinct steps of the same individual, something crucial in the evaluation of the movement patterns in some settings. In addition, the mean curve can be a summary measure of the information that is very sensitive to outliers, something that is frequent in biomechanical data. This is particularly true in the measurement of movements performed at high or low speed, where sensor and/or human variability often increase. Moreover, we often need to compare the effect of an intervention along the different training sessions on different days, and for this, we have several repeated measures per individual in different periods. For all this, a more natural analysis is to exploit the advantages of multi-level functional models that allow the analysis of several hierarchy levels. With these models, we can incorporate in a natural way, into the same statistical models, biomechanical patterns using a significant fraction of training session, test or event data from a complete session. These methods also allow the capture the variations in different periods at a intra-inter individual level.
Surprisingly, there is little use of these FDA techniques within the literature, either in sport or indeed other clinical areas [Ullah and Finch, 2013]. Both of these areas would benefit from larger data sets, be that longitudinally or more strides or conditions. The objective of this paper is two-fold. First, we will introduce the analysis of multi-level FDA from the methodological point of view. After, we will illustrate from an applied point of view that these biomechanical methods analyze several exciting research questions. For this purpose, we use a sample composed of 19 athletes during two different training sessions, one moderate and one high speed, in a controlled laboratory environment. During these training sessions we measured knee patterns with a tridimensional sensor of 20 strides during the stance phase.
The structure of the paper is as follows. First, we introduce multi-level models with FDA, as review the literature. We then describe the study sample; we do a more in-depth analysis of tridimensional changes to acquire new biomechanic knowledge. Finally, we discuss the results and future challenges in multi-level models in biomechanics and other sport biosensor data.
Given the double target audience of this paper, to maintain the interest for biomechanical practitioners that do not have a specific interest in the mathematical details, we a illustrate multi-level methodology to show biomechanical applications to our data analysis example: • What are the correlations between knee functional running patterns during a HIIT training session and the loss of force production in training?
• What is the reliability of the functional running parameters in two independent HIIT training seasons?
• Are there differences in knee angles between a continuous running session and HIIT training?
• Are functional biomechanics patterns very individual between runners?
• Is it appropriate to use p-value to detect biomechanical changes in the practice? 4 2 Multi-level functional data analysis 2.1 General overview of state-of-art of multi-level models Functional data analysis [Febrero Bande and Oviedo de la Fuente, 2012, Cuevas, 2014, Wang et al., 2016 with a multi-level structure and repeated measurements is a field that has received substantial attention in recent years (see for example [Lee et al., 2018, Li et al., 2020b). In the statistical community, it appeared in the literature as an essential new methodology to statistical practitioners. These techniques have been applied successfully to answer central scientific questions in such heterogeneous domains, to study the variability between subjects, days, tests, physical activity patterns, speech, or sleep quality monitoring [Xiao et al., 2015, Huang et al., 2019, Park et al., 2018, Martinez et al., 2013, Pouplier et al., 2017, Di et al., 2009. Probably the first work that addressed mixed functional data problem was back in 2003-in this work [Morris et al., 2003], using wavelets and a Bayesian estimation procedure, the effect of type f dietary fat on O6-methylguanine-DNA-methyltransferase (MGMT), an important biomarker in early colon carcinogenesis, was explored. Since then, different models have appeared in the functional data analysis literature that have modeled different hierarchy levels, including nested and crossed structures-using diverse estimation strategies adapted to the nature of the real problem that the authors treat to solve. From a general point of view, the different data characteristics involved in the creation of new mixed functional data models are the number of data recorded [Zipunnikov et al., 2011], the density of the functional data [Di et al., 2014], the number of replicates in each unit of the hierarchical structure [Zipunnikov et al., 2011], the structure of dependence between levels of the hierarchy or replicates [Staicu et al., 2010, Staicu et al., 2012, or the dependence between covariates in multidimensional functional problems [Volkmann et al., 2021]. For example, in some relevant applications, such as the analysis of longitudinal data obtained from medical images using nuclear magnetic resonance, the high computational demands of the estimation of image correlation operators and the calculation of projections between subjects and visits have resulted in a series of papers using new computationally efficient multi-level methods that scale well in problems involving millions of data [Zipunnikov et al., 2011, Zipunnikov et al., 2014. In a similar way, essential progress has been made in recent years in the smoothing of correlation operators, whereby at the moment, in problems with half a million covariates, the smoothing can be done in a few seconds [Xiao et al., 2016, Cederbaum et al., 2018. In the reverse situation with sparse functional data, efficient estimation methods have also been proposed [Di et al., 2014, Xiao et al., 2018, Li et al., 2020a, but from a statistical perspective, data variability increases due to the low density of functional data in this framework, and the estimation problems are magnified. Using this approach, several works have proposed different inferential contributions such as re-sampling bootstrap methods to perform inferential tasks such as calculating confidence intervals or comparing the equality of means between groups of subjects by exploiting the rich source of information several measures of the same individual in biological problems [Crainiceanu et al., 2012, Goldsmith et al., 2013, Park et al., 2018. Also, some of the previous multilevel models have been generalized to introduce the impact of specific covariates on the levels of variability of the different levels of hierarchy and subjects so in su-pervised and unsupervised problems [Crainiceanu et al., 2009, Gertheiss et al., 2013, Xiao et al., 2015, Scheipl et al., 2015. Furthermore, new methods have been proposed in a more general set up as other complex objects such as functional matrix structure data [Huang et al., 2019]. This technique, has facilitated opportunities to study the variations of physical activity patterns in a group of subjects with cardiac pathology over several weeks, highlighting the power of these models to solve real complex problems.

Mathematical foundations of standard functional principal component analysis
Let X(t), t ∈ [0, 1], be a random function with mean µ(t) = E(X(t)) and covariance function Σ(t, s) = E((X(t) − µ(t)(X(s) − µ(s)) for all t, s ∈ [0, 1]. The heart of many functional data analysis models is based on calculating modes of variability of the random function X(t) based on the spectral decomposition of the covariance operator Σ(·, ·) in a set of eigen-functions {e i (·)} ∞ i=1 and eigenvalues {λ i } ∞ i=1 , and where, we suppose that λ 1 ≥ λ 2 ≥ · · · . Thus, from the decomposition of Karhunen-Loève we have where, c k = 1 0 (X(t) − µ(t))e k (t)dt being c k 's incorrelated random variables with mean zero, and variance λ k . These variables are usually known as scores or loading variables.
In the real-world setting, we have n realizations generally independent of the process X(·), X 1 (·), . . . , X n (·), but only we observe a sample of n vectors of length m, X 1 , . . . , X n , in a grid {0 ≤ t 1 < · · · < t m ≤ 1}, and where X i j = X i (t j ) for all i = 1, . . . , n, j = 1, . . . , m. Next, by simplicity, to refer X i j , we use X i (t j ).
The simpler estimator of Σ is, whereμ(t j ) = 1 n n i=1 X i (t j ), for all j = 1, . . . , m. A major step here, in many applications where observations are subjected to a large measurement error, is the smoothing process, to ensure the optimal performance of the empirical estimatorΣ. Three different strategies have generally been used in the literature [Shang, 2014, Cederbaum et al., 2018: i) Smoothing of the original functional data; ii) Introduction of a regularization term in the estimation ofΣ; and iii) direct application of a smoothing procedure in the raw estimation ofΣ. Subsequently, we denote byΣ Smooth the smoothed version of Σ by any of the three previous procedures. The next step is to calculate the auto-vectors and auto-values ofΣ Smooth , according to the spectral theory of linear algebra, as happens in the classical context of principal component analysis in multivariate statistics. After performing this procedure, and selecting the first K < m auto-vectors {ê i } K i=1 and auto-values {λ i } K i=1 we obtain the following decomposition: kêi k (i = 1, · · · , n; j = 1, · · · , m), whereĉ i = X i −μ,ê i , denoting with , the usual scalar-product and beinĝ e i j the j component of the auto-vectorê i . More details of this procedure can be found in the following reviews and general books of functional data analysis [Horváth and Kokoszka, 2012, Shang, 2014, Kokoszka and Reimherr, 2017, where different estimation procedures of the number K, of components are established [Li et al., 2013]. For more theoretical aspects of the estimators such as asymptotic properties we refer the reader to [Hall and Hosseini-Nasab, 2006].

Introduction of functional multilevel models
In the previous Section, we have seen how to carry out a principal component analysis when we observed n independent functional data. In practice, in biomedical and sports applications, when patients or athletes are analyzed at different moments in time, for example, the training load throughout a season, it is common to have several repeated and correlated measurements. Therefore, the previous procedure may be inadequate. Before starting, we introduce some extra notation to describe this scenario.
To start, let-consider the two-way functional ANOVA model, whose structure is introduced below: where µ(t) is the mean global, ν j (t), is the mean of measure j, Z i (t) is the subject-specific deviation from the measure-specific mean function, and W i,j (t) is the residual subject-and measure-specific deviation from the subject-specific mean. In this framework, µ(t) and ν j (t) are treated as fixed functions, while Z i (t) and W i,j (t) are treated as random function of mean zero. Moreover, with the proposal of identification correctly the model, we assume that Z i (t), and W i,j (t) are random uncorrelated functions. In many applications, we note ν j (t) could be set to zero when functional responses are interchangeable within different measures, and the model becomes a one-way functional ANOVA.
In the literature of multi-level models, the functions Z i (t)'s are known as the 1-level of functions, while W i,j (t)'s functions composed the 2-level.
Again, the foremost step in a multi-level functional component analysis model is to rely on the Karhunen-Loève decomposition. For example, in the model defined by the equation 4, we have where {e k } ∞ k=1 are the auto-functions related to the random functions of the levels 1, 2, respectively, while {c i k } ∞ k=1 , {d i,j k } ∞ k=1 are the scores o loading variables for all i = 1, . . . , n; j = 1, . . . , J.
(6) Importantly, in this models, the functions {e k } ∞ k=1 are orthonormal basis in the space of square functions, but in general the functions that compose each function bases are not orthogonal with each other, which implies that the estimation of scores is not simple in practice, a topic that we discuss later. Moreover, score's {c i k } ∞ k=1 and {d i,j k } ∞ k=1 are random variables of mean zero and with variance are given by covariance funtions of stochastic process Z i (t)'s and W i,j (t)'s.
Below, we explain how to calculate the auto-functions and auto-values of the model defined in the equation 6. Let Σ T (s, t) = Cov(X i,j (s), X i,j (t)) be the overall covariance function and Σ B (s, t) = Cov(X i,j (s), X i,k (t)) the covariance function between the units of the second level setting the effect of the first level. Applying Mercer's theorem, we can see where the indices T ,B, and W are used to refer to the "total," "between," and "within" subject covariances.
As in the previous Section, the curves are observed uniquely, in a grid of points {0 ≤ t 1 < · · · < t m ≤ 1} and, in this situation we have to perform the empirical estimatorsΣ T (t s , t k ),Σ W (t s , t k ),Σ B (t s , t k ). Unlike in the preceding Section, we only observe directly information from the process X i,j (t), and, it is possible to estimate the covariance matrix,Σ T (t s , t k ), according the usual empirical estimator, that is, To estimate the covariance operator Σ B (s, t), it is enough to appeal to the method of moments or in an equivalent way to the covariance estimator through a U -statistic estimator, is not necessarily a defined positive matrix in the sample context, we have to trim the eigenvalue-eigenvector pair where the eigenvalue is negative.
Finally, for the calculation of the scores {c i k } ∞ k=1 y {d i,j k } ∞ k=1 (for all i = 1, . . . , n; j = 1, . . . , J). Different estimation strategies can be used. We highlight computationally intensive methods such as Markov chain Monte Carlo (MCMC) [Di et al., 2009], projection algorithms designed for this problem, or the most computationally efficient and used in practice method: the best linear unbiased prediction estimator for mixed models (BLUP) [Robinson et al., 1991], which is detailed for the multilevel problem defined in this Section, in the following reference [Di et al., 2014].

More general extensions
Different levels of hierarchy may appear in real problems that can be nested as in the previous Section or crossed. Following [Shou et al., 2015], the different situations that usually occur are listed in Table 1. All these models have the same structure The latent processes are assumed to be zero-mean and square-integrable so that they are identifiable, and the standard statistical assumptions for scalar outcomes can be generalized to functional data. In this way, the total variability of a functional outcome is decomposed into a sum of process-specific variations plus σ 2 .
Both nested and crossover models can be used to employ a general estimation strategies. Below we summarize the steps necessary to do so, which are analogous to those explained in the previous Section: 1. Estimate the means and covariance functions involved in the differents models via moment methods.
2. With the estimated covariance functions, calculate an appropriate number of K-auto-values and auto-vectors along the different levels of hierarchy that collect the different modes of variability in a precise way so that the problem we want to address. In the step 1, to estimate the covariance functions in the different models mentioned above, a general estimation strategy proposed in [Koch, 1968] can be used. For example, following the notation and problem (N 2) defined in the previous Section,Σ T ,Σ W ,Σ B , can be expressed with the following sandwich structure:Σ where X is a matrix of size (nJ) × m that records the different observations of all the individuals and levels of hierarchy, while G T , G W and G B are designspecific matrices of dimension m × m.
In particular, the usual co-variance matrixΣ T is written asΣ T = XG T X T , where G T = 1 nJ (I − 11 T ) where I denotes the identity matrix, and with 1, we denote m length vector with all ones. More details about these procedures, as well as about the selection of the components and the score estimation, can be found in the following references [Di et al., 2009, Shou et al., 2015.

Intra-class correlation coefficient (ICC)
A significant problem when several repeated measurements are collected from a subject over different days or other periods is to determine how much variability is explained by the subjects' effect and how much by making different measurements over different levels of the hierarchy. This issue is in the literature and is known as the process of estimating the coefficient of intra-class correlation (ICC) [Müller and Büttner, 1994] that pursues to estimate the variability of measuring a subject in conditions that are assumed to be standardized across different tests. The estimation of ICC is crucial, for example, in the field of clinical laboratory testing, where we want to use clinical variables for the monitoring and diagnosis of patients that are not modified abruptly between days by a problem of error of measurement of the device, and by the intra-day variability of individuals, see an example of the above in diabetes in [Selvin et al., 2007]. In biomechanics and exercise sciences, the ICC's quantification is also critical in searching for objective criterium to assess performance and control the individual's degree of fatigue [Van Gheluwe et al., 2002, Koldenhoven andHertel, 2018]. Although a variable may have a high variability, it can be a very useful criterion for decisionmaking. In this case, it is necessary to make several measurements to capture that variable accurately. The ICC can also quantify how many measures we have to make to capture variable distribution with enough accuracy. The first model where the ICC was estimated is (N 2) of the Table 1. In this scenario, we have Fixed t ∈ [0, 1], by analogy with a univariate non-functional case, the proportion of the total explained variable by the effect of the subjects at that point, is given by being ρ(t), the intra-class correlation coefficient in the point, t.
In a straightforward way, the ICC can be generalized as a global measure at the functional level, see for example ([Shou et al., 2013]) comparing the total variability collected by the involved co-variance operators with the variability modes K A , and K W and the global white noise (t) (covariance functions at the levels 1 and 2). Thus, the global ICC, denoted as ρ, is wherewith tr(·), we denote the trace operation, and where we are using the notation of the Section 2.2.2. We have to point out that the homoscedastic errorterm (t) has been included according to the convention followed in the Section 2.2.3-a more general setting that source of random variable is decomposed into an independent term.
The ICC can be calculated in more complex multi-level models. Suppose that we wish to use model (N3) and then, we have three levels. To develop such a task, it is enough to divide the source of variability generated by the hierarchy associated with subjects by all variability sources, that is: Recently, the intra-class correlation coefficient has been extended for objects that live in complex spaces where similarity between objects can be computed by the particular distance [Xu et al., 2020].

Hypothesis testing between different levels
Consider the model (N3) specified on Table 1 : where i = 1, 2, . . . , n; j = 1, 2, . . . , n i ; k = 1, 2, . . . , K ij . Without a loss of generality, suppose that the first level is the individual, the second is the test performed (HIIT or CTR), and the last level is the stride number. Comparisons between the differences in HIIT run and CTR are of considerable interest in biomechanical studies. To do this, we need to compare the difference between 2−levels effect functions In practice, fixed k, we can test univariate distribution changes with the estimated score of the second level composed of HIIT and CRT runs effects respectively, For this purpose, we can use the rich family test that provides energy distance methodology [Rizzo and Székely, 2016], or with classical tests such as Kolmogorov-Smirnov or Crammer-Von-Misses. As we applied univariate-test m times and obtained m marginal p-values (for each score), we must apply false discovery rate [Benjamini and Hochberg, 1995] or other criteriums to performer corrections for multiple comparisons to control type 1 error under the null hypothesis. Finally, we return as global p-value min{p * 1 , · · · , p * m }, where p * i denotes the adjusted p-value for score i. A similar methodology introduced here, was used in the standard set-up of hypothesis testing with functional data [Pomann et al., 2013] out of a multilevel data framework.

Summary of functional multi-level models. What is
the reason that this models are so important?
The increasing ability to store different profiles and functions of different variables that measure individuals' health from a broad spectrum of perspectives at different time scales provides several methodological challenges of statistical analysis that multilevel models can solve. In particular: • We can obtain a vectorial representation for each individual that captures the differences between individuals in a context of repeated and longitudinal measures.
• We can obtain the same representation for each individual in different hierarchical levels, for example, in a specific run and specific step recorded.
• For a specific individual, we can estimate the differences between different hierarchical levels. In addition, we can quantify intra and inter individuals' variability in all model levels. With this model, we can see under specific conditions, the specific modes of variability and compare with other conditions.
• We can obtain reliability measures as ICC or compare through hypothesis testing changes along a group of individuals or test conditions with paired and repeated measures. We can do this with the methodology previously established or following [Crainiceanu et al., 2012].

General description of study and variables
In order to assess biomechanical changes in typical training sessions in recreational runners on an equal level, 20 participants (10 women and 10 men) were initially selected to complete four typical training sessions sufficiently spaced in time. Two were high-intensity interval training, and the rest were continuous training. In the first case, athletes ran 6 × 800(m) intervals at 1 km/h below their maximum aerobic speed with 1 : 1 recovery. While in the second one, the athletes completed a continuous at a speed below maximum steady stater.
The duration of the continuous run was individualised to the same estimated energy expenditure as the interval training. In addition, the training sessions were conducted at the same time of day to avoid possible daytime fluctuations.
Running kinematics was measure with three-Dimensional motion analysis system that collect data at frequency of 500 Hz. All sessions were performed in an environmentally controlled laboratory setting, the athletes all used the same treadmill. Isometric strength was for various actions were recorded pre and post run.
The participants' basic characteristics can be found in Table 2. The strength changes in the last 800-m interval in Hip Abduction, Hip Adduction, Knee Extension are shown in Figure 1.
In this paper, we analyze 20 cycles of the stance phase for each run over 19 participants. For security reasons, we have excluded one of the participants, due to the presence of some outliers and missing data in some part of strides.
In our analysis, we have only focused on what happens in the three-dimensional knee segments: Knee-X, Knee-Y and Knee-Z. Additional details, about the study design and how the measurements were made can be found for example in [Riazati et al., 2020].

Aims of the analysis, clinical implications and statistical analysis
Knee injuries are the most common injury among runners of all levels. Therefore analyzing changes in posture and differences during typical training sessions has high clinical value in acquiring new epidemiological knowledge related to the causes of injuries. Examining the reliability between two interval training sessions is of fundamental importance. This will enable us to know how much information needs to be captured in order to characterize an athlete's biomechanical profile in a training session . In addition, changes in biomechanical patterns are very individualized and variable between individuals, therefore, using statistical tools that put the focus on the average subject in the study population rather than on individual variations can be very misleading, particularly in studies where the sample size is minimal. Finally, it is of interest to know if the information registered through the functional profiles and analyzed by the multilevel models provided information of the changes that occur during the athletes' training sessions.
In order to answer some of the questions mentioned above, we have divided the statistical analysis made with the three-dimensional information of the knee, in the following items: 1. Examine the correlation between the scores obtained after applying the functional component analysis with the force production changes in the training session.
2. To estimate the multilevel functional intraclass correlation coefficient to measure the reliability between two interval training sessions using the 20-step information.
3. To establish if statistically significant differences exist between a continuous and an interval training session by means of a hypothesis test exploiting the representation constructed from a functional multilevel model.
In particular, we have selected the three-level nested model (N3) from the Table 1 to carry out the model of all the previous issues.
In the analyses outlined, the results will be accompanied by graphs that help us understand and discuss how individualized the biomechanical changes are and how it is a useful hypothesis test to infer conclusions in this context. 14 Figures 3 and 4 contain the information about 20 strides per individual in different two HIIT runs in Knee-X, Knee-Y and Knee-Z. We can see that there are subjects in whom there are hardly any differences in their biomechanical profiles between the two runs. However, in others, differences seem to be present. In addition, we can also see that the patterns between the two runs are quite individual; no common pattern in angle values exists across the all the runners examined. Figure 2 show the bivariate association between each functional scores after applied multilevel principal component analysis and changes in the strength production between athletes. The results show that in some scores, there is a significant correlation against changes in strength production. However, in other cases, the correction is very poor. Notwithstanding this, we examined the marginal association and probability can be the interaction in more complex models between scores, but the limited sample size of this study remains fit more complex model.
The functional ICC for Knee-X is 0.55, Knee-Y 0.54, and Knee-Z 0.61. Likewise, at a significance level of 5%, no statistical differences were found between the biomechanical patterns during interval training and continuous running, with the p-values for Knee-X, Knee-Y, and Knee-Z respectively of 0.17, 0.12 and 0.4. Figures 5, 6 shows the measured curves of each group taken with the average of the 20 steps, along with the biomechanical profiles of two athletes we consider representative. They show that the biomechanical changes between an interval run and a continuous run among the athletes are very changeable. In some individuals, there is no biomechanical changes in running style, while in others, the biomechanic profiles are very different. However, the mean values are not significantly different.

Discussion
Knee injuries are one of the most frequent problems faced by recreational runners [Van Gent et al., 2007]. An accurate characterization of the biomechanical changes that occur in typical training sessions can be critical in identifying the etiology of injuries [Donoghue et al., 2008] and developing predictive models to detect injury risk [Ceyssens et al., 2019]. Here, we have illustrated how to exploit the functional information of different steps during different training sessions from multilevel models to: i) examine the correlation between knee angles and changes in force production in the same training session; ii) measure the reliability between two training sessions: iii) see that there are no statistically significant differences between a continuous run and an interval training with the same energy expenditure, although remarkable differences exist if we visually analyze some individuals.
The complete analysis of each cycle through functional analysis techniques that analyze the curve in its totality has lead to more nuanced findings [Donoghue et al., 2008]. Traditional techniques that analyze either fixed angles, the average angle, the range of movement or other measures summarized, result in the loss of information that its use entails. Complementary, interesting problems can be identified when using more informative gait points. Recent statistical methodologies can be used to address this problem [Berrendero et al., 2016, Poß et al., 2020 Functional multilevel models are an essential weapon in the challenge to exploit information from monitoring athletes or patients, to optimize decision making using different sources of information and measurements, made at different resolution levels. These tools can help integrate and analyze the information together, obtain a representation of the individuals along with different levels of hierarchy, and establish the different forms of variability in the different levels considered. These tools are remarkable if we want to analyze all training records or physiological variables of a group of athletes over a season or different micro-macro-cyc*les [Lambert andBorresen, 2010, Halson, 2014]. For example, there is not yet a sufficiently good methodology to represent the information inherently as proposed by these models [Matabuena and Rodríguez-López, 2019, Piatrikova et al., 2021, Kalkhoven et al., 2021. Despite being an exciting research topic with high relevance, we believe that there are not many methodologies to address relevant problems in biomechanics to date. For example, a specific need of this field could be to build a multilevel model that considers the different time length of the step, and not lose information of the step geometry with the standardization of all the strides to the [0, 1] interval.
The multilevel models have allowed us to calculate the intraclass correlation coefficient between the two interval training sessions taking into account the 20 steps recorded in each session. To the best of our knowledge, this is a novel approach in this area since the traditional approaches previously used to measure reliability rely on the compression of information in the average curve and only beetween two conditions [Pini et al., 2019]. At the same time, we have introduced a new hypothesis test to test the statistical differences between continuous and multilevel running, taking advantage of the representation we obtained with the multilevel model at the second level of the hierarchy. This also represents an advance, since with the inclusion of the 20 steps in the model in each test, we have more information, and with the new procedure, we can see if there are statistical differences between the different levels of hierarchy or groups of patients/athletes taking into account the differences in the study design.
An important aspect to consider in analyzing the results is that the individuals' movement patterns seem unique. This is not new, and several papers have exempted the individuality of human walking and running [Horst et al., 2019]. In this sense, since the biomechanical patterns are probably grouped in clusters [Phinyomark et al., 2015, Jauhiainen et al., 2020, standard hypothesis tests applied to the whole sample are not the best way to establish biomechanical differences. There are some discrepancies between studies when examining these issues. Also, in the biomechanics literature, as in other biomedical literature areas, there is some controversy about the use of p-value [Benjamin et al., 2018], and the use of other approaches such as effect size [Browne, 2010] or e-values [Vovk and Wang, 2019] may be recommended.
A limitation of this study is the sample size, together with the fact that we are analyzing the biomechanical variations of the knee, without taking into account the possible multivariate structure of knee movement. However, due to the reduced number of data, we can gain a greater interpretation in this type of study of a more exploratory character with this procedure. Moreover, this work's main objective is to illustrate the use of multilevel models with biomechanical data.
The rise of biosensors [Ferber et al., 2016, Phinyomark et al., 2018 in the area of biomechanics and medicine is causing an unprecedented revolution in the evaluation of athletes and patients care. It is likely that in the coming years, many of the clinical decisions will also be supported by the values predicted from the algorithms in many contexts, such as the prediction of injuries [Clermont et al., 2020, Van Hooren et al., 2020 or optimal surgery recovery [Karas et al., 2020, Kowalski et al., 2021 so in sport and general populations. Undoubtedly, the introduction of the data analysis techniques discussed here will help practitioners analyze objects that vary in a continuum repeatedly and that appear more and more frequently in biomedical data [Dunn et al., 2018].