Testing Linear Operator Constraints in Functional Response Regression with Incomplete Response Functions

Hypothesis testing procedures are developed to assess linear operator constraints in function-on-scalar regression when incomplete functional responses are observed. The approach enables statistical inferences about the shape and other aspects of the functional regression coefficients within a unified framework encompassing three incomplete sampling scenarios: (i) partially observed response functions as curve segments over random sub-intervals of the domain; (ii) discretely observed functional responses with additive measurement errors; and (iii) the composition of former two scenarios, where partially observed response segments are observed discretely with measurement error. The latter scenario has been little explored to date, although such structured data is increasingly common in applications. For statistical inference, deviations from the constraint space are measured via integrated $L^2$-distance between the model estimates from the constrained and unconstrained model spaces. Large sample properties of the proposed test procedure are established, including the consistency, asymptotic distribution and local power of the test statistic. Finite sample power and level of the proposed test are investigated in a simulation study covering a variety of scenarios. The proposed methodologies are illustrated by applications to U.S. obesity prevalence data, analyzing the functional shape of its trends over time, and motion analysis in a study of automotive ergonomics.


Introduction
We develop a new scope of inferential procedures for testing the shape constraints on the regression coefficients in function-on-scalar regression models through the linear operator representation, where functional responses are incompletely observed.We assume that the functional response Y i (t) is available for t ∈ I i , where I i ⊂ [0, 1] is an individual-specific random subset independent of the stochastic mechanism that generates the complete functional response Y i for i = 1, . . ., n.We allow I i for a union of sub-intervals, discrete subsets, or the composition of the two scenarios.The functional response fully available on the domain is a special case by simply letting I i = [0, 1].Recently [1], [2], [3] and [4] studied functional data analysis of partially observed curves, including principal component analysis, mean and covariance functions estimation, and optimal reconstruction of individual curves, but hypothesis testing problem has been less developed.
We consider a class of testing composite null hypotheses on the functional regression model (1) of the form equivalently H 0 : β ∈ ker(C), where C is a linear operator that maps vector functions to the function space of inferential interest, and ker(C) is the kernel space of C. In this study we focus on testing a dual formation of the null hypothesis (2) expressed by where V (r) = {v l ∈ L 2 : l = 1, . . ., r} is an orthonormal basis that specifies the parametric family of ker(C).It is worth mentioning that (3) is a generalization of the classical linear contrast hypothesis.For example, let C ∈ R d×p be of full rank d.The null hypothesis H 0 : Cβ(t) = 0 studied in [14] identifies β(t) = u 0 (t) + d l=1 b l u l for some vector function u 0 (t) = (u 0,1 (t), . . ., u 0,p (t)) satisfying Cu 0 (t) = 0 and b = (b 1 , . . ., b d ) ∈ R d , where U (d) = {u l ∈ R p : l = 1, . . ., d} is the orthonormal basis of ker(C) in R p .Here we extend the theory from the finite dimensional constraints such as H 0 to the potentially infinite dimensional linear operator constraints in (2).
An important class of the null hypothesis (2) includes testing the shape of regression functions.For example, [15] evaluated a physical mechanism for the conjectured linear trend in the Northern hemisphere cooling analysis.[16] also tested linear or (piecewise) cubic time-course variations in gene expression experiments.The functional trends can be evaluated by the coefficient function associated with the constant covariate X = 1 and its shape constraint H 0 : β ∈ span{V }, where V is a set of L 2 -functions or a basis that specifies the functional trend of interest.In this case, the space of shape-constrained regression functions is expressed by (2), where the kernel space of C is spanned by V .Previously [17] studied similar topics by testing individual probes in the form of c, β j = 0 for a fixed known L 2 -function c as a special case of (2).Later, [18] proposed a residual-based permutation test for performing a hypothesis test on the shape of a mean function, although the large sample properties and the power behaviors of the proposed method were not investigated.Related work also includes [19], [20], [21], [22], [23], [24], [25], [26], [27], [28], [29], [30].Recently, [31] and [32] developed goodness-of-fit tests for functional models evaluated by empirical processes, and the significance of the family of models against general alternatives is tested by wild bootstrap resampling.But their applications to statistical inference are mainly aligned with validating functional linear models against a general class of non-structured functional models.Moreover, the extension of the existing methods to incomplete functional data has not been investigated.
The main contributions of our study are as follows.We extend the the goodness-of-fit test to the general testing framework (3), applicable to the model with incomplete functional response data.Our framework includes three scenarios as can often occur in practice; (i) the partially observed functional responses with random missing segments, where we have access to observations only for individual-specific sub-interval of the domain, but observation is not available on its complement, (ii) the functional responses observed with measurement errors on randomly spaced discrete evaluation points asynchronous across subjects, and (iii) a more challenging composite case, where individual curves are discretely observable over random sub-intervals of the domain.We especially investigate the theoretical property of the composition of the sub-interval censoring and discrete sampling of functional responses, where the proposed test procedure is applicable to a wide class of incomplete functional data.The asymptotic null distribution of the test statistic is also derived together with the consistency of the test with local alternatives H 1n : β = β 0 + n −τ /2 ∆, where τ ∈ [0, 1], for some β 0 = (β 0,1 , . . ., β 0,p ) and ∆ = (∆ 1 , . . ., ∆ p ) satisfying Cβ 0 = 0 and C∆ = 0, respectively.
The methodology and basic theory of the proposed test procedures under incomplete functions responses are developed in Section 2. In Section 3, we present numerical simulations, where the finite sample performance of the proposed test is evaluated in several scenarios.We also illustrate two applications from an obesity prevalence study and an automotive ergonomic experiment in Section 4. Our concluding discussion is in Section 5. Technical details, including the numerical implementation steps and theoretical proofs, are relegated to the Appendix.

Partially sampled functional responses
We first formulate the partially observed functional data as proposed in [1].Let δ 1 , . . ., δ n be a random sample of a stochastic process, defined on [0, 1], satisfying the following conditions.
The partially observed functional responses are defined by {Y i (t) : t ∈ I i , i = 1, . . ., n}, where Various types of incomplete functional data structures satisfy conditions C1-C4, including dense functional snippets [33], fragmented functional data [3], or functional data with single or multiple random missing intervals.More examples can be found in [34].Although C3 does not allow a sparse irregular sampling scheme, we consider the discretized noisy collection of partial data under the unified framework in Section 2.3.

Estimation of functional regression coefficients and asymptotics
leads to αw (t; -and (n×q)-design matrices of full rank, respectively.Substituting αw (t; β(t)) for α(t) in (4), we obtain ( where I = diag(1 n ) and P = Z(Z Z) −1 Z are the projection matrices that only depend on 1 n and Z.It follows that is the weighted least-squares estimator of β(t), where X = (I − P)X is the design matrix orthogonal to Z.This enables testing the hypothesis for β(t) while the nuisance regression coefficients related to Z are unspecified.Indeed, βw (t) represents a pointwise least-square estimator calculated based on a subset of samples, where its response information is available at given location t.It also follows from αw The expression (6) will be used in the next subsection to define the model space.
Theorem 2.1.Under tr(γ) < ∞ and conditions C1-C4, where Theorem 2.1 implies that pointwise βw (t) uniformly converges to β(t) over t ∈ [0, 1] and further follows asymptotic Gaussian process with root-n rates of convergence even under partial sampling structure.We also note that the condition on covariance function tr(γ) = 1 0 γ(t, t) dt < ∞ is commonly adopted in developing asymptotic theories on regression coefficient estimators under fully observed functional response.In practice, if we observe an undefined β w (t) at a certain range of the domain under a finite sample size, it can be estimated using interpolation or smoothing methods when the smoothness and continuity of β(t) is assumed.

The test statistic
To test the appropriateness of the shape-constrained null hypothesis (2) or equivalent (3), we compare the model estimates from the unrestricted space M = {µ = Xβ + Zη : β j ∈ L 2 [0, 1], j = 1, . . ., p} and the reduced space M 0 = {µ 0 = Xβ 0 + Zη : β 0,j ∈ span{V (r)}, j = 1, . . ., p}.To this end, we construct a test statistic which is based on the L 2 -distance between μ and μ0 defined by μ = argmin where • denotes the standard 2 -norm in R n .The objective functions with the weight matrix W(t) imply the pointwise optimization under the partially sampled responses.It can be verified that μ(t) = Xβ w (t) + Zη w (t) as in (6).Next, we define a linear operator L : L 2 [0, 1] → span{V (r)} as where f, g = 1 0 f (t)g(t) dt, and let L denote the multivariate operator that applies L in an element-wise fashion.We then get μ0 (t) = Xβ w 0 (t)+Zη w (t), where βw 0 = L βw .Even though we have partial response information for each observation, the uniformly consistent estimator βw provides the consistent model estimates μ(t) and μ0 (t) over t ∈ [0, 1].We next use them to propose a test statistic Note that T n is the integrated squared distance between the model fits obtained under M and M 0 , respectively, and we reject the null hypothesis if T n is large.Under the orthogonality between X and Z, distance between μ and μ0 is translated to the weighted distance between two coefficient estimates βw and βw 0 .While similar types of the L 2 -norm based test-statistic have been employed in [35], [14], and [22] for conventional hypothesis testing, such as testing the nullity of functional coefficients, our study considers a more general scope of the null hypothesis, using linear operator constraints, and the scope of response function sampling, allowing for partially observed functional response data.

Asymptotics and power considerations
In this section, we derive the limit law of the proposed test statistic T n under the the null and local alternative hypotheses using the asymptotic Gaussianity of βw shown in Theorem 2.1.Since a Gaussian process is closed under a linear operator and X X/n converges to Ψ in probability, under the null hypothesis (3), we can derive where C = I − L for I element-wisely operating the identity map I and We then consider sequences of local alternatives of the form where τ ∈ [0, 1], and ∆(t) = (∆ 1 (t), . . ., ∆ p (t)) represents a normalized functional deviation from the null hypothesis, independent of n.Then the asymptotic distribution of the test statistic is derived as the theorem below.
, where φ m , m = 1, 2, . .., are eigenfunctions of θ(s, t).Then, the test statistic T n converges to T ∆ in probability, defined as where λ m are decreasing-ordered eigenvalues of θ(s, t), corresponding to eigenfunctions φ m , and Based on Theorem 2.2, we obtain the null distribution of the test statistic T n and asymptotic power derivations as follows.
Corollary 2.3.Assume the same conditions as in Theorem 2.2.
(i) Under the null hypothesis, i.e., C∆ = 0, Theorem 2.2 implies that the null distribution of the test statistic T n converges to T 0 in distribution, where T 0 = ∞ m=1 λ m A m with λ m , decreasing-ordered eigenvalues of θ(s, t), and (ii) Suppose that C∆ = 0, that is, the local alternative, and ).Then, Theorem 2.2 yields the asymptotic power of the test as; lim n→∞ P (T n ≥ t α |H 1n ) = 1, where t α is the upper-α quantile of the null distribution T 0 in the case (i).
As we can see in the proof of the Corollary 2.3 in the Appendix, the asymptotic power of the test goes 1 under H 1n of (13) with any τ ∈ [0, 1) and non-zero ∆, which is desirable.When considering τ = 1, where the local alternative tends to the null with root-n rate, the non-trivial asymptotic power goes to 1 when ∞ m=1 π 2 m = ∞.In Section 3 of simulation studies, we consider different magnitudes of null-deviated signals π 2 m under τ = 1 to investigate the power in the practical setting.

Discrete observations with measurement errors
In this section, we extend the proposed test to the case where functional responses are observed with measurement errors over finite discrete points in their domains, and possibly sampled asynchronously across subjects.Being different from the partially observed sampling scheme with continuum measurement over the subset of [0, 1], we consider the case that functional measurements are collected on a discrete subset of the functional domain with additive measurement errors.Specifically, let {(Y * i , T i , X i , Z i ) : i = 1, . . ., n} be a random sample of (Y * , T, X, Z), where ) is the finite observations of the i-th subject associated with evaluation points The sampling design differs from the ones considered in the previous subsections as functional outcomes are prone to measurement errors, denoted by ε i,m , and finite observations are only available.We note that ε i,m = 0 is a special case that follows the same model (1).For statistical analysis, we assume that ε i,1 , . . ., ε i,Ni are i.i.d. as ε such that E(ε|X, Z) = 0 and E|ε| k < ∞ for some k > 2. The finite evaluation points T i,1 , . . ., T i,Ni are randomly generated by a probability density function λ(t) bounded away from zero and infinity whose derivative also exists and is bounded.We also assume that N 1 , . . ., N n are i.i.d. as an independent random integer N ≥ 1 that asymptotically increases as the sample size n becomes large.This sampling framework is similar to those considered by [36,37,38], and [39].We refer to the theorem and remark below for technical details.However, it is infeasible to apply the same procedure demonstrated in Section 2.1 because functional responses are only available at discrete evaluation points.Unlike the partially sample functional responses, we lose functional continuum in outcome variables which bases point-wise estimates (5) to calculate the test statistic (10).More importantly, the magnitude of false signals is not ignorable in the presence of measurement errors.This means that, even though infinitely many evaluation points are available, the coefficient function estimates will be biased as we may not achieve consistency.As a result, Corollary 2.3 may not serve as a reference distribution for testing (3).
To tackle the bottleneck, we employ kernel smoothing to recover the unobserved functional responses, where the false signals produced by measurement errors are mitigated, and substitute the estimated curves for the true functional responses to perform the test demonstrated in Section 2.1.Formally, as a twostep procedure, we first kernel smooth discrete observations for each subject as the Nadaraya-Watson kernel estimator of for each i = 1, . . ., n. Ỹ * i (t) is , where h > 0 is a bandwidth and K h (t) = K(t/h)/h is the scaled kernel of a symmetric density function K.Then, we define a kernel-smoothed test statistic as where β * , where t α is the level-α critical value for T 0 in Corollary 2.3.
In the pre-smoothing approach, it is critical to recover individual curves with a uniform rate of convergence on the entire domain [0, 1] since the kernelsmoothed test statistic T * n is defined as a weighted L 2 norm of C β * , while β * (t) is given by a point-wise estimate.The local constant smoothing, also known as Nadaraya-Watson type estimation, is easy to implement, but it is less preferred when the reconstruction of individual curves is of main interest in functional data analysis because the asymptotic bias near the boundary of the domain may vary with individual smoothing.However, Theorem 2.4 below shows that we can still attain the consistency of the test procedure with local constant smoothing.
, where we define T n in (10) with δ i = 1 for all i = 1, . . ., n. Remark 1.For each i-th subject, the optimal rate of univariate bandwidth for kernel estimation is typically given by h N −1/5 i [40,41].Since N i ≥ n θ for all 1 ≤ i ≤ n with probability tending to 1 (Lemma A.2 in the Appendix), the use of a common rate h n −θ/5 in Theorem 2.4 allows us to employ the existing bandwidth selectors [42,43].
However, we note that the classical pre-smoothing approach such as [17] requires densely observed functional responses over the entire domain for all subjects.In practice, this requirement is implausible when the observations are relatively sparse.In the following subsection, we introduce a new scope of partially observed functional data to ease the limitation.

Composition of partial filtering and discrete sampling
In Section 2.1, partially observed data were assumed to be evaluated over continuous subsets of the functional domain.If such data are observed discretely rather than continuously, then the observation framework reduces to that of discretely observed functional response data, and the smoothing approach of Section 2.2 may be applied.In this case, we assume that the complete observations for responses are given by {Y δ i : i = 1, . . ., n}.For random evaluation points T i = (T i,1 , . . ., T i,Ni ) and the indicator process δ i , we define a random subset We assume that T i and δ i are independent.
The corresponding discrete functional observations are given by for some j = j m ∈ I * i , which can be viewed as the discrete sampling of Y i composed with the partial filtering process δ i .Then, we define for each i = 1, . . ., n.Also, we define a kernel-smoothed test statistic as where To investigate the theoretical property of the proposed method, we assume that for some p > 2. Also, suppose that there exists an absolute constant C > 0 satisfying The reciprocal moment condition (20) implies that that the length of the random sub-interval dv is positive (a.s.).Hence, together with (21), the composition sampling has discrete observations densely available on each subinterval, but not necessarily over the entire domain.We also refer to Remark 2 below for the equivalent expression of (21).
The above observation implies that, even if discrete observations are sampled from the random segments of functional responses, the proposed method works for this case if we impose additional conditions on the filtering process δ i so that the density of T * i,m given by satisfies the key design condition of Theorem 2.5.
For the boundedness of λ * , we note that conditions (C1)-(C4) and the assumptions on λ imply the uniform lower bound, where b 0 = inf t b(t) and λ 0 = inf t λ(t).Also, (20) gives the uniform upper bound, For the smoothness of λ * , we note that where ) since λ has a bounded derivative.The moment condition (20) and the dominated convergence theorem give To analyze A i2 (s, t), using Hölder's inequality with p −1 +q −1 = 1 for p, q > 1, we have . Doing the same with A i3 , we claim that lim sup Indeed, it follows from (21) that provided that p > 2 > p p−1 = q.Therefore, combining ( 23), (24), and ( 25), we conclude that the derivative of λ * is given by .
The boundedness of (λ * ) can also be shown similarly as (22).
Remark 2. The condition (21) can also equivalently understood as To see this, we note that Similarly, we have

It follows that
Similarly, we have Therefore, (21) and (26) are equivalent because Remark 3. We provide one simple example of δ that satisfies the conditions (20) and (21).Suppose that U (1) < • • • < U (2p+k+2) be order statistics of a Uniform(0, 1) random sample of size To verify the condition (21), let s < t without loss of generality.We note that 3 dudv satisfying g(s, s) = 0, the Leibniz rule and integration by parts give for some non-zero constants c 0 , . . ., c 3 .Therefore, it follows from the mean value theorem that g (0,1) (s, u) where C * p,k = C p,k max c .The case for P δ(s) = 0, δ(t) = 1 can also be verified similarly, and we get the condition (21).

Simulation studies
In this section, we study the finite sample performance of the proposed testing procedure in terms of size control and powers under various settings.The performances under incomplete functional response models are compared to the benchmark performance, where functional responses are fully observed without measurement errors.
For each scenario, we apply three incomplete sampling schemes.First, we consider the partially observed functional responses with the random missing period M i , on which functional values on the ith trajectory are removed.By following a part of the setting in Remark 3, we generate where U (1) < • • • < U (2p+k+2) are order statistics of independent random samples of a size (2p + k + 2) from Uniform(0,1).We note that 1 − δ(t) in Remark 3 is set as our indicator process, where employing reversed indicator process does not affect the remarked conclusion.We here set constant parameters p, k, as p = k = 3.On average, for each simulation set, 30.4 % of each trajectory is removed by missing interval M i .Second, we consider functional responses irregularly collected over 80 asynchronous grid points with i.i.d.measurement errors following N (0, 0.5 2 ) added to each Y i (T i,m ), m = 1, . . ., 80.The locations of 80 grid points are uniformly sampled among 100 grids from each observation.

Empirical size and power
We examine the empirical sizes and powers of the proposed procedures for models from fully, partially, irregular, and partially irregular error-prone functional response data using their corresponding test statistics, denoted as T Full n , T n , T * n , and T * * n respectively.Practical implementation steps for each test statistic are provided in the Appendix.All simulation results below were based on 5,000 simulation replicates, and the critical value of the test was estimated by 5,000 bootstrap samples in each simulation run.To calculate the test statistic T * n and T * * n involving kernel smoothing, we chose a common bandwidth that minimizes the leave-one-out cross-validation [44,45] across all subjects in each simulation sample.
Table 1 summarizes results for hypothesis (31) at 5% nominal level un-    der scenario A from test statistics from corresponding functional response data structures, for τ = 1, 0.8, 0.67.It can be seen that the empirical sizes are reasonably controlled around the nominal level 0.05.Although the sizes under error-prone partially observed structure, corresponding to the test statistic T * * n , show slightly larger values around 0.07, and it is due to loss of original information with missing intervals and additive noise.In terms of power, we investigate the results depending on τ , which regulates the rate that the null-deviated model approaches the null model.As expected, the empirical power increases as τ decreases or as d A increases.In addition, the power reasonably approaches to 1 under all settings.Especially for τ = 1 of T Full n , the power approaches 1 even with moderate magnitudes of the null-deviated signals, indicating that the condition of 3 is not restrictive in practical application.The relatively deflated powers from T * n might be due to some loss of the null-deviated signal after applying the smoothing process to noisy data.We observe that the power from T n goes to 1 with a reasonable but slightly slower rate than T FULL n shows, and it is from the smaller effective sample sizes at each grid due to partial sampling.Although the lowest powers are observed from T * * n under all settings due to most significant loss of original information with missing periods and noisy discretized measurements, we still see the power gradually increases towards 1.Indeed, our extra simulations considering larger values of d A show that powers under τ = 1 from T * n and T * * n become 1 when d A = 13 and 17, respectively.
The simulation results from scenario B are illustrated in Figure 3.The results under n = 100 and n = 200 are represented by full and dotted lines, respectively.We observe a very similar pattern to that under scenario A with reasonable size controlling at the 0.05 nominal levels and with the behaviors of the power for d B > 0. We again confirm that power approaches to 1 when τ = 1 under the moderate magnitudes of the null-deviated signals.The power tends to 1 with relatively slower but reasonable rates with an increase of d B for T n and T * n due to the same reasons described in results from scenario A. Again, we observe the lowest powers achieved from T * * n under d B > 0, they gradually approaches towards 1.We note that extra simulations considering larger values of d B show that powers under τ = 1 for T * and T * * are attained as 1 when d B = 2.1 and 3, respectively.It implies empirically consistent properties of our proposed tests.
Although we have only illustrated the simulation result for investigating the finite sample performance of T * n with N i = 80 in Table 1 and Figure 3, we observed that the power and size of the proposed test are also well achieved with N i = 60, where relatively rich response information is available over the domain.However, under the sparse setting, N i = 10 or 30, we observed relatively unsatisfactory results with the finite sample analysis.We note that, in our simulation setting, the null-deviated signals visualized in Figure 1 are quite subtle, with a delicate difference in the trend and visually detectable discrepancies only at boundaries.Hence, we report a limitation of the proposed method for T * n such that the estimated regression coefficients calculated from noisy functional responses collected over sparse grids may not be able to effectively detect subtle trend differences near the boundary.

Real data application 4.1 The obesity prevalence trend change
We illustrate the practical application of the proposed testing procedure through an analysis of the U.S. overweight and obesity prevalence data from 2011 to 2020.It is a part of the data of the U.S. residents regarding their healthrelated risk behaviors and chronic health conditions, collected by Behavioral Risk Factors Surveillance System (BRFSS) through the state-based telephone interview survey in cooperation with the Centers for Disease Control and Prevention (CDC).The dataset consists of percentages (%) of adults aged 20 and over populations with the weight status of obese, overweight, normal weight, and underweight from 50 states.Along with weight status, socioeconomic status is also measured through educational and income levels of samples.In terms of income, each survey sample is classified into one of five categories; less than $15,000, $15,000-$24,999, $25,000-$34,999, $35,000-$49,999, and over $50,000.The full dataset can be found at: https://chronicdata.cdc.gov/Behavioral-Risk-Despite growing recognition of the problem, the obesity epidemic continues in the U.S. with steadily rising obesity rates.For example, 1999-2000 through 2017-2018, U.S. obesity prevalence increased from 30.5% to 42.4%.Figure 4 illustrates such trends during recent 10 years from 50 states for obese, overweight, and normal weight groups.The bold lines represent the sample mean trajectories of each group, where its calculation is specified later with the model specification (32).With the rising obesity rates, we observe decreasing proportions of normal weight population along with seemingly constant rates of overweight population.We apply the proposed methods to identify shape of the tendency on prevalence rates for each group of the weight status.Furthermore, the gap of obesity prevalence between low and high income groups changes during this time is also examined.
We first investigate the shape of overall prevalence trend for each weight group.Since data is collected over regular grids for all states with a few missing values, we adopt the test statistic T n for partially observed functional data.Let Y i (t m ) denote the observed prevalence rate for given weight status group from ith state.We formulate the intercept-only model with Z = 0 in (1) and rescaled discrete time points 2011, . . ., 2020 to equally spaced t m ∈ [0, 1], Based on it, we obtain the least square estimates β(t m ) = 50 i=1 Y i (t m )/50 as the sample trajectories of each weight group, illustrated with bold lines in Figure 4. To identify its shape, we consider the null hypotheses for the constant and linear spaces, corresponding to H 0,c : β(t) ∈ span{V (1)} and H 0,l : β(t) ∈ span{V (2)}, respectively.Here, V (r) = {v l (t); t ∈ [0, 1]} r l=1 is an orthonormal set we obtain by applying the Gram-Schmidt process to the polynomial basis Table 2 shows calculated test statistic T n and corresponding calculated p-values for each null hypothesis from each weight group.Calculation details and numerical implementation steps are provided in the Appendix.In Table 2, we reject the null hypothesis of constant space H 0,c for the obese and normal groups, but not H 0,l .That is, at a significance level less than 0.001, obesity prevalence has linearly risen over time, while rates of normal weight population is linearly decreased.On the other hand, we could not find any significant trend as we retain the constant shape hypothesis H 0,c at level 0.1 for the null hypothesis H 0,c .
We next investigate the obesity prevalence over time associated with income levels.In recent literature, statistical analyses on the association between in-come levels and obesity rates have repeatedly reported that obesity prevalence has been significantly increased at a faster rate mostly in relatively low-income levels [46,47,48].Figure 5 (a) illustrates obesity prevalence rates for five income levels and their mean trajectories.While all five income levels present increasing obesity prevalence over time, the group of income less than $15, 000 shows the highest rates while the groups of income over $50, 000 illustrates the lowest rates.We first apply the functional ANOVA to this data, a special case of our proposed testing procedures corresponding to a part of [49].To do this, we formulate the model based on (4) by setting (250 × 4) matrix for X = diag{1 50 , . . ., 1 50 } and (250×1) vector of 1's for Z.The null hypothesis for fANOVA corresponds to (3), where V = {0}; i.e., H 0 : β j (t) = 0, for t ∈ [0, 1], and j = 1, . . ., 4. By applying the proposed testing procedure, we obtain p-value < 0.001 and conclude that significant differences on obesity rates among different income groups exist.We then apply a type of post hoc test, specifically to examine how the gap of prevalence among lowest highest income group changes over time.Figure 5 (b) shows mean trajectory of gap in obesity rates between the lowest and highest income groups.It is observed that this gap tends to decrease over time and fitted linear and quadratic trending lines are illustrated, respectively.We also consider the fit with the piecewise linear bases, where its fitting details and testing results are specified later.We apply the proposed procedure based on the test statistic T n to identify the shape of this gap.Let where Y level1 (t m ) is a vector of length 50 with the elements of obesity rates for the lowest income group from 50 states at mth year.Similarly, Y level5 (t m ) denotes a vector for the highest income group.We then specify the model based on (4) with X = (1 50 , 0 50 ) and Z the length 100 vector of 1's.Under given model formulation, β(t) represents the difference between two groups means.Then two null hypotheses of linear and quadratic functional spaces are considered, H 0,l : β(t) ∈ span{V (2)} and H 0,q : β(t) ∈ span{V (3)}, respectively.By applying the proposed testing procedures; for the test under H 0,l , we obtain T n = 10.75 with the p-value 0.02; and for the test under H 0,q , T n = 5.03 and p-value is 0.23.Under significance level α = 0.05, we fail to reject the quadratic null space and conclude that the gap of obesity prevalence between lowest and highest income groups is significantly decreasing with the quadratic shape.To demonstrate the further application of our method with null hypothesis with other types of bases, we try the hypothesis testing for H 0,pl : β(t) ∈ span{U (3)}, where U (3) represent a set of three orthonormal B-spline bases derived from the piecewise linear functions with knots at 0, 0.5, and 1, where the internal knot 0.5 is chosen by the estimated peak from the quadratic fit.Under this null hypothesis, we obtain T n = 4.70 and p-value 0.28.By comparing obtained p-value 0.28 with the p-value 0.23 derived from the null hypothesis H 0,q , we observe slightly stronger statistical evidence on the conclusion for the piecewise linear shape on gap between two groups under given sample sizes.We note that results from smoothed trajectories through the test statistic T * n leads the same inferential conclusions for H 0,l , H 0,q , and H 0,pl , although they are not presented here.It empirically demonstrates the performance of our proposed method in detecting significant functional shape even under non-smoothed raw trajectories.states and (b) mean of differences of the obesity rates between group of income less than $15, 000 and group of income over $50, 000, with fitted lines using linear bases ( ), using quadratic bases ( ) and piecewise linear bases ( ).

Human motion analysis in ergonomics
We illustrate another data example in automotive ergonomics, previously analyzed by [50], [35], [22], [32], and among others.The Center for Ergonomics at the University of Michigan collected data on body motions of an automobile driver.As part of the project, the right elbow angles of the test driver were captured as time-varying responses when the driver's hand leaves the steering wheel until reaching 20 different locations in the car.There were 3 repeated reaches to each of the different targets located near the glove compartment, headliner, radio panel, and gear shifter.We associate observed discrete trajectories of elbow angles R ij (t ij,m ) with with the (x, y, z)-coordinate of a reaching target with extra variables as for i = 1, . . ., 20, j = 1, 2, 3, and k = 1, . . ., N i , where (c i1 , c i2 , c i3 ) represents the (x, y, z)-coordinate of a target location with its origin at the initial hand posture on the steering wheel and d ik 's are 0-1 dummy variables indicating four nominal areas of different targets.Specifically, d i1 = 1 if the target is located near the headliner, d i2 = 1 if the radio, d i3 = 1 if the gear shifter, and zeros otherwise so that we set the glove compartment for the baseline location.By adding the nominal target information to the conventional model, we are able to statistically compare the changes of elbow angles from different experimental conditions.Among 60 experiments, we drop one trial which has been excluded in the literature, where the researchers revealed that the driver's motion was mistaken while reaching the target.See [50] for more details about the experimental settings.
Since observed discrete trajectories of elbow angles reveal some noises due to the measurement errors, [50] applied the smoothing splines to raw data to respect the smoothness of human motion and obtained pre-smoothed angle random curve.We denote it as R * ij (t), where the tracking time points t ij,1 , . . ., t ij,Nij were re-scaled to [0, 1] for each of 60 reaches.The pre-smoothed random sample { R * ij : i = 1, . . ., 20, j = 1, 2, 3} can be analyzed by the standard one-way functional ANOVA.We note that the model ( 33) turns out to be adequate for the data as the bootstrap-based test [50] does not reject the lack of fit compared with the functional ANOVA model (p-value = 0.436).[51,52] also considered similar approaches to the driver's motion prediction in a larger dataset by adding extra variables to statistically control different experimental conditions.
In this example, we aim to analyze the shape of the time-varying motion changes rather than find a predictive model for an arbitrary target location.Based on the asymptotic equivalence between splines and certain class of kernel estimates [53,54], we apply the proposed method to the pre-smoothed random sample { R * ij : i = 1, . . ., 20, j = 1, 2, 3} to test the null hypothesis H α 0 : {α 1 , α 2 , α 3 } ∈ span{V (2)} with the same V (2) defined in Section 4.1.We perform inference using T * n from Section 2.2 and find that the null hypothesis cannot be rejected (T * n = 3.16, p-value = 0.660).This result together with Figure 6 shows that, compared to the glove reaching experiment, the driver stretched their elbow less and moved slower at a constant relative angular velocity when reaching different area.We also individually test several hypotheses such as H We close this section by reporting that all twelve hypotheses we have tested were still not rejected after applying the multiple comparison adjustment, both the Bonferroni and Benjamini-Hochberg corrections, at 5% significance level, implying statistically significant linear trends on them.

Discussion
We have presented a statistical procedure for testing shape-constrained hypotheses on regression coefficients in function-on-scalar regression models, generalizing existing methods such as fANOVA that consider nullity hypotheses only.The approach presented here enables inferences about temporal/spatially vary- ing coefficient effects as well.The large sample properties of the proposed test were investigated by deriving the asymptotic null distribution of the test statistic and consistency of the test against local alternatives.The methodology was demonstrated under three incomplete sampling situations; (i) partially observed, (ii) irregularly observed error-prone, and (iii) composition of former two incomplete functional response data.A few studies have recently illustrated goodness-of-fit tests for functional linear models under fully observed responses, but handling incomplete sampling designs was not studied either in theory or practice.Furthermore, the critical value in our methodology can be approximated with the spectral decomposition of covariance function, for which one can easily exploit the existing methods in the recent developments in functional data analysis.A key aspect of the methodology developed here is the specification of a relevant shape hypothesis of interest.Ideally, the application defines the relevant shape space.Otherwise, we can use standard curve-fitting hypotheses defined by, for example, polynomial basis functions, exponential functions, or periodic functions cycling at different frequencies.
In Section 2.2, we considered functional data, where each sample path is observed on randomly spaced discrete points of size N i .Assuming that N i 's are increasing as the sample size increase, which is often called "densely observed" functional data, we adopted the individual smoothing strategy as interpolation.Another interesting and challenging situation is when the functional data are so sparsely observed that the individual smoothing strategy employed here is not effective.In this case, one may consider a functional principal components based approach to reconstruct individual curves.Recently, [4] proposed optimal reconstruction of individual curves in which each of the incomplete n functions is observed at discrete points considerably smaller than n in finite sample analysis.They showed that the functional principal components based approach can provide better rates of convergence than conventional smoothing methods, where min i {N i } n θ as n → ∞ for some θ > 0. However, hypothesis testing under the functional principal component analysis framework remains to be developed.
Z Y(t m ).We then derive the discretized γ(s, t) following the formula in (12), denoted as Γ, by calculating Γ = Γ − Γ (c) − Γ (r) + Γ (c,r) , where Γ (c) = (γ (c)1 , . . ., γ(c)N ) is the matrix of the fitted multi-response regression values based on N separate regressions, with each column of Γ as the response, and the r-columns of V as covariates.Simiarly, Γ (r) is the matrix of fitted multi-response regression values, with each row of Γ as the response, and the r-columns of V as covariates.Lastly, Γ (c,r) is the matrix of fitted values by applying previous two steps to Γ subsequently.We then generate a large bootstrap samples of T0 = p , where K denotes the number of positive eigenvalues of Γ, and use its (1 − α)% quantile as the critical value.
To perform the proposed hypothesis testing under irregularly collected response data with additive measurement errors, we replace Y i (t) by kernel smooth estimates of Ỹ * i (t), obtained by ( 16), and apply the procedures described for the fully observed response data.To find the optimal smooth parameter in kernel estimation, one may adopt the leave-one-out cross-validation.
For the application of the proposed testing procedure to partially observed functional response data, we calculate (p × N ) matrix D w = βw − L βw , where D w = (D w 1 , . . ., D w N ) and βw j = ( βw j (t 1 ), . . ., βw j (t N )) , and approximate To find the empirical critical value for the level α test, we estimate Γ by employing nonparametric covariance surface estimation method applicable to sparse functional data, available through the function GetCovSurface in R package 'fdapace'.The optimal bandwidth for the surface estimation can be found through the cross-validation steps.Next, we calculate v(t m , t m ) = Lastly, we can calculate T * * n by combining two previous steps, smoothing process over observed trajectories and calculation of test-statistic under partial structures.

A.2 Technical Details for Section 2.1
Proof of Theorem 2.1 Suppose E( X|Z) = 0 without loss of generality.Under the null hypothesis, and let Z n (t) = √ n( βw (t) − β 0 (t)).Then we can write Let X = ( X1 , . . ., Xn ) , where Xi = (x i1 , . . ., xip ) , and let be the p-variate random functions with the zero mean, corresponding to the third term in (34).Its j-th element is specifically written as , where T Q = {t q ∈ [0, 1] : q = 1, . . ., Q} is a finite collection of any Q time points, for Q ≥ 1.By the multivariate CLT and the mutual independence among xij , δ i , and i , we have where Ξ = ϑ qq 1≤q,q ≤Q is the Q × Q covariance matrix with ϑ qq = γ(t q , t q )v(t q , t q )b(t q ) −1 b(t q ) −1 , Ψ = Ψ jj 1≤j,j ≤p is the p × p matrix with Ψ = E(Var(X|Z)), and the Kronecker product of Ξ and Ψ is given by We specifically derive Ξ ⊗ Ψ as follows.For p-variate random variable V n (t q ), the diagonal of its asymptotic covariance matrix, i.e., (j, j)-th element of the matrix, is derived as γ(t q , t q )b(t q ) −1 E(x 2 ij ) = ϑ qq Var(x ij ), and the (j, j )-th element of the covariance matrix, for j = j , is γ(t q , t q )b(t q ) −1 E(x ij xij ) = ϑ qq Cov(x ij , xij ).That is, the block diagonal covariance matrix of vec(V n ) is ϑ(t q , t q )Ψ.We then examine the block off-diagonal covariance matrix of vec(V n ) by calculating the covariance between V n (t q ) and V n (t q ), for q = q .we can show that the diagonal (j, j)-th element of the covariance matrix is ϑ qq Var(x ij ) and the (j, j )-th element of the matrix, for j = j , is ϑ qq Cov(x ij , xij ).That is, p × p off-diagonal block covariance matrix of vec(V n ) is written as ϑ qq Ψ.By [55] and [56], the multivariate process {V n (t) : t ∈ [0, 1]} converges to the multivariate Gaussian process in distribution as where the finite-dimensional restrictions of ϑ is given by the covariance matrix Ξ. Next, we can show that the (p × p) matrix X W(t) X/ n i=1 δ i (t) in the second term of (34) converges to Ψ in probability, under the conditions C2 and C4.
where e k ∈ R p be a unit vector whose k-th component is 1.We note that β * j (t) = βj (t) + e j ( X X) −1 X Ỹ * (t) − Y(t) for each j = 1, . . ., p, where βj (t) = e j ( X X) −1 X Y(t) is the least-squares estimator of β j (t) with fully observed data.The large sample property of β(t) = ( β1 (t), . . ., βp (t)) follows from Theorem 2.1 by letting I i = [0, 1] for all i = 1, . . ., n.Since I −L op ≤ 1, it follows from the Cauchy-Schwarz inequality and Lemma A.2 that The above result is analogous to the proof of theorems in [14].On the other hand, (11) gives
Figure 2 (a) illustrates a randomly selected set of fully observed response trajectories, and three other sets of trajectories in Figure 2 (b), (c), (d) display partially observed response trajectories filtered by missing random intervals, noisy responses generated over irregular grid points with additive measurement errors, and noisy partially observed responses over irregular grids with additive measurement errors, respectively.

Figure 2 :
Figure 2: Randomly selected six simulated trajectories of (a) fully observed response data, (b) partially observed response data filtered by independent missing intervals, (c) irregularly observed data with added measurement, and (d) partially observed data over irregular grids with measureme errors.

Figure 4 :
Figure 4: Percentages (%) of (a) obese, (b) overweight, and (c) normal weight adults among the U.S. adults aged 20 and over populations from 2011 to 2020, from 50 states (gray lines) and sample means (solid lines)

Figure 5 :
Figure 5: (a) Percentages (%) of obese prevalence by five income levels from 50states and (b) mean of differences of the obesity rates between group of income less than $15, 000 and group of income over $50, 000, with fitted lines using linear bases ( ), using quadratic bases ( ) and piecewise linear bases ().

Figure 6 :
Figure 6: The regression coefficient estimates of the model (33) are depicted.The solid and dot-dashed lines are coefficient estimates and their 95% confidence bands, respectively.The long-dashed lines show the estimates under the null hypotheses H 0,l .

n
i=1 δ i (t m )δ i (t m )/n, and b(t m ) = n i=1 δ i (t m )/n.Then (N × N ) matrix Ξ, discretized version of ϑ(s, t) in Theorem 2.1, can be derived, where its (m, m )-th element is calculated as Ξ mm = Γ mm Π mm .Here, (N × N ) matrix Π has its (m, m )-th element asΠ mm = v(t m , t m ) b(t m ) −1b (t m ) −1 .We next calculate Ξ = Ξ − Ξ (c) − Ξ (r) + Ξ (c,r) , where Ξ (c) and Ξ (r) are obtained by following definitions of each term, described in the implementation for the fully observed data.In practical application, we adopt the standardized test statistic Tn = N −1 N m=1 Dw m ( X X) Dw m , where Dw m = D w m b(t m )v(t m , t m ) −1/2 and obtain an approximate critical value from Ξ * = Ξ * −Ξ * (c) −Ξ * (r) +Ξ * (c,r) , where Ξ * mm = Γ mm Π * mm with Π * representing the standardized matrix of Π having unit variance for diagonals.The standardized testing procedure empirically shows the improved performance in size controlling in simulation studies.

Table 1 :
Empirical size and power at the 5% nominal level for testing H