Efficient Estimation for Varying-Coefficient Mixed Effects Models with Functional Response Data

In this article, we focus on the estimation of varying-coefficient mixed effects models for longitudinal and sparse functional response data, by using the generalized least squares method coupling a modified local kernel smoothing technique. This approach provides a useful framework that simultaneously takes into account the within-subject covariance and all observation information in the estimation to improve efficiency. We establish both uniform consistency and pointwise asymptotic normality for the proposed estimators of varying-coefficient functions. Numerical studies are carried out to illustrate the finite sample performance of the proposed procedure. An application to the white matter tract dataset obtained from Alzheimer’s Disease Neuroimaging Initiative study is also provided.


Introduction
Longitudinal and functional data are frequently encountered in biomedicine, epidemiology, and other fields of natural and social sciences. See, Hand and Crowder (1996), Ramsay and Silverman (2005) and Ferraty (2011) for abundant interesting examples. In recent years, various functional regression models are proposed to conduct a better description and analysis for this complicated data, since it allows the responses or the regressors, or both, are curves. According to the difference in the concrete forms of the responses and regressors, functional regression models are usually subdivided into three categories, including scalar-on-function regression (see, e.g., Cai and Hall 2006;Chen, Hall and Müller 2011;Sang et al. 2018), function-on-function regression (Yao et al. 2005a;He et al. 2010;Wang et al. 2018) and function-on-scalar regression (Wu et al. 2013;Zhang and Wang 2015). For more details about the functional regression models and its applications, one can see Morris (2015), Wang et al. (2016) and Reiss et al. (2017) for an overview.
Motivated by the regression analysis of imaging (or functional) data, such as diffusion tensor imaging (DTI) and positron emission tomography (PET) (Zhu et al. 2007;Friston 2009), many nonparametric and semi-parametric regressions for functional responses have been proposed, see, for example, Luo et al. (2016) and Li et al. (2017). A functional varying-coefficient mixed effects model, which can be used for modeling the dynamic impacts of the covariates of interest on the response, is one of the most powerful tools in addressing these scientific problems for imaging data. Let Y i (t) be the functional response variable for subject i, i = 1, . . . , n, and x i = (x i1 , . . . , x i p ) T be the associated p-dimensional covariate vector, then the functional varying-coefficient mixed effects model takes the form where α(t) = (α 1 (t), . . . , α p (t)) T is the p-dimensional unknown but smooth coefficient function vector, ε i (t) are measurement errors, and η i (t) is the random effect function that characterizes individual curve variation (subject-effect) from x T i α(t). Without loss of generality, throughout this paper we assume that t ∈ T = [0, 1]. Moreover, {η i (t) : t ∈ T } is a stochastic process indexed by t ∈ T that captures the within-curve dependence. We assume that η i (t) are independent copies of η(t) ∼ SP(0, ), where SP(μ, ) denotes a stochastic process with mean function μ(t) and covariance function (s, t).
Model (1), as a classical function-on-scalar regression has been widely applied for functional data analysis and extensively studied in the existing literature. For instance, for model (1) without involving the random effect item, Wu and Chiang (2000) proposed two kernel estimators based on componentwise local least-squares criteria to estimate the varying-coefficient functions; Chiang et al. (2001) employed a componentwise smoothing spline method to estimate the time-varying coefficients and studied the corresponding asymptotic properties. Zhang and Chen (2007) adopted the idea of "smoothing first, then estimation" to construct the local polynomial kernel reconstruction-based estimator for mean, covariance and error variance functions for model (1) and derived their asymptotics. Zhu et al. (2012) extended model (1) to a multivariate case and systematically studied the theoretical properties of the resulting estimators obtained by local linear smoother. Zhu et al. (2014) further proposed a spatially varying coefficient model (SVCM) for neuroimaging study and investigated the asymptotic properties of the parameter estimators. Other more studies on varying-coefficient models and mixed effects models for longitudinal and functional data, readers may refer to Fan and Zhang (2008) and Chen and Wang (2011) for comprehensive discussions.
However, these aforementioned estimation approaches for model (1) are mainly based on kernel smoothing and smoothing spline, both of which ignore the within-subject correlation information in the estimation procedures. For longitudinal and functional data, it is well known that there exists a high correlation between the observations on the grid of each datum. As a result, these existing methods might lead to a loss of efficiency. Additionally, the kernel smoothing method utilizes only observations in the neighborhood of the target point for estimation. It is noticeable that these neighborhoods usually contain only a few observations in longitudinal studies, and the within-subject correlation may causes observations out of the neighborhoods carrying significant information. In view of this, the local kernel method might result in inefficient estimators. Therefore, some new techniques improving the efficiency of the conventional local kernel method are highly desirable in longitudinal and functional data analysis. Chen et al. (2010) proposed a global partial likelihood method to study the nonparametric proportional hazards model. The main idea of their global approach is to use all observations instead of only local observations near the target point for estimation to improve the efficiency. Zhou et al. (2018) employed that global smoothing idea combining with the full quasi-likelihood method to conduct the estimation of mean and covariance functions for sparse functional data. Motivated by the idea of global smoothing, we present an efficient estimation procedure for estimating the varying-coefficient functions in model (1) for longitudinal and sparse functional response data. This procedure combines the generalized least squares method and a modified local kernel smoothing criterion which uses the observations on the whole grid of each subject in estimation. An iterative but simple algorithm is proposed to compute the estimator.
The proposed method overcomes the potential drawbacks of existing local kernel methods and improves the efficiency of the resulting estimators. Specifically, there are several advantages for the new method: (i) the proposed estimation procedure involves the global smoothing idea that not only absorbs the advantages of the local linear smoother by using the observations near the target point for estimation, but also utilizes the information available at data outside the neighborhood, which is not used in the local kernel smoothing principle; (ii) the within-subject covariance structure is also incorporated into the estimation procedure to improve the estimators; and (iii) the iterative procedure is easy to implement so that it can be applied to solve various application problems. We derive both uniform consistency and asymptotic normality for the estimated varying coefficient functions. Simulation studies are conducted to compare the performances of our method and the local linear method in terms of root mean squared error, and the results reveal that the proposed method has a smaller root mean squared error than its competitor. This discovery further indicates that the results obtained by our method can more accurately characterize the dynamic relationship between a set of covariates and the functional response variable in application problems. Finally, a real data example is presented to illustrate the newly proposed approach.
The rest of this article is organized as follows. Section 2 introduces the estimation procedure. The theoretical properties including consistency and asymptotic normality of the estimators are established in Sect. 3. Section 4 presents the simulation studies which are performed to evaluate the performance of the method. Section 5 is devoted to applying this method to analyze the white matter tract dataset obtained from Alzheimer's Disease Neuroimaging Initiative (ADNI). Section 6 concludes the paper and discusses future work. Technical proofs are given in the "Appendix".

Methods
In this paper, we focus on sparse data situations. Suppose that we observe the random functions Y i at discrete points t i j , that is, where the observation locations t i j are independent and identically distributed random variables with density function f (·). Throughout this paper, we shall assume that x i , t i j , η(·) and ε(·) are all independent and the noise terms, ε i j , are independent and identically distributed random errors with mean zero and variance σ 2 . Suppose that 2 ≤ m i ≤ M, for some finite constant M, for all i = 1, . . . , n. Our proposed estimation procedure contains two key steps: (I) calculate the initial estimateα(t) of α(t) and construct the nonparametric covariance estimation for (s, t); (II) use the estimated covariance function and the initial estimateα(t) to develop a refined estimate for α(t), denoted by α(t).

Initial estimation
We adopt the local linear scatter-plot smoother (Fan and Gijbels 1996) to estimate the varying coefficient function α(t) and the covariance function (s, t). For any T in a small neighborhood of t, applying the Taylor's expansion, we have . . , a p ) T and A = (a 1 , . . . , a p , b 1 , . . . , b p ) T . Let K (·) be a kernel function and K h (·) = h −1 K (·/h) be a scaled kernel with bandwidth h. Then the local linear estimate of the regression coefficients, α(t), is simplyα(t) = A 0 where ) be the "raw" covariances, wherẽ α(·) is the estimated varying coefficient function obtained from (3). Note that is an indicator function. Therefore the local linear surface smoother for covariance function (s, t) can be achieved by the raw covariances with their diagonal elements being removed, and this yields that (s, t) = β 0 , where β = ( β 0 , β 1 , β 2 ) T and In order to estimate σ 2 , we first estimate the variance The estimate for the variance σ 2 is then, where the notation |T | denotes the length of the interval T . The more detailed process for the estimation of mean and covariance functions for longitudinal and sparse functional data, one can refer to Yao et al. (2005b) and Li and Hsing (2010).

Refined estimation for varying-coefficient functions
We know that the local kernel smoothing uses only local observations near the target point for estimation. However, the within-subject correlation, presented in the longitudinal and functional data, may result in observations outside the neighborhood of the target point carrying relevant information. Therefore, the local linear method used in Sect. 2.1 for estimating the varying-coefficient functions may reduce to loss of efficiency. In this subsection, we employ the generalized least squares method incorporating a modification of local linear smoothing to conduct a refined estimator for α(·). This new method not only adopts the idea that uses all observations instead of the local observations near the target point for estimation but also incorporates the estimated covariance structure which can be used to account for the within-subject correlation into the estimation process to improve the resulting estimators. T and I m i is the identity matrix of order m i for i = 1, . . . , n. For convenience, let V i = i + σ 2 I m i . Then the generalized least squares estimator of α(t) is obtained by minimizing the follow objective function Firstly, we replace V i with its estimates. Since we do not assume any structure on the covariance function of Y (t) other than that it is a smooth function, the local linear estimates presented in Sect. 2.1 is recommended. Let i be the estimated covariance matrix of i and σ 2 be the estimated variance of measurement errors obtained from (4), (5) and (6). Then V i is substituted by V i = i + σ 2 I m i . Secondly, we adopt a similar strategy in Chen et al. (2010) and Zhou et al. (2018) to approximate α(t). Specifically, let B n (t) be a neighborhood of t, it is easy to see that Based on the idea of local linear smoother, we use linear function ζ +ζ (t i j − t) to approximate α(t i j ) in the first term of (8) and replace I (t i j ∈ B n (t)) with a kernel function K {(t i j − t)/h}. Then for t i j in a neighborhood of t, we have This expansion is of great use in longitudinal and functional data analysis by the reason of it not only absorbs the advantage of local linear smoother but also utilizes the information available at data outside the neighborhood region B n (t). Substituting this approximation and the estimated covariance matrix into (7), we can obtain the refined estimator for ϑ = (ζ T ,ζ T ) T by solving the following where Note that the refined estimator can not be directly obtained by solving the the Eq. (9) as there still exists an unknown function vector α(t i ). Therefore, we present an iterative procedure by substituting α(t i ) with its previous estimator to address this problem. Let t 1 < t 2 < · · · < t M be the distinct points of t i j for i = 1, . . . , n and j = 1, . . . , m i . The iterative algorithm is described as follows.
• Step 0 Taking the estimator fulfilled by (3) as the initial values of α(t), denoted by Step 1 For each t = t 1 , . . . , t M , solve the Eq. (9) in which α(t) is replaced by its previous estimator, that is, (10) yields the estimator ϑ(t), given by Step 1 until convergence.
The selection of the bandwidth h for the estimation procedure presented above is achieved by the leave-one-curve-out cross-validation (Rice and Silverman 1991). Specifically, let α (−i) (·) be the estimated varying-coefficient function of α(·) after removing the ith subject from the data. Then the optimal bandwidth h can be chosen by minimizing the cross-validation score based on the squared prediction error given by

Remark 1
The proposed iterative estimation procedure can be readily implemented as there is a closed-form at each step. Compared to the local linear method used in Zhu et al. (2012), the newly presented estimation approach may be more efficient since it takes the within-subject correlation structure and the observations on the whole grid of each subject into consideration, two of which are significantly important in longitudinal and functional data analysis. It is not hard to see that the estimation in (3) is, in essence, a local marginal estimation, assuming that the within-subject covariance matrices are identities. This leads to a working independence estimator for α(t) which ignores the correlation among the observations on the grid of each subject, and further causes a loss of efficiency.

Remark 2 Replacing α(t)
with its estimate to model (1), we also can conduct the estimation of each random effect function η i (t). Many existing non-parametric methods such as the local linear regression have been developed to directly estimate all individual functions η i (t), i = 1, . . . , n, (see, e.g., Zhu et al. 2012;Li et al. 2017). However, these non-parametric methods may be not suitable for our sparse data situations when the number of observation points for each subject is small. Functional principle component analysis (FPCA) as a key technique in characterizing the features of unknown functions in functional data analysis, is another feasible method for the estimation of η i (t). We recommend using the version of FPCA developed by Yao et al. (2005b) which they refer to as principal components analysis through conditional expectation (PACE) for longitudinal data to estimate each random effect function η i (t). The PACE approach is appealing since it is implemented mainly based on the estimated population covariance function of η i (t), which utilizes the whole sample information and can be easily obtained by using the method presented in Sect. 2.1. Comprehensive surveys on FPCA for longitudinal or sparse functional data were developed in Yao et al. (2005b) and Hall et al. (2006).

Asymptotic properties
In this section, we study the asymptotic properties of our estimators. Some regularity conditions needed for the theoretical developments are stated in "Appendix A.1". All these conditions are very mild and most of these are also adopted by Chen et al. (2010) and Zhou et al. (2018). The detailed proofs for the following theorems are given in "Appendix A.2-A.3".
To start with, we give some notations which will be used in the theorems. Let ν r = u r K (u)du and υ r = u r K 2 (u)du for notational simplicity. For a vec- Let I be the identity operator, and A be the linear operator satisfying

Theorem 1 Under the regularity Conditions
Theorem 2 Suppose that conditions of Theorem 1 hold. Then, for any t ∈ (0, 1), we have and ϕ is a random vector following the standard normal distribution.
Theorem 2 implies the following Corollary 3.

Corollary 3 Under the Conditions (C1)-(C8) in the "Appendix
Remark 3 Theorem 1 shows the proposed estimator α j (·), j = 1, . . . , p, is uniformly convergent. Theorem 2 shows the asymptotic bias of ( , and the asymptotic variance is of order O((nh) −1 ). As a consequence, the theoretical optimal bandwidth is of order n −1/5 . These results are similar to that of Zhou et al. (2018). The asymptotic properties presented above can be viewed as a generalization of the results in Zhou et al. (2018) from mean function estimation to varying-coefficient function estimation. These results will be more widely used in solving various practical problems than that of Zhou et al. (2018). Especially, it can provide theoretical guidance for analyzing the dynamic relationship between a set of covariates and the functional response variable, rather than only exploring the feature of functional variables.

Simulation studies
In this section, we are going to perform some simulation studies to illustrate the finite sample properties of the proposed estimation method. In consideration of the estimation procedure we proposed is aimed at estimating the varying-coefficient functions, we will use the accuracy of the varying-coefficient function estimators to demonstrate how well the procedure works. The performances are evaluated via the bias, standard deviation (SD), and root mean squared error (RMSE), at given grid points, which also are used by Zhou et al. (2018). Specifically, the bias, SD and RMSE of an estimator α(·) are respectively defined as is estimated by the sample mean based on the N simulated data sets, and t j ( j = 1, . . . , n g ) are the grid points at which the function α(·) is evaluated. In addition, comparisons are also made between our method and the local linear method employed in Zhu et al.
(2012) (Local hereinafter). In the following simulated examples, the Epanechnikov kernel is used for nonparametric regression. We set grid points n g = 100 and compute the results based on N = 500 replications.
Example 4 Generate data from the following model and ξ i2 are independent random variables. We set (σ, λ 1 , λ 2 ) = (0.3, 1, 0.5), and the varying-coefficient functions, and eigenfunctions as follows: We first conducted simulation studies under the following two regular time points settings.
• Setting 2: n = 150 with various time points: m i = 7 and m i = 13.
Setting 1, two different sample size n = 200, 400, and a fixed time points m i = 5 are designed. Setting 2, sample size n = 150 is fixed with various time points settings m i = 7 and m i = 13, which guarantees the total number of observations is comparable to that in Setting 1, by decreasing the number of subjects and increasing the observation number from the same individual. Tables 1 and 2 report the numerical results for the proposed method and Local method under Setting 1 and Setting 2, respectively. Table 1 The simulation results for Setting 1 in Example 4. Bias, SD and RMSE of α 1 (·), α 2 (·) and α 3 (·) using the proposed method and the Local method Prop.  Table 2 The simulation results for Setting 2 in Example 4. Bias, SD and RMSE of α 1 (·), α 2 (·) and α 3 (·) using the proposed method and the Local method Prop. Example 5 To illustrate that the proposed procedure is applicable even to the case that the observation times for each trajectory are randomly missed, we replicated Example 4 with an irregular time points design used in Fan, Huang and Li (2007). Specifically, we will generate the observation times randomly, assuming that each individual has a set of "scheduled" observation points, {0, 1, 2, . . . , 11}, and each scheduled time, except time 0, has a 20% probability of being skipped, then a uniform [0, 1] random variable is added to a nonskipped scheduled time so that the actual observation time is in fact a random perturbation of a scheduled time. Table 3 presents the numerical results for this case with sample size n = 200 and 400.
Tables 1, 2 and 3 summarize the numerical performance of the proposed method and the Local method through the Bias, SD and RMSE. It is shown from Tables 1, 2 and 3 that the RMSE of the proposed method is consistently smaller than that of the Local method for all varying-coefficient function estimators α 1 (t), α 2 (t) and α 3 (t), which indicates that the proposed method performs better than the Local method. Moreover, the proposed method has less bias and variance than the Local method, mainly because of the incorporation of the correlation information and the modified local linear technique. Tables 1, 2 and 3 also show that the Bias, SD and RMSE values decrease as the sample size increases, which corroborates with our consistency and asymptotic normality established in Sect. 3. From Tables 1 and 2, we find that the resulting estimators obtained in Example 4 with Setting 1 have smaller RMSE than those obtained in Example 4 with Setting 2 when the total number of observations for these two settings is close. This finding implies that the sparsity is beneficial to the performance of the proposed method and the estimators of the varying-coefficient functions can be improved with increasing the number of subjects. For all simulation cases, we chose h 0 = 0.12 as the bandwidth for the initial estimates of the varyingcoefficient functions and the estimates of the covariance functions. We discover that the proposed method is not sensitive to the initial values. Furthermore, the optimal bandwidth for the proposed method is slightly larger than that for the Local method. This is likely because the correlation matrix is incorporated into the estimation process Table 3 The simulation results for Example 5. Bias, SD and RMSE of α 1 (·), α 2 (·) and α 3 (·) using the proposed method and the Local method Prop.  of the proposed method, which can be treated as the weight as well as the kernel matrix, and further lower data usage. Therefore, the larger bandwidth for the proposed method is needed to enlarge the included range of data. We also apply the PACE method to Example 5 to conduct the estimation of η i (t). A comparison is made between the PACE method and the local linear approach. The PACE method is implemented with the FPCA function available in the package fdapace (Dai et al. 2019). We choose the first four functional principal components to construct the estimators of η i (t) (i = 1, . . . , n) for the PACE method. The bandwidth for the local linear method is selected by generalized cross-validation (GCV). We evaluate the performances through the root average square errors (RASE) given by where t * i j ( j = 1, . . . , n * g ) are the grid points at which the function η i (·) is evaluated. Since the observation points for each subject are sparse, we set n * g = 50 to calculate the RASE. Table 4 presents the mean and standard deviation (SD) of RASE, both of which are computed based on 200 simulation replications. We can see from Table 4 that all the mean and SD of RASE for the PACE method are significantly smaller than that of the local linear method. These results imply that the PACE method is superior to the local linear method for the estimation of η i (t). Moreover, we observe that the mean value of RASE for the PACE method decreases as the sample size increases, while the value for the local linear method has hardly changed. This is likely due to the fact that the PACE method for constructing the estimators η i (t) (i = 1, . . . , n) is mainly based on the estimated population covariance function of η i (t), which applies the whole sample information in estimation. Nevertheless, the local linear method for directly estimating each η i (t) only utilizes the observations on the grid of its corresponding curve, therefore, the resulting estimators could not be improved by increasing the size of samples.
To have a more visual display about how well the proposed estimation procedure works, we choose Example 5 with sample size n = 200 as an example. Figure 1 exhibits the averaged estimated varying-coefficient functions based on 50 replications and the boxplot of RASE of the varying-coefficient functions over 500 simulations. It can be seen from Fig. 1a-c that the estimated varying-coefficient curves (broken lines), obtained by the proposed method, closely resemble the corresponding true curves (solid lines), reflecting the good performance of the proposed procedure. Figure 1d shows the proposed method performs better than the Local method in terms of RASE.

Real data analysis
We apply model (2) and the newly proposed procedure to analyze a real diffusion tensor imaging (DTI) dataset with n = 213 subjects that was collected from NIH Alzheimer's Disease Neuroimaging Initiative (ADNI) study. One goal of the NIH ADNI is to test whether genetic, structural and functional neuroimaging, and clinical data can be integrated to evaluate the progression of mild cognitive impairment (MCI) and early Alzheimer's disease (AD). For more details on this dataset, see the ADNI publicly available database (http://adni.loni.usc.edu/). In diffusion data analysis, fractional anisotropy (FA), which quantifies the directional strength of white matter tract structure, at a particular location, is one of the most used measures, and has been widely applied to statistical analyses in many imaging studies. During the stage of data processing, we calculate the FA curve for each subject using a Tract-Based Spatial Statistics (TBSS) pipeline analysis (Smith et al. 2006) and the FMRIB software library tool. The more details for this process, one can refer to Li et al. (2017).
In this study, 213 subjects with age of 48.4-73.2 years old are sampled to assess the impact of some indicators which are commonly occurred in diffusion data, such as gender (123 male and 90 female), handiness (193 righthand and 20 left-hand), education level (years, ranges from 9 to 20 (years), mean 15.91), Alzheimer's disease (AD) status (19.6%), mild cognitive impairment (MCI) status (55.1%) and Mini-Mental State Exam (MMSE) of ADNI on the FA. We aim to evaluate the effects of these indicators on the FA using model (2). Specifically, the response Y i (t) represents the FA curve and the covariate vector x i includes an intercept term, the gender variable coded by a dummy variable indicating for male, the age of the subject, a handiness indicator coded by a dummy variable indicating for left-hand, the education level, an indicator for AD status, an indicator for MCI status and MMSE. We randomly selected 15 observations from the 83 grid points of each FA curve, so that we can compare the results obtained from the sparse data and the complete data and further validate the new approach. We standardize all variables to be mean zero and variance one, and apply the proposed procedure established in Sect. 2 to estimate the corresponding coefficients. Figure 2 shows the estimated varying-coefficient functions corresponding to intercept, gender, age, education, AD status, MCI status and MMSE indicator and the associated 95% pointwise confidence intervals for both sparse and complete data. The pointwise confidence intervals were obtained with the standard error based on 200 bootstrap resampling. The dashed lines are the estimates fulfilled by the proposed method for sparse data and the solid lines are estimated by the local linear method for complete data. Figure 2 indicates that the proposed estimation method works well even though data are sparse; it also reveals that age, handiness, education, AD and MCI have more effects on FA than that of gender and MMSE. The intercept function characterizes the nonlinear overall trend of FA. It can be seen from Fig. 2 that the coefficient curve for gender is positive and that for MMSE is near to the zero line, while the coefficients for the other five covariates of interest are negative, at most of the grid points. These findings may indicate that males tend to have higher FA values than females, MMSE has no remarkable effects on FA, AD and MCI patients tend to have lower FA values and being older, left-hand, or more educated may lead to smaller FA values. These discoveries are also consistent with the previous analysis in Li et al. (2017). Finally, we calculated the squared prediction error (i.e., i is the predicted trajectory for the ith subject by using the data with the ith subject being pulled out.) both for the proposed method and the local linear method. The corresponding SPE values for these two approaches are respectively 18.34 and 18.42, which implies that the proposed method has a smaller prediction error than that of the local linear method for analyzing this data. Figure 3 presents the estimated eigenvalues and eigenfunctions for the covariance function (s, t) with sparse data. The top four nonzero eigenvalues are 0.0028, 0.00047, 0.0003 and 0.00016 respectively. The first four eigenvalues explain about 94% of the total variability, while, it can be seen from Fig. 3 that the remaining eigenvalues rapidly drop to zero. The first eigenfunction, with a dominant eigenvalue accounting for 70% of the total variation, varies simply in structure. While the other eigenfunctions also have quite simple structures and most of them are roughly sinusoidal. The estimated varying coefficients for intercept, sex, age, handiness, education, AD, MCI and MMSE and their 95% pointwise confidence intervals for the ADNI data. The dashed lines are the estimates fulfilled by our proposed estimation procedure with sparse data and the solid lines are estimated by the local linear method with complete data

Conclusions
In this paper, we used a functional varying-coefficient mixed effects model to investigate the dynamic association between a set of covariates of interest and the functional response. An iterative estimation procedure was developed to estimate the varyingcoefficient functions in this model for longitudinal or sparse functional data. The proposed estimation procedure not only employs the idea that uses all observations instead of the local observations near the target point for estimation but also incorporates the estimated covariance structure into the estimation process to improve the resulting estimators. The asymptotic properties including uniform consistency and  Fig. 3 The estimated eigenvalues in decreasing order (left), the estimated eigenfunctions (middle) and the corresponding cumulate FVE (right) for ADNI data analysis asymptotic normality for the proposed estimators were also established. Finally, we compared the proposed procedure with the local linear method numerically through several simulation studies. The simulation results reveal that our method performed better than the local linear method in terms of the root mean squared error. A real data example is also provided to illustrate our approach.
Our future research will be focused on addressing several important issues. First, it would be interesting to gain a more efficient estimator for the covariance function by using the newly proposed estimation procedure. Second, applying the proposed estimation procedure to investigate other more complex regression models, such as functional index models (Jiang and Wang 2011), single-index varying coefficient (SIVC) models (Luo et al. 2016) and functional additive models (Zhang, Park and Wang 2013), is also a meaningful issue. Third, it is also worthwhile to conduct statistical inferences such as hypothesis tests for longitudinal and functional data.
proposed estimators, are imposed to facilitate the proof and are similar to Zhou et al. (2018). They are not the weakest possible conditions, but mainly for mathematical simplicity and may be modified if necessary.
(C1) The kernel K (·) is bounded symmetric probability density function with the compact support [−1, 1] and satisfy ν r = u r K (u)du < ∞, υ r = u r K 2 (u)du < ∞. (C2) max 1<i<n m i ≤ M, where M is a constant, and the dimension of covariate vector is bounded. (C3) t i j , i = 1, . . . , n, j = 1, . . . , m i are independent and identically distributed copies of a random variable T defined on [0,1]. Its density function f (·) is positive and has continuous second derivatives. (C4) The covariate vectors x i 's are independent and identically distributed with and γ 2 (t, z) are defined in the following.
Conditions (C1)-(C5) are mild assumptions for the kernel function, covariates and the function of interest, which are commonly used in many longitudinal literature. Condition (C6) guarantees the convergence of the initial estimators which can be easily achieved by the exiting kernel and spline estimation method. Condition (C7) is the regular condition for the bandwidth, which also has been used in Chen et al. (2010) and Zhou et al. (2018). Condition (C8) ensures identifiable.
Lemma 1 Suppose that Conditions (C1)-(C7) hold, ζ l (·) ∈ F 1 is any bounded function having continuous third derivatives on its compact support, where l = 1, . . . , p. Then, for j = 1, . . . , p, The proof of this lemma is similar to that of Theorem 1 in Zhou et al. (2018) and follows from Theorem 37 and Example 38 in Chapter 2 of Pollard (1984).