Basis expansion approaches for functional analysis of variance with repeated measures

The methodological contribution in this paper is motivated by biomechanical studies where data characterizing human movement are waveform curves representing joint measures such as flexion angles, velocity, acceleration, and so on. In many cases the aim consists of detecting differences in gait patterns when several independent samples of subjects walk or run under different conditions (repeated measures). Classic kinematic studies often analyse discrete summaries of the sample curves discarding important information and providing biased results. As the sample data are obviously curves, a Functional Data Analysis approach is proposed to solve the problem of testing the equality of the mean curves of a functional variable observed on several independent groups under different treatments or time periods. A novel approach for Functional Analysis of Variance (FANOVA) for repeated measures that takes into account the complete curves is introduced. By assuming a basis expansion for each sample curve, two-way FANOVA problem is reduced to Multivariate ANOVA for the multivariate response of basis coefficients. Then, two different approaches for MANOVA with repeated measures are considered. Besides, an extensive simulation study is developed to check their performance. Finally, two applications with gait data are developed.


Introduction
The well known Analysis of Variance (ANOVA) methodology aims at comparing more than two groups and/or treatments with respect to an scalar response variable.This comparison is based on the total variability decomposition of an experiment in independent components that are attributed to different reasons.In general terms, it determines if the discrepancy between the averages of the treatments is greater than what would be expected within the treatments.Thus, ANOVA can be seen as the natural extension of two-sample classical statistical tests to the case of more than two populations.From its formulation, ANOVA has been constantly object of study to adapt it for different scenarios: single or multiple factors, parametric-non parametric cases, repeated measures frameworks or longitudinal data, among others.The more complex problem of testing the equality of a large number of populations has been recently studied in Jiménez-Gamero et al. (2022).Motivated by biomechanical studies where the aim is to test the differences between the means of the human movement curves observed on independent samples of subjects under different treatments, the current manuscript is focused on addressing this tool from a Functional Data Analysis (FDA) perspective.
FDA is a branch of statistics devoted to analyse the information of functions (usually curves) that evolve over time, space or another continuous argument.FDA is able to explore the curves in all the domain.This fact avoids the loss of important information usually produced when functional data are analysed through multivariate approaches from discrete measures.The great computational advances experimented by technology in last years make possible to model the complete curve and retain its main features.Key text books in FDA and related topics offer a broad vision of the most general methodologies, applications and computational aspects of this field, see e.g.Ramsay andSilverman (2002, 2005); Ramsay et al. (2009); Ferraty and Vieu (2006); Horvath and Kokoszka (2012).The majority of the classical multivariate statistical techniques have been extended for functional data: principal and independent component analysis (Aguilera and Aguilera- (Araki et al. 2009;Aguilera-Morillo and Aguilera 2020), among other recent papers.ANOVA problem for functional data (FANOVA) has also been considered in the literature.There are available different parametric (Cuevas et al. 2004; Cuesta-Albertos and Febrero-Bande 2010) and nonparametric (Delicado 2007;Hall and Van Keilegom 2007;Jiménez-Gamero and Franco-Pereira 2021) approaches to tackle the traditional m-sample problem.A deep review of the most important aspects for FANOVA problem is developed in the book by Zhang (2014).An interesting comparison of multiple tests based on the idea of B-spline tests (Shen and Faraway 2004) can be seen in Górecki and Smaga (2015).More recently, a novel approach based on functional principal component analysis was introduced in Aguilera et al. (2021a).The purpose is to reduce the dimension of the problem and conduct a multivariate ANOVA on the vector of the most explicative principal component scores.Additionally, several authors have also focused their efforts on providing some tools to deal with the multivariate ANOVA problem for functional data (Górecki and Smaga 2017;Acal et al. 2021).
Despite its noticeable interest in applications with real data, the repeated measures setting for functional data (FANOVA-RM) is barely considered.This theoretical framework deals with the situation where a functional response variable is observed under different conditions (also called treatments) or time periods for each subject.The most of the works available in the literature only concerns the paired sample layout.The first statistic to solve this problem was introduced by Martínez- Camblor and Corral (2011).This statistic is given by the integral of the difference between the sample mean functions, whose null distribution was approximated by considering multiple bootstrap and permu-tation methods.Additionally, Smaga (2019) developed a new way to approximate the distribution of the same statistic by means of Box-type approximation.The performance of these approximations was very similar in relation to the sample size and power, but the Box-type approximation proved to be faster from a computational viewpoint.However, this statistic only takes into account the between group variability.In order to solve the lack of information about the within group variability, two new statistics adapted from the classical paired t-test were introduced in Smaga (2020).These statistics were extended by considering the basis expansion of the sample curves with the aim of detecting changes in air pollution during the COVID-19 pandemic in Acal et al. (2021).In the more general context of testing homogeneity of paired functional data, a Cramér-von-Mises type test statistic is introduced in Ditzhaus and Gaigall (2021) with application to stock market returns.
The current work is focused on two-way FANOVA problem, in which the subjects are classified in independent groups and the response variable is observed under different conditions for each individual.Thus, one factor represents the repeated measures effect (treatments) and the second one denotes the group contribution.This is the case of one of the aplications developed in this paper (see Section 4.2) where the aim is to detect if there are significant differences in the kinetic curves (angle of knee) recorded under three different velocities (repeated measures) on two independent samples (age groups).A simple FANOVA model can be expressed as follows in terms of the global mean function, main-effect functions and i.i.d.errors: where x ijk (t) is the observed value of the response variable for the kth subject in the jth age group, measured under the ith running velocity at moment t in a continuous time interval T, (i = 1, 2, 3; j = 1, 2; k = 1, . . ., n j ) with n j being the sample size for each age group.
The main idea is to assume the basis expansion of the sample curves in order to turn the FANOVA-RM into a multivariate ANOVA with repeated measures, where the vector of basis coefficients of the sample curves would be the dependent multivariate variable.As far as we know, this theoretical setting has not been ever addressed in the literature, only in Martínez-Camblor and Corral (2011) it is briefly indicated how the tests could be generalized for the case of more than two samples, but without further details about the case of independent groups.
In addition to this introduction, the manuscript has four additional sections.The theoretical aspects related to the proposed methodology are in Section 2.
An extensive Monte Carlo simulation study to show the good performance of two considered FANOVA-RM approaches is developed in Section 3. The applications with biomechanics data that motivate this study are presented in Section 4. Finally, the most important conclusions from this research are summarized in Section 5.

Model set up
As natural extension of the multivariate case, the aim of Functional Analysis of Variance is to estimate the effect of one or more factors on a functional response variable.In this paper, two factors are considered (twoway FANOVA model) which usually represent groups and treatment conditions.Two different basis expansion approaches are developed in what follows.

Two-way FANOVA model
Let {x ijk (t) : i = 1, 2, . . ., m; j = 1, 2, . . ., g; k = 1, 2, . . ., n ij ; t ∈ T } denote g × m independent samples of curves defined on a continuous interval T. That is, x ijk (t) represents the response variable of the kth subject in the jth group under ith treatment at instant t.Note that each sample contains n ij observations and the total sample size is n = m i=1 g j=1 n ij .Sample curves can be seen as observations of i.i.d.stochastic processes (functional variables) {X ijk (t) : i = 1, 2, . . ., m; j = 1, 2, . . ., g; k = 1, 2, . . ., n ij ; t ∈ T } with distribution SP (µ ij (t), γ(t, s)), with µ ij (t) being the mean function and γ(t, s) being the common covariance function associated with each of the g × m stochastic processes.It is supposed that these stochastic processes are second order, continuous in quadratic mean and with sample functions in the Hilbert space L 2 [T ] of squared integrable functions with the usual inner product In Two-way FANOVA model, functional data verify the functional linear model given by where µ(t) is the overall mean function; α i (t) and β j (t) are the ith and jth main-effect functions of treatments and groups, respectively; and ǫ ijk (t) are i.i.d.errors with distribution SP (0, γ(s, t)) ∀i = 1, 2, . . ., m; j = 1, 2, . . ., g; k = 1, 2, . . ., n ij .This model is the generalization of the model proposed by Zhang (2014) when In model (1), possible interactions between groups and treatments are assumed to be ignorable.This means that the effect of a certain factor's level is the same for each level of any other factor.In our case the effect of the treatments on the functional response would be the same for all groups.The interaction functional parameter can be added in model (1) as follows It is well known that FANOVA model is not identifiable (the functional parameters are not uniquely defined).This problem is overcome by applying certain constraints.By extending the usual constraints in the balanced multivariate case (the cell sizes n ij are equal), the main effects and interaction effects functions sum to zero.An appropriate sequence of positive weights must be considered to define the constraints in the unbalanced design (Zhang 2014).From now on, we will assume the constraints (3) Under these constraints, we have that where µ i. (t) and µ .j(t) are the marginal mean functions of the functional response for each treatment and each group, respectively.The most interesting hypothesis tests associated with two-way FANOVA model are given by the following null hypotheses against the alternative, in each case, that its negation holds: -Testing if the main-effects of treatments are statistically significant (equality of the unknown treatments mean functions) -Testing if the main-effects of groups are statistically significant (equality of the unknown groups mean functions) -Testing if the main-effects of the groups and treatments are simultaneously null -Testing if there is a significant interaction-effect between groups and treatments Different test statistics have been proposed in literature for FANOVA testing problems.A detailed study of point-wise F-test, L 2 -norm-based test and functional F-type test under Gaussian assumption, together with χ 2 and bootstrap approaches for non Gaussian samples can be seen in Zhang (2014).Exhaustive simulation studies to compare multiple existing tests for one-way FANOVA testing are presented in Górecki and Smaga (2015).Different basis expansion and functional principal component approaches are proposed in Aguilera et al. (2021a) and Aguilera et al. (2021b) with applications in electronic and Google Trends, respectively.
Under the constraints defined in (3), the usual least squares estimators of the functional parameters in model ( 2) are obtained by minimizing

Estimation and computation from basis expansions
In practice, functional data are not continuously observed over time but only discrete observations are available for each sample curve.What is more, the number of observations and the location of observed time points could be different for each curve.Because of this issue, the first step in FDA is to reconstruct the original functional form of the data by using some functional projection approach.
Let us suppose that {φ h (t)} h=1,...,∞ is a basis of the functional space L 2 [T ] the curves belong to.Then, each curve admits an expansion into this basis as follows where the basis coefficients y ijkh are generated by random variables with finite variance.An approximated representation is usually obtained by truncating this basis expansion in terms of a number p of basis functions sufficiently large to assure an accurate representation of each curve.
From now on, it will be assumed that sample curves belong to the space generated by the basis {φ 1 (t), . . ., φ p (t)} .In vectorial form with y ijk = (y ijk1 , . . ., y ijkp ) ′ and Φ(t) = (φ 1 (t), . . ., φ p (t)) ′ .This way, the initial curves x ijk are replaced with their vectors of basis coefficients y ijk .The main advantage of this approach is that the dimension of the data depends only on the number of curves and on the order p of the expansion.Due to the fact that the curves are observed with error, some smoothing approach (e.g.least squares approximation) is usually performed to estimate the basis coefficients.In relation to the choice of a suitable basis, there are multiple options depending on the characteristics of the sample curves but the most common are Fourier, B-Splines or wavelets bases.Another key point is to select the dimension of the basis, which can be worked out, for instance, by generalized cross validation (Craven and Wahba 1978).
Hence, by assuming the basis expansion in (9), the estimators of the functional parameters can be expressed in terms of basis functions as follows where, x ijk (t) x ijk (t), i = 1, . . ., m; j = 1, . . ., g x ijk (t), j = 1, . . ., g with y ... , y i.. , y .j. and y ij.being the grand mean vector, the treatment mean vector, the group mean vector and the interaction mean vector, respectively, associated with the vector of sample curves basis coefficients y ijk .
The results above proves that Two-Way FANOVA model turns into Two-Way multivariate ANOVA model for the p-dimensional response variable (Y 1 , Y 2 , . . ., Y p ) that generates the basis coefficients of the functional variable X.

Repeated measures approaches
Two-Way FANOVA model presented above corresponds with independent samples.Nevertheless, in many fields such as medicine, social sciences, education or psychology, among others, it is very common to deal with a repeated measures design in which measurements on one or more response variables are conducted at several occasions (longitudinal data) or under different treatment conditions on the same subject.When a single response variable is observed, the design is called univariate repeated measures design.In this document, this approach is extended to test the effect of two factors (treatment and group) on a functional random variable.This new approach is known as Two-Way FANOVA-RM.Obviously, the simplest approach would be the model where the treatment factor is only considered.This model can be derived through the model with two factors.
Let us now suppose that we have a repeated measures design with g independent samples of curves (one per group) so that the response functional variable X is repeatedly measured on each subject at m different time periods (longitudinal functional data) or under m different treatment conditions.
Let {x ijk (t) : i = 1, 2, . . ., m; j = 1, 2, . . ., g; k = 1, 2, . . ., n j ; t ∈ T } denote g independent samples of curves defined on a continuous interval T. That is, x ijk (t) is the response of the kth subject in the jth group under the ith treatment.Note that the response of each subject is observed m times so that we have n = g j=1 n j subjects and n × m sample curves.It is assumed that each treatment is applied to all subjects (balanced design).This fact will be an essential aspect later.
The main objective of this manuscript is to adapt Two-Way FANOVA model to the case of repeated measures by taking the intra-subject variability into account.We propose to perform a multivariate analysis of variance with repeated measures on the multivariate response defined by the random basis coefficients of the functional variable.
As far as the authors know, there are two different models to include the intra-subject effect in the analysis: Doubly Multivariate Model (DMM) and Mixed Multivariate Model (MMM).Both assume the multivariate normality hypothesis and homogeneity of covariance matrices.The difference between them arises in the assumptions on the covariance matrix.DMM only assumes that the covariance matrix is positive definite, whereas MMM requires the multivariate sphericity condition.This restrictive condition is not verified in many real situations, so that DMM is more frequently used.However, if the sphericity condition is satisfied, MMM should be employed because it is more powerful (Bock 1975).Several reviews, new results and comparisons of both models by standing out their principal characteristics and behaviour on applications were developed in Timm (1980) and Boik (1988Boik ( , 1991).An application with data from an educational survey can be seen in Filiz (2003).Three different methods to solve the lack of variance homogeneity are studied in Lix and Lloyd (2007).Finally, a new statistic based on DMM is developed in Hirunkasi and Chongcharoen (2011) for the tricky scenario where the dimension of the response variable is greater than the number of observations.

Doubly Multivariate Model
In our functional data context, FANOVA-RM model can be written as a MANOVA-RM model for the basis coefficients of the sample curves as where the response Y is the matrix n × (p × m) whose rows contains the p-dimensional basis coefficients of the functional response variable X for the n subjects (distributed amongst g independent samples) examined under each of the m treatments (p × m dimensional response ordered within each row according to treatment an within treatment according to the basis coefficients).
Les us observe that in this model we have p×m response variables that represent the p-dimensional vector of basis coefficients for each treatment.In addition, X is the between group design matrix n × g and B is the unknown parameter matrix g × (p × m), with E being the error matrix n × (p × m) whose rows E i are i.i.d.
where I n is the identity matrix and The hypotheses tests of interest related with the statistical significance of treatment, group and interaction effects, i.e. ( 5), ( 6) and ( 8), respectively, can be expressed in terms of basis functions and formulated through the following general linear hypothesis where 0 is a matrix of zeros with appropriate order, G is a matrix g × s (rank s) which contains the coefficients for between group tests and T is a matrix m × q (rank q) which contains the coefficients for within treatments tests.The columns of the matrix G are composed by the coefficients of s estimable between group functions and the columns of the matrix T are the coefficients of q linear functions of the m treatments.Without loss of generality, it is assumed that T is chosen to be orthonormal T ′ T = I.Depending on the type of contrast and the objective of the study, matrices G and T will have different expressions.Deep studies related with these topics were developed in Timm (1980), Thomas (1983), Hand andTaylor (1987) andTimm (2002).
A DMM testing problem is worked out by means of the usual MANOVA statistics, e.g., Wilks's lambda (W), Lawley-Hotelling's trace (LH), Pillai's trace (P) or Roy's maximum root (R), associated with the following reduced qp-dimensional multivariate linear model These MANOVA statistics are based on the relation between the sum of squares and cross product matrices corresponding to error and hypothesis obtained by where B is the maximum likelihood estimator of B given by B = (X ′ X ) − X ′ Y, with (X ′ X ) − being any generalized inverse of X ′ X .In practice, W, LH, P and R are usually approximated by F-tests statistics through Rao's approximation (Rencher and Christensen 2012).Finally, it is important to keep in mind that DMM can only be used when n > p × m, since otherwise the matrix S e would be singular.
The multivariate mixed-effect model is a generalization of Scheffé's Univariate Mixed Model (Scheffé 1956).In our FANOVA-RM approach, the model for the p-dimensional response of basis coefficients can be expressed as where π k are i.i.d.subject effects with distribution N (0, Σ π ) and ǫ ijk are i.i.d.errors with distribution N (0, Σ ǫ ).In addition, the covariance matrix of the pdimensional response is given by Σ = Σ π + Σ ǫ because π k (t) and ǫ ijk (t) are mutually independent.MMM can be expressed in terms of the linear model for multivariate repeated measures defined in ( 12) by rearranging the data matrix Y as follows.
Let us denote by y i the ith row of the response matrix Y in model ( 12) that represents the pmdimensional vector of responses values for the ith sample subject.Then, the vector y i is rearranged to obtain a m × p matrix Y * i such that vec((Y * i ) ′ ) = y i .The vec()-operator stacks the columns of a matrix.Thus, the rows and columns of Y * i correspond with the treatments and the dependent variables, respectively.Considering this transformation, the rearranged response matrix for MMM analysis is If B and E are rearranged in the same way, a (g × m)×p matrix of unknown parameters and a (n×m)×p matrix of random errors are obtained.These rearranged matrices verify that vec(Y After this transformation, ( 12) and ( 15) can be written as Let us observe that the rows of the error matrix E * are normally distributed but not independent because the rows corresponding to the observation of the response variable on the same individual under the different treatments conditions could be correlated.Therefore, ( 16) is a mixed model that keep the influence of the individuals on the response variable in mind (individual random effects).Now, the general null hypothesis of interest is given by that is the same than ( 14).Then, the matrices corresponding to error and hypothesis for testing (17) are For a MMM to be valid it must be verified that the S e and S h matrices have to be independently distributed as Wishart matrices.Independence is derived from multivariate normality and homogeneity of the errors given in ( 13) but a new assumption on the covariance structure, called multivariate sphericity, is a necessary and suffcient condition for these matrices to be Wishart.Note that sphericity is a situation more general of the composed symmetry.A likelihood ratio test for checking the multivariate sphericity is derived in Thomas (1983) but the asymptotic distribution may be a poor approximation when the sample size is moderate.An approximation that solves the lack of power for moderate sample sizes is developed in Boik (1988) by applying Box's expansion of the characteristic function (Box 1949).For the univariate case, Box (1954) proposed a factor of correction with the goal of giving a solution when the sphericity is not verified.This method consisted of disrupting the degrees of freedom of F-statistic.In this line, Boik (1988) formulated an analogous approach for the multivariate case.
From a theoretical viewpoint, by considering the restricted data matrix Y(T ⊗ I p ) that examines q functions of the treatments, multivariate sphericity is a condition for the structure of the covariance matrix Ω of Y(T ⊗ I p ).Thus, Ω = Cov(Y(T ⊗ I p )) = (T ′ ⊗I p )Σ(T⊗I p ).In particular, this condition assumes that Ω = I q ⊗ Γ with Γ being a p × p positive definite matrix of covariances among the p response variables.The variation of the combinations of treatment levels is captured in the (q × p) × p diagonal blocks of Ω.In this sense, multivariate sphericity gives raise that all these blocks are identical and that the q subvectors of each row of restricted data matrix are independent.Thus, multivariate sphericity can be seen as a condition about the variation of the dissimilarities between treatment modalities.
Finally, let us observe that MMM only can be used when n × m > p, although this assumption is almost always fulfilled in practice.

Simulation study
In this section, an extensive Monte Carlo simulation study is carried out in order to test the performance of Two-Way FANOVA-RM methods proposed in previous section.Specifically, four different simulation studies are developed with the purpose of evaluating the behaviour of the tests by considering different type of errors, sample sizes and shapes of the functions that provide the curves.Sample curves are generated artificially according to the following functional mixed ANOVA model where α i (t) and β j (t) are the ith and jth main-effect functions of treatments and groups, respectively; θ ij (t) is the (i, j)th interaction-effect between treatments and groups; γ k sin(πt) represents the subject-effect and ǫ ijk (t) is the error function.
In the four studies, functional data have been generated in discretized versions x ijk (t r ) for r = 1, ..., 101 with t 1 , ..., t 101 being chosen equidistant in the interval [0,1].Least squares approximation in terms of a basis of cubic B-splines with 14 functions was employed in all cases in order to reconstruct the functional form of sample curves.Therefore, the sample of each of the three treatments (m = 3) is represented by a vector of 14 dependent variables (p = 14).Besides, it is considered n 1 = n 2 = n with n = 50, 100 to check the power of the tests when n > (≈) p × m and n ≫ p × m.Inspired by Durban et al. (2005), the random subjecteffect in all models is given by γ The tests were replicated 500 times for each of the scenarios to be specified below.Significance level was established as α = 0.05.Finally, Wilks' Lambda statistics was conducted both DMM and MMM for testing if the profiles for each variable are parallel and whether there are differences in treatments and in groups.For G and T from ( 14), the following matrices were employed

First scenario (M1)
Three different forms are assumed for the main-effect functions of treatments and groups, and two for the interaction-effect functions.Thus, 18 different models are obtained by making all possible combinations among them.The selection of these functions was inspired by Cuevas et al. (2004), Górecki and Smaga (2015) and Górecki and Smaga (2017) and they are given by M1.A1: The null hypothesis of interest holds for models with M1.A1, M1.B1 and M1.I1, whereas the opposite happens for the remainder.In order to evaluate the results, it is important to keep in mind that the interaction function M1.I2 changes for each level of treatment and group.The main differences between M1.A2 and M1.A3 (the same for M1.B2 and M1.B3) are that the main effects in M1.A2 and M1.B2 are quite separated so that the testing problem should be less harder.Finally, ǫ ijk (t r ) are i.i.d.random variables N (0, σ ǫ ) with σ ǫ = 0.10, 0.20, 0.40.The latter value, i.e. σ ǫ = 0.40, is introduced for checking the performance of the tests under extreme situations.Figure 1 displays the differences among main effects functions when the treatment, group and interaction functions are considered to be different.Table 1 shows the acceptance proportions for each scenario.The obtained outcomes can be summed up in the following commentaries: 1.For the cases where the three null hypothesis are true, the tests reach good results (the error rates are lower than 7% in all cases).Furthermore, for the rest of cases, the decision of the test when σ ǫ = 0.10 is really satisfactory even in the more extreme situations.2. MMM gets better results than DMM, being the differences in some occasions almost of 20%.In this context, as long as the MMM's conditions are satisfied is better to use this approach.3. The sample size and the dispersion parameter play an important role in the analysis.As the value of σ ǫ increases the error rate turn into larger, especially when the sample size is similar to the dimension of response variables.4. Regarding the check of differences among treatments, the tests provide adequate outcomes.It is only found out an error rate a little higher when jointly n = 50 and σ ǫ = 0.40 for M1.A3 (the acceptance proportions vary between 0.340-0.392and 0.244-0.260 in DMM and MMM, respectively).Remind that in M1.A3 the differences between the treatments functions are smaller.Secondly, in relation to the differences among groups, the performance of the tests is suitable except for σ ǫ = 0.40.In this case good results are appreciated only when n = 100 and the similarities between the groups are quite separated (M1.B2).It is also discovered a great error rate when n = 50, σ ǫ = 0.20 and M1.B3 were considered together.Finally, the power of the tests for the difference of the interactions between groups and treatment is too small when, on the one hand, σ ǫ = 0.40, and on the other hand, σ ǫ = 0.20 and n = 50 are considered at the same time.It is worrisome the frequency of a correct decision in these situations.

Second scenario (M2)
This second study (M2) is motivated to analyse the quality of the proposed methods when others settings for error functions are used.In particular, it is assumed that ǫ ijk (t) = 20 −1 B(t), where B(t) is a standard brownian process with dispersion parameter σ ǫ .M2 has been inspired by Martínez-Camblor and Corral (2011).The form for the functions of treatments, groups and interactions are the same than in M1.In Table 2 appears the achieved results.The conclusions done above about the treatments are maintained.However, the outcomes for the difference among groups and about the interactions are much better than in M1 when the corresponding H 0 is false.The error rate for σ ǫ = 0.20 does not exceed the 8% in any case.For its part, when jointly σ ǫ = 0.40 and n = 50 (because if n = 100 the behaviour of the tests is really good) the acceptance proportion varies between 0.184-0.254and 0.292-0.358 in MMM and DMM, respectively, for the case of parallelism.Likewise, for the effect of the groups, the proportion changes between 0.244-0.280and 0.440-0.494 in MMM and DMM, respectively, when M2.B3 was applied.Therefore, the type of error is another key point in this kind of analysis.

Third scenario (M3)
M3 is carried out under the same conditions than in M1 except the form of the interaction functions which were modified.The reason of contemplating this scenario is due to the fact that the obtained results for testing the hypothesis of parallelism were unsatisfactory in M1.Hence, with this scenario it is intended to study the impact of changing the form of the functions in the power of the tests.Consequently, the particular forms for interaction functions are the following M3.I1: θ ij (t) = sin(πt) 13 , i = 1, 2, 3; j = 1, 2; M3.I2: θ ij (t) = sin(πt) 21−2ij , i = 1, 2, 3; j = 1, 2.
Fig. 1: Main-effect (treatments and groups) and interaction-effect functions (first and second scenarios).
Table 3 contains the outcomes of this study.The conclusions about the hypothesis of the treatment and group effects are similar to that given in the first study.Nevertheless, we notice an important improvement in relation to the hypothesis of parallelism.In particular, when null hypothesis is false and σ ǫ = 0.40, the error rate is lower than 9.2% for n = 50 and being 0% for n = 100.Thus, this study lay bare another interesting factor that influences in the performance of the tests, i.e., the quality of the test relies on the shapes considered for the curves.

Fourth scenario (M4)
This study is worked out to corroborate the last affirmation made in previous section.For that purpose, M4 presents the same characteristics than M1 but interchanging the functions of the groups and interactions, that is, now we have M4.B1: β j (t) = [sin(2πt 2 )] 5 , j = 1, 2; M4.B2: The obtained results in M4 (see Table 4) confirm the suspicions about the importance of the curves form in the performance of the tests.The investigation brings excellent outcomes until σ ǫ = 0.20.It is only appreciated an error rate slightly high for testing the parallelism of the profiles in DMM when n = 50 and M4.I2 are considered at the same time (as maximum, the acceptance proportions reaches 10.8%).On the other hand, when σ ǫ = 0.40 two different behaviours are detected: 1.For all cases where M4.I1 is considered, the behaviour of the tests is really acceptable for the group effect (there are only little deviations when DMM is employed for n = 50 y M4.B3).Regarding the treatments, the remarks are similar to the rest of the previously applied models.Not a single problem was discovered for testing the hypothesis about interactions.As it has already been commented, the power of the tests is actually satisfactory when H 0 is true for any effect.2. For all cases where M4.I2 is considered, the results are deficient for checking the parallelism.Besides, there are some problems for testing differences among groups when n = 50 and M4.B1 are assumed, being the first time that it happens during the simulation.This should not occur because, although M4.B1 represents the case of no difference among groups, it is considered that the interaction depends on each level of the treatments and groups.
To sum up, this exhaustive simulation study displays sufficient evidences for concluding that this new methodologies, based on basis expansion of sample curves, is a suitable candidate to deal with FANOVA-RM problem.We are only slightly concerned about the variability of the power of the tests when the dispersion parameter is great.However, it is important to keep in mind that the subject effect also plays a fundamental role in the analysis and here it is assumed a high value for σ k .This fact produces important noise in the curves and the tests could convert into less conservative.Hence, this is another reason that endorses the goodness of the approaches presented in this paper, since even when the variability among subjects is large, the tests works very well in general terms.Finally, although the assumptions are satisfied by the own construction of the models, the normality, variance homogeneity and the sphericity (for MMM) are checked with a ratio of rejection lower than 5% of the cases.The developed simulation and the results shown in this section have been computationally implemented in R-cran.

Applications in Biomechanics
In this section, two different applications with real data associated with human body movement are developed.

Human activity data
The data of this application correspond with a research carried out by Anguita et al. (2013) where the human activity recognition over 30 subjects was analysed.In the current manuscript, the study is focused on the variable called linear acceleration (metre per second squared) which is measured on the axis X.The information about this variable is recorded under three different treatments (walking, walking upstairs and walking downstairs) in 128 equidistant knots at the interval [0,2.56].A subject was removed for being considered as outlier.This dataset was also used in Aguilera-Morillo and Aguilera (2020) for a functional linear discriminant analysis approach to classify a set of kinematic data.In the line of this work, sample curves were reconstructed through a cubic B-spline basis of dimension 27 with 25 equally spaced knots in the interval [0,2.56].The smoothed curves for each treatment are displayed in Figure 2 together with the sample mean function of each group (bottom-right).Based off this graph, it seems reasonable to think that there are significant differences among the three stimulus, i.e., the linear acceleration on the axis X is affected by the type of movement.This fact is numerically corroborated by means of the FANOVA-RM approaches considered in this paper.
FANOVA-RM analysis conducted in this manuscript is summarised as follows.On the one hand, there is a single available group (g = 1).Then, it is only possible to check whether there are differences between the treatments.On the other hand, DMM can not be considered because the sample space is smaller than the space of the variables (n < pm, being n = 29, p = 27 and m = 3).Therefore, MMM is employed in order to test the differences aforementioned.Finally, the p-value is calculated under the following conditions: 1. MMM is applied by assuming that the normality and the sphericity are satisfied.2. Due to the fact that the multivariate sphericity assumption is rejected (the likelihood ratio test provided a p-value lower than α = 0.05), an adjusted MMM by correcting the degrees of freedom of Fstatistic is performed.3. Given that the normality assumption is also in question, the permutation testing procedure proposed by Górecki and Smaga (2015) is adapted for the repeated measures scenario.The steps of this procedure are described below: (A) Calculate the value of the test statistic S 0 for the original sample data.(B) For each subject, it is necessary to permute randomly without replacement its observed values on the treatments.If there were more than one group, once permuted the values of all subjects, the following step would be to join all subjects in 'a single group'.Later, choose randomly without replacement n 1 observations for the first new then from the remainder of the observations draw randomly without replacement n 2 observations for the second new group, and repeat this process up to complete the g groups.
The p-values obtained after applying the different methods described above figure in Table 5.As a result, we can conclude that there are significant differences among the three types of stimulus on the linear acceleration on the axis X.

Dataset of running biomechanics
The public dataset available in Fukuchi et al. (2017) contains the biomechanical information of 28 regular runners.Concretely, lower-extremity kinematics and kinetics were registered meanwhile the subjects ran at different velocities (2.5 m/s, 3.5 m/s and 4.5 m/s) on an instrumental treadmill.Other relevant data such as demographics information, running-training characteristics or previous injuries were also collected.The current application is focused on analysing the angle (in degrees) of right knee on the axis X, which has been recorded over 101 time points.In particular, the aim is to detect if there are significant differences in this functional variable among the different velocities (repeated measures) according to the age.The variable age has been discretized in two independent groups.The first group is formed by 14 runners with less than 35 years old and the second one by the rest (≥ 35 years).The functional reconstruction of the curves by means of a cubic B-spline basis of dimension 18 can be seen in Figure 3.The mean curves of each age group under the different conditions are displayed in Figure 4. We observe certain differences in the angle regarding the velocity, being greater for the velocity equals to 2m/s.Nevertheless, it is less clear whether the angle depends on the age group; the angle of the runners older than 35 year old is higher in all the domain, especially at the end of the cycle, but the discrepancies are not very notice-   able.In order to validate statistically these assertions, the new methodology is applied.
In particular, this FANOVA-RM analysis contains the following characteristics: there are two independent groups of runners classified by the age, that is, g = 2.
The number of treatments is m = 3, whereas p = 18 Parallelism Differences-groups Differences-treatments 0.3852 0.7685 0.0020 Table 6: P-values after applying the MMM approach for FANOVA-RM by means of permutation testing procedure.
because of the dimension of the functional reconstruction.In these conditions, only MMM can be considered.
As the model assumptions are not verified, the permutation testing procedure aforementioned will be applied making use of the Pillai's trace statistic which is more robust than the other statistics in relation to the violation of model assumptions (Olson 1974).The results of this analysis can be seen in Table 6.The results show that the effect of the treatments does not depends on the levels of the age (there is not interaction).Additionally, it is corroborated that the differences between the age group are not significant, while the running speed plays an important role in the angle of the right knee on the axis X.

Conclusions
Functional analysis of variance with repeated measures aims at checking if the mean functions of a functional response variable observed in different time periods or treatments are equal or not.In spite of its great interest in practice, only a few works related to this topic are available in the literature.The current manuscript is focused on addressing Two-Way FAVOVA-RM problem.The first factor represents the multiple levels in which the functional variable is observed (repeated measures), while the second one constitutes the independent groups in which the subjects of the sample are distributed (independent measures).Under this scenario, it is necessary to study both the between-group and intra-subject variability.As far as we know, this theoretical setting has not been ever dealt yet from a functional data analysis viewpoint.Hence, a new approach based on basis expansion of sample curves is introduced in order to solve this problem.In particular, we prove that Two-Way FANOVA-RM model turns into Two-Way multivariate ANOVA-RM model for the multivariate response defined by the basis coefficients of the functional variable.In this point, mixed multivariate model or doubly multivariate model can be conducted to take the intra-subject variability into account in the analysis.An extensive simulation study has shown that the multivariate mixed model approach has a better performance than the doubly multivariate model, although both approaches provide good results in general terms.Only in extreme situations (small differences among functions, small sample size or great dispersion) the tests become conservative.The new methodology has also been applied to two real biomechanical datasets.In the first application, the objective is to evaluate how three type of stimulus affect on the linear acceleration of human movement, whereas the second study is focused on analysing the influence of age and speed in the knee flexion angle while running.

Fig. 2 :
Fig. 2: Sample group mean functions (bottom-right) and all the B-spline smoothed registered curves for each stimulus.
(C) Compute the value of the test statistic for the new sample generated in previous step.(D) Repeat steps (B)-(C) F times, being F a large number.Each achieved value in (C) will be denoted by S f with f = 1, . . ., F .

Fig. 3 :
Fig. 3: B-spline smoothed registered curves for each velocity according to the age group.

Fig. 4 :
Fig. 4: Sample mean functions according to the age group and velocity.

Table 1 :
Observed acceptance proportions for each scenario at a significance level 0.05 in the case of simulation study M1.

Table 2 :
Observed acceptance proportions for each scenario at a significance level 0.05 in the case of simulation study M2.

Table 3 :
Observed acceptance proportions for each scenario at a significance level 0.05 in the case of simulation study M3.

Table 4 :
Observed acceptance proportions for each scenario at a significance level 0.05 in the case of simulation study M4.