# Quadratic inference with dense functional responses

Pratim Guha Niyogi  Ping-Shou Zhong  Department of Biostatistics, Johns
Hopkins Bloomberg School of Public Health, Maryland, USA  Department of
Mathematics, Statistics and Computer Science, University of Illinois at
Chicago, Illinois, USA

###### Abstract

We address the challenge of estimation in the context of constant linear
effect models with dense functional responses. In this framework, the
conditional expectation of the response curve is represented by a linear
combination of functional covariates with constant regression parameters. In
this paper, we present an alternative solution by employing the quadratic
inference approach, a well-established method for analyzing correlated data,
to estimate the regression coefficients. Our approach leverages non-
parametrically estimated basis functions, eliminating the need for choosing
working correlation structures. Furthermore, we demonstrate that our method
achieves a parametric $\sqrt{n}$-convergence rate, contingent on an
appropriate choice of bandwidth. This convergence is observed when the number
of repeated measurements per trajectory exceeds a certain threshold,
specifically, when it surpasses $n^{a_{0}}$, with $n$ representing the number
of trajectories. Additionally, we establish the asymptotic normality of the
resulting estimator. The performance of the proposed method is compared with
that of existing methods through extensive simulation studies, where our
proposed method outperforms. Real data analysis is also conducted to
demonstrate the proposed method.

###### keywords:

Constant Linear-Effect Models, Functional Principal Component Analysis,
Quadratic Inference, Semi-parametric Functional Regression.

##  1 Introduction

Longitudinal data analysis (LDA) involves tracking repeated measurements on
the same individuals over time, allowing us to study changes over time and
identify influencing factors. Unlike cross-sectional studies, which capture
only “between-individual” responses at a single time point, LDA can capture
“within-individual” changes through repeated measurements. Longitudinal data
is often observed in clusters, with each cluster representing measurements
from one individual. As data complexity increases with advancing technology,
functional data analysis (FDA) has become a vital tool, extending our
understanding from finite to infinite dimensions. In LDA, the data are
generally observed with noise for measurements at each time-point [6, 14].
Moreover, a few repeated measurements are required in LDA, and the data are
observed sparsely with noise. On the other hand, in FDA, data are densely
observed as a continuous-time stochastic process without or with noise [34,
20]. Often, the sampling plan can have an effect on the performance of the
estimation procedures and inference [12]. In some situations, data are
typically functions by nature and are observed densely over time. Chiou et al.
[4] proposed a class of semi-parametric functional regression models to
describe the influence of vector-valued covariates on a sample of the response
curve. When data collection leads to experimental error, smoothing is
performed at closely spaced time-points in order to reduce the effect of
noise. The current developments of functional regression techniques have been
rigorously studied in Fan et al. [8], Hall et al. [11], Chen et al. [3].

The generalized estimating equation (GEE) technique proposed by Liang and
Zeger [21] has been extensively used in LDA for estimation of parameters.
Although an efficient technique, the GEE may not be efficient when the
correlation matrix is not correctly specified. Hence, without requiring the
estimation of the correlation parameters, the quadratic inference function
(QIF) approach proposed by Qu et al. [26] is useful for parameter estimation
in longitudinal studies [6] and cluster randomized trials [28]. By
representing the inverse of the working correlation matrix in terms of linear
combinations of the basis matrices and involving multiple sets of score
functions, the QIF approach has improved efficiency over GEE when the working
correlation matrix is not specified correctly. Although it maintains the same
efficiency as in the situation where the working correlation matrix is
specified correctly, the QIF method is not independent of the choice of the
working correlation matrix. A QIF method-based approach to varying-coefficient
models for longitudinal data was proposed by Qu and Li [25]. The related work
of Bai et al. [1] is an extension of QIF for the partial linear model. An
alternative method was presented in Yu et al. [31] where each set of score
equations was solved separately and their solutions were combined afterward;
thereby providing results on inference for an optimally weighted estimator and
extending those insights to the general setting with over-identified
estimating equations.

The fundamental limitations that all the above-mentioned powerful techniques
suffer from are: (1) all the above methodologies require prior information on
the working correlation structure or choose appropriate basis matrices, and
(2) the performance of the classical QIF approach is unknown for dense
functional data. Our study is motivated by problems from multiple real data
applications that involve dense functional data when information on the
working correlation structure is lacking. Let us discuss two motivating
examples that we will use to illustrate the proposed method in this paper (see
Section 5 for details). In the Beijing2017-data example, particulate matter
(PM) with a diameter of less than 2.5 micrometers is collected over different
time-points in different locations in China. Scientists are interested in
knowing the linear dependence of the pollution factor $\text{PM}_{2.5}$ with
other atmospheric chemicals [22]. Figure 1 (left) in Section 5 pictorially
demonstrates the readings of $\text{PM}_{2.5}$ for the given locations over
several hourly time-points; therefore, dense functional data analysis can be
implemented. In another example from a neuroimaging study Apnea-data,
scientists are interested in modeling the change of white matter structure
among voxels in each region of interest (ROI) of the human brain. Xiong et al.
[29] investigates white matter structural alterations using diffusion tensor
imaging (DTI) in obstructive sleep apnea (OSA) patients. Here, the change of
DTI parameters such as fractional anisotropy (FA) is investigated by using
constant linear effect models with the interaction of the count of lapses
obtained from the Psychomotor Vigilance Task and voxel locations as
predictors. We applied this model to each ROI and compared the results
obtained from all the ROIs.

We propose a data-driven way to select the working covariance matrix and
express the inverse of the covariance function in terms of the empirical
eigen-functions of the covariance operator. The covariance operator can be
estimated as in Hsing and Eubank [16] and other related methods based on
functional principal component analysis (FPCA) as found in Dauxois et al. [5],
Yao et al. [30], Hall and Hosseini-Nasab [12], Hall et al. [11], Li and Hsing
[19], Zhang and Li [32]. Note that the estimation of the eigen-functions is
nonparametric and it introduces some errors in the proposed estimation method.
In this article, we answer the following question: while we estimate the
eigen-functions nonparametrically from data, is the estimation of coefficient
vectors in a semi-parametric problem $\sqrt{n}$\- consistent in dense
functional data, and can we achieve asymptotic normality? The advantages of
our proposed method are the following: First, our method preserves the good
properties of the QIF method and is easier to implement since the eigen-
functions can be estimated using the existing packages in statistical software
such as R. Second, under some mild conditions, our proposed estimator can
obtain the optimal convergence rate and is asymptotically normally distributed
with less variance as compared to the classical QIF methods. Third, asymptotic
results show the estimation accuracy of the coefficients in semi-parametric
functional model, therefore, making the influence of the dimension reduction
step using FPCA redundant. The error in the estimation of the eigen-functions
contributes to the error in the estimation of the parameters. Under some mild
bandwidth conditions, the above-mentioned error contribution is of the same
order of magnitude as an error in parameter estimation if eigen-functions are
known in advance.

The rest of the paper is organized as follows. In Section 2, we introduce the
basic concept of QIF along with our proposed method. The asymptotic results
for the proposed estimator are presented in Section 3. In Section 4, we
demonstrate the performance for finite samples. We also apply the proposed
method to real data-sets in Section 5. We conclude with some remarks in
Section 6. All technical proofs are presented in the Appendix of the article.

##  2 Functional response model and estimation procedure

###  2.1 Basic model

To analyze longitudinal data, a straightforward application of a generalized
linear model [24] for single response variables is not applicable due to the
lack of independence between repeated measures. To account for the high
correlation in the longitudinal data, some special techniques are required. A
seminal work by Liang and Zeger [21] proposed the use of GLM for the analysis
of longitudinal data. The model we consider in this article is commonly
observed in spatial modeling, where associations among variables do not change
over the functional domain (see Zhang and Banerjee [33] and references
therein); which is termed as constant linear effects model. In this paper, the
variable time is used as a functional domain variable.

Let $y(t)$ be the response variable at time-point $t$ and $\mathbf{x}(t)$ be
$p$-dimensional covariates observed at time $t\in\mathcal{T}$ where without
loss of generality we assume $\mathcal{T}=\left[0,1\right]$, is the spectrum
of the time-points. The stochastic process $y(t)$ is square-integrable with
marginal mean $\mathbb{E}\\{y(t)\mid\mathbf{x}(t)\\}$ and finite covariance
function; the regression parameter $\boldsymbol{\beta}$ is unknown and is to
be efficiently estimated. Thus, linear models with longitudinal data have the
following expression.

$y(t)=\mathbf{x}(t)^{\text{T}}\boldsymbol{\beta}+e(t),$ (1)

where the stochastic process $y(t)$ is decomposed into two parts: one is the
mean function $\mu(t)=\mathbf{x}(t)^{\text{T}}\boldsymbol{\beta}$ that depends
on time-varying covariates and vector-valued coefficient vector
$\boldsymbol{\beta}$, and other is the random error part $e(t)$ where
$\mathbb{E}\\{e(t)\\}=0$ and has finite second-order covariance $R(s,t)$. Let
$y_{i}(t)$ be i.i.d. copies of the stochastic process and for each individual,
the measurements are taken on $m_{i}$ discrete time-points $T_{ij}$ for
$j=1,\cdots,m_{i};i=1,\cdots,n$. Therefore, at time $T_{ij}$, we observe a
$m_{i}\times 1$ response vector $y_{i}(T_{ij})$ and corresponding covariates
$\mathbf{x}_{i}(T_{ij})$ for the $i$-th subject. We assume that $m_{i}$’s are
all of the same order as $m=n^{a}$ for some $a\geq 0$, thus $m_{i}/m$ is
bounded below and above by some constants. Functional data are considered to
be sparse depending on the choice of $a$ [12]. Data with bounded $m$ or $a=0$
are called sparse functional data and if $a\geq a_{0}$ where $a_{0}$ is a
transition point, are called dense functional data. Moreover, the regions
$(0,a_{0})$ are sometimes referred to as moderately dense. Furthermore, we
denote $y_{ij}$ and $\mathbf{x}_{ij}$ as $y_{i}(T_{ij})$ and
$\mathbf{x}_{i}(T_{ij})$ respectively. $(y_{i1},\cdots,y_{im_{i}})^{\text{T}}$
and $(\mu_{i1},\cdots,\mu_{im_{i}})^{\text{T}}$ are $m_{i}$ component vectors,
denoted as $\mathbf{y}_{i}$ and $\boldsymbol{\mu}_{i}$ respectively. The
derivative of $\boldsymbol{\mu}$, denoted as $\dot{\boldsymbol{\mu}}$, is a
$m_{i}\times p$ matrix.

In the classical problem of GEE, we estimate $\boldsymbol{\beta}$ by solving
the quasi-likelihood equations [21]:

$\sum_{i=1}^{n}\dot{\boldsymbol{\mu}}_{i}^{\text{T}}\mathbf{V}_{i}^{-1}(\mathbf{y}_{i}-\boldsymbol{\mu}_{i})=0.$
(2)

We denote
$\mathbf{V}_{i}=\nu\mathbf{A}_{i}^{1/2}\mathbf{R}_{i}(\rho)\mathbf{A}_{i}^{1/2}$
where $\mathbf{R}_{i}(\rho)$ is the working correlation matrix, $\nu$ is an
over-dispersion parameter and $\mathbf{A}_{i}$ is a diagonal matrix where
entries are marginal variances $\text{Var}(y_{i1})$, $\cdots$,
$\text{Var}(y_{im_{i}})$. In this article, we simply set $\nu=1$ while the
extension to a general $\nu$ is straightforward. In practice, the prior
knowledge of the working correlation matrix is unknown, and the estimation of
the coefficient is influenced by its choice. Therefore, Qu et al. [26]
suggested an expansion of the inverse of the working correlation matrix as
$\mathbf{R}(\rho)^{-1}=\sum_{k=1}^{\kappa_{0}}a_{k}(\rho)\mathbf{M}_{k}$ where
$\mathbf{M}_{k}$ are some basis matrices. Zhou and Qu [35] modified linear
representation by grouping the basis matrices into an identity matrix and some
symmetric basis matrices. For example, if the working correlation matrix is
exchangeable/compound symmetric,
$\mathbf{R}(\rho)^{-1}=c_{1}\mathbf{I}_{m}+c_{2}\mathbf{J}_{m}$ where
$\mathbf{I}_{m}$ is the $m\times m$ identity matrix and $\mathbf{J}_{m}$ is
the $m\times m$ matrix such that $0$ is in diagonal and $1$ is in off-diagonal
positions. On the other hand, for first-order auto-regressive correlation
matrix,
$\mathbf{R}(\rho)^{-1}=c_{1}\mathbf{I}_{m}+c_{2}\mathbf{J}_{m}^{(1)}+c_{3}\mathbf{J}^{(2)}_{m}$
where $\mathbf{J}^{(1)}_{m}$ is a matrix with $1$ on the two main off-diagonal
positions and $0$ otherwise, $\mathbf{J}^{(2)}_{m}$ is a matrix such that $1$
is in the corner positions, viz. $(1,1)$ and $(m,m)$ and $0$ elsewhere. Here,
$c_{k}$s are real constants that depend on the nuisance parameter $\rho$.
Therefore, Equation (2) reduces to the linear combination of the score
vectors:

$\displaystyle\bar{\mathbf{g}}(\beta)$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\mathbf{g}_{i}(\beta)=\begin{bmatrix}n^{-1}\sum_{i=1}^{n}\dot{\boldsymbol{\mu}}^{\text{T}}_{i}\mathbf{A}_{i}^{-1/2}\mathbf{M}_{1}\mathbf{A}_{i}^{-1/2}(\mathbf{y}_{i}-\boldsymbol{\mu}_{i})\\\
\vdots\\\
n^{-1}\sum_{i=1}^{n}\dot{\boldsymbol{\mu}}^{\text{T}}_{i}\mathbf{A}_{i}^{-1/2}\mathbf{M}_{\kappa_{0}}\mathbf{A}_{i}^{-1/2}(\mathbf{y}_{i}-\boldsymbol{\mu}_{i})\\\
\end{bmatrix}.$ (3)

Due to the higher dimension of $\overline{\mathbf{g}}$, Qu et al. [26] used
the generalized method of moments (GMM) [15] for which the method of
estimation boils down to minimization of the quadratic inference function
$\mathcal{Q}(\boldsymbol{\beta})=n\bar{\mathbf{g}}(\boldsymbol{\beta})^{\text{T}}\widehat{\mathbf{C}}(\boldsymbol{\beta})^{-1}\bar{\mathbf{g}}(\boldsymbol{\beta})$
where
$\widehat{\mathbf{C}}(\boldsymbol{\beta})=n^{-1}\sum_{i=1}^{n}\bar{\mathbf{g}}_{i}(\boldsymbol{\beta})\bar{\mathbf{g}}_{i}(\boldsymbol{\beta})^{\text{T}}$
is the sample covariance matrix of Equation (3). To obtain the solution of
$\boldsymbol{\beta}$, the Newton-Raphson method is used which iteratively
updates the value of $\boldsymbol{\beta}$.

###  2.2 Incorporating eigen-functions in QIF

Now, due to standard Karhunen-Loève expansions of
$e_{i}(t)=y_{i}(t)-\mu_{i}(t)$ [18, 23]

$e_{i}(t)=\sum_{r=1}^{\infty}\xi_{ir}\phi_{r}(t),$ (4)

with probability 1, where independently distributed random variables
$\xi_{ir}=<e_{i},\phi_{r}>\sim(0,\lambda_{r})$ for ordered eigen-values
$\lambda_{r}$ such that $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq 0$ and
$\phi_{r}$s are orthonormal eigen-functions such that
$\int\phi_{r}(t)\phi_{l}(t)=\textbf{1}(r=l)$. Due to Mercer’s theorem [17], a
symmetric continuous non-negative definite kernel function $R$ has the
representation $R(s,t)=\sum_{r=1}^{\infty}\lambda_{r}\phi_{r}(s)\phi_{r}(t)$
where the sum is absolutely and uniformly continuous. We extract the main
directions of the variation of the response variables using FPCA. In this
situation, we take the first $\kappa_{0}$ terms which provide a good
approximation of the infinite sum in Equation (4) by considering that the
majority of the variations in the data are contained in the subspace spanned
by few eigen-functions [3]. For finite $\kappa_{0}\geq 1$, we therefore
consider the rank-$\kappa_{0}$ FPCA model,

$\mathbb{E}\\{y(t)|\mathbf{x}(t)\\}=\mu(t)+\sum_{r=1}^{\kappa_{0}}\mathbb{E}\\{\xi_{r}|\mathbf{x}(t)\\}\phi_{r}(t).$
(5)

Unlike existing models in Chiou et al. [4] where the main interest of the work
is to study how the vector covariates influence the whole response curve, we
study point-wise effect of covariates on functional response the influence of
covariates vectors indirectly on the functional response through the
regression of scores on covariates, viz., $\mathbb{E}\\{\xi_{r}|x(t)\\}$. An
analogue of the truncated empirical version of Equation (4) and Mercer’s
representation can be provided easily, and we discuss the proposed method
based on this truncated version. Moreover, we discuss how to choose
$\kappa_{0}$ in our situation in Section 3 in detail.

In this chapter, we propose a data-driven way to compute the basis matrices to
obtain the approximate inverse of $\mathbf{V}$ as discussed earlier. In this
approach, it is enough to find the eigen-functions to construct a GEE. Let us
define

$\bar{\mathbf{g}}(\boldsymbol{\beta})=\begin{bmatrix}n^{-1}\sum_{i=1}^{n}\dot{\boldsymbol{\mu}}^{\text{T}}_{i}\widehat{\boldsymbol{\Phi}}_{i1}(\mathbf{y}_{i}-\boldsymbol{\mu}_{i})\\\
\vdots\\\
n^{-1}\sum_{i=1}^{n}\dot{\mu}^{\text{T}}_{i}\widehat{\boldsymbol{\Phi}}_{i\kappa_{0}}(\mathbf{y}_{i}-\boldsymbol{\mu}_{i})\\\
\end{bmatrix},$ (6)

where for $k=1,\cdots,\kappa_{0}$, we define
$\widehat{\boldsymbol{\Phi}}_{ik}=\left(m_{i}^{-2}\widehat{\phi}_{k}(T_{ij})\widehat{\phi}_{k}(T_{ij^{\prime}})\right)_{j,j^{\prime}=1,\cdots,m_{i}}$.
Since the dimension of $\overline{\mathbf{g}}$ in Equation (6) is greater than
the number of parameters to estimate, instead of setting
$\overline{\mathbf{g}}$ to zero, we minimize quadratic function.
$\widehat{\boldsymbol{\beta}}=\arg\min_{\boldsymbol{\beta}}\mathcal{Q}(\boldsymbol{\beta})~{}\text{where}~{}~{}\mathcal{Q}(\boldsymbol{\beta})=n\bar{\mathbf{g}}(\boldsymbol{\beta})^{\text{T}}\widehat{\mathbf{C}}(\boldsymbol{\beta})^{-1}\bar{\mathbf{g}}(\boldsymbol{\beta})$
where,
$\widehat{\mathbf{C}}(\boldsymbol{\beta})=n^{-1}\sum_{i=1}^{n}\mathbf{g}_{i}(\boldsymbol{\beta})\mathbf{g}_{i}(\boldsymbol{\beta})^{\text{T}}$.
For the existence of $\widehat{\mathbf{C}}^{-1}$ we need the additional
restriction: $n\geq dim(\mathbf{g}_{i})=p\times\kappa_{0}$ where $\kappa_{0}$
is the number of eigen-functions. Under the given set-up, by Equation (8) in
Qu et al. [26] the estimating equation for $\boldsymbol{\beta}$ will be
$\dot{\mathcal{Q}}(\boldsymbol{\beta})\approx
2\dot{\bar{\mathbf{g}}}(\boldsymbol{\beta})^{\text{T}}\widehat{\mathbf{C}}(\boldsymbol{\beta})^{-1}\bar{\mathbf{g}}(\boldsymbol{\beta})$.
To obtain the solution of the above equation, we use a Newton-like method. In
practice, the standard Newton method does not lead to a decrease in the
objective function, that is, at each step of the iteration, there is no
guarantee that
$\mathcal{Q}(\boldsymbol{\beta}_{s+1})<\mathcal{Q}(\boldsymbol{\beta}_{s})$.
Therefore, we use the Algorithm 1 to estimate $\boldsymbol{\beta}$ using the
Quasi-Newton method with halving [10].

Algorithm 1 Estimation of $\boldsymbol{\beta}$ using the Quasi-Newton method
with halving.

Set $\widehat{\boldsymbol{\beta}}_{0}$ (initial estimates)

Set $\epsilon_{0}$ (threshold, a small number)

Set max.count (maximum number of repetition)

Calculate
$\mathcal{Q}(\widehat{\boldsymbol{\beta}_{0}}),\dot{\mathcal{Q}}(\widehat{\boldsymbol{\beta}_{0}}),\ddot{\mathcal{Q}}(\widehat{\boldsymbol{\beta}_{0}})$
respectively.

while $\text{Error}>\epsilon_{0}$ do

Calculate $\dot{\mathcal{Q}}(\widehat{\boldsymbol{\beta}}_{1})$ and
$\ddot{\mathcal{Q}}(\widehat{\boldsymbol{\beta}}_{1})$ based on
$\widehat{\boldsymbol{\beta}}_{1}$

Initialise $r_{0}=1$

$\widehat{\boldsymbol{\beta}}_{2}\leftarrow\widehat{\boldsymbol{\beta}}_{1}-r_{0}\ddot{\mathcal{Q}}(\widehat{\boldsymbol{\beta}}_{1})^{-1}\dot{\mathcal{Q}}(\widehat{\boldsymbol{\beta}}_{1})$

Calculate ${\mathcal{Q}}(\widehat{\boldsymbol{\beta}}_{1})$ and
${\mathcal{Q}}(\widehat{\boldsymbol{\beta}}_{2})$ based on
$\widehat{\boldsymbol{\beta}}_{1}$ and $\widehat{\boldsymbol{\beta}}_{2}$
respectively using proposed method

while
$\mathcal{Q}(\widehat{\boldsymbol{\beta}}_{2})>\mathcal{Q}(\widehat{\boldsymbol{\beta}}_{1})$
do

$r_{0}\leftarrow r_{0}/2$

$\widehat{\boldsymbol{\beta}}_{2}\leftarrow\widehat{\boldsymbol{\beta}}_{1}-r_{0}\ddot{\mathcal{Q}}(\widehat{\boldsymbol{\beta}}_{1})^{-1}\dot{\mathcal{Q}}(\widehat{\boldsymbol{\beta}}_{1})$

Calculate $\mathcal{Q}(\widehat{\boldsymbol{\beta}}_{1})$ and
$\mathcal{Q}(\widehat{\boldsymbol{\beta}}_{2})$ based on
$\widehat{\boldsymbol{\beta}}_{1}$ and $\widehat{\boldsymbol{\beta}}_{2}$
respectively using proposed method

end while

Calculate
$\text{Error}=\|\widehat{\boldsymbol{\beta}}_{2}-\widehat{\boldsymbol{\beta}}_{1}\|^{2}$

$\widehat{\boldsymbol{\beta}}_{0}\leftarrow\widehat{\boldsymbol{\beta}}_{1}$

$\widehat{\boldsymbol{\beta}}_{1}\leftarrow\widehat{\boldsymbol{\beta}}_{2}$

end while

Output $\boldsymbol{\beta}$

###  2.3 Estimation of eigen-functions

Estimation of eigen-functions is an important step in our proposed quadratic
inference technique. In general, FPCA plays an important role as a dimension-
reduction technique in functional data analysis. Some important theories on
FPCA have been developed in recent years. In particular, Hall and Hosseini-
Nasab [12] proved various asymptotic expressions for FPCA for densely observed
functional data. Later, Hall and Hosseini-Nasab [13] showed more common
theoretical arguments, including the effect of the gap between eigen-value
(a.k.a., spacing) on the property of eigen-value estimators. In Li and Hsing
[19], uniform rates of convergence of the mean and covariance functions are
given, which are equipped for all possible choices/scenarios of $m_{i}$s.

Note that the error process $e(t)$ has mean zero, defined on compact set
$\mathcal{T}=[0,1]$ satisfying
$\int_{\mathcal{T}}\mathbb{E}\\{e^{2}\\}<\infty$. The functional principal
components can be constructed via the covariance function $R(s,t)$ (induces
the kernel operator $\mathcal{F}$) defined as
$R(s,t)=\mathbb{E}\\{e(s)e(t)\\}$ which is assumed to be square-integrable. An
empirical analogue of the spectral decomposition of $R$ can be obtained
$\widehat{R}(s,t)=\sum_{r=1}^{\infty}\widehat{\lambda}_{r}\widehat{\phi}_{r}(s)\widehat{\phi}_{r}(t)$
where the random variables
$\widehat{\lambda}_{1}\geq\widehat{\lambda}_{2}\geq\cdots\geq 0$ are the
eigen-values of the estimated operator $\widehat{\mathcal{F}}$ and the
corresponding sequence of eigen-functions are
$\widehat{\phi}_{1},\widehat{\phi}_{2},\cdots$. Further, assume that
$\int_{\mathcal{T}}\phi_{r}\widehat{\phi}_{r}\geq 0$ to avoid the issue
regarding change of sign [12] for practical comparison of eigen-functions,
otherwise there is no impact on the convergence rate of eigen-functions and
hence the proposed estimators.

Suppose that $T_{ij}$ are observational points with a positive density
function $f_{T}(\cdot)$. Assume $m_{i}\geq 2$ and define
$N=\sum_{i=1}^{n}N_{i}$ where $N_{i}=m_{i}(m_{i}-1)$. This approach is based
on local linear smoother which is popular in functional data analysis,
including Fan and Gijbels [7], Li and Hsing [19] among many others. Let
$K(\cdot)$ be a symmetric probability density function on $[-1,1]$, which is
used as kernel and $h>0$ be bandwidth, thus re-scaled kernel function is
defined as $K_{h}(\cdot)=h^{-1}K(\cdot)$. Therefore, for given
$s,t\in\mathcal{T}$, choose
$(\widehat{a}_{0},\widehat{b}_{1},\widehat{b}_{2})$ be the minimizer of the
following equation.

$\displaystyle n^{-1}\sum_{i=1}^{n}{N_{i}^{-1}}\underset{j_{1}\neq
j_{2}}{\sum_{j_{1}=1}^{m_{i}}\sum_{j_{2}=1}^{m_{i}}}\\{e_{i}(T_{ij_{1}})e_{i}(T_{ij_{2}})-a_{0}-b_{1}(T_{ij_{1}}-s)-b_{2}(T_{ij_{2}}-t)\\}^{2}$
$\displaystyle\qquad\times
K_{h}\left({T_{ij_{1}}-s}\right)K_{h}\left({T_{ij_{2}}-t}\right).$ (7)

Thus we estimate $R(s,t)=\mathbb{E}\\{e(s)e(t)\\}$ using the quantity
$\widehat{a}_{0}$, viz., $\widehat{R}(s,t)=\widehat{a}_{0}$. The operator
$\widehat{\mathcal{F}}$ is in general positive semi-definite and the estimated
eigen-values $\widehat{\lambda}_{r}$ are non-negative; indeed, $\widehat{R}$
is symmetric.

##  3 Asymptotic properties

In this section, we study the asymptotic properties of the proposed estimator.
Let us introduce some notations. Assume that $m_{i}$s are all of the same
order, viz, $m\equiv m(n)=n^{a}$ for some $a\geq 0$. Define,
$d_{n1}(h)=h^{2}+h\overline{m}/m$ and
$d_{n2}(h)=h^{4}+h^{3}\overline{m}/m+h^{2}\overline{\overline{m}}/m^{2}$ where
$\overline{m}=\lim_{n\rightarrow\infty}n^{-1}\sum_{i=1}^{n}m/m_{i}$ and
$\overline{\overline{m}}=\lim_{n\rightarrow\infty}n^{-1}\sum_{i=1}^{n}(m/m_{i})^{2}$.
Denote $\delta_{n1}(h)=\left\\{d_{n1}(h)\log n/(nh^{2})\right\\}^{1/2}$ and
$\delta_{n2}(h)=\left\\{d_{n2}(h)\log n/(nh^{4})\right\\}^{1/2}$. Further,
$\nu_{a,b}=\int t^{a}K^{b}(t)dt$. Define,
$\mathbf{W}=(\boldsymbol{\phi}(t_{1})^{\text{T}},\cdots,\boldsymbol{\phi}(t_{m})^{\text{T}})^{\text{T}}$
is a matrix of order $m\times\kappa_{0}$ obtaining after stacking all
$\boldsymbol{\phi}_{k}$s and random components
$\xi_{i}=(\xi_{i1},\cdots,\xi_{i\kappa_{0}})^{\text{T}}$. Further,
$\boldsymbol{\xi}$ has mean zero and variance $\boldsymbol{\Lambda}$ which is
a diagonal matrix with components $\lambda_{1},\cdots,\lambda_{\kappa_{0}}$.
The sign $``\lesssim"$ indicates that the left-hand side of the inequality is
bounded by the right-hand upto a multiplicative positive constant, i.e. for
two positive variables $f_{1}$ and $f_{2}$ we define $f_{1}\lesssim f_{2}$ as
$f_{1}\leq Cf_{2}$ where $C$ is a positive constant not involving $n$. The
following conditions are needed for further discussion of the asymptotic
properties.

  1. (C1)

Kernel function $K(\cdot)$ is a symmetric density function defined on bounded
support $[-1,1]$.

  2. (C2)

Density function $f_{T}$ of $T$ is bounded above and away from infinity. Also,
the density function is bounded below away from zero. Moreover, $f$ is
differentiable and the derivative is continuous.

  3. (C3)

$R(s,t)$ is twice differentiable, and all second-order partial derivatives are
bounded on $[0,1]^{2}$.

  4. (C4)

$\mathbb{E}\\{\sup_{t\in[0,1]}|e(t)|^{\gamma}\\}<\infty$ and
$\mathbb{E}\\{\sup_{t\in[0,1]}|\mathbf{x}_{i}(t)|^{2\gamma}\\}<\infty$ for
some $\gamma\in(4,\infty)$.

  5. (C5)

$h\rightarrow 0$ as $n\rightarrow\infty$ such that $d_{n1}^{-1}(\log
n/n)^{1-2/\gamma}\rightarrow 0$ and $d_{n2}^{-1}(\log
n/n)^{1-4/\gamma}\rightarrow 0$ for $\gamma\in(4,\infty)$.

  6. (C6)

Condition for eigen-components.

    1. (a)

for each $1\leq k<r<\infty$
${\max\left\\{\lambda_{k},\lambda_{r}\right\\}}/{|\lambda_{k}-\lambda_{r}|}\leq
C_{0}{\max\left\\{k,r\right\\}}/{|k-r|}$ for non-zero finite generic constant
$C_{0}$,.

    2. (b)

For some $\alpha>0$, with the condition
$V_{r}\lambda_{r}^{-2}r^{1+\alpha}\rightarrow 0$ as $r\rightarrow\infty$ where
$V_{r}=\mathbb{E}\\{\int\dot{\mu}(t)\phi_{r}(t)dt\\}^{2}$.

    3. The above two conditions hold if $\lambda_{r}=r^{-\tau_{1}}\Lambda(r)$ and $V_{r}=r^{-\tau_{2}}\Gamma(r)$ for slowly varying functions $\Lambda$ and $\Gamma$ where $\tau_{2}>1+2\tau_{1}>3$ and $\tau=\alpha+\tau_{1}$.

    4. (c)

$\int\phi_{k}^{4}(t)dt$ and $\int\dot{\mu}^{2}(t)\phi_{k}^{2}(t)dt$ are finite
for all $k\geq 1$. $\int\phi_{k}^{2}(t)\phi_{l}^{2}(t)dt$ is finite.

  7. (C7)

$\widehat{\mathbf{C}}(\boldsymbol{\beta})$ converges almost surely to an
invertible matrix
$\mathbf{C}_{0}=\mathbb{E}\\{\mathbf{g}(\boldsymbol{\beta}_{0})\mathbf{g}(\boldsymbol{\beta}_{0})^{-1}\\}$.

  8. (C8)

Conditions for $h$ and $\kappa_{0}$.

    1. (a)

If $a>1/4$, $\kappa_{0}=O(n^{1/(3-\tau)})$ and $n^{-1/4}\lesssim h\lesssim
n^{-(a+1)/5}$

    2. (b)

If $a\leq 1/4$, $\kappa_{0}=O(n^{4(1+a)/5(3-\tau)})$ and $h\lesssim n^{-1/4}$

######  Remark 1.

Condition (C1) is commonly used in non-parametric regression. The bound
condition for the density function of time-points has the standard Condition
(C2) for random design. Similar results can be obtained for a fixed design
where the grid-points are pre-fixed according to the design density
$\int_{0}^{T_{j}}f(t)dt=j/m$ for $j=1,\cdots,m$, for $m\geq 1$. Furthermore,
it is important to note that this approach does not involve the requirement to
obtain sample path differentiation when we invoke the estimation of
eigenfunctions from Li and Hsing [19]. Therefore the method could be
applicable for Brownian motion which has a continuously non-differentiable
sample path. To expand in Taylor series, Condition (C3) is required, and is
also common in non-parametric regression. Condition (C4) is required for a
uniform bound for certain higher-order expectations to show uniform
convergence. This is a similar condition adopted from Li and Hsing [19].
Smoothness conditions in (C5) and (C8) is common in kernel smoothing and
functional data to control bias and variance. The first condition for tuning
the parameters mentioned in (C5) is similar to Li and Hsing [19]. The required
spacing assumptions for eigen-values in Conditions (C6)a and (C6)b are similar
as in Hall and Hosseini-Nasab [13]. Condition (C6)c is a trivial assumption
that frequently arises in functional data analysis literature. In most of the
situations, this condition automatically holds. Using the weak law of large
numbers, Condition (C7) holds for large $n$. Similar kind of conditions can be
invoked such as convexity assumption, i.e.,
$\lambda_{r}-\lambda_{r+1}\leq\lambda_{r-1}-\lambda_{r}$ for all $r\geq 2$.
Condition (C8) is determined to control the rate of the number of repeated
measurements.

We present the following theorem to provide the asymptotic expansion and
consistency of the proposed estimator for $\widehat{\boldsymbol{\beta}}$.

######  Theorem 1.

Let $\boldsymbol{\beta}_{0}$ be the true value of $\boldsymbol{\beta}$. Under
the Conditions (C1)-(C6), for $k=1,\cdots,\kappa_{0}$, we have the asymptotic
mean square error for $\overline{\mathbf{g}}_{k}(\boldsymbol{\beta}_{0})$ as

$\text{AMSE}\\{\overline{\mathbf{g}}_{k}(\boldsymbol{\beta}_{0})\\}=O\left(n^{-1}+n^{-1}\kappa_{0}^{3-\tau}R_{n}(h)\right)\qquad\text{almost
surely}$ (8)

where
$R_{n}(h)=\left\\{h^{4}+(1/n)+(1/{nmh})+(1/{n^{2}m^{2}h^{2}})+(1/{n^{2}m^{4}h^{4}})+(1/{n^{2}mh})+(1/{n^{2}m^{3}h^{3}})\right\\}$.
Moreover, under Condition (C8),
$\text{AMSE}\\{\widehat{\mathbf{g}}_{k}(\boldsymbol{\beta}_{0})\\}=O(n^{-1})$.
Therefore, if in addition, Condition (C7) holds, as $n\rightarrow\infty$,
$\|\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}_{0}\|=O(n^{-1/2})$ in
probability.

The following theorem states the results of asymptotic normality.

######  Theorem 2.

Define
$\mathcal{C}_{i}=\sum_{k_{1}=1}^{\kappa_{0}}\sum_{k_{2}=1}^{\kappa_{0}}\boldsymbol{\Phi}_{k_{1}}\mathbf{X}_{i}\mathbf{C}_{k_{1},k_{2}}^{-1}\mathbf{X}_{i}^{\text{T}}\boldsymbol{\Phi}_{k_{2}}$,
where $\mathbf{C}_{k_{1},k_{2}}^{-1}$ is a $(k_{1},k_{2})$ block of
$\mathbf{C}_{0}^{-1}$ with
$\mathbf{C}_{0}=\mathbb{E}\left\\{\mathbf{g}(\boldsymbol{\beta}_{0})\mathbf{g}_{i}^{\text{T}}(\boldsymbol{\beta}_{0})\right\\}$.
Assume that the conditions for Theorem 1 holds. Then
$\sqrt{n}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}_{0})\xrightarrow{d}N(0,\boldsymbol{\Sigma})$
where $\boldsymbol{\Sigma}=\mathbf{B}^{-1}\mathbf{A}\mathbf{B}^{-1}$. where
the quantities $\mathbf{A}$ and $\mathbf{B}$ are, respectively, limits of
$\widehat{\mathbf{A}}=n^{-1}\sum_{i=1}^{n}\mathbf{X}_{i}^{\text{T}}\widehat{\mathcal{C}}_{i}\widehat{\mathbf{e}}_{i}\widehat{\mathbf{e}}_{i}^{\text{T}}\widehat{\mathcal{C}}_{i}\mathbf{X}_{i}$
and
$\widehat{\mathbf{B}}=n^{-1}\sum_{i=1}^{n}\mathbf{X}_{i}^{\text{T}}\widehat{\mathcal{C}}_{i}\mathbf{X}_{i}$,
and $``\xrightarrow{d}{}"$ denotes the convergence in distribution.

######  Remark 2.

Here, the selection of the bandwidth only affects the second-order term of the
MSE of $\widehat{\boldsymbol{\beta}}$ and has no effect on the asymptotic
result of normality as long as $h$ satisfies the Conditions (C5) and (C8)
along with some restrictions on $\kappa_{0}$.

Outline of the proofs of the proposed theorems is discussed in the Appendix
with additional technical details.

##  4 Simulation studies

We conduct extensive numerical studies to compare the finite sample
performance of our proposed method to that of Qu et al. [26] under different
correlation structures.

###  4.1 Simulation set-up

Consider the normal response model

$y_{i}(T_{ij})=\mathbf{x}_{i}(T_{ij})^{\text{T}}\boldsymbol{\beta}+e_{i}(T_{ij}).$
(9)

For $p=2$, we set coefficient vectors,
$\boldsymbol{\beta}=(\beta_{1},\beta_{2})^{\text{T}}$ where $\beta_{1}=1$ and
$\beta_{2}=0.5$. Consider the covariates
$x_{ik}(t)=\chi_{i1}^{(k)}+\chi_{i2}^{(k)}\sqrt{2}\sin\left(\pi
t\right)+\chi_{i3}^{(k)}\sqrt{2}\cos\left(\pi t\right)$ where the coefficients
$\chi_{i1}^{(k)}\sim N(0,(2^{-0.5(k-1)})^{2})$, $\chi_{i2}^{(k)}\sim
N(0,(0.85\times 2^{-0.5(k-1)})^{2})$, $\chi_{i3}^{(k)}\sim N(0,(0.7\times
2^{-0.5(k-1)})^{2})$ and $\chi_{ij}$s are mutually independent for each
trajectories $i$ and each $j$. In fixed design situations, associated
observational times are fixed. Sample trajectories are observed at $m=100$
equidistant time-points $\\{t_{1},\cdots,t_{m}\\}$ on $[0,1]$. Set number of
trajectories $n\in\\{100,300,500\\}$. The residual process $e_{i}(t)$ is a
smoothed function with mean zero and unknown covariance function, where each
$e_{i}$ is distributed as $e_{i}=\sum_{k\geq 1}\xi_{i}\phi_{i}$ and $\xi_{k}$s
are independent normal random variables with mean zero and respective
variances $\lambda_{k}$. For numerical computation, we truncate the finite
series at $k=3$ in Karhunen-Loève expansions for Situations (a), (b), and (c)
as described below. In Situations (d) and (e), the error process is generated
from given covariance functions.

  1. (a)

Brownian motion. The covariance function for the Brownian motion is
$\min(s,t)$, $\lambda_{k}=4/{\pi^{2}(2k-1)^{2}}$ and
$\phi_{k}(t)=\sqrt{2}\sin(t/\sqrt{\lambda_{k}})$.

  2. (b)

Linear process. Consider the eigen-values be $\lambda_{k}=k^{-2l_{0}}$ and
$\phi_{k}(t)=\sqrt{2}\cos(k\pi t)$. We fix $l_{0}\in\\{1,2,3\\}$.

  3. (c)

Ornstein Uhlenbeck (OU) process. For positive constants $\mu_{0}$ and
$\rho_{0}$, we have a stochastic differential equation for $e(t)$ as $\partial
e(t)=-\mu_{0}e(t)\partial t+\rho_{0}\partial w(t)$ for the Brownian motion
$w(t)$. It can be shown that
$\text{cov}\\{e(t),e(s)\\}=c\exp\\{-\mu_{0}|t-s|\\}$ where
$c=\rho_{0}^{2}/2\mu_{0}$. Here we assume $c=1$. Thus, by solving the integral
equation we have $\phi_{k}(t)=A_{k}\cos(\omega_{k}t)+B_{k}\sin(\omega_{k}t)$
and $\lambda_{k}={2\mu_{0}}/{(\omega_{k}^{2}+\mu_{0}^{2})}$ where $\omega$ is
solution of $\cot(\omega)={(\omega^{2}-\mu_{0}^{2})}/{2\mu_{0}\omega}$. The
constants $A_{k}$ and $B_{k}$ are defined as $B_{k}=\mu_{0}A_{k}/\omega_{k}$
where
$A_{k}=\sqrt{\\{}{2\omega_{k}^{2}}/{(2\mu_{0}+\mu_{0}^{2}+\omega_{k}^{2})}\\}$.
Here $\mu_{0}$ is chosen to be $1$ or $3$.

  4. (d)

Power exponential. $R(s,t)=\exp\\{(-|s-t|/a_{0})^{b_{0}}\\}$ where scale
parameter $a_{0}=1$ and shape parameter $b_{0}\in\\{1,2,5\\}$.

  5. (e)

Rational quadratic.
$R(s,t)=\left\\{1+{(s-t)^{2}}/{a_{0}^{2}}\right\\}^{-b_{0}}$ where scale
parameter $a_{0}=1$ and shape parameter $b_{0}\in\\{1,2,5\\}$.

###  4.2 Comparison and evaluation

For each of the situations, we perform 500 simulation replicates. To execute
Qu et al. [26]’s approach, we construct the scores using basis matrices as
described in Example 1 (approximation of the compound symmetric correlation
structure, denoted as ldaCS in the tables) and Example 2 (for the first-order
autoregressive correlation structure, denoted as ldaAR in the tables) in their
paper. Ordinary least squares estimate (init) is taken as the initial estimate
of $\boldsymbol{\beta}$ for both ours and Qu et al. [26]’s approach. We
indicate fda-k as our proposed method with the number of basis functions k.
The estimation procedure in the iterative algorithm converges when the square
difference between the estimated values of two consecutive steps is bounded by
a small number $10^{-10}$ or the maximum number of steps crosses 500,
whichever happens earlier. To make theoretical results and numerical examples
consistent, we use FPCA function in R which is available in fdapace packages
[9]. Bandwidths are selected using generalized cross-validation and the
Epanechnikov kernel $K(x)=0.75(1-x^{2})_{+}$ is used for estimation where
$(a)_{+}=\max(a,0)$.

The means and standard deviations (SD) of the regression coefficients based on
500 simulations are given as summary measures. We calculate the standard
deviation mentioned in the tables based on 500 estimates from 500 replications
that can be viewed as the true standard error. Moreover, we also compute
absolute bias (defined as
$\text{AB}=\sum_{b=1}^{500}|\widehat{\boldsymbol{\beta}}_{b}-\boldsymbol{\beta}|/500$)
and mean square error (defined as
$\text{MSE}=\sum_{b=1}^{500}(\widehat{\boldsymbol{\beta}}_{b}-\boldsymbol{\beta})^{2}/500$)
to compare the performance of estimation, where for $b$-th replication
$\widehat{\boldsymbol{\beta}}_{b}$ be the estimated value for
$\boldsymbol{\beta}$, MSEs are reported in the order of $10^{-2}$. The number
of selected eigen-functions plays a critical role in the proposed method. We
choose $\kappa_{0}$ based on a scree plot where the elbow of the graph is
found and the components to the left are considered as significant.

###  4.3 Simulation results

Simulation results associated with the Brownian motion are shown in Table 1.
In this situation, we observe that our approach produces better results in
terms of the dispersion measures. Tables LABEL:Chapter2-qif-Table:decay1,
LABEL:Chapter2-qif-Table:decay2 and LABEL:Chapter2-qif-Table:decay3 show
results for linear processes, our proposed method performs better in
situations with a working correlation matrix as AR but is comparable for
exchangeable structure for $l_{0}=1,2,3$. Moreover, in our proposed method, as
$l$ increases all dispersion measures such as MSE decrease. Results based on
the OU process are documented in Tables LABEL:Chapter2-qif-Table:ou1 and
LABEL:Chapter2-qif-Table:ou3. Our method outperforms the existing methods for
both situations and as $\mu_{0}$ increases, the MSE decreases. For three
different parameter choices of the power exponential and rational quadratic
covariance structure, the numerical results are presented in Tables
LABEL:Chapter2-qif-Table:pe1, LABEL:Chapter2-qif-Table:pe2,
LABEL:Chapter2-qif-Table:pe5 and LABEL:Chapter2-qif-Table:quad1,
LABEL:Chapter2-qif-Table:quad2, LABEL:Chapter2-qif-Table:quad5, respectively.
As before, we observe that our proposed method is finer than the existing ones
in all sub-cases; but interestingly, when $b_{0}$ increases, the MSE decreases
for the power exponential, whereas it increases for the rational quadratic
covariance structure as expected due to the covariance structure. Overall, we
observe that for all the above situations, as sample size increases, the
dispersion measures, for example, SD and MSE decrease. It establishes that as
sample size increases, the parameter estimates get closer and closer to the
true parameters. In each of the above situations, the SDs for the proposed
methods decrease as we increase $\kappa_{0}$ and stabilize after some value of
$\kappa_{0}$ where the fraction of variance (FVE) is approximately 100%.

Table 1: Performance of the estimation procedure where the residuals are generated from Brownian motion. Mean of the estimated coefficients, standard deviation, absolute bias, mean square error $(\times 100)$, and FVE in percentage are summarized. Method | $\beta_{1}$ | $\beta_{2}$ | FVE %-age
---|---|---|---
| Mean | SD | AB | MSE | Mean | SD | AB | MSE |
| $n=100$
init | 0.9999 | 0.0373 | 0.0297 | 0.1391 | 0.4995 | 0.0486 | 0.0384 | 0.2354 |
ldaAR | 0.9991 | 0.0331 | 0.0265 | 0.1095 | 0.5004 | 0.0445 | 0.0353 | 0.1972 |
ldaCS | 0.9987 | 0.0316 | 0.0253 | 0.1000 | 0.4997 | 0.0411 | 0.0322 | 0.1685 |
fda-1 | 0.9998 | 0.0564 | 0.0447 | 0.3180 | 0.5006 | 0.0743 | 0.0587 | 0.5516 | 86.2672
fda-2 | 1.0001 | 0.0269 | 0.0213 | 0.0723 | 0.4971 | 0.0362 | 0.0290 | 0.1314 | 96.3746
fda-3 | 0.9998 | 0.0231 | 0.0181 | 0.0532 | 0.4978 | 0.0317 | 0.0251 | 0.1010 | 99.9220
fda-4 | 1.0004 | 0.0052 | 0.0014 | 0.0028 | 0.4994 | 0.0092 | 0.0022 | 0.0085 | 99.9979
fda-5 | 0.9999 | 0.0021 | 0.0008 | 0.0004 | 0.4999 | 0.0051 | 0.0012 | 0.0026 | 100.0000
| $n=300$
init | 1.0002 | 0.0200 | 0.0162 | 0.0401 | 0.5003 | 0.0288 | 0.0226 | 0.0825 |
ldaAR | 1.0003 | 0.0184 | 0.0147 | 0.0336 | 0.5000 | 0.0259 | 0.0203 | 0.0670 |
ldaCS | 1.0007 | 0.0170 | 0.0134 | 0.0288 | 0.4995 | 0.0242 | 0.0190 | 0.0583 |
fda-1 | 1.0002 | 0.0309 | 0.0251 | 0.0955 | 0.5008 | 0.0443 | 0.0350 | 0.1962 | 86.7578
fda-2 | 1.0002 | 0.0142 | 0.0114 | 0.0202 | 0.4998 | 0.0213 | 0.0169 | 0.0451 | 96.4747
fda-3 | 1.0002 | 0.0122 | 0.0098 | 0.0150 | 0.4992 | 0.0179 | 0.0144 | 0.0321 | 99.9745
fda-4 | 1.0002 | 0.0021 | 0.0003 | 0.0004 | 0.4998 | 0.0032 | 0.0005 | 0.0010 | 99.9993
fda-5 | 1.0000 | 0.0002 | 0.0001 | 0.0000 | 0.5000 | 0.0003 | 0.0002 | 0.0000 | 100.0000
| $n=500$
init | 1.0002 | 0.0148 | 0.0117 | 0.0219 | 0.5006 | 0.0223 | 0.0177 | 0.0497 |
ldaAR | 1.0005 | 0.0138 | 0.0111 | 0.0189 | 0.5000 | 0.0206 | 0.0162 | 0.0422 |
ldaCS | 1.0000 | 0.0128 | 0.0102 | 0.0164 | 0.4992 | 0.0184 | 0.0146 | 0.0340 |
fda-1 | 1.0007 | 0.0234 | 0.0185 | 0.0545 | 0.5012 | 0.0348 | 0.0277 | 0.1213 | 86.7520
fda-2 | 0.9996 | 0.0105 | 0.0083 | 0.0110 | 0.5002 | 0.0157 | 0.0126 | 0.0247 | 96.5174
fda-3 | 0.9991 | 0.0091 | 0.0074 | 0.0084 | 0.4999 | 0.0133 | 0.0107 | 0.0177 | 99.9851
fda-4 | 1.0000 | 0.0002 | 0.0001 | 0.0000 | 0.5000 | 0.0003 | 0.0002 | 0.0000 | 99.9996
fda-5 | 1.0000 | 0.0002 | 0.0001 | 0.0000 | 0.5000 | 0.0003 | 0.0002 | 0.0000 | 100.0000

##  5 Real data analysis

In this section, we apply our proposed method to motivating examples in two
different data-sets.

###  5.1 Beijing’s $\text{PM}_{2.5}$ pollution study

In the atmosphere, suspended microscopic particles of solid and liquid matter
are commonly known as particulates or particulate matter (PM). Such
particulates often have a strong noxious impact toward human health, climate,
and visibility. One such common and fine type of atmospheric particle is
$\text{PM}_{2.5}$ with a diameter less than $2.5$ micrometers. Many developed
and developing cities across the world are experiencing chronic air pollution,
with major pollutants being $\text{PM}_{2.5}$; Beijing and a substantial part
of China are among such places. Some studies show that there are many non-
ignorable sources of variability in the distribution and transmission pattern
of $\text{PM}_{2.5}$, which are confounded with secondary chemical generation.
The atmospheric $\text{PM}_{2.5}$ data used in our analysis were collected
from the UCI machine learning repository
https://archive.ics.uci.edu/ml/datasets/Beijing+Multi-Site+Air-Quality+Data
[22]. The data-set includes daily measurements of $\text{PM}_{2.5}$ and
associated covariates at twelve different locations in China, namely
Aotizhongxin, Changping, Dingling, Dongsi, Guanyuan, Gucheng, Huairou,
Nongzhanguan, Shunyi, Tiantan, Wanliu, and Wanshouxigong during January 2017.
After excluding missing data, there were 608 hourly data points in
Beijing2017-data. We assume that the atmospheric measurements are independent
since they are located quite apart. The objective of our analysis is to
describe the trend of functional response $\text{PM}_{2.5}$ (as shown in
Figure 1 (left)) and evaluate the effect of covariates including the chemical
compounds such as sulfur dioxide $(\text{SO}_{2})$, nitrogen dioxide
$(\text{NO}_{2})$, carbon monoxide $(\text{CO})$ and ozone $(\text{O}_{3})$
over time. We smooth the covariates and responses to reduce variability and
center them. Subsequently, we consider the following model

$Y_{i}(t)=\beta_{0}+\text{SO}_{2}(t)\beta_{1}+\text{NO}_{2}(t)\beta_{2}+\text{CO}(t)\beta_{3}+\text{O}_{3}(t)\beta_{4}+e_{i}(t).$
(10)

We use Algorithm 1 to estimate the coefficients of the regression model
mentioned above. Through the simulation results, we observe that if the values
of $\kappa_{0}$ increase, the standard deviation of the coefficients
decreases. For small FVE such as $50\%$, the corresponding $\kappa_{0}=1$ and
the estimation procedure performs poorly; whereas for large FVE percentages,
the estimation procedure has adequately improved in terms of standard error.
The estimated values for $\beta_{0},\beta_{1},\beta_{2},\beta_{3},$ and
$\beta_{4}$ produce similar results across different choices of $\kappa_{0}$.
From the estimated standard error, using scree-plot of standard error, we
conclude that the suitable choice of $\kappa_{0}$ is approximately $10$. The
estimated scaled eigen-values are provided in Figure 1 (right) which clearly
shows their decay rate. The estimated coefficients with standard error are
0.0009 (1.1644), 0.0829 (0.2584), 0.9503 (0.1586), 0.0196 (0.0037) and 1.1523
(0.1198) respectively.

Fig. 1: Beijing2017-data: (left) Reading of hourly $\text{PM}_{2.5}$ measures
for twelve different locations over 608 hourly time-points during January
2017. (right) Scree plots of a fraction of variance explained (FVE).

###  5.2 DTI study for sleep apnea patients

The diffusion tensor imaging data (DTI) to understand the white matter
structural alteration using diffusion tensor imaging (DTI) are used to
illustrate the application of the proposed method and its estimation procedure
[29]. MRI is a powerful technique for investigating the structural and
functional changes in the brain during pathological and neuro-psychological
processes. Due to the advancement in diffusion tensor imaging (DTI), several
studies on white matter alterations associated with clinical variables can be
found in recent literature. For our analysis, we use Apnea-data obtained from
one such study on obstructive sleep apnea (OSA) patients [29]. The data
consists of 29 male patients between the ages 30-55 years who underwent a
study for the diagnosis of continuous positive airway pressure (CPAP) therapy.
Among those who have sleep disorders other than OSA, night-shift workers,
patients with psychiatric disorders, hypertension, diabetes, and other
neurological disorders were excluded. In this study, the psychomotor vigilance
task (PVT) was performed in which a light was randomly switched-on on a screen
for several seconds in a certain interval of time, and subjects were asked to
press a button as soon as they saw the light appear on screen; such an
experiment provides a numerical measure of sleepiness by counting the number
of “lapses” for each individual. Psychomotor vigilance task was performed
after continuous positive airway pressure treatment which counts the number of
lapses in the attention tracking test. In such a test, subjects were asked to
press a button as soon as the light appeared on the screen which was turned on
randomly for several seconds in a certain interval of time. The interesting
measure is sustained attention and provides a numerical measure of sleepiness
by counting the number of “lapses” for each individual

DTI was performed at 3T MRI scanner using a commercial 32-channel head coil,
followed by the analysis using tract-based spatial statistics to investigate
the difference in fractional anisotropy (FA) and other DTI parameters. The
image acquisition is as follows. An axial T1-weighted image of the brain
(3D-BRAVO) is collected with repetition time (TR) = 12ms, echo time (TE) =
5,2ms, flip angle = $13^{\circ}$, inversion time = 450 ms, matrix = $384\times
256$, voxel size = $1.2\times 0.57\times 0.69$mm and scan time = 2 min 54 sec.
DTI are obtained in the axial plane using a spin-echo echo planner imaging
sequence with TR = 4500ms, TE = 89.4ms, field of view = $20\times 20$cm2,
matrix size = $160\times 132$, slice thickness = 3mm, slice spacing = 1mm,
b-values = 0, 1000 s/mm2.

Our objective is to investigate the structural alteration of white matter
using DTI in patients with OSA over each voxel at various regions of the brain
(called ROIs). Thus, our response variable is one of the DTI parameters, viz.,
fractional anisotropy (FA) and we are interested in studying the effect of the
changes of FA over continuous domains such as voxels with the interaction of
the lapses and the voxel locations in each ROIs. We consider the following
model for each ROI.

$\text{FA}_{i}(s)=\beta_{0}+\beta_{1}\text{lapses}_{i}\times s+e_{i}(s)$ (11)

where $s\in\mathcal{S}$, a set of voxels in the considered ROIs. Using the
Algorithm 1, we estimate the coefficients $\beta_{1}$ and $\beta_{2}$ as
mentioned in the model 11 and the results are presented in Table 2. We find
that the coefficient estimates are close enough to their initial estimates and
the estimated standard error is smaller for the coefficients based on the
proposed method. Here $\kappa_{0}$ (i.e., the number of eigen-functions) are
determined by FVE for simplicity.

Table 2:  Estimated values and associated standard errors for the regression coefficients are provided upto four decimal places based on the existing and proposed methods. First line corresponding to each ROI shows results based on initial estimates and the second line corresponds to that of proposed estimates. |  | $\beta_{0}$ | $\beta_{1}$
---|---|---|---
region | # functional points | Estimate | Std. Error $(\times 100)$ | Estimate $(\times 100)$ | Std. Error $(\times 100)$
ROI.6 | 659 | 0.4512 | 0.1343 | -0.0606 | 0.0130
|  | 0.4512 | 0.0983 | -0.0605 | 0.0028
ROI.7 | 1362 | 0.5048 | 0.0628 | 0.0309 | 0.0061
|  | 0.5050 | 0.0681 | 0.0342 | 0.0007
ROI.8 | 1370 | 0.5256 | 0.0586 | -0.0667 | 0.0057
|  | 0.5271 | 0.0346 | -0.0733 | 0.0006
ROI.9 | 690 | 0.4951 | 0.0910 | 0.2904 | 0.0088
|  | 0.5443 | 0.0874 | 0.1660 | 0.0014
ROI.10 | 699 | 0.4951 | 0.0892 | 0.3314 | 0.0086
|  | 0.5262 | 0.1398 | 0.4231 | 0.0014
ROI.11 | 968 | 0.4372 | 0.0979 | 0.1323 | 0.0095
|  | 0.4380 | 0.0637 | 0.1311 | 0.0009
ROI.12 | 968 | 0.4529 | 0.0948 | 0.0965 | 0.0092
|  | 0.4664 | 0.0750 | 0.0504 | 0.0013
ROI.13 | 992 | 0.5448 | 0.1060 | 0.3453 | 0.0103
|  | 0.5449 | 0.0856 | 0.3559 | 0.0011
ROI.14 | 992 | 0.5435 | 0.1068 | 0.3432 | 0.0104
|  | 0.5436 | 0.0754 | 0.3437 | 0.0003
ROI.37 | 1236 | 0.3695 | 0.0779 | -0.1126 | 0.0076
|  | 0.3713 | 0.0669 | -0.1175 | 0.0017
ROI.38 | 1155 | 0.3564 | 0.0819 | -0.1356 | 0.0079
|  | 0.3578 | 0.0420 | -0.1356 | 0.0009
ROI.39 | 1124 | 0.4618 | 0.0760 | 0.1972 | 0.0074
|  | 0.4621 | 0.0615 | 0.1996 | 0.0007
ROI.40 | 1125 | 0.4786 | 0.0658 | 0.0953 | 0.0064
|  | 0.4780 | 0.0369 | 0.1016 | 0.0005
ROI.45 | 380 | 0.4189 | 0.1071 | 0.1647 | 0.0104
|  | 0.4190 | 0.0175 | 0.1648 | 0.0001
ROI.46 | 376 | 0.4074 | 0.1033 | 0.1988 | 0.0100
|  | 0.4074 | 0.0159 | 0.1994 | 0.0002
ROI.47 | 596 | 0.4596 | 0.0932 | 0.1304 | 0.0090
|  | 0.4594 | 0.0191 | 0.1349 | 0.0001
ROI.48 | 600 | 0.4045 | 0.0868 | 0.1100 | 0.0084
|  | 0.4036 | 0.0644 | 0.1067 | 0.0006

##  6 Discussion

In this article, we propose an estimation procedure for the constant linear
effects model, which is commonly used in statistics [33] especially in spatial
modeling. One of the key factors of this estimation procedure is the fact that
it is based on the quadratic inference methodology that has served a huge role
in the analysis of correlated data since it was discovered by Qu et al. [26].
In contrast with the existing method, our approach allows the number of
repeated measurements to grow with sample size; therefore, the trajectories of
individuals can be observed on a dense grid of a continuum domain. Instead of
assuming a working correlation structure, we propose a data-driven way by
estimating the eigen-functions that are obtained by functional principal
component analysis. Here, we achieve $\sqrt{n}-$consistency of the parametric
estimates in the regression model, even though the eigen-functions are
estimated non-parametrically.

Additionally, our method is easy to implement in a wide range of applications.
The applicability of the proposed method is illustrated by extensive
simulation studies. Moreover, two real-data applications in different
scientific domains are provided which confirm the efficacy of the proposed
method.

## Acknowledgement

We would like to thank Dr. Xiaohong Joe Zhou of University of Illinois at
Chicago for providing the Apnea-data used in Section 5.2.

##  Appendix A Some preliminary definitions and concepts of operators

Consider the standard $\mathcal{L}^{2}[0,1]$ space that defines the set of
square-integrable functions defined on the closed set $[0,1]$ that takes
values on the real line. The space $\mathcal{L}^{2}[0,1]$ is equipped with an
inner product and is defined as $\left<f,g\right>=\int_{0}^{1}f(t)g(t)dt$ for
$f$ and $g$ in that space, and forms a Hilbert space. Moreover, we denote the
norm $\|\cdot\|_{2}$ in $\mathcal{L}^{2}$ which is defined as
$\|f\|_{2}=\left\\{\int f^{2}(u)du\right\\}^{1/2}$. Define $\mathcal{F}$ be an
operator that assigns an element $f$ in $\mathcal{L}^{2}[0,1]$ to a new
element $\mathcal{F}f$ in $\mathcal{L}^{2}[0,1]$ moreover, $\mathcal{F}$ is
linear and bounded. A linear mapping $\mathcal{F}f(\cdot)=\int
R(\cdot,u)f(u)du$ for any function $f\in\mathcal{L}^{2}[0,1]$ and for some
integrable function $R(\cdot,\cdot)$ on $[0,1]\times[0,1]$. This function is
preferably known as integral operator and the bivariate function $R$ is known
as a kernel in statistics and functional analysis literature. Furthermore,
under the assumption that $\int\int R^{2}(u,v)dudv<\infty$. It is easy to see
that $\mathcal{F}f(\cdot)$ is uniformly continuous and compact for a non-
negative definite symmetric kernel $R$. For some $\lambda$, in Fredholm
integral equation , $\mathcal{F}\phi=\lambda\phi$ has non-zero solution $\phi$
then we call $\lambda$ as eigen-value of $\mathcal{F}$ and the solution of the
eigen-equation is called eigen-functions, altogether, the pair of eigen-values
and eigen-function, viz., $(\lambda,\phi)$ are called eigen-elements. Due to
non-negative definiteness of $\mathcal{F}$, the eigen-values are ordered as
$\lambda_{1}\geq\lambda_{2}\geq\cdots\geq 0$ .

Suppose for self-adjoint compact operator on Hilbert space $\mathcal{H}$
consider two operators $\mathcal{F}$ and $\mathcal{G}$, define perturbation
operator $\Delta=\mathcal{G}-\mathcal{F}$ such that
$\mathcal{G}=\mathcal{F}+\Delta$ where $\mathcal{G}$ is an approximation to
$\mathcal{F}$ where $\Delta$ amount of error is occurred. Let $\mathcal{F}$
and $\mathcal{G}$ have kernels $F$ and $G$ respectively with eigen-elements
$(\theta_{r},\psi_{r})$ and $(\lambda_{r},\phi_{r})$. For simplicity, we
assume that the eigen-values are distinct. Then the following Lemma provides
perturbation of the eigen-functions.

######  Lemma 1 (Theorem 5.1.8 in Hsing and Eubank [16]).

Let $(\lambda,\phi)$ be the eigen-components of $\mathcal{F}$ and
$(\theta,\psi)$ be that of $\mathcal{G}$ with multiplicity of all eigen-values
are restricted to be 1. Define $\eta_{k}=\min_{r\neq
k}|\lambda_{r}-\lambda_{k}|$. Assume $\left<\phi_{r},\psi_{r}\right>\geq 0$
and $\eta_{k}>0$. Then

$\psi_{k}-\phi_{k}=\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\theta_{k}-\lambda_{r})^{-1}\mathcal{P}_{r}\Delta\psi_{k}+\mathcal{P}_{k}(\psi_{k}-\phi_{k}).$
(12)

The above equation follows

$\psi_{k}-\phi_{k}=\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\theta_{k}-\lambda_{r})^{-1}\mathcal{P}_{r}\Delta\psi_{k}+O(\|\Delta\|^{2}).$
(13)

######  Remark 3.

Equation (13) plays an important role in finding the bound of the proposed
estimator introduced in Section 2.2. Note that $\sup_{r\geq
1}|\theta_{r}-\lambda_{r}|\leq\|\Delta\|\leq\inf_{r\neq
k}|\lambda_{k}-\lambda_{r}|$ (see Theorem 4.2.8 in Hsing and Eubank [16] for
proof). Thus, it is easy to see,
$|\theta_{r}-\lambda_{r}|\leq|\lambda_{k}-\lambda_{r}|$ which implies from
Equation (12)

$\displaystyle\psi_{k}-\phi_{k}=\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}\sum_{s=0}^{\infty}\left\\{{(\lambda_{k}-\theta_{r})}/{(\lambda_{k}-\lambda_{r})}\right\\}^{s}\mathcal{P}_{r}\Delta\\{\phi_{k}+(\psi_{k}-\phi_{k})\\}+\mathcal{P}_{k}(\psi_{k}-\phi_{k})$
$\displaystyle=\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}\mathcal{P}_{r}\Delta\phi_{k}+\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}\mathcal{P}_{r}\Delta(\psi_{k}-\phi_{k})$
$\displaystyle\qquad+\underset{r\neq
k}{\sum_{r=1}^{\infty}}\sum_{s=1}^{\infty}\left\\{{(\lambda_{k}-\lambda_{s})^{s}}/{(\lambda_{k}-\lambda_{r})^{s+1}}\right\\}\mathcal{P}_{r}\Delta\psi_{k}+\mathcal{P}_{k}(\psi_{k}-\phi_{k}).$
(14)

Moreover, using Bessel’s inequality, we can bound last three terms in the
above equation by $\|\Delta\|^{2}$.

##  Appendix B Some useful lemmas

In this section, we present some useful lemmas. For convenience, let us recall
the notation. Assume that $m_{i}$s are all of the same order, viz, $m\equiv
m(n)$. Define, $d_{n1}(h)=h^{2}+h\overline{m}/m$ and
$d_{n2}(h)=h^{4}+h^{3}\overline{m}/m+h^{2}\overline{\overline{m}}/m^{2}$ where
$\overline{m}=\lim_{n\rightarrow\infty}n^{-1}\sum_{i=1}^{n}m/m_{i}$ and
$\overline{\overline{m}}=\lim_{n\rightarrow\infty}n^{-1}\sum_{i=1}^{n}(m/m_{i})^{2}$.
Denote $\delta_{n1}(h)=\left\\{d_{n1}(h)\log n/(nh^{2})\right\\}^{1/2}$,
$\delta_{n2}(h)=\left\\{d_{n2}(h)\log n/(nh^{4})\right\\}^{1/2}$ and
$\overline{\delta}_{n}(h)=h^{2}+\delta_{n1}(h)+\delta_{n2}^{2}(h)$. Further,
$\nu_{a,b}=\int t^{a}K^{b}(t)dt$. Define,
$\mathbf{W}=(\boldsymbol{\phi}(t_{1})^{\text{T}},\cdots,\boldsymbol{\phi}(t_{m})^{\text{T}})^{\text{T}}$
be matrix of order $m\times\kappa_{0}$ obtained after stacking all
$\boldsymbol{\phi}_{k}$s and random components
$\boldsymbol{\xi}_{i}=(\xi_{i1},\cdots,\xi_{i\kappa_{0}})^{\text{T}}$.
Further, $\boldsymbol{\xi}$ has mean zero and variance $\boldsymbol{\Lambda}$
which is a diagonal matrix with components
$\lambda_{1},\cdots,\lambda_{k_{0}}$. The sign ${}^{\prime}\lesssim^{\prime}$
indicates that the left-hand side of the inequality is bounded by the right-
hand upto a multiplicative positive constant and the constant does not
involves $n$ or function of $n$.

######  Lemma 2.

Consider $Z_{1},\cdots,Z_{n}$ be independent and identically distributed
random variables with mean zero and finite variance. Suppose that there exists
an $M$ such that $P(|Z_{i}|\leq M)=1$ for all $i=1,\cdots,n$. Let
$T_{n}=n^{-1}\sum_{i=1}^{n}Z_{i}$. then, $n^{-1}\sum_{i=1}^{n}Z_{i}=O((\log
n/n)^{1/2}$ almost surely. If $\sqrt{Var(T_{n})}=O((\log n/n)^{1/2})$ then
$T_{n}=O(\log n/n)$ almost surely.

###### Proof.

Bernstein’s inequality states that if $Z_{1},\cdots,Z_{n}$ be centered
independent bounded random variables with probability 1. Let
$T_{n}=n^{-1}\sum_{i=1}^{n}Z_{i}$, then let
$\text{Var}\\{T_{n}\\}=\sigma_{n}^{2}$. Then for any positive real number $u$,
we have $P(|T_{n}|\geq u)\leq\exp\\{-{nu^{2}}/(2\sigma^{2}_{n}+2Mu/3)\\}$
where $M$ is such that $P(|Z_{i}|\leq M)=1$. Moreover, if $T_{n}$ converges to
its limit in probability fast enough, then it converge almost surely in the
limit, i.e., if for any $u>0$, $\sum_{n=1}^{\infty}P(|T_{n}|\geq u)<\infty$
them $T_{n}$ converges to zero almost surely. Now, choose
$u=\sqrt{\\{}{4\sigma^{2}_{n}\log n}/{n}\\}+{4M\log n}/{3n}$. Thus,
$\sum_{n=1}^{\infty}P(|T_{n}|\geq u)<\sum_{n=1}^{\infty}1/n^{2}$ which is
finite. Therefore, $T_{n}=O(u)$ almost surely. Now let
$\sigma_{n}\leq\sqrt{\\{}4M^{2}\log n/9n\\}$, we have, $T_{n}=O(\log n/n)$ and
if $\sigma_{n}=O(1)$ then $T_{n}=O((\log n/n)^{1/2})$ almost surely. ∎

######  Lemma 3.

Suppose $T_{ij}$ are i.i.d. with density $f_{T}$. Then for fixed
$i=1,\cdots,n$, any $k$ and $l\geq 1$, under assumptions (C2) and (C6)c,
$m_{i}^{-1}\sum_{j=1}^{m_{i}}\phi_{k}(T_{ij})\phi_{l}(T_{ij})=\textbf{1}(k=l)+O(\left({\log
m_{i}}/{m_{i}}\right)^{1/2})$ almost surely where 1 is the indicator function.

###### Proof.

Observe,
$\mathbb{E}\left\\{m_{i}^{-1}\sum_{j=1}^{m_{i}}\phi_{k}(T_{ij})\phi_{l}(T_{ij})\right\\}=\int\phi_{k}(t)\phi_{l}(t)dt=\textbf{1}(k=l)$
and,

$\displaystyle\text{Var}\left\\{m_{i}^{-1}\sum_{j=1}^{m_{i}}\phi_{k}(T_{ij})\phi_{l}(T_{ij})\right\\}=\mathbb{E}\left\\{m_{i}^{-1}\sum_{j=1}^{m_{i}}\phi_{k}(T_{ij})\phi_{l}(T_{ij})\right\\}^{2}-\textbf{1}(k=1)$
$\displaystyle=m_{i}^{-2}\sum_{j=1}^{m_{i}}\mathbb{E}\\{\phi_{k}^{2}(T_{ij})\phi_{l}^{2}(T_{ij})\\}+m_{i}^{-2}\underset{j_{1}\neq
j_{2}}{\sum_{j_{1}=1}^{m_{i}}\sum_{j_{2}=1}^{m_{i}}}\mathbb{E}\\{\phi_{k}(T_{ij_{1}})\phi_{k}(T_{ij_{2}})\phi_{l}(T_{ij_{1}})\phi_{l}(T_{ij_{2}})\\}$
$\displaystyle\qquad-\textbf{1}(k=l)$
$\displaystyle=m_{i}^{-2}\sum_{j=1}^{m_{i}}\mathbb{E}\\{\phi_{k}^{2}(T_{ij})\phi_{l}^{2}(T_{ij})\\}+m_{i}^{-2}\underset{j_{1}\neq
j_{2}}{\sum_{j_{1}=1}^{m_{i}}\sum_{j_{2}=1}^{m_{i}}}\mathbb{E}\\{\phi_{k}(T_{ij_{1}})\phi_{l}(T_{ij_{1}})\\}\mathbb{E}\\{\phi_{k}(T_{ij_{2}})\phi_{l}(T_{ij_{2}})\\}$
$\displaystyle\qquad-\textbf{1}(k=l)$
$\displaystyle=\begin{cases}m_{i}^{-1}\int\phi_{k}^{4}(t)dt+(m_{i}-1)/m_{i}(\int\phi_{k}^{2}(t)dt)^{2}-1&\text{if
}k=l\\\
m_{i}^{-1}\int\phi_{k}^{2}(t)\phi_{l}^{2}(t)dt+(m_{i}-1)/m_{i}(\int\int\phi_{k}(t)\phi_{l}(t)\phi_{k}(t^{\prime})\phi_{l}(t^{\prime})dtdt^{\prime})&\text{if
}k\neq l\end{cases}$ $\displaystyle=O(1/m_{i}).$ (15)

Therefore, applying the Lemma 2, the result is immediate. ∎

######  Lemma 4.

Suppose $T_{ij}$ are i.i.d with density $f_{T}$. Then for fixed
$i=1,\cdots,n$, for any $k\geq 1$, under assumptions (C2), (C6)c,
$m_{i}^{-1}\sum_{j=1}^{m_{i}}\dot{\mu}_{i}(T_{ij})\phi_{k}(T_{ij})=\int\dot{\mu}_{i}(t)\phi_{k}(t)dt+O\left((\log
m_{i}/m_{i})^{1/2}\right)$ almost surely.

###### Proof.

Observe,
$\mathbb{E}\left\\{m_{i}^{-1}\sum_{j=1}^{m_{i}}\dot{\mu}_{i}(T_{ij})\phi_{k}(T_{ij})\right\\}=\int\dot{\mu}_{i}(t)\phi_{k}(t)dt$
and,

$\displaystyle\text{Var}\left\\{m_{i}^{-1}\sum_{j=1}^{m_{i}}\dot{\mu}_{i}(T_{ij})\phi_{k}(T_{ij})\right\\}\leq\mathbb{E}\left\\{m_{i}^{-1}\sum_{j=1}^{m_{i}}\dot{\mu}_{i}(T_{ij})\phi_{k}(T_{ij})\right\\}^{2}$
$\displaystyle=m_{i}^{-2}\sum_{j=1}^{m_{i}}\mathbb{E}\left\\{\dot{\mu}_{i}^{2}(T_{ij})\phi_{k}(T_{ij})^{2}\right\\}+m_{i}^{-2}\underset{j_{1}\neq
j_{2}}{\sum_{j_{1}=1}^{m_{i}}\sum_{j_{2}=1}^{m_{i}}}\mathbb{E}\\{\dot{\mu}_{i}(T_{ij_{1}})\dot{\mu}_{i}(T_{ij_{2}})\phi_{k}(T_{ij_{1}})\phi_{k}(T_{ij_{2}})\\}$
$\displaystyle=O(1/m_{i}),\qquad\text{since
}\int\dot{\mu}^{2}(t)\phi_{k}^{2}(t)dt<\infty.$ (16)

Therefore, applying the Lemma 2, the result is immediate. ∎

######  Lemma 5.

Define, $\mathcal{M}_{ir}=\int\dot{\mu}_{i}(t)\phi_{r}(t)dt$ and
$V_{r}=\mathbb{E}\\{\int\dot{\mu}(t)\phi_{r}(t)dt\\}^{2}$ for $r\geq 1$. Then
under Conditions (C6)a and (C6)b, for some $\alpha>0$ such that
$V_{r}\lambda_{r}^{-2}r^{1+\alpha}\rightarrow 0$ as $r\rightarrow\infty$ (due
to Condition (C6)b, $\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}n^{-1}\sum_{i=1}^{n}\mathcal{M}_{ir}\xi_{ik}=O\left((\log
n/n)^{1/2}\lambda_{k}^{1/2}k^{(1-\alpha)/2}\right)$ almost surely.

###### Proof.

It is easy to see that, $\mathbb{E}\Big{\\{}\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}n^{-1}\sum_{i=1}^{n}\mathcal{M}_{ir}\xi_{ik}\Big{\\}}=0$.
Using the spacing condition among the eigen-values in (C6)a, for each $1\leq
k<r<\infty$ and for nonzero finite generic constant $C_{0}$,

${\max\left\\{\lambda_{k},\lambda_{r}\right\\}}/{|\lambda_{k}-\lambda_{r}|}\leq
C_{0}{\max\left\\{k,r\right\\}}/{|k-r|}.$ (17)

Similar kind of conditions can be invoked such as the convexity assumption,
i.e. $\lambda_{r}-\lambda_{r+1}\leq\lambda_{r-1}-\lambda_{r}$ for all $r\geq
2$. Thus, using Inequality (17), for some $\alpha>0$, with condition
$V_{r}\lambda_{r}^{-2}r^{1+\alpha}\rightarrow 0$ as $r\rightarrow\infty$, we
can write

$\displaystyle\underset{r\neq
k}{\sum_{r=1}^{\infty}}V_{r}(\lambda_{k}-\lambda_{r})^{-2}\lesssim\underset{r\neq
k}{\sum_{r=1}^{\infty}}V_{r}\left\\{{\max(k,r)}/{|k-r|\max(\lambda_{k},\lambda_{r})}\right\\}^{2}$
$\displaystyle=\sum_{r\leq
k/2}V_{r}\lambda_{r}^{-2}{k^{2}}/{(k-r)^{2}}+\sum_{r>2k}V_{r}\lambda_{k}^{-2}{r^{2}}/{(k-r)^{2}}$
$\displaystyle\qquad+\sum_{k/2<r<k}V_{r}\lambda_{r}^{-2}{k^{2}}/{(k-r)^{2}}+\sum_{k<r<2k}V_{r}\lambda_{k}^{-2}{r^{2}}/{(k-r)^{2}}$
$\displaystyle\lesssim\sum_{r\leq
k/2,r>2k}V_{r}\lambda_{r}^{-2}+k^{2}\sum_{k/2<r<2k}V_{r}\lambda_{r}^{-2}(k-r)^{-2}$
$\displaystyle\lesssim 1+k^{1-\alpha}\sum_{k/2<r<2k}(k-r)^{-2}\lesssim
k^{1-\alpha}.$ (18)

This follows the line of proofs in Hall and Hosseini-Nasab [13] in different
contexts. Thus, using the Inequality (B), it follows that

$\displaystyle\text{Var}\Big{\\{}\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}n^{-1}\sum_{i=1}^{n}\mathcal{M}_{ir}\xi_{ik}\Big{\\}}$
$\displaystyle=n^{-1}\lambda_{k}\underset{r\neq
k}{\sum_{r=1}^{\infty}}V_{r}(\lambda_{k}-\lambda_{r})^{-2}=O\Big{(}n^{-1}\lambda_{k}k^{(1-\alpha)}\Big{)}.$
(19)

Therefore, applying Lemma 2, the proof is immediate. ∎

######  Lemma 6.

For $\mathcal{M}_{ir}=\int\dot{\mu}(t)\phi_{r}(t)dt$ and $\eta_{k}=\min_{r\neq
k}|\lambda_{k}-\lambda_{r}|>0$, under Conditions (C6)a and (C6)b,
$\underset{r_{1}\neq k}{\sum_{r_{1}\neq 1}^{\infty}}\underset{r_{2}\neq
k}{\sum_{r_{2}\neq
1}^{\kappa_{0}}}(\lambda_{k}-\lambda_{r_{1}})^{-1}(\lambda_{k}-\lambda_{r_{2}})^{-1}n^{-1}\sum_{i=1}^{n}\mathcal{M}_{ir_{1}}\xi_{ir_{2}}=O\left((\log
n/n)^{1/2}\kappa_{0}^{(3-\alpha)/2}\lambda_{\kappa_{0}}^{-1}\left\\{\sum_{r=1}^{\kappa_{0}}\lambda_{r}\right\\}^{1/2}\right)$
almost surely.

###### Proof.

It is not difficult to see that,

$\mathbb{E}\Big{\\{}\underset{r_{1}\neq
k}{\sum_{r_{1}=1}^{\infty}}\underset{r_{2}\neq
k}{\sum_{r_{2}=1}^{\kappa_{0}}}(\lambda_{k}-\lambda_{r_{1}})^{-1}(\lambda_{k}-\lambda_{r_{2}})^{-1}n^{-1}\sum_{i=1}^{n}\mathcal{M}_{ir_{1}}\xi_{ir_{2}}\Big{\\}}=0.$

Moreover, using the spacing condition mentioned in (C6)b, one can derive the
upper bound of $\eta_{k}^{-1}$ by

$\displaystyle\eta_{k}^{-1}$ $\displaystyle=\left\\{\min_{r\neq
k}|\lambda_{k}-\lambda_{r}|\right\\}^{-1}=\max_{r\neq
k}|\lambda_{k}-\lambda_{r}|^{-1}$ $\displaystyle\lesssim\max_{r\neq
k}\left\\{{\max(k,r)}/{|k-r|\max(\lambda_{k},\lambda_{r})}\right\\}\leq\lambda_{k}^{-1}k.$
(20)

Due to the monotonic decreasing property of eigen-values, for fixed
$k=1,\cdots,\kappa_{0}$, we have

$\displaystyle\underset{r\neq
k}{\sum_{r=1}^{\kappa_{0}}}\lambda_{r}(\lambda_{k}-\lambda_{r})^{-2}$
$\displaystyle\lesssim\eta_{k}^{-2}\sum_{r=1}^{\kappa_{0}}\lambda_{r}\lesssim\lambda_{k}^{-2}k^{2}\sum_{r=1}^{\kappa_{0}}\lambda_{r}\lesssim\lambda_{\kappa_{0}}^{-2}\kappa_{0}^{2}\sum_{r=1}^{\kappa_{0}}\lambda_{r}.$
(21)

Therefore, the following holds under similar conditions to obtain the
Inequality (B),

$\displaystyle\mathbb{E}\left\\{\underset{r_{1}\neq
k}{\sum_{r_{1}=1}^{\infty}}\underset{r_{2}\neq
k}{\sum_{r_{2}=1}^{\kappa_{0}}}(\lambda_{k}-\lambda_{r_{1}})^{-2}(\lambda_{k}-\lambda_{r_{2}})^{-2}\left(n^{-1}\sum_{i=1}^{n}\mathcal{M}_{ir_{1}}\xi_{ir_{2}}\right)^{2}\right\\}$
$\displaystyle=n^{-1}\underset{r_{1}\neq
k}{\sum_{r_{1}=1}^{\infty}}\underset{r_{1}\neq
k}{\sum_{r_{1}=1}^{\kappa_{0}}}(\lambda_{k}-\lambda_{r_{1}})^{-2}(\lambda_{k}-\lambda_{r_{2}})^{-2}V_{r_{1}}\lambda_{r_{2}}$
$\displaystyle\lesssim n^{-1}k^{1-\alpha}\underset{r\neq
k}{\sum_{r=1}^{\kappa_{0}}}\lambda_{r}(\lambda_{k}-\lambda_{r})^{-2}\lesssim
n^{-1}\kappa_{0}^{3-\alpha}\lambda_{\kappa_{0}}^{-2}\sum_{r=1}^{\kappa_{0}}\lambda_{r}.$
(22)

Therefore, applying the Lemma 2, the result is immediate. ∎

######  Remark 4.

Define,
$d_{ik}(T_{ij_{1}},T_{ij_{2}}):=\widehat{\phi}_{k}(T_{ij_{1}})\widehat{\phi}_{k}(T_{ij_{2}})-\phi_{k}(T_{ij_{1}})\phi_{k}(T_{ij_{2}})$.
Now replace $\mathcal{G}$, $\theta$ and $\psi$ by $\widehat{\mathcal{F}}$,
$\widehat{\lambda}$ and $\widehat{\phi}$ respectively since
$\widehat{\mathcal{F}}$ be the approximation of $\mathcal{F}$ and $\Delta$ is
the corresponding perturbation operator in Equation (12). Therefore, Lemma 1
immediately implies the following expansion, which is the key fact to
represent the objective function in QIF.

$\widehat{\phi}_{k}-\phi_{k}=\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}\left<\phi_{r},\Delta\phi_{k}\right>\phi_{r}+O(\|\Delta\|^{2})\qquad\text{almost
surely},$ (23)

where $\Delta$ be the integral operator with kernel $\widehat{R}-R$.

##  Appendix C Proof of Theorem 1

For the $k-$th element of $\overline{\mathbf{g}}_{n}(\boldsymbol{\beta}_{0})$
$(1\leq k\leq\kappa_{0})$,

$\displaystyle\overline{\mathbf{g}}_{n,k}(\boldsymbol{\beta}_{0})$
$\displaystyle=n^{-1}\sum_{i=1}^{n}m_{i}^{-2}\dot{\boldsymbol{\mu}}_{i}^{\text{T}}\widehat{\boldsymbol{\Phi}}_{k}(\mathbf{y}_{i}-\boldsymbol{\mu}_{i})=n^{-1}\sum_{i=1}^{n}m_{i}^{-2}\dot{\boldsymbol{\mu}}_{i}^{\text{T}}\widehat{\boldsymbol{\Phi}}_{k}\mathbf{W}\boldsymbol{\xi}_{i}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}m_{i}^{-2}\dot{\boldsymbol{\mu}}_{i}^{\text{T}}\boldsymbol{\Phi}_{k}\mathbf{W}\boldsymbol{\xi}_{i}+n^{-1}\sum_{i=1}^{n}m_{i}^{-2}\dot{\boldsymbol{\mu}}_{i}^{\text{T}}(\widehat{\boldsymbol{\Phi}}_{k}-\boldsymbol{\Phi}_{k})\mathbf{W}\boldsymbol{\xi}_{i}:=J_{k}^{n1}+J_{k}^{n2}.$
(24)

Now, using Lemmas 3 and 4, the first part of the expression of
$\overline{\mathbf{g}}_{n,k}(\boldsymbol{\beta}_{0})$ becomes

$\displaystyle J_{k}^{n1}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}m_{i}^{-2}\dot{\boldsymbol{\mu}}_{i}^{\text{T}}\boldsymbol{\Phi}_{k}\mathbf{W}\boldsymbol{\xi}_{i}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}m_{i}^{-2}\sum_{j_{1}=1}^{m_{i}}\sum_{j_{2}=1}^{m_{i}}\sum_{l=1}^{\kappa_{0}}\dot{\mu}_{i}(T_{ij_{1}})\phi_{k}(T_{ij_{1}})\phi_{k}(T_{ij_{2}})\phi_{l}(T_{ij_{2}})\xi_{il}$
$\displaystyle\lesssim n^{-1}\sum_{i=1}^{n}\left\\{\mathcal{M}_{ik}+O((\log
m/m)^{1/2})\right\\}\left\\{1+O((\log m/m)^{1/2})\right\\}\xi_{ik},$
$\displaystyle\qquad\text{where
}\mathcal{M}_{ik}=\int\dot{\mu}_{i}(t)\phi_{k}(t)dt$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\mathcal{M}_{ik}\xi_{ik}\left\\{1+O\left((\log
m/m)^{1/2}\right)\right\\}=O\left((\log n/n)^{1/2}\left\\{1+(\log
m/m)^{1/2}\right\\}\right)$ $\displaystyle\qquad\text{almost surely}.$ (25)

On the other hand, the last part of
$\overline{\mathbf{g}}_{n,k}(\boldsymbol{\beta})$ can be expressed as

$\displaystyle J_{k}^{n2}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}{m_{i}^{-2}}\dot{\boldsymbol{\mu}}_{i}^{\text{T}}(\widehat{\boldsymbol{\Phi}}_{k}-\boldsymbol{\Phi}_{k})\mathbf{W}\boldsymbol{\xi}_{i}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}{m_{i}^{-2}}\sum_{j_{1}=1}^{m_{i}}\sum_{j_{2}=1}^{m_{i}}\sum_{l=1}^{\kappa_{0}}\dot{\mu}_{i}(T_{ij_{1}})d_{ik}(T_{ij_{1}},T_{ij_{2}})\phi_{l}(T_{ij_{2}})\xi_{il},$
(26)

where
$d_{ik}(T_{ij_{1}},T_{ij_{2}}):=\widehat{\phi}_{k}(T_{ij_{1}})\widehat{\phi}_{k}(T_{ij_{2}})-\phi_{k}(T_{ij_{1}})\phi_{k}(T_{ij_{2}})$
as defined in Remark 4. Therefore, using the discussion in Remark 4, we can
obtain the following expression almost surely.

$\displaystyle
d_{ik}(T_{ij_{1}},T_{ij_{2}}):=\widehat{\phi}_{k}(T_{ij_{1}})\widehat{\phi}_{k}(T_{ij_{2}})-\phi_{k}(T_{ij_{1}})\phi_{k}(T_{ij_{2}})$
$\displaystyle=\left\\{\phi_{k}(T_{ij_{1}})+\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}\left<\phi_{r},\Delta\phi_{k}\right>\phi_{r}(T_{ij_{1}})+O(\|\Delta\|^{2})\right\\}$
$\displaystyle\qquad\times\left\\{\phi_{k}(T_{ij_{2}})+\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}\left<\phi_{r},\Delta\phi_{k}\right>\phi_{r}(T_{ij_{2}})+O(\|\Delta\|^{2})\right\\}-\phi_{k}(T_{ij_{1}})\phi_{k}(T_{ij_{2}})$
$\displaystyle=\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}\left<\phi_{r},\Delta\phi_{k}\right>\Big{\\{}\phi_{r}(T_{ij_{1}})\phi_{k}(T_{ij_{2}})+\phi_{k}(T_{ij_{1}})\phi_{r}(T_{ij_{2}})\Big{\\}}$
$\displaystyle\qquad+\sum_{r_{1}\neq k}\sum_{r_{2}\neq
k}(\lambda_{k}-\lambda_{r_{1}})^{-1}(\lambda_{k}-\lambda_{r_{2}})^{-1}\left<\phi_{r_{1}},\Delta\phi_{k}\right>\left<\phi_{r_{2}},\Delta\phi_{k}\right>\phi_{r_{1}}(T_{ij_{1}})\phi_{r_{2}}(T_{ij_{2}})$
$\displaystyle\qquad+O(\|\Delta\|^{2})$
$\displaystyle:=I_{ik}^{n1}(T_{ij_{1}},T_{ij_{2}})+I_{ik}^{n2}(T_{ij_{1}},T_{ij_{2}})+O(\|\Delta\|^{2}).$
(27)

Thus, almost surely, we can write

$\displaystyle
n^{-1}\sum_{i=1}^{n}m_{i}^{-2}\dot{\boldsymbol{\mu}}_{i}^{\text{T}}(\widehat{\boldsymbol{\Phi}}_{k}-\boldsymbol{\Phi}_{k})\mathbf{W}\boldsymbol{\xi}_{i}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}m_{i}^{-2}\sum_{j_{1}=1}^{m_{i}}\sum_{j_{2}=1}^{m_{i}}\sum_{l=1}^{\kappa_{0}}\dot{\mu}_{i}(T_{ij_{1}})d_{ik}(T_{ij_{1}},T_{ij_{2}})\phi_{l}(T_{ij_{2}})\xi_{il}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}m_{i}^{-2}\sum_{j_{1}=1}^{m_{i}}\sum_{j_{2}=1}^{m_{i}}\sum_{l=1}^{\kappa_{0}}\dot{\mu}_{i}(T_{ij_{1}})d_{ik}(T_{ij_{1}},T_{ij_{2}})\phi_{l}(T_{ij_{2}})\xi_{il}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}m_{i}^{-2}\sum_{j_{1}=1}^{m_{i}}\sum_{j_{2}=1}^{m_{i}}\sum_{l=1}^{\kappa_{0}}\dot{\mu}_{i}(T_{ij_{1}})\Big{\\{}I_{ik}^{n1}(T_{ij_{1}},T_{ij_{2}})+I_{ik}^{n2}(T_{ij_{1}},T_{ij_{2}})\Big{\\}}\phi_{l}(T_{ij_{2}})\xi_{il}+O(\|\Delta\|^{2})$
$\displaystyle:=J_{k1}^{n2}+J_{k2}^{n2}+O(\|\Delta\|^{2}).$ (28)

Under assumptions (C1)-(C5), by using Theorem 3.3 in Li and Hsing [19],
$\|\Delta\|^{2}=O(h^{4}+\delta_{n2}^{2}(h))$ almost surely. Now observe that

$\displaystyle
J_{k1}^{n2}=n^{-1}\sum_{i=1}^{n}m_{i}^{-2}\sum_{j_{1}=1}^{m_{i}}\sum_{j_{2}=1}^{m_{i}}\sum_{l=1}^{\kappa_{0}}\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}\dot{\mu}_{i}(T_{ij_{1}})\Big{\\{}\phi_{r}(T_{ij_{1}})\phi_{k}(T_{ij_{2}})+\phi_{k}(T_{ij_{1}})\phi_{r}(T_{ij_{2}})\Big{\\}}$
$\displaystyle\qquad\times\phi_{l}(T_{ij_{2}})\left<\phi_{r},\Delta\phi_{k}\right>\xi_{il}$
$\displaystyle\lesssim
n^{-1}\sum_{i=1}^{n}m_{i}^{-1}\sum_{j_{1}=1}^{m_{i}}\underset{r\neq
k}{\sum_{r=1}^{\infty}}\sum_{l=1}^{\kappa_{0}}(\lambda_{k}-\lambda_{r})^{-1}\dot{\mu}_{i}(T_{ij_{1}})\phi_{r}(T_{ij_{1}})\left\\{\textbf{1}(l=k)+O((\log
m/m)^{1/2})\right\\}\left<\phi_{r},\Delta\phi_{k}\right>\xi_{il}$
$\displaystyle\qquad+n^{-1}\sum_{i=1}^{n}m_{i}^{-1}\sum_{j_{1}=1}^{m_{i}}\underset{r\neq
k}{\sum_{r=1}^{\infty}}\sum_{l=1}^{\kappa_{0}}(\lambda_{k}-\lambda_{r})^{-1}\dot{\mu}_{i}(T_{ij_{1}})\phi_{k}(T_{ij_{1}})\left\\{\textbf{1}(r=l)+O((\log
m/m)^{1/2})\right\\}\left<\phi_{r},\Delta\phi_{k}\right>\xi_{il}$
$\displaystyle\lesssim
n^{-1}\sum_{i=1}^{n}m_{i}^{-1}\sum_{j_{1}=1}^{m_{i}}\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}\dot{\mu}_{i}(T_{ij_{1}})\phi_{r}(T_{ij_{1}})\left<\phi_{r},\Delta\phi_{k}\right>\xi_{ik}\left\\{1+O((\log
m/m)^{1/2})\right\\}$
$\displaystyle\qquad+n^{-1}\sum_{i=1}^{n}m_{i}^{-1}\sum_{j_{1}=1}^{m_{i}}\underset{r\neq
k}{\sum_{r=1}^{\kappa_{0}}}(\lambda_{k}-\lambda_{r})^{-1}\dot{\mu}_{i}(T_{ij_{1}})\phi_{k}(T_{ij_{1}})\left<\phi_{r},\Delta\phi_{k}\right>\xi_{ir}\left\\{1+O((\log
m/m)^{1/2})\right\\}$ $\displaystyle:=(U_{k1}^{n}+U_{k2}^{n})\left\\{1+O((\log
m/m)^{1/2})\right\\}\qquad\text{almost surely}.$ (29)

Then applying the triangle inequality, we have

$\displaystyle U_{k1}^{n}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}m_{i}^{-1}\sum_{j=1}^{m_{i}}\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}\dot{\mu}_{i}(T_{ij})\phi_{r}(T_{ij})\left<\Delta\phi_{k},\phi_{r}\right>\xi_{ik}$
$\displaystyle\lesssim\|\Delta\phi_{k}\|\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}n^{-1}\sum_{i=1}^{n}\left\\{\mathcal{M}_{ir}+O((\log
m/m)^{1/2})\right\\}\xi_{ik}$ $\displaystyle=\|\Delta\phi_{k}\|\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}n^{-1}\sum_{i=1}^{n}\mathcal{M}_{ir}\xi_{ik}\left\\{1+O((\log
m/m)^{1/2})\right\\},$ (30)

where $\mathcal{M}_{ir}=\int\dot{\mu}_{i}(t)\phi_{r}(t)dt$. By Lemma 6 of Li
and Hsing [19], under conditions (C1)-(C5), for any measurable bounded
function $e$ on $[0,1]$,
$\|\Delta\phi_{k}\|=O(h^{2}+\delta_{n1}(h)+\delta^{2}_{n2}(h))\equiv
O(\overline{\delta}_{n}(h))$ almost surely, where
$\overline{\delta}_{n}(h)=h^{2}+\delta_{n1}(h)+\delta^{2}_{n2}(h)$. Thus, in
addition with Inequalities (19) in Lemma 5, we obtain
$U_{k1}^{n}=O\left(\overline{\delta}_{n}(h)(\log
n/n)^{1/2}\lambda_{k}^{1/2}k^{(1-\alpha)/2}\\{1+(\log m/m)^{1/2}\\}\right)$
almost surely. Next, under the spacing condition mentioned earlier and in
assumption (C6)a, using the Inequality (B), recall
$\eta_{k}^{-1}\lesssim\lambda_{k}^{-1}k$. Thus, observe that

$\displaystyle U_{k2}^{n}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}m_{i}^{-1}\sum_{j=1}^{m_{i}}\underset{r\neq
k}{\sum_{r=1}^{\kappa_{0}}}(\lambda_{k}-\lambda_{r})^{-1}\dot{\mu}_{i}(T_{ij})\phi_{k}(T_{i_{j}})\left<\Delta\phi_{k},\phi_{r}\right>\xi_{ir}$
$\displaystyle\lesssim\|\Delta\phi_{k}\|\underset{r\neq
k}{\sum_{r=1}^{\kappa_{0}}}(\lambda_{k}-\lambda_{r})^{-1}\left\\{n^{-1}\sum_{i=1}^{n}\mathcal{M}_{ik}\xi_{ir}\right\\}\left\\{1+O((\log
m/m)^{1/2})\right\\}$
$\displaystyle\lesssim\|\Delta\phi_{k}\|\eta_{k}^{-1}\underset{r\neq
k}{\sum_{r=1}^{\kappa_{0}}}\left\\{n^{-1}\sum_{i=1}^{n}\mathcal{M}_{ik}\xi_{ir}\right\\}\left\\{1+O((\log
m/m)^{1/2})\right\\}.$ (31)

Using condition (C6)b, we also have $V_{k}^{1/2}\eta_{k}^{-1}\lesssim
V_{k}^{1/2}\lambda_{k}^{-1}k=O(k^{(1-\alpha)/2})$. Finally, combining with the
bounds for $U_{k1}^{n},U_{k2}^{n}$, we have, almost surely,

$\displaystyle J_{k1}^{n2}$ $\displaystyle=O\Big{\\{}(\log
n/n)^{1/2}\overline{\delta}_{n}(h)\Big{(}\lambda_{k}^{1/2}k^{(1-\alpha)/2}+\eta_{k}^{-1}V_{k}^{1/2}\underset{r\neq
k}{\sum_{r=1}^{\kappa_{0}}}\lambda_{r}^{1/2}\Big{)}\left(1+(\log
m/m)^{1/2}\right)\Big{\\}}$ $\displaystyle:=O(\omega_{k1}(n,h)),$ (32)

where $\omega_{k1}(n,h)=(\log
n/n)^{1/2}\overline{\delta}_{n}(h)k^{(1-\alpha)/2}\sum_{k=1}^{\kappa_{0}}\lambda_{r}^{1/2}\left\\{1+(\log
m/m)^{1/2}\right\\}$.

It is easy to see that
$\sum_{r=1}^{\kappa_{0}}\lambda_{r}^{1/2}\sim\kappa_{0}^{-{\tau_{1}}/{2}+1}$.
Therefore, for $\tau=\alpha+\tau_{1}$, $\omega_{k1}(n,h)\sim(\log
n/n)^{1/2}\overline{\delta}_{n}(h)\kappa_{0}^{(3-\tau)/2}\left\\{1+(\log
m/m)^{1/2}\right\\}$. Similarly, to the derivation of the bound for
$J_{k1}^{n2}$, we can write

$\displaystyle J_{k2}^{n2}$
$\displaystyle={n}^{-1}\sum_{i=1}^{n}{m_{i}^{-2}}\sum_{j_{1}=1}^{m_{i}}\sum_{j_{2}=1}^{m_{i}}\sum_{l=1}^{\kappa_{0}}\dot{\mu}_{i}(T_{ij_{1}})I_{ik}^{n2}(T_{ij_{1}},T_{ij_{2}})\phi_{l}(T_{ij_{2}})\xi_{il}$
$\displaystyle={n}^{-1}\sum_{i=1}^{n}{m_{i}^{-2}}\sum_{j_{1}=1}^{m_{i}}\sum_{j_{2}=1}^{m_{i}}\sum_{l=1}^{\kappa_{0}}\underset{r_{1}\neq
k}{\sum_{r_{1}=1}^{\infty}}\underset{r_{2}\neq
k}{\sum_{r_{1}=1}^{\infty}}(\lambda_{k}-\lambda_{r_{1}})^{-1}(\lambda_{k}-\lambda_{r_{2}})^{-1}\left<\phi_{r_{1}},\Delta\phi_{k}\right>\left<\phi_{r_{2}},\Delta\phi_{k}\right>$
$\displaystyle\qquad\times\dot{\mu}_{i}(T_{ij_{1}})\phi_{r_{1}}(T_{ij_{1}})\phi_{r_{2}}(T_{ij_{2}})\phi_{l}(T_{ij_{2}})\xi_{il}$
$\displaystyle\lesssim{n}^{-1}\sum_{i=1}^{n}{m_{i}}^{-1}\sum_{j_{1}=1}^{m_{i}}\sum_{l=1}^{\kappa_{0}}\underset{r_{1}\neq
k}{\sum_{r_{1}=1}^{\infty}}\underset{r_{2}\neq
k}{\sum_{r_{2}=1}^{\infty}}(\lambda_{k}-\lambda_{r_{1}})^{-1}(\lambda_{k}-\lambda_{r_{2}})^{-1}\left<\phi_{r_{1}},\Delta\phi_{k}\right>\left<\phi_{r_{2}},\Delta\phi_{k}\right>$
$\displaystyle\qquad\times\dot{\mu}_{i}(T_{ij_{1}})\phi_{r_{1}}(T_{ij_{1}})\left\\{\textbf{1}(r_{2}=l)+O((\log
m/m)^{1/2})\right\\}\xi_{il}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}{m_{i}}^{-1}\sum_{j_{1}=1}^{m_{i}}\underset{r_{1}\neq
k}{\sum_{r_{1}=1}^{\infty}}\underset{r_{2}\neq
k}{\sum_{r_{2}=1}^{\kappa_{0}}}(\lambda_{k}-\lambda_{r_{1}})^{-1}(\lambda_{k}-\lambda_{r_{2}})^{-1}\left<\phi_{r_{1}},\Delta\phi_{k}\right>\left<\phi_{r_{2}},\Delta\phi_{k}\right>$
$\displaystyle\qquad\times\dot{\mu}_{i}(T_{ij_{1}})\phi_{r_{1}}(T_{ij_{1}})\xi_{ir_{2}}\left\\{1+O((\log
m/m)^{1/2})\right\\}$
$\displaystyle\lesssim\|\Delta\phi_{k}\|^{2}\underset{r_{1}\neq
k}{\sum_{r_{1}=1}^{\infty}}\underset{r_{2}\neq
k}{\sum_{r_{2}=1}^{\kappa_{0}}}(\lambda_{k}-\lambda_{r_{1}})^{-1}(\lambda_{k}-\lambda_{r_{2}})^{-1}\left({n}^{-1}\sum_{i=1}^{n}\left\\{\mathcal{M}_{ir_{1}}+O((\log
m/m)^{1/2})\right\\}\xi_{ir_{2}}\right)$
$\displaystyle\qquad\times\left\\{1+O((\log m/m)^{1/2})\right\\}$
$\displaystyle=\|\Delta\phi_{k}\|^{2}\underset{r_{1}\neq
k}{\sum_{r_{1}=1}^{\infty}}\underset{r_{2}\neq
k}{\sum_{r_{2}=1}^{\kappa_{0}}}(\lambda_{k}-\lambda_{r_{1}})^{-1}(\lambda_{k}-\lambda_{r_{2}})^{-1}\times\left\\{{n}^{-1}\sum_{i=1}^{n}\mathcal{M}_{ir_{1}}\xi_{ir_{2}}\right\\}\left\\{1+O((\log
m/m)^{1/2})\right\\}.$ (33)

Therefore, Inequality (C) immediately follows using Lemma 6,

$\displaystyle J_{k2}^{n2}$ $\displaystyle=O\left((\log
n/n)^{1/2}\overline{\delta}^{2}_{n}(h)\kappa_{0}^{(3-\alpha)/2}\lambda_{\kappa_{0}}^{-1}\\{\sum_{r=1}^{\kappa_{0}}\lambda_{r}\\}^{1/2}\left\\{1+(\log
m/m)^{1/2}\right\\}\right)$
$\displaystyle:=O(\omega_{k2}(n,h)),\qquad\text{almost surely}.$

where

$\omega_{k2}(n,h)=(\log
n/n)^{1/2}\overline{\delta}^{2}_{n}(h)\kappa_{0}^{(3-\alpha)/2}\lambda_{\kappa_{0}}^{-1}\sum_{r=1}^{\kappa_{0}}\lambda_{r}\left\\{1+(\log
m/m)^{1/2}\right\\}.$

We observe that
$\sum_{r=1}^{\kappa_{0}}\lambda_{r}\sim\kappa_{0}^{-\tau_{1}+1}$ under
assumption (C6)b. Thus, $\omega_{k2}(n,h)\sim(\log
n/n)\overline{\delta}^{2}(h)\kappa_{0}^{(4-\tau)/2}\\{1+(\log m/m)^{1/2}\\}$.
Since $\omega_{n2}=O(\omega_{n1})$ and
$\overline{\delta}_{n}^{2}(h)=O(\omega_{n1})$, in summary, for each
$k=1,\cdots,\kappa_{0}$,

$\overline{\mathbf{g}}_{k}(\boldsymbol{\beta})=O\left((\log
n/n)^{1/2}\left\\{1+(\log m/m)^{1/2}\right\\}+\omega_{k1}(n,h)\right)$

almost surely. Since $1/(nm)=O(1/n)$,

$\text{AMSE}\\{\overline{\mathbf{g}}_{k}(\boldsymbol{\beta}_{0})\\}=O\left(n^{-1}+n^{-1}\kappa_{0}^{3-\tau}R_{n}(h)\right),$

where
$R_{n}(h)=\left\\{h^{4}+\frac{1}{n}+\frac{1}{nmh}+\frac{1}{n^{2}m^{2}h^{2}}+\frac{1}{n^{2}m^{4}h^{4}}+\frac{1}{n^{2}mh}+\frac{1}{n^{2}m^{3}h^{3}}\right\\}$.
Combining the above conditions, we find that if $a>1/4$,
$\kappa_{0}=O(n^{1/(3-\tau)})$ and $n^{-1/4}\lesssim h\lesssim n^{-(a+1)/5}$
then
$\text{AMSE}\\{\overline{\mathbf{g}_{k}}(\boldsymbol{\beta}_{0})\\}=O(1/n)$.
On the other hand, if $a\leq 1/4$, $\kappa_{0}=O(n^{4(1+a)/5(3-\tau)})$ and
$h\lesssim n^{-1/4}$
$\text{AMSE}\\{\overline{\mathbf{g}_{k}}(\boldsymbol{\beta}_{0})\\}=O(1/n)$.

Note that for a three-dimensional array $(\partial
C/\partial\beta_{1},\cdots,\partial C/\partial\beta_{p})$ such that the
following is a $p\times 1$ vector.

$\overline{\mathbf{g}}(\boldsymbol{\beta}_{0})^{\text{T}}\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\dot{\mathbf{C}}(\boldsymbol{\beta}_{0})\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\overline{\mathbf{g}}(\boldsymbol{\beta}_{0})$
(34)

Therefore,

$n^{-1}\dot{\mathcal{Q}}(\boldsymbol{\beta}_{0})=2\dot{\overline{\mathbf{g}}}(\boldsymbol{\beta}_{0})^{\text{T}}\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\overline{\mathbf{g}}(\boldsymbol{\beta}_{0})-\overline{\mathbf{g}}(\boldsymbol{\beta}_{0})^{\text{T}}\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\dot{\mathbf{C}}(\boldsymbol{\beta}_{0})\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\overline{\mathbf{g}}(\boldsymbol{\beta}_{0})$
(35)

and

$n^{-1}\ddot{\mathcal{Q}}(\boldsymbol{\beta}_{0})=2\dot{\overline{\mathbf{g}}}(\boldsymbol{\beta}_{0})^{\text{T}}\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\dot{\overline{\mathbf{g}}}(\boldsymbol{\beta}_{0})+r_{n1}+r_{n2}+r_{n3}+r_{n4}$
(36)

where

$\begin{split}r_{n1}&=2\ddot{\overline{\mathbf{g}}}(\boldsymbol{\beta}_{0})\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\overline{\mathbf{g}}(\boldsymbol{\beta}_{0})\\\
r_{n2}&=-4\dot{\overline{\mathbf{g}}}(\boldsymbol{\beta}_{0})^{\text{T}}\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\dot{\mathbf{C}}(\boldsymbol{\beta}_{0})\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\overline{\mathbf{g}}(\boldsymbol{\beta}_{0})\\\
r_{n3}&=2\dot{\overline{\mathbf{g}}}(\boldsymbol{\beta}_{0})\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\dot{\mathbf{C}}(\boldsymbol{\beta}_{0})\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\overline{\mathbf{g}}(\boldsymbol{\beta}_{0})\\\
r_{n4}&=-\overline{\mathbf{g}}(\boldsymbol{\beta}_{0})^{\text{T}}\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\ddot{\mathbf{C}}(\boldsymbol{\beta}_{0})\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\overline{\mathbf{g}}(\boldsymbol{\beta}_{0})\end{split}$
(37)

Since $\overline{\mathbf{g}}(\boldsymbol{\beta}_{0})=O_{P}(n^{-1/2})$ and the
weight matrix converges almost surely to an invertible matrix,
$\overline{\mathbf{g}}(\boldsymbol{\beta}_{0})^{\text{T}}\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\dot{\mathbf{C}}(\boldsymbol{\beta}_{0})\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\overline{\mathbf{g}}(\boldsymbol{\beta}_{0})=o(n^{-1})$
almost surely. Furthermore, $r_{n1}=O(n^{-1/2})$, $r_{n2}=o(n^{-1/2})$,
$r_{n3}=O(n^{-1/2})$, and $r_{n4}=O(n^{-1})$ almost surely. Combining these
bounds, we have $r_{n}=o(1)$ almost surely. Therefore,
$\|n^{-1}\dot{\mathcal{Q}}-2\dot{\overline{\mathbf{g}}}(\boldsymbol{\beta}_{0})^{\text{T}}\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\overline{\mathbf{g}}(\boldsymbol{\beta}_{0})\|=o_{P}(n^{-1})$
and
$\|n^{-1}\ddot{\mathcal{Q}}-2\dot{\overline{\mathbf{g}}}(\boldsymbol{\beta}_{0})^{\text{T}}\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\dot{\overline{\mathbf{g}}}(\boldsymbol{\beta}_{0})\|=o_{P}(1)$.

The following lines are based on common steps in the GEE literature that
includes McCullagh and Nelder [24], Balan et al. [2], Tian et al. [27] among
many others. Let
$\boldsymbol{\beta}_{n}=\boldsymbol{\beta}_{0}+\delta\mathbf{d}$ where set
$\delta=n^{-1/2}$. We have to show that for any $\epsilon>0$ there exists a
large constant $c$ such that

$P\\{\inf_{\|d\|=c}\mathcal{Q}(\boldsymbol{\beta}_{n})\geq\mathcal{Q}(\boldsymbol{\beta}_{0})\\}>1-\epsilon$
(38)

Note that the above statement is always true if $\epsilon\geq 1$. Thus, we
assume that $\epsilon\in(0,1)$. Due to Taylor series expansion,

$\mathcal{Q}(\boldsymbol{\beta}_{n})=\mathcal{Q}(\boldsymbol{\beta}_{0}+\delta\mathbf{d})=\mathcal{Q}(\boldsymbol{\beta}_{0})+\delta\mathbf{d}^{\text{T}}\dot{\mathcal{Q}}(\boldsymbol{\beta}_{0})+0.5\delta\mathbf{d}^{\text{T}}\ddot{\mathcal{Q}}(\boldsymbol{\beta}_{0})\mathbf{d}+\|\mathbf{d}\|^{2}o_{P}(1)$
(39)

Now, observe that, using Equation (35),

$\delta\mathbf{d}^{\text{T}}\dot{\mathcal{Q}}(\boldsymbol{\beta}_{0})=\|\mathbf{d}\|O_{P}(\sqrt{n}\delta)+\|\mathbf{d}\|O_{P}(\delta)$
(40)

and

$0.5\delta^{2}\mathbf{d}^{\text{T}}\ddot{\mathcal{Q}}(\boldsymbol{\beta}^{*})\mathbf{d}=n\delta^{2}\mathbf{d}^{\text{T}}\dot{\overline{\mathbf{g}}}(\boldsymbol{\beta}_{0})\mathbf{C}^{-1}(\boldsymbol{\beta}_{0})\dot{\mathbf{g}}(\boldsymbol{\beta}_{0})\mathbf{d}+n\delta^{2}\|\mathbf{d}\|^{2}o_{P}(1)$
(41)

Therefore, for given $\epsilon>0$, there exists a large enough $c$ such that
the above equation (38) holds. This implies that there exists a
$\hat{\boldsymbol{\beta}}$ that satisfies
$\|\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}_{0}\|=O_{P}(\delta)$. Thus,
for large $n$, with probability 1, $\mathcal{Q}(\boldsymbol{\beta})$ attains
the minimal value at $\hat{\boldsymbol{\beta}}$ and therefore,
$\dot{\mathcal{Q}}=0$.

##  Appendix D Proof of Theorem 2

Recall,
$\mathcal{C}_{i}=\sum_{k_{1}=1}^{\kappa_{0}}\sum_{k_{2}=1}^{\kappa_{0}}\boldsymbol{\Phi}_{k_{1}}\mathbf{X}_{i}\mathbf{C}_{k_{1},k_{2}}^{-1}\mathbf{X}_{i}^{\text{T}}\boldsymbol{\Phi}_{k_{2}}$,
where $\mathbf{C}_{k_{1},k_{2}}^{-1}$ is the $(k_{1},k_{2})$ block of
$\mathbf{C}_{0}^{-1}$. Similarly, we can define
$\widehat{\mathcal{C}}_{i}=\sum_{k_{1}=1}^{\kappa_{0}}\sum_{k_{2}=1}^{\kappa_{0}}\widehat{\boldsymbol{\Phi}}_{k_{1}}\mathbf{X}_{i}\mathbf{C}_{k_{1},k_{2}}^{-1}\mathbf{X}_{i}^{\text{T}}\widehat{\boldsymbol{\Phi}}_{k_{2}}$.
It is easy to observe that
$\widehat{\mathcal{C}}_{i}=\mathcal{C}_{i}+\sum_{k_{1}=1}^{\kappa_{0}}\sum_{k_{2}=1}^{\kappa_{0}}(\widehat{\boldsymbol{\Phi}}_{k_{1}}-\boldsymbol{\Phi}_{k_{1}})\mathbf{X}_{i}\mathbf{C}_{k_{1},k_{2}}^{-1}\mathbf{X}_{i}^{\text{T}}(\widehat{\boldsymbol{\Phi}}_{k_{2}}-\boldsymbol{\Phi}_{k_{2}})+2\sum_{k_{1}=1}^{\kappa_{0}}\sum_{k_{2}=1}^{\kappa_{0}}\boldsymbol{\Phi}_{k_{1}}\mathbf{X}_{i}\mathbf{C}_{k_{1},k_{2}}^{-1}\mathbf{X}_{i}^{\text{T}}(\widehat{\boldsymbol{\Phi}}_{k_{2}}-\boldsymbol{\Phi}_{k_{2}})$.
Therefore,
$n^{-1}\sum_{i=1}^{n}\dot{\boldsymbol{\mu}}_{i}^{\text{T}}(\widehat{\mathcal{C}}_{i}-\mathcal{C})\mathbf{X}_{i}=n^{-1}\sum_{i=1}^{n}\sum_{k_{1}=1}^{\kappa_{0}}\sum_{k_{2}=1}^{\kappa_{0}}\mathbf{P}_{ik_{1}}\mathbf{C}_{k_{1},k_{2}}^{-1}\mathbf{P}_{ik_{2}}$
and
$n^{-1}\sum_{i=1}^{n}\dot{\boldsymbol{\mu}}_{i}^{\text{T}}(\widehat{\mathcal{C}}_{i}-\mathcal{C})\mathbf{y}_{i}=n^{-1}\sum_{i=1}^{n}\sum_{k_{1}=1}^{\kappa_{0}}\sum_{k_{2}=1}^{\kappa_{0}}\mathbf{P}_{ik_{1}}\mathbf{C}_{k_{1},k_{2}}^{-1}\mathbf{Q}_{ik_{2}}$
where
$\mathbf{P}_{i,k}=\dot{\boldsymbol{\mu}}_{i}^{\text{T}}\mathbf{D}_{ik}\mathbf{X}_{i}$
and
$\mathbf{Q}_{ik}=\dot{\boldsymbol{\mu}}_{i}^{\text{T}}\mathbf{D}_{ik}\mathbf{y}_{i}$
with $\mathbf{D}_{ik}$ be the difference matrix with $(j_{1},j_{2})$-th
element is $d_{i}(T_{ij_{1}},T_{ij_{2}})$. Thus, note that, almsot surely we
have the following relation,

$\displaystyle\mathbf{P}_{ik}$
$\displaystyle=m_{i}^{-2}\sum_{j_{1}=1}^{m_{i}}\sum_{j_{2}=1}^{m_{i}}\dot{\mu}_{i}(T_{ij_{1}})d_{i}(T_{ij_{1}},T_{ij_{2}})x_{i}(T_{ij_{2}})$
$\displaystyle\lesssim
m_{i}^{-2}\sum_{j_{1}=1}^{m_{i}}\sum_{j_{2}=1}^{m_{i}}\dot{\mu}_{i}(T_{ij_{1}})\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}\left<\phi_{r},\Delta\phi_{k}\right>\left\\{\phi_{r}(T_{ij_{1}})\phi_{k}(T_{ij_{2}})+\phi_{k}(T_{ij_{1}})\phi_{r}(T_{ij_{2}})\right\\}$
$\displaystyle\qquad+O(\|\Delta\|^{2})$ $\displaystyle\lesssim\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}\|\Delta\phi_{k}\|+O(\|\Delta\|^{2})=O(\varpi),$
since $m_{i}^{-1}\sum_{j=1}^{m_{i}}\dot{\mu}(T_{ij})\phi_{k}(T_{ij})$ and
$m_{i}^{-1}\sum_{j=1}^{m_{i}}x_{i}(T_{ij})\phi_{k}(T_{ij})$ are finite

where $\varpi=\underset{r\neq
k}{\sum_{r=1}^{\infty}}(\lambda_{k}-\lambda_{r})^{-1}\overline{\delta}_{n}(h)+h^{2}+\delta_{n2}^{2}(h)$.
A similar result can be obtained for $\mathbf{Q}_{ik}$. Combining such
results, in summary, we have
$-{2}{n^{-1}}\sum_{i=1}^{n}\mathbf{X}_{i}^{\text{T}}\widehat{\mathcal{C}}_{i}(\mathbf{y}_{i}-\mathbf{X}_{i}^{\text{T}}\widehat{\boldsymbol{\beta}})=-{2}{n^{-1}}\sum_{i=1}^{n}\mathbf{X}_{i}^{\text{T}}\mathcal{C}_{i}(\mathbf{y}_{i}-\mathbf{X}_{i}^{\text{T}}\widehat{\boldsymbol{\beta}})+O(\varpi_{n})$.
Since, for $n\rightarrow\infty$, $\mathcal{Q}(\boldsymbol{\beta})$ attains a
minimal value at $\boldsymbol{\beta}=\widehat{\boldsymbol{\beta}}$, we
therefore have $\dot{\mathcal{Q}}(\widehat{\boldsymbol{\beta}})=0$. Thus,

$\displaystyle\dot{\mathcal{Q}}(\widehat{\boldsymbol{\beta}})=-{2}{n^{-1}}\sum_{i=1}^{n}\mathbf{X}_{i}^{\text{T}}\widehat{\mathcal{C}}_{i}(\mathbf{y}_{i}-\mathbf{X}_{i}^{\text{T}}\widehat{\boldsymbol{\beta}})$
$\displaystyle=-{2}{n^{-1}}\sum_{i=1}^{n}\mathbf{X}_{i}^{\text{T}}(\widehat{\mathcal{C}}_{i}-\mathcal{C}_{i})(\mathbf{y}_{i}-\mathbf{X}_{i}^{\text{T}}\widehat{\boldsymbol{\beta}})-{2}{n^{-1}}\sum_{i=1}^{n}\mathbf{X}_{i}^{\text{T}}{\mathcal{C}}_{i}(\mathbf{y}_{i}-\mathbf{X}_{i}^{\text{T}}\widehat{\boldsymbol{\beta}})=0.$

Therefore, almost surely, we have,

$\displaystyle-{2}{n^{-1}}\sum_{i=1}^{n}\mathbf{X}_{i}^{\text{T}}\mathcal{C}_{i}(\mathbf{y}_{i}-\mathbf{X}_{i}^{\text{T}}\widehat{\boldsymbol{\beta}})+O(\varpi_{n})$
$\displaystyle=0$
$\displaystyle-{2}{n^{-1}}\sum_{i=1}^{n}\mathbf{X}_{i}^{\text{T}}\mathcal{C}_{i}(\mathbf{X}_{i}^{\text{T}}\boldsymbol{\beta}_{0}+\mathbf{e}_{i}-\mathbf{X}_{i}^{\text{T}}\widehat{\boldsymbol{\beta}})+O(\varpi_{n})$
$\displaystyle=0$
$\displaystyle\sqrt{n}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}_{0})\left\\{n^{-1}\sum_{i=1}^{n}\mathbf{X}_{i}^{\text{T}}\mathcal{C}_{i}\mathbf{X}_{i}\right\\}$
$\displaystyle={n^{-1}}\sum_{i=1}^{n}\mathbf{X}_{i}^{\text{T}}\mathcal{C}_{i}\mathbf{e}_{i}$
(44)

Now, using the central limit theorem, we can obtain the following.

$\sum_{i=1}^{n}\mathbf{X}_{i}^{\text{T}}\mathcal{C}_{i}\mathbf{e}_{i}/{\sqrt{n}}\xrightarrow[]{d}N(0,\mathbf{A})$
(45)

In addition, by the law of large numbers
$n^{-1}\sum_{i=1}^{n}\mathbf{X}_{i}^{\text{T}}\mathcal{C}_{i}\mathbf{X}_{i}\rightarrow\mathbf{B}$
in probability. Therefore, using the Slutsky theorem, we complete the proof of
Theorem 2.

## References

  * Bai et al. [2008] Y. Bai, Z. Zhu, W. K. Fung, Partial linear models for longitudinal data based on quadratic inference functions, Scandinavian Journal of Statistics 35 (2008) 104–118.
  * Balan et al. [2005] R. M. Balan, I. Schiopu-Kratina, et al., Asymptotic results with generalized estimating equations for longitudinal data, The Annals of Statistics 33 (2005) 522–541.
  * Chen et al. [2019] X. Chen, H. Li, H. Liang, H. Lin, Functional response regression analysis, Journal of Multivariate Analysis 169 (2019) 218–233.
  * Chiou et al. [2003] J.-M. Chiou, H.-G. Müller, J.-L. Wang, Functional quasi-likelihood regression models with smooth random effects, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 65 (2003) 405–423.
  * Dauxois et al. [1982] J. Dauxois, A. Pousse, Y. Romain, Asymptotic theory for the principal component analysis of a vector random function: some applications to statistical inference, Journal of multivariate analysis 12 (1982) 136–154.
  * Diggle et al. [2002] P. Diggle, P. J. Diggle, P. Heagerty, K.-Y. Liang, P. J. Heagerty, S. Zeger, et al., Analysis of longitudinal data, Oxford University Press, 2002.
  * Fan and Gijbels [1996] J. Fan, I. Gijbels, Local polynomial modelling and its applications, Chapman & Hall/CRC, 1996.
  * Fan et al. [1999] J. Fan, W. Zhang, et al., Statistical estimation in varying coefficient models, The annals of Statistics 27 (1999) 1491–1518.
  * Gajardo et al. [2021] A. Gajardo, C. Carroll, Y. Chen, X. Dai, J. Fan, P. Z. Hadjipantelis, K. Han, H. Ji, H.-G. Mueller, J.-L. Wang, fdapace: Functional Data Analysis and Empirical Dynamics, 2021. R package version 0.5.7.
  * Givens and Hoeting [2012] G. H. Givens, J. A. Hoeting, Computational statistics, volume 703, John Wiley & Sons, 2012.
  * Hall et al. [2007] P. Hall, J. L. Horowitz, et al., Methodology and convergence rates for functional linear regression, The Annals of Statistics 35 (2007) 70–91.
  * Hall and Hosseini-Nasab [2006] P. Hall, M. Hosseini-Nasab, On properties of functional principal components analysis, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68 (2006) 109–126.
  * Hall and Hosseini-Nasab [2009] P. Hall, M. Hosseini-Nasab, Theory for high-order bounds in functional principal components analysis, in: Mathematical Proceedings of the Cambridge Philosophical Society, volume 146, Cambridge University Press, p. 225.
  * Hand and Crowder [2017] D. Hand, M. Crowder, Practical longitudinal data analysis, Routledge, 2017.
  * Hansen [1982] L. P. Hansen, Large sample properties of generalized method of moments estimators, Econometrica: Journal of the Econometric Society (1982) 1029–1054.
  * Hsing and Eubank [2015] T. Hsing, R. Eubank, Theoretical foundations of functional data analysis, with an introduction to linear operators, volume 997, John Wiley & Sons, 2015.
  * J Mercer [1909] B. J Mercer, Xvi. functions of positive and negative type, and their connection the theory of integral equations, Phil. Trans. R. Soc. Lond. A 209 (1909) 415–446.
  * Karhunen [1946] K. Karhunen, Zur spektraltheorie stochastischer prozesse, Ann. Acad. Sci. Fennicae, AI 34 (1946).
  * Li and Hsing [2010] Y. Li, T. Hsing, Uniform convergence rates for nonparametric regression and principal component analysis in functional/longitudinal data, The Annals of Statistics 38 (2010) 3321–3351.
  * Li et al. [2022] Y. Li, Y. Qiu, Y. Xu, From multivariate to functional data analysis: Fundamentals, recent developments, and emerging areas, Journal of Multivariate Analysis 188 (2022) 104806.
  * Liang and Zeger [1986] K.-Y. Liang, S. L. Zeger, Longitudinal data analysis using generalized linear models, Biometrika 73 (1986) 13–22.
  * Liang et al. [2015] X. Liang, T. Zou, B. Guo, S. Li, H. Zhang, S. Zhang, H. Huang, S. X. Chen, Assessing beijing’s pm2. 5 pollution: severity, weather impact, apec and winter heating, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 471 (2015) 20150257.
  * Loève [1946] M. Loève, Functions aleatoire de second ordre, Revue science 84 (1946) 195–206.
  * McCullagh and Nelder [1989] P. McCullagh, J. Nelder, Generalized linear models ii, 1989.
  * Qu and Li [2006] A. Qu, R. Li, Quadratic inference functions for varying-coefficient models with longitudinal data, Biometrics 62 (2006) 379–391.
  * Qu et al. [2000] A. Qu, B. G. Lindsay, B. Li, Improving generalised estimating equations using quadratic inference functions, Biometrika 87 (2000) 823–836.
  * Tian et al. [2014] R. Tian, L. Xue, C. Liu, Penalized quadratic inference functions for semiparametric varying coefficient partially linear models with longitudinal data, Journal of Multivariate Analysis 132 (2014) 94–110.
  * Turner et al. [2017] E. L. Turner, F. Li, J. A. Gallis, M. Prague, D. M. Murray, Review of recent methodological developments in group-randomized trials: part 1—design, American journal of public health 107 (2017) 907–915.
  * Xiong et al. [2017] Y. Xiong, X. J. Zhou, R. A. Nisi, K. R. Martin, M. M. Karaman, K. Cai, T. E. Weaver, Brain white matter changes in cpap-treated obstructive sleep apnea patients with residual sleepiness, Journal of Magnetic Resonance Imaging 45 (2017) 1371–1378.
  * Yao et al. [2005] F. Yao, H.-G. Müller, J.-L. Wang, Functional linear regression analysis for longitudinal data, The Annals of Statistics (2005) 2873–2903.
  * Yu et al. [2020] H. Yu, G. Tong, F. Li, A note on the estimation and inference with quadratic inference functions for correlated outcomes, Communications in Statistics-Simulation and Computation (2020) 1–12.
  * Zhang and Li [2022] H. Zhang, Y. Li, Unified principal component analysis for sparse and dense functional data under spatial dependency, Journal of Business & Economic Statistics 40 (2022) 1523–1537.
  * Zhang and Banerjee [2021] L. Zhang, S. Banerjee, Spatial factor modeling: A bayesian matrix-normal approach for misaligned data, Biometrics (2021).
  * Zhang and Wang [2016] X. Zhang, J.-L. Wang, From sparse to dense functional data and beyond, The Annals of Statistics 44 (2016) 2281–2321.
  * Zhou and Qu [2012] J. Zhou, A. Qu, Informative estimation and selection of correlation structure for longitudinal data, Journal of the American Statistical Association 107 (2012) 701–710.