{'fraction_non_alphanumeric': 0.04914899835870225, 'fraction_numerical': 0.017334802257346495, 'mean_word_length': 3.977394807520143, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
PROLIFIC: Projection-based Test for Lack of Importance of Smooth Functional Effect in Crossover Design

Wearable devices for continuous monitoring of electronic health increased attention due to their richness in information. Often, inference is drawn from features that quantify some summary of the data, leading to a loss of information that could be useful when one utilizes the functional nature of the response. When functional trajectories are observed repeated over time, it is termed longitudinal functional data. This work is motivated by the interest to assess the efficacy of a noninflammatory medication, meloxicam, on the daily activity levels of household cats with a pre-existing condition of osteoarthritis under a crossover design. These activity profiles are recorded at a minute level by accelerometer over the entire study period. To this aspect, we propose an orthogonal projection-based test pseudo generalized F test for significance of the functional treatment effect under a functional additive crossover model after adjusting for the carryover effect and other baseline covariates. Under mild conditions, we derive the asymptotic null distribution of the test statistic when the projection function for the underlying Hilbert space is estimated from the data. In finite sample numerical studies, the proposed test maintains the size, is powerful to detect the significance of the smooth effect of meloxicam, and is very efficient compared to bootstrap-based alternatives.


Introduction
The prospect of monitoring health electronically has led to rapid use of wearable devices that are capable of collecting abundant amount of health related information continuously over time. Accelerometer is one of the most popular wearable devices that can objectively measure the physical activity (PA) count as densely as at a minute level (Bussmann et al., 2001). Most literature on this topic, summarize the massive activity data using various summary measures in pursuit of explaining the association between activities and health outcomes (Reider et al., 2020). However, such summaries-based approaches completely eliminate the variation of the PA over time. Functional data analysis (FDA)-based approaches view the daily PA counts as realization of some latent stochastic process and focus on modeling, prediction and studying the association between response with various important covariates such as age, gender, etc. Modelling daily PA involves handling complex functional dependence that inherently exists in the data. When the daily PA activities are recorded over multiple days, Goldsmith et al. (2015) proposed a multilevel functional data method to analyze PA profiles, which was extended by Xiao et al. (2015) to account for subject specific covariates. However, in many applications response trajectories are not observed daily for each subject, and the days over which they are observed are different for each subject. In our motivating study, household cats were subjected to a placebo-controlled four period crossover design (Ratkowsky et al., 1992) with intermediate washout period where daily PA profiles were measured at every minute level for each subject through accelerometer over 12 weeks (Table 1) [ Table 1 about here.] The objective of the study is to formally test the efficacy of an active drug meloxicam on the cats with degenerative joint disease (DJD) in the form of an improved daily physical activity counts. Let {τ (·, d)} be the mean response that is specific to the meloxicam treatment during the dth day, since the beginning of the treatment period. We want to formally assess whether, H 0 : τ (·, d) = 0 ∀ d vs H A : τ (·, d) = 0 for some d.
Most often, especially in a two period two treatment crossover design (chapter 2 of Jones and Kenward, 2014), the effect of treatment at period 2 is confounded by the residual effect of the treatments applied in period 1. This residual effect is termed as "carryover" effect of treatment (Cochran et al., 1941). Efficient estimation of direct treatment effect after removing the inherent high between-subject variability, even in the presence of significant carryover, can be done by choosing a suitable crossover design with more than two periods and/or more than two groups, which is the case in our meloxicam study. A significant amount of research has been done in the late twentieth century to address the estimation of direct treatment effect in presence of carryover (Hills and Armitage, 1979).
The literature of testing procedures for the bivariate mean function or bivariate smooth effect of a predictor involved in the context of LFD is sparse. Much of the existing methodologies on LFD concentrates on modelling through functional principal component analysis (FPCA) (Park and Staicu, 2015;Scheffler et al., 2020). To the best of our knowledge, the only functional testing procedure developed for LFD tests for invariance of the smooth bivariate mean µ(s, d) along the longitudinal component d (Park et al., 2018;Koner et al., 2021). In this paper, we develop a pseudo Generalized F (pGF) test for the significance of smooth bivariate effect of treatment on the PA in a hierarchically structured longitudinal functional crossover design. Our methodology relies on projecting τ (·, d) onto a set of orthogonal basis functions {φ k (·)} k 1 and testing the dependence of coefficient functions over d by extending the generalized F test developed by Wang and Chen (2012) in a more complex dependence structure. This allows us to transform the global null hypothesis into a number of simpler hypotheses along the longitudinal direction d. Compared to the pseudo likelihood ratio test (pLRT) developed in Koner et al. (2021), our pGF test is flexible, in the sense that it allows testing for the significance of a smooth bivariate effect in the presence other smooth effects in the model. This is important in our crossover design because we have to test for direct treatment effect in the presence of overall mean effect and the confounding carryover effect.
Testing for the direct effect of the treatment in the presence of carryover effect has its own challenges. Denoting λ(·, d) the carryover effect of the treatment that influences the response in the dth day of the washout period. Because the carryover is the residual effect of treatment, under the H 0 , τ (·, d) = 0 implies λ(·, d) = 0. On the other hand, under H A the objective is to test for the direct effect of treatment, not the carryover λ(·, d). This makes the hypotheses non-trivial in the sense that the H 0 and H 1 combined do not span the entire parameter space. As a result, a generalized F based statistic considering the null model where both the treatment and the carryover effect are zero, and the full model where both of them are present, will fail to maintain the empirical type 1 error. To mitigate this problem, we propose a two-stage procedure to carryout the test, where we first test for significance for the carryover effect in the presence of treatment, and in the second stage, we test for treatment effect under a model where carryover effect is present or absent depending on the conclusion of the test of carryover conducted at the first stage. The two-stage procedure maintains the type 1 error, has an excellent power in small samples and computationally efficient.
The rest of the article is organized as follows. In Section 2 we formulate the problem and introduce the model framework. The testing procedure along with the theoretical results related to the asymptotic null distribution are described in Section 3. Numerical studies are presented in Section 4 to demonstrate the finite sample performance of the test. Section 5 summarizes the findings on the efficacy of the active treatment meloxicam on physical activity, based on the conclusion drawn from the test. Assumptions related to the main theorem of the paper are in Appendix. Detailed proofs of the theorems as well as additional results related to the real data applications are provided in the supplementary material.

Model Framework
where Y ipj (·) is one-dimensional response trajectory (PA) observed for ith subject at the jth day during period p, d ipj for j = 1, . . . , m ip , p = 1, . . . , 4, along with a set of subjectspecific baseline covariates C i , = 1, . . . , L. The response curve is observed over a fine grid s ipjr , r = 1, . . . , R ipj with R ipj large. We assume that for every i, p and j, the set {s ipjr : r = 1, . . . , R ipj } is dense in compact set S. Without loss of generality, we assume that R ipj = R and s ipjr = s r for all r, i, p and j and use the index s instead of s r to denote a typical observation in the entire trajectory. It is assumed that the number of curves in each period, m ip , is small for each i but the collection {d ipj : j = 1, . . . , m ip , p = 1, . . . , 4, i = 1, . . . , n}, over all the subjects is dense in a compact set D. Let g i ∈ {1, 2} be the group identifier for each subject, i.e. g i = 1 if the subject is in group 1 and g i = 2 otherwise. As depicted in Table 1, in a crossover design, the treatment regime applied at a period p is identified by the group indicator g i . Let I ip,τ be the variable indicating whether the active drug meloxicam is applied on subject i at period p, i.e. I ip,τ = 1 if g i = 1, p = 1 or g i = 2, p = 3 and I ip,τ = 0 otherwise. Similarly, let I ip,λ be the indicator variable for the carryover effect in the washout period i.e., I ip,λ = 1 only if g i = 1, p = 2 or g i = 2, p = 4. We model the response trajectory using a functional additive crossover model (FACM) as, where µ(·, d) is the population mean, τ (·, d) is the direct effect of the treatment, and λ(·, d) is the carryover effect for the day d since the beginning of the period. Additionally, β (s) quantifies the smooth effect of the baseline covariate on the response curve. It is assumed that all these population level effects are smooth uni/bivariate functions. Finally, i (s, d ipj ) is a mean zero random deviation independent and identically distributed across all subjects i. The error process is meant to capture the variability in the response trajectory along with the variation in the response across different days in the period and the measurement error.
Under the FACM in (1) the hypothesis for significance of treatment effect translates to, Under the null hypothesis of no treatment effect, the carryover effect is constrained to be zero, whereas under the alternative of significant treatment effect, the carryover effect can either be zero or non-zero. Testing problem of this kind is atypical in the literature as the nuisance parameter λ(s, d) is dependent on the actual parameter of interest τ (s, d). Since the carryover effect size is informative of the treatment effect, we propose to pursue it sequentially, by first testing for carryover and then for treatment effect, as we will describe in the next section.

Alternate formulation of original hypothesis
Let {φ k (s) : s ∈ S} k 1 be a set of orthonormal basis functions in L 2 (S), with S φ k (s)φ k (s)ds = I(k = k ) for k, k 1. Then, the continuous function τ (s, d) can be represented uniquely Using the same set of orthogonal basis, we can expand the carryover effect as λ(s, d) Then, the original null hypothesis in (2) is equivalent to testing a series of simpler hypotheses, We discuss the procedure in detail in Section 3. Prior to that, we set up the framework for testing H 01,k and H 02,k in the FACM (1) below.

Testing framework under the projected model
Given the set orthogonal basis functions {φ k (s)}, consider the projected data, Y ipj,k = S Y ipj (s)φ k (s)ds, for k 1. The integral can be computed numerically with a very high precision since the grid at which the functional trajectories are observed is dense. The FACM in (1) for the projected response transforms to, where the components in the projected model are obtained by projecting each term of the original model onto φ k (·), i.e. µ k (d) := S µ(s, d)φ k (s)ds, β ,k = S β (s)φ k (s)ds, and In the context of the meloxicam study, the response trajectories are observed over 20 days, one can approach the hypothesis problem in (3) as testing of 20-dimensional vector (τ k (d 1 ), . . . , τ k (d 20 )) = 0. However, when the time points at which trajectories are observed are very sparse and different for each subject so that the set D can not be construed as a finite set with small number of elements, or when the variation of τ k (d) is smooth over time, standard analysis of variance (ANOVA) type approaches are not powerful.
Hypothesis testing of smooth effects carried out by a mixed model as above has been discussed in nonparameteric regression literature; Crainiceanu and Ruppert (2004)  and Chen assumed that the error e k are independently distributed across all i, p, and j. In our setup, the covariance structure of e i,k in (5) is non-trivial. Assuming a completely unknown dependence structure, let σ 2 k Σ ik be the m i × m i covariance matrix of e i,k in the linear where γ k (·, ·) is a continuous covariance function and σ 2 k > 0. Under the independence of the subjects, Cov(e k ) := σ 2 k Σ k = σ 2 k diag(Σ 1,k , . . . , Σ n,k ). Defining, π k := σ 2 µ,k /σ 2 k , The above testing framework requires a specified set of orthogonal basis system {φ k (s)} k 1 for the space L 2 (S) to compute the projected response and test H 0k . We take the eigenfunctions {φ k (s)} from the spectral decomposition of the so-called "marginal covariance" Ξ(s, s ) of the error process i (s, d ipj ), as our choice of orthogonal bases. A detailed description on the choice of orthogonal bases is provided in Section S8 of the supplement. However, the projected response Y ipj,k in model (4) is unobserved since the true eigenfunctions {φ k (·)} k 1 of the marginal covariance Ξ(s, s ) are unknown. However, we can compute the "quasi projections" W ipj,k := S Y ipj (s) φ k (s)ds as a proxy to the unobserved Y ipj.k . The rate of accuracy in the estimation of the eigenfunction φ k (·) ensures that the quasi projections W ipj,k are sufficiently close to Y ipj,k . For each k = 1, 2, . . . , K, stack the W ipj,k s for each subject i to construct to be the subset of X k by removing the columns corresponding to α λ,k . Similarly, define X −τ,k := ( X b,k , X λ,k ) and the same for Z −λ,k and Z −τ,k . Denote r := L + h µ + h τ + h λ + 3 as the rank of X k and P k := X k X k X k −1 X k as the projection matrix onto the column space to be the generalized projection matrix onto the column space of X k under the full model (5). With this setup, we layout the two stages of our testing procedure below.
Stage 1: Testing for carryover under the full model First, we test the significance of the carryover under model (5). The residual sum of squares (RSS) for the full model (5) using the quasi-projections is qRSS k (π k , η k , γ k ) : Under the null hypothesis H 01,k , α λ,k = 0 and γ k = 0.
The RSS under H 01,k simplifies to qRSS S1 0,k (π k , η k ) : The pseudo quasi GF (pqGF) statistic for testing H 01,k can be constructed as, where π k , η k , and γ k are estimated via restricted maximum likelihood (REML) under the full model (5). We call it as a pqGF statistic because it is constructed via quasi projections and the components of the model (5) are scaled by a pseudo estimator of the true covariance.
The test statistic has a similar form to what Wang and Chen (2012) considered. However, the fundamental difference between the RSS of Wang and Chen (2012) and qRSS k is that the latter is constructed through the quasi-projections W ipj,k , not Y ipj,k . This makes the derivation of the null distribution significantly more challenging as W ipj,k are no longer independent across i, because {φ k (·)} k 1 are obtained from the full data. Nonetheless, the next theorem states that if the eigenfunctions are estimated with high accuracy, then the null distribution changes by a minimal amount.
and suppose that Assumptions 1-2 hold for the true model (1). Assume that sup s∈S | φ k (s) − φ k (s)| → 0 in probability as n → ∞ and the Assumptions 3-4 hold for the elements of the projected model in (5) Then, for every k = 1, . . . , K, under the null hypothesis (2) the test statistic has an approx- where u s iid ∼ N(0, 1) and independently distributed with χ 2 h λ +1 and χ 2 N −r−Q λ . The quantities π k , η k and γ k are the minimizer of the spectral decomposition of the negative log-profile restricted likelihood under the alternative, A detailed proof is provided in Section S12 of the supplement. The asymptotic null distribution in (7) (2019), which can be directly applied to our case. In that sense, uniform convergence of eigenfunctions is viable in the context of our model (1).
Fixing an α ∈ (0, 1), let pqGF S1 ∞,k,α be the 100(1 − α)% percentile of the distribution of the random variable on the right hand side of (7). Then a α-level test for the null hypothesis H 01,k has the rejection region, where Assume that H 01,k was tested at Stage 1 and the decision was to reject. Evidenced by the significance of the carryover effect, using the full model (5) construct the RSS under H 02,k as qRSS S2a under H 02,k . Then, the pqGF statistic for testing H 02,k effect under the full model (5) can be constructed as, pqGF S2a where π k , η k , and γ k are estimated via REML under the full model. Fixing an α ∈ (0, 1), let pqGF S2a ∞,k,α be the 100(1 − α)% percentile of the distribution of the random variable on the right hand side of (14) in Section S9 of the supplementary material. A α-level test for testing the null hypothesis H 02,k under the model (5) has rejection region Stage 2b: Testing for treatment in the absence of carryover Consider the situation that H 01,k was tested and the results indicated lack of evidence of carryover effect. Then we test H 02,k under a simpler model of (5), omitting the terms due to the carryover effect as, The RSS under the model (10) is qRSS S1 0,k (π k , η k ), defined in the context of Stage 1. Additionally when H 02,k is true, the RSS simplifies to qRSS S2b and η k = 0. Then the pqGF statistic for testing H 02,k under the model (10) can be constructed as pqGF S2b N,k := N {qRSS S2b 0,k ( π k ) − qRSS S1 0,k ( π k , η k )}/qRSS S1 0,k ( π k , η k )s. Fixing an α ∈ (0, 1), let pqGF S2b ∞,k,α be the 100(1 − α)% percentile of the distribution of the random variable on the right hand side of (15) in Section S9 of the supplementary material. A α-level test for the null hypothesis H 02,k under the model (10) has rejection region

Two-stage test rule
We are now ready to present the proposed two-stage test along direction k. Fix a level of significance α 1 ∈ (0, 1) for testing the carryover at Stage 1. For any α ∈ (0, 1), a level α-test for testing the null hypothesis H 0k has the rejection region, The above test along the direction of φ k (·) is the key ingredient of the PROLIFIC, presented in (13).
The proof follows by the fact that both the tests in (9) and (11) reject the null hypothesis H 0 with a probability at most α when H 0 is true. It is important that both the tests at the second stage for the treatment effect are conducted at the same level, α, in order to ensure that the overall two-stage procedure has an overall size α. Moreover, the size of the overall test in (12) is independent of the level (α 1 ) at which the carryover is tested at Stage 1. This is because under the Gaussian data generating assumption the contrasts contributing to the direct treatment effect and that for the carryover effect are asymptotically independent. Put it differently, the inference drawn from test of carryover effect at Stage 1 does not influence the conclusion from the tests at Stage 2a and 2b, which is also evidenced by the numerical results presented in Section 4. This is the main reason behind the widely discussed criticism against the two-stage procedure for the typical AB/BA design (Senn, 2002, chapter 3) as the two tests at the Stage 1 and the Stage 2 are not independent (Freeman, 1989).
Under the null hypothesis H 02,k , both the tests with rejection regions R S2a α,k and R S2b α,k are equally capable of making a correct decision about the significance of treatment effect, up to an error level α. Thus testing for the carryover effect at the first stage does not have any direct implication when H 0 is true. However, the impact of testing for the carryover at the first stage will be profound on the power to detect departure from null. To justify this, suppose that in truth, both the treatment and the carryover effects are significant. Then the power to detect the significant treatment effect using the test in (11) will be much smaller than using solely of (9), as the later is conducted on a correct model. Therefore, the test for the carryover effect at the first stage provides a tool against possible model misspecification, while testing for treatment effect.
Finally, our projection based test, PROLIFIC, formally assesses the global null hypothesis H 0 by simultaneous testing of H 0k along the K directions and using a Bonferroni multiple testing correction to control for the family-wise error rate. Fix the nominal level α 1 ∈ (0, 1) for the test of carryover at Stage 1. Then, for every K 1 and pre-specified nominal level α, the rejection region for PROLIFIC is, Corollary 2: Assume the setup and the conditions of the Theorem 1. Furthermore, assume that the null hypothesis in (2) is true. Then, for any K 1, for every α, α 1 ∈ (0, 1).
The choice of the truncation parameter K does not affect the size the PROLIFIC. However, it affects the power to detect departure from null. In an hypothetical situation, when the treatment effect τ (s, d) is not different zero along the direction of the eigenfunctions {φ k (s)} K k=1 , but it is significantly different zero along a direction φ j (s) for some j > K, then the test does not have any power. On the other hand, choosing a large value of K will make the level for the individual hypothesis testing very small, α/K and α 1 /K, leading to a loss of power. In numerical results, we see that pre-specifying the percentage of variation explained (PVE) to 90%, PROLIFIC has desirable size and strong power performance.
The choice of nominal level α 1 for the test for the carryover at Stage 1 should be determined based on the implication of finding a significant carryover effect. A test for the carryover effect may be important to assess the usefulness of the washout period. A very small value of α 1 will lead to poor identification of the carryover, and as a result we might end up testing for the direct treatment effect in a wrong model. In practice, we recommend choosing a slightly higher value α 1 (say 10%) compared to α (say 5%). We conclude this section by describing the steps associated to implement PROLIFIC in Algorithm 3.1.
1 Construct a smooth estimator of the components of mean model as µ Obtain a smooth estimator of Ξ(s, s ) using the demeaned responses; Compute the p-value p k of test in (12) with a specified level of significance α 1 /K for the test of carryover at Stage 1; 8 end 9 Reject H 0 in (2) if min{p k : k = 1, . . . , K} < α/K at some level of significance α; Algorithm 3.1: PROLIFIC

Data generation
To assess the performance of PROLIFIC, we generate synthetic data for sample size n varying from 100 to 300. As described, we consider a crossover design with 4 periods and within each period the response profiles are observed sparsely over m ip time points. The number of profiles in each period, m ip , is generated randomly from {8, 9, . . . , 12} (low sparsity level).
are observed over a dense grid R = 101 points equally spaced over S = [0, 1]. With the above simulation design, the data is generated from the model The residual term in the model is generated as is mean zero subject specific random deviation that influences the response trajectories at every time point d ipj , along with a smooth random variation ε sm ipj (s) that is presumed to capture the additional variability at that specific time point and white noise process ε wn ipj (s). The random components of the model are generated from the following mechanism: On the other hand, the parameter γ controls the magnitude of carryover relative to the treatment. Setting δ > 0 and γ = 0, we can enforce absence of carryover even when the treatment is significant. In our simulation study, we take the shape parameters of the beta density as a = 2, b = 4. Further computational details are in Section S10 of the supplement.

Remarks on competing methods
One can adopt other procedures suitable to test for significance of unknown smooth function in univariate functional data and apply them to test for the projection of treatment effect on φ k (s), i.e. to test for H 02,k : τ k (d) = 0 under the projected model (4). In this regard,  Remarkably, the test fails to maintain the nominal level by large margin, at least for sample size up to 100. It is possible that a sample of size is n = 100 is not large enough to fairly approximate the asymptotic null distribution. On the other hand, when the null distribution is approximated by bootstrap, the test exhibits a rather conservative type 1 error.
Power: Fix the level of significance α = 0.05. The empirical power of PROLIFIC is plotted as a function of δ for small and medium sample size in Figure 1, based on 1000 simulations.
The left column pertains to the situation when both the carryover and the treatment effect are significant, and the right column when carryover is absent and only the direct treatment effect is significant. We do not present the power of the Ad-ZC when the null distribution is approximated by mixture of chi-squares, because it fails to maintain the size. As expected, in both the cases, we see that the power of the test increases rapidly with the increment in the sample size and as δ goes away from zero. The upper panel demonstrates that PROLIFIC is more powerful than Ad-ZC to detect departure from null, irrespective whether carryover is zero or not. The overlapping plots in the lower panel illustrate that the power of PROLIFIC is not affected by the level (α 1 ) at which the carryover is tested at the first stage.
The strength of PROLIFIC in detecting a very slight departure from null even with very small sample size can be attributed to the fact that both the contrasts for the treatment and the carryover are estimable after removing the variation due to subject and that for every subject, we observe the functional observations over four periods. Overall, the numerical results testify for the effectiveness of the two-stage procedure to detect the significant direct treatment effect in a crossover design, when both the treatment and carryover contrasts are separately estimable, in contrast to the widely criticized lack of power of the two-stage procedure in the case of AB/BA crossover design (Senn, 2002, chapter 3).

Meloxicam study of cats with osteoarthritis
The data originates from the meloxicam study of 58 household cats with existing condition of osteoarthritis. These cats were enrolled in a completely randomized double masked placebocontrolled crossover trial conducted at the College of Veterinary Medicine of North Carolina State University. The subjects were randomized into two groups. As described in Table 1, the first group received the single dose of active drug meloxicam for the 20 days in the first period, followed by placebo during the last three periods. Whereas, group 2 received the drug at period 3 and received placebo at all the remaining three periods. The objective of study is to understand the efficacy of an active drug meloxicam on the joint pain as reflected by an improved PA counts, measured at every minute level during the day by an activity monitor. See Gruen et al. (2015) for a complete details of the study.

Data preprocessing
Figure 3(a) presents the daily raw activity counts recorded by Accelerometer for a randomly selected cat over 5 days in every period. Since the cats in general stays in a resting state for a long period of time, followed by a sudden jump due to some external factors, the raw activity profiles are condensed by a lot of zeros between two high peaks. To reduce the large scale of variation in the activity counts, we add them by 1 and take the logarithm, i.e. x → log(x + 1). Let us call these log transformed PA counts as logPA. As a part of the preprocessing step, we take the cumulative average of the logPA at every minute in the day as Y ipj (s r ) := r −1 r =1 log(1 + PA ipj (s )), where r = 1, . . . , 1440 denotes the minutes in a day with r = 1 referring to the midnight (12:00:00 AM) and PA ipj (s ) is the PA counts for the ith at the th minute of the jth day in the pth period. We focus on the time of the day from 5AM to 10PM, when the owners are more likely to be awake; Figure 4 and 5 show the cumulative average of the logPA for four randomly selected cats during this time over some days in all the four periods. As the profiles are relatively smooth, we work with Y ipj (s r ) between 5AM to 10PM (i.e. r = 300 to 1320) as our response profile.
There are several baseline covariates collected at the beginning of the study, notable of them are age (in days), weight (WT), a numeric radiologist evaluated disease severity score called as DJD score. The number of PA profiles in each period (m ip ) varies across the subjects and a frequency distribution of {m ip } n i=1 for all the four periods is provided in Figure 6(a).
The age of cats in the study varies between 6 years to 21 years with median age of 12 years.
Based on a simple boxplot analysis (Figure 6(b)), we removed the cats with 6 years (cat number 14) and 21 years (cat number 15) of age from further analysis.

Data analysis
To test for the significance of the treatment effect, we posit the FACM, where Age ipj is the age of the i-th cat at the j-th day of the p-th period, WE ipj is the weekend indicator, i.e. it takes values 1 if the j-th day in the p-th period is a weekend, otherwise it is zero. All the other components in the model are defined previously. The quantity Age ipj can be easily computed by adding the baseline age of the cats with the total number days spared in the study period. Instead of including age of the cat at the baseline in an additive manner, we consider that the mean of the response evolves as a smooth function of age, which allows us to the model the effect of age more generally. The DJD score is a factor that is expected to affect the PA. The activities of the cats are also expected to be different over weekdays or weekends, as their owners stay at home and spend more time with them.
[ . The growing correlation along the center is a direct consequence of the cumulative average of activities, described in the preprocessing step. Using the { φ k (·)}, we obtain the projected response W ipj,k = R r=1 Y ipj (s r ) φ k (s r ) and consider the projected model The framework of PROLIFIC allows us to test for the H 0,k : τ k (d) = 0, λ k (d) = 0 vs

Appendix
The assumptions on which Theorem 1 relies are, Employing the smoothness of τ k (d), we expand τ k (d) as a truncated polynomial basis as propriately placed knots (Ruppert et al., 2003). Similarly, we can expand the other smooth the vector of spline coefficients for the smooth treatment effect τ k (d); and α λ,k , b λ,k as the vector of polynomial basis coefficients and the spline coefficients respectively for the carryover effect λ k (d). As it is common in the literature we treat the coefficients of the polynomial terms as fixed but unknown parameters and the coefficients of the non-polynomial terms as random.
Using the mixed model representation we can write µ k (d ipj ) = X ipj,µ α µ,k + Z ipj,µ b µ,k , where hµ + ) and b qµ,k 's, are assumed to be iid with mean zero and variance σ 2 µ,k for q = 1, . . . , Q µ . For the treatment effect we can Let X i,µ be the m i × (h µ + 1) fixed design matrix constructed by row-stacking X ipj,µ over p and j and Z i,µ be the m i × Q µ random design matrix obtained by row-stacking Z ipj,µ , where Similarly, construct X i,τ and Z i,τ for the treatment effect and X i,λ and Z i,λ for the carryover effect respectively. Further construct a m i ×L matrix X i,w corresponding to the baseline covariates C i in the model (4), i.e. X i,w = 1 m i ⊗[C i1 , . . . , C iL ] and column stack it with X i,µ to construct where 1 m i is the m i -length column vector 1's and ⊗ denotes the Kronecker product. Denote by N := n i=1 m i the total number of curves for all the subjects, by X b = X 1,b , . . . , X n,b the N × (L + h µ + 1) matrix of X i,b 's, by X τ = X 1,τ , . . . , X n,τ the N × (h τ + 1) matrix of fixed effect for τ k (d), by X λ = X 1,λ , . . . , X n,λ the N × (h λ + 1) matrix of fixed effect for λ k (d), by Z τ = Z 1,τ , . . . , Z n,τ the N × Q τ matrix of Z i,τ 's and by Z λ = Z 1,λ , . . . , Z n,λ the N × Q λ matrix of Z i,λ 's. Furthermore, let Y k = Y 1,k , . . . , Y n,k with Y i,k be the columns vector of the projected responses for ith subject by stacking Y ipj,k over p and j's, and the residual vector e k = e 1,k , . . . , e n,k with e i,k is constructed by stacking i,k (d ipj ) over all p and j.

S8. Selection of the orthogonal basis
Testing the original hypothesis problem in (2) reduces to simultaneous sequential testing of H 01,k and H 02,k , k = 1, . . . , K for a large value of K. Moreover, the above testing framework requires a specified set of orthogonal basis system {φ k (s)} k 1 for the space L 2 (S) to compute the projected response and test H 0k under the projected model (4). Theoretically, any known preset orthogonal basis function such Fourier basis, wavelets or Legendre basis will work.
However the selection of truncation parameter K becomes difficult, and typically that will require to test a very large number of simpler hypotheses of the form H 0,k . To avoid this, we choose a set of data-driven eigenbases from an appropriate covariance function, as adopted by Koner et al. (2021). Specifically, define the so-called "marginal covariance" Ξ(s, s ) of the error process i (s, d ipj ) in model (1) by marginalizing over the the sampling distribution of the design points d ipj 's. As Ξ(s, s ) is guaranteed to be a proper covariance function (Park and Staicu, 2015), we extract the eigenfunctions {φ k (s)} from the spectral decomposition of Ξ(s, s ) (Mercer, 1909). Finally, the truncation parameter K is chosen as the minimum value of K such that K k=1 λ k / tr(Ξ) PVE, where λ 1 λ 2 · · · 0 are the eigenvalues, tr(Ξ) := Ξ(s, s)ds is the trace of the covariance Ξ(s, s ) and the PVE (typically 95%) is some pre-specified threshold indicating the percentage of variation explained.
Although using a set of eigenfunctions identifies principal sources of variations in the data and provides an objective framework for choosing K parsimoniously, in practice these eigenfunctions are unknown. As a result, the projected response Y ipj,k can not be computed unless the eigenfunction are estimated with high accuracy. A detailed description of the estimation of the eigenfunctions from the marginal covariance is laid out in Koner et al.
(2021), we omit it here to avoid redundancy. We develop the testing procedure using these estimated set of eigenfunctions φ k (·) : k = 1, . . . , K as our choice of orthogonal basis functions. Furthermore, we derive the null distribution of the proposed test statistic; the results rely on the uniform convergence of the eigenfunctions estimators (Theorem 1).

S10. Details of estimation of FACM
To obtain the smooth estimates of bivariate mean function µ(s, d), treatment effect τ (s, d) and carryover effect λ(s, d) we fit the FACM using gam() function in R package mgcv (Wood, 2004). Using the residulals, we estimate the marginal covariance using the bivariate sandwich smoother by Xiao et al. (2013) implemented in the fpca.face() in R package refund (Goldsmith et al., 2020). After estimating K and the eigenfunctions φ k (s) with a PVE of 90%, we project the response onto the direction of the eigenfunctions. Next, we fit the smooth components of the projected model (4) using a truncated linear basis (h µ = h τ = h λ = 1) by placing the knots κ 1 , . . . , κ Q at a equally spaced quantile levels of the observed visit times {d ipj } m ip j=1 : i, p with a number of knots Q = max{20, min(0.25 × number of unique d ipj , 40)} (Ruppert et al., 2003), which is also the same for µ, τ and λ.
For each k, the covariance function Σ k of i,k (d ipj ) is estimated nonparametrically using fpca.sc() function in refund package. The number of eigenfunction is chosen with PVE of 90%. After denoising the quasi projections W k,ij with the inverse square root of estimated covariance matrix Σ W,k , we fit the LMM in (5) using the lme() function in nlme (Pinheiro et al., 2021) package and conduct the two-stage test in (12) by simulating from null distribution of the test statistics in (7), (14) and (15) implementing algorithm B of Wang and Chen (2012). Finally the overall conclusion for the hypothesis (2) is drawn combining the results of each of the k = 1, . . . , K tests as per the rule (13).

S11. Discussion of assumptions
Assumption 1 ensures that the number of curves for all subjects are finite. The moment condition in Assumption 2 is very common in FDA literature (Zhang et al., 2016;Xiao et al., 2020). It relates to the continuity of the sample paths of the error process (s, d).
The condition ensures that the projection of the response trajectory onto the eigenfunction φ k (·), Y ipj,k = Y ipj (s)φ k (s)ds is consistently defined and that the projected response has finite second moment. Assumption 3 states that the unobserved projected response Y k in model (5) is multivariate Gaussian and is the key ingredient to derive the null distribution of the test statistic in Theorem 1. Gaussianity of Y ipj,k follows if the original response Y ipj (·) is distributed as a Gaussian process with continuous sample paths and the eigenfunctions {φ k (·)} k 1 are continuous (Taylor and Karlin, 2014, chapter 8). We want to point out that we do not make any distributional assumption for the quasi-projections W ipj,k , which are based on the eigenfunctions φ k (·) that are estimated from the full data. Assumption 4, inspired by Staicu et al. (2014), is crucial to justify that the approximation error by pre-whitening the response with the estimator of the true covariance function Σ k goes away as n → ∞.
When the eigenfunctions {φ k (·)} k 1 are estimated consistently at a certain uniform rate of convergence, the accuracy in the estimation of Σ k through the quasi projections W ipj,k does not degrade compared to when Σ k is estimated through the unobserved projected response Y ipj,k . A mathematical justification of this is provided in the appendix of Koner et al. (2021).

S12. Proof of theorems and corollaries
The proof of Corollary 3 and 4 goes exactly in the same way as the proof of Theorem 1. We provide a detailed proof of the Theorem 1 and omit the proofs of the corollaries to avoid redundancy.
We first derive the asymptotic null distribution of pqGF S1 N,k , assuming that π k and η k are known. The linear mixed effect model in (5) can be equivalently written as where α −λ,k := (α b,k , α τ,k ) and under assumption 3, E k follows a Normal distribution with Denote the scaled version data and the design matrices of the model by the inverse square root of the true covariance of the data under the null model as Y 0 k := Then the linear mixed effect model can be conveniently written as For notational simplicity, denote the covariance matrix in the statement of the theorem, 0,k . Further, define P 0,k := X 0,k ( X 0,k X 0,k ) −1 X 0,k be the projection matrix onto the column space of X 0,k . Further define, P −λ,0,k := X −λ,0,k ( X −λ,0,k X −λ,0,k ) −1 X −λ,0,k be the projection matrix onto the column space of X −λ,0,k , where X −λ,0,k is the submatrix of X 0,k after removing the columns corresponding to α λ,k . Since I N − P 0,k is a projection matrix with rank N − r, there exists a N × (N − r) matrix G 0,k with orthonormal columns such that, where D k := diag( ζ 1,k , . . . , ζ N −r,k ) with ζ s being the sth eigenvalue of the (N − r) × (N − r) matrix G 0,k Z λ,0,k Z λ,0,k G 0,k . See the supplementary material of Wang and Chen (2012) for the construction of this matrix G 0,k . Moreover, since P 0,k − P −λ,0,k is also a projection matrix with rank (h λ + 1), there exists a N × (h λ + 1) matrix with orthonormal columns G λ,0,k such that G λ,0,k G λ,0,k = P 0,k − P −λ,0,k .
The most fascinating thing in the proof is that even after using the quasi-projection W ipj,k instead of the unobserved projection Y ipj,k we get null distribution that is same to what Oh et al. (2019) obtained using a known eigenfunction, upto a a remainder term that sharply goes to zero in probability at a rate that is dependent on the convergence rate of the eigenfunctions φ k (s) and accuracy rate of estimation of Σ k . By the definition of W ipj,k and Y ipj,k , the difference between these two quantities are, Note that under the null, δ ipj,k : are stacked version of δ i,k 's for all subjects. Similarly define δ k := (δ 1,k , . . . , δ n,k ) and δ k := (δ 1,k , . . . ,δ n,k ) . Further define the scaled version of these as δ 0,k := V −1/2 W,kδ k . This implies W 0,k = Y 0,k + δ 0,k and δ 0,k = δ 0,k + δ 0,k , which we will used in the later part of the proof. With all the notations defined, we can now move onto proving the null distributions of pqGF S1 N,k . Omitting the dependence on (π k , η k , γ k ) for brevity, the RSS under the null model can be written as, Similarly, the RSS for the full model can be expressed as, Now we will work with the quantity in the center of the quadratic form above. Note that by Woodbury matrix inversion identity, Using this and by one more application of Woodbury identity, This implies, and, Combining equation (20) and (21) with (16), we get, Putting W 0,k = Y 0,k + δ 0,k , in (19) and (22) we obtain, Although, ζ s,k : s = 1, . . . , N − r are the eigenvalues of the matrix G 0,k Z λ,0,k Z λ,0,k G 0,k , it is not of full rank if Q λ < N − r. In fact, the non-zero eigenvalues of this matrix coincide with the eigenvalues of the Q λ × Q λ matrix Z λ,0,k G 0,k G 0,k Z λ,0,k = Z λ,k V k (π k , η k , 0) −1 (I N − H k (π k , η k , 0)) Z λ,k as in Theorem 1. This means D k = diag(ξ λ,1,k , . . . , ξ λ,Q λ ,k , 0, . . . , 0). Now, separate the terms involving δ 0,k in equation (23) and (24) and define Λ k := σ −1 k G 0,k E 0,k and θ k := σ −1 k G λ,0,k E 0,k . Under the null hypothesis, Y 0,k = X −λ,0,k α −λ,k + E 0,k and thus, G 0,k X −λ,0,k = 0, G λ,0,k X −λ,0,k = 0, and G 0,k δ 0,k = 0, G λ,0,k δ 0,k = 0. Therefore, the RSS under the null and alternate hypotheses reduce to, where, This implies that under the null the test-statistic is equal to, By Assumption 4 and an application of continuous mapping theorem (see Staicu et al., 2014, lemma 7.1), Λ s,k and sth element of σ −1 k G 0 k E 0 k are asymptotically equivalent. Moreover, as Gaussian random variable with mean 0 and variance σ 2 k . Thus Λ s,k are asymptotically independent standard Gaussian. By the same argument, θ s,k are also asymptotically independent are also asymptotically independent of . Thus, it remains to show that T S1 0,k −T k and T k /N converges to zero in probability as n → ∞. Defining ω k := G 0,k δ 0,k and ν k := G λ,0,k δ 0,k , From equation (27) and (28) To show that where the first line is due to independence of the error across i, third line follows by Hölder's inequality. The proof is now complete by the uniform convergence of the eigenfunctions φ k (s), Assumption 1 and 2. When the variance parameters π k and η k are unknown but consistently estimated, V k ( π k , η k , 0) −1/2 V 0,k V k ( π k , η k , 0) −1/2 = I N + o p (1). Then the variance of Y 0,k under the null hypothesis will be approximately σ 2 I N and under the full model will be approximately V 1,k . By replacing the π k and η k with π k and η k , we obtain the spectral decomposition of the test statistic. The decomposition of the restricted likelihood under the full model can be derived in the same way as done in the online appendix of Wang and Chen (2012).

S13. Additional results relevant to the simulation study
[ Table 3