# Large-scale multiple testing of cross-covariance functions with applications
to functional network models

Qin Fang  University of Sydney Business School, Sydney, Australia Qing Jiang
Center for Statistics and Data Science, Beijing Normal University, Zhuhai,
China  Xinghao Qiao  Faculty of Business and Economics, The University of Hong
Kong, Hong Kong

###### Abstract

The estimation of functional networks through functional covariance and
graphical models have recently attracted increasing attention in settings with
high dimensional functional data, where the number of functional variables $p$
is comparable to, and maybe larger than, the number of subjects. In this
paper, we first reframe the functional covariance model estimation as a
tuning-free problem of simultaneously testing $p(p-1)/2$ hypotheses for cross-
covariance functions. Our procedure begins by constructing a Hilbert–Schmidt-
norm-based test statistic for each pair, and employs normal quantile
transformations for all test statistics, upon which a multiple testing step is
proposed. We then explore the multiple testing procedure under a general
error-contamination framework and establish that our procedure can control
false discoveries asymptotically. Additionally, we demonstrate that our
proposed methods for two concrete examples: the functional covariance model
with partial observations and, importantly, the more challenging functional
graphical model, can be seamlessly integrated into the general error-
contamination framework, and, with verifiable conditions, achieve theoretical
guarantees on effective false discovery control. Finally, we showcase the
superiority of our proposals through extensive simulations and functional
connectivity analysis of two neuroimaging datasets.

Keywords:  False discovery control; Functional covariance model; Functional
graphical model; High-dimensional functional data; Partially observed
functional data; Power.

##  1 Introduction

Recent advances in information technology have led to the growing prevalence
of multivariate or even high-dimensional functional datasets across various
applications. Examples include time-course gene expression data in genomics,
and various types of brain imaging data in neuroscience, such as
electrocorticography, magnetoencephalography and functional magnetic resonance
imaging scans, among others. The brain signals are typically recorded in the
form of a location-by-time matrix for each individual subject, where the rows
and columns correspond to a set of brain regions and a number of observational
time points spanning over minutes, respectively. To capture the non-stationary
and dynamic patterns, recent proposals involve modeling brain signals as
multivariate random functions, treating the time course data of each region as
a random function (e.g., Zapata et al., 2022; Lee et al., 2023; Xue et al.,
2023). Such high-dimensional functional data can be represented as ${\bf
X}_{i}(\cdot)=\\{X_{i1}(\cdot),\dots,X_{ip}(\cdot)\\}^{\mathrm{\scriptscriptstyle
T}}$ defined on a compact interval ${\cal U}$ with marginal- and cross-
covariance functions, which together form the covariance function matrix

$\boldsymbol{\Sigma}(\cdot,\cdot)=\\{\Sigma_{jk}(\cdot,\cdot)\\}_{p\times
p},~{}~{}~{}\Sigma_{jk}(u,v)=\text{Cov}\\{X_{ij}(u),X_{ik}(v)\\}~{}~{}\text{for}~{}~{}(u,v)\in{\cal
U}^{2}.$

We observe ${\bf X}_{i}(\cdot)$ for $i=1,\dots,n,$ where the dimension $p$ is
large relative to, and maybe greater than the number of subjects $n.$

Our motivation lies in the functional connectivity analysis based on brain
imaging data. The literature primarily focuses on two types of multivariate
functional data analysis approaches for estimating brain connectivity
networks, consisting of $p$ nodes. The first method, proposed by Fang et al.
(2023), considers a functional covariance model depicting the marginal
correlation information in ${\bf X}_{i}(\cdot).$ This approach aims to
identify the functional sparsity structure in
$\boldsymbol{\Sigma}(\cdot,\cdot),$ i.e., discovering edges $(j,k)$’s such
that $\Sigma_{jk}(u,v)\neq 0$ for some $(u,v)\in{\cal U}^{2}.$ By applying
adaptive functional thresholding to the entries of the sample version of
$\boldsymbol{\Sigma}(\cdot,\cdot),$ they achieve estimation and support
recovery consistencies. The second method frames the network estimation as a
functional graphical modelling problem Zhu et al. (2016); Qiao et al. (2019).
This model characterizes the conditional dependence structure of $p$ Gaussian
random functions, i.e., nodes $j$ and $k$ are connected by an edge if and only
if $X_{ij}(\cdot)$ and $X_{ik}(\cdot)$ are dependent, conditional on the
remaining $p-2$ functions. This line of work has witnessed numerous
advancements, see, e.g., Li and Solea (2018); Zapata et al. (2022); Solea and
Li (2022); Lee et al. (2023) and Tsai et al. (2024). Both types of models
involve developing the regularized estimation of functional network structures
that leverage the functional sparsity information. However, the selection of
practically reasonable regularization parameters poses a challenging task.
While large regularization parameters result in sparse networks and may not be
effective in discovering edges with small weights, small regularization
parameters can produce excessive false edges, leading to high false discovery
rates, as evidenced by simulation results in Section 5.

In this paper, we reframe the functional covariance model estimation as a
tuning-free problem of simultaneously testing $p(p-1)/2$ hypotheses for cross-
covariance functions

$H_{0,jk}:\Sigma_{jk}(u,v)=0\text{ for any }(u,v)\in{\cal
U}^{2}~{}~{}\text{vs}~{}~{}H_{1,jk}:\Sigma_{jk}(u,v)\neq 0\text{ for some
}(u,v)\in{\cal U}^{2},$ (1)

where $1\leq j<k\leq p$ and we identify a significant edge between nodes $j$
and $k$ if and only if $H_{0,jk}$ is rejected. Compared to the multiple
testing for the local covariance structures of high-dimensional vectors Cai
(2017), the infinite-dimensional nature of functional data introduces
additional methodological and theoretical complexities that arise from
constructing the test statistic under a functional norm and addressing various
technical obstacles within an abstract Hilbert space. Furthermore, when
handling the practical scenario for brain imaging data with each trajectory
$X_{ij}(\cdot)$ being observed over a densely sampled grid, nonparametric
smoothing is frequently employed to obtain estimated curves
$\widehat{X}_{ij}(\cdot),$ which requires considering functional error-
contaminated versions of $X_{ij}(\cdot)$ satisfying:

$\widehat{X}_{ij}(\cdot)=X_{ij}(\cdot)+e_{ij}(\cdot),~{}~{}i=1,\dots,n,~{}j=1,\dots,p.$
(2)

Additionally, we present a novel procedure based on nodewise functional
regressions that converts the functional graphical model estimation into the
multiple testing task for the cross-covariance structures of functional
regression errors. As a result, we need to deal with estimated functional
regression errors (i.e., functional residuals) instead of true ones, following
a similar form as (2). In both scenarios, accounting for the functional
errors, such as $e_{ij}(\cdot)$’s in (2), introduces an extra layer of
complexity to our theoretical analysis.

Our paper aims to develop a general framework for large-scale multiple testing
of cross-covariance functions from both methodological and theoretical
perspectives. We begin by constructing a Hilbert–Schmidt-norm-based test
statistic for each pair $j<k,$ which is shown to have asymptotically the same
limiting null distribution as an infinite mixture of chi-squares. We then
employ normal quantile transformations for all test statistics, upon which a
multiple testing procedure is proposed to account for the multiplicity and
dependence among the transformed test statistics. Theoretically, we establish
that our procedure can control false discoveries asymptotically under both
fully observed and error-contaminated functional scenarios, and, furthermore,
be applied to two concrete examples: the functional covariance model with
partial observations and the more challenging functional graphical model.
Specifically, we demonstrate the seamless integration of our proposed methods
for both examples into the general error-contamination framework, and, with
verifiable conditions, achieve theoretical guarantees on false discovery
control. Empirically, we conduct simulations to showcase the uniform
superiority of our proposals over competitors in terms of false discovery
control and power for both fully and partially observed functional data within
both functional covariance and graphical models. We also apply our method to
identify network structures using two brain imaging datasets, and observe
scientifically interpretable patterns.

The main contributions of our paper are threefold.

  * •

First, we make the first attempt in the literature of functional data analysis
and multiple testing to develop a general procedure with theoretical
guarantees for the simultaneous testing of a large collection of hypotheses
for the functional sparse covariance structures. Our approach is fully
functional in the sense it does not rely on dimension reduction techniques,
thereby avoiding any incurred information loss.

  * •

Second, we extend our method and theory to the more general error-
contamination setting (2), encompassing the common practical scenario of
partially observed functional data as a special case. On the method front, our
procedure remains valid by simply replacing each $X_{ij}(\cdot)$ with its
estimated surrogate $\widehat{X}_{ij}(\cdot).$ Theoretically, we demonstrate
that, under an additional condition, our proposal still ensures control over
false discoveries. Such condition can be verified under mild circumstances for
partially observed functional data.

  * •

Third, we propose a novel method to formulate the functional graphical model
estimation as a multiple testing problem for the cross-covariance structures
between $p(p-1)/2$ pairs of functional regression errors formed by
respectively regressing each pair $X_{ij}(\cdot)$ and $X_{ik}(\cdot)$ ($1\leq
j<k\leq p$) on the remaining $p-2$ functional variables. Employing the
penalized functional regression technique to estimate functional coefficients,
we obtain $p(p-1)/2$ pairs of functional residuals, which can be nicely
integrated into a generalized version of our error-contamination framework. By
deriving relevant convergence rates to validate the extra condition and
applying our established theory, we show that such proposal combined with our
multiple test procedure achieves effective false discovery control.

Our work lies in the intersection of high-dimensional functional data analysis
and large-scale multiple testing. In addition to aforementioned functional
covariance and graphical models, various functional sparsity assumptions are
commonly imposed on the model parameter space for high-dimensional functional
data, and different regularized estimation procedures have been developed for
the respective learning tasks including functional additive regressions Fan et
al. (2015); Kong et al. (2016); Luo and Qi (2017); Chang et al. (2024),
functional linear discriminant analysis Xue et al. (2023) and sparse
functional principal components analysis Hu and Yao (2022). There is also a
wealth of literature focused on the simultaneous testing of a large collection
of hypotheses with false discovery control for the local covariance
structures. Examples include, but are not limited to, studies on correlations
Cai and Liu (2016), conditional dependence between variables or subvectors in
Gaussian graphical models Liu (2013); Xia, Cai and Cai (2018), regression
coefficients over multiple responses Xia, Cai and Li (2018), cross-covariance
matrices between subvectors in multimodal integrative analysis Xia et al.
(2020) and more complex dependence structures of multimodal imaging data Chang
et al. (2023).

The rest of the paper is organized as follows. Section 2 presents the test for
a given pair of cross-covariance functions and a multiple testing procedure
for all pairs, along with the established theoretical properties. Section 3
explores the multiple testing under a general error-contamination framework
and presents the supporting theory. Section 4 illustrates the proposed methods
with theoretical guarantees using two functional network models. We
demonstrate the superior finite-sample performance of the proposed methods
through extensive simulations in Section 5 and the analysis of two brain
imaging datasets in Section 6. All technical proofs are relegated to the
Supplementary Material.

We summarize some notation to be used throughout the paper. Denote by
$I(\cdot)$ the indicator function. Let $\overset{d}{=}$ denote the equality in
distribution. For a positive integer $q,$ we write $[q]=\\{1,\ldots,q\\}$. For
any $x,y\in\mathbb{R}$, let $x\vee y=\max(x,y)$ and $(x)_{+}=x\vee 0$. For any
set ${\cal A}$, we denote its cardinality by $|{\cal A}|.$ For two sequences
of positive numbers $\\{a_{n}\\}$ and $\\{b_{n}\\}$, we write $a_{n}\lesssim
b_{n}$ or $b_{n}\gtrsim a_{n}$ if there exists some positive constant $c$ such
that $\lim\sup_{n\to\infty}a_{n}/b_{n}\leq c$, and $a_{n}\ll b_{n}$ if
$\lim\sup_{n\rightarrow\infty}a_{n}/b_{n}=0$. We write $a_{n}\asymp b_{n}$ if
and only if $a_{n}\lesssim b_{n}$ and $b_{n}\lesssim a_{n}$ hold
simultaneously. For a vector ${\bf
b}=(b_{1},\ldots,b_{p})^{{\mathrm{\scriptscriptstyle T}}}$, we denote its
$\ell_{2}$ norm by $\|{\bf b}\|=(\sum_{j=1}^{p}|b_{j}|^{2})^{1/2}$. For a
matrix ${\bf B}=(B_{ij})_{p\times q}$, we denote its Frobenius norm by $||{\bf
B}||_{\text{F}}=\big{(}\sum_{i,j}\text{B}_{ij}^{2}\big{)}^{1/2}.$ Let ${\cal
L}_{2}({\cal U})$ be a Hilbert space of square-integrable functions on a
compact interval ${\cal U}$ equipped with the inner-product $\langle
f,g\rangle=\int_{{\cal U}}f(u)g(u)\,{\rm d}u$ for $f,g\in{\cal L}_{2}({\cal
U})$ and the induced norm $\|\cdot\|=\langle\cdot,\cdot\rangle^{1/2}$. For a
Hilbert space $\mathbb{H}\subseteq L_{2}({\cal U}),$ we denote the $p$-fold
Cartesian product by $\mathbb{H}^{p}=\mathbb{H}\times\cdots\times\mathbb{H}$.
For any ${\cal B}\in\mathbb{S}\equiv L_{2}({\cal U}\times{\cal U})$, we denote
its Hilbert–Schmidt norm by $\|{\cal B}\|_{{\cal S}}=\\{\iint{\cal
B}(u,v)^{2}\,{\rm d}u{\rm d}v\\}^{1/2}$. For any $\boldsymbol{\cal B}=({\cal
B}_{1},\ldots,{\cal B}_{p})^{{\mathrm{\scriptscriptstyle T}}}$ with each
${\cal B}_{j}\in\mathbb{S}$, we denote its functional versions of $\ell_{1}$
and $\ell_{\infty}$ norms by $\|\boldsymbol{\cal B}\|_{{\cal
S},1}=\sum_{j=1}^{p}\|{\cal B}_{j}\|_{{\cal S}}$ and $\|\boldsymbol{\cal
B}\|_{{\cal S},\max}=\max_{j\in[p]}\|{\cal B}_{j}\|_{{\cal S}},$ respectively.

##  2 Test with fully observed functional data

This section considers the setting where each trajectory $X_{ij}(\cdot)$ for
$i\in[n]$ and $j\in[p]$ is fully observed without any error. In Section 2.1,
we propose a marginal test statistic for (1) with a given pair $1\leq j\neq
k\leq p,$ and study its theoretical properties. In Section 2.2, we develop a
procedure for simultaneously testing multiple hypotheses of (1) for $1\leq
j<k\leq p$ with theoretical guarantees on false discovery control.

###  2.1 Testing a cross-covariance function

Based on $n$ i.i.d. functional observations ${\bf X}_{1}(\cdot),\dots,{\bf
X}_{n}(\cdot),$ we obtain the sample covariance function matrix
$\widehat{\boldsymbol{\Sigma}}(\cdot,\cdot)=\\{\widehat{\Sigma}_{jk}(\cdot,\cdot)\\}_{p\times
p}$ with its $(j,k)$th entry given by

$\widehat{\Sigma}_{jk}(u,v)=\frac{1}{n-1}\sum_{i=1}^{n}\\{X_{ij}(u)-\widebar{X}_{j}(u)\\}\\{X_{ik}(v)-\widebar{X}_{k}(v)\\}~{}\text{for}~{}(u,v)\in{\cal
U}^{2}~{}\text{and}~{}j,k\in[p]$ (3)

and $\widebar{X}_{j}(\cdot)=n^{-1}\sum_{i=1}^{n}X_{ij}(\cdot)$ for $j\in[p].$

We propose the following Hilbert–Schmidt-norm-based test statistic for testing
$H_{0,jk}$ with a given pair $j\neq k,$ i.e.,

$T_{n,jk}=n\|\widehat{\Sigma}_{jk}\|_{{\cal S}}^{2}\,.$ (4)

###### Remark 1

Note that the test statistic $T_{n,jk}$ measures the discrepancy between the
sample cross-covariance function $\widehat{\Sigma}_{jk}(u,v)$ and the true
function using an $L_{2}$-type norm, which tends to be small under the null
hypothesis $H_{0,jk}.$ For functional data with local spikes, one may consider
a supermum-norm-based test statistic, upon which our multiple testing
procedure in Section 2.2 remains valid under the associated technical
analysis. See more detailed discussions in Section 7. By comparison, our
$L_{2}$ type-based test performs better when a (eventually not so strong)
signal in the cross-covariance function is spread out over the entire domain
${\cal U}^{2},$ making it well-suited for accommodating the brain imaging
data. See Figures S1 and S2 of the Supplementary Material for the supporting
empirical evidence.

Denote the $j$th subject-effect function
$Z_{ij}(\cdot)=X_{ij}(\cdot)-\mu_{j}(\cdot)$ with the mean function
$\mu_{j}(\cdot)=\mathbb{E}\\{X_{ij}(\cdot)\\}$ for $j\in[p].$ Denote
$\Gamma_{jk}(u_{1},v_{1},u_{2},v_{2})=\text{Cov}\\{Z_{j}(u_{1})Z_{k}(v_{1}),Z_{j}(u_{2})Z_{k}(v_{2})\\}$
for $(u_{1},v_{1}),(u_{2},v_{2})\in{\cal U}^{2}$, which can be shown to be the
asymptotic covariance function of $\sqrt{n}\widehat{\Sigma}_{jk}(u,v)$ for
$(u,v)\in{\cal U}^{2}.$ By performing an eigenanalysis of
$\Gamma_{jk}(\cdot,\cdot,\cdot,\cdot)$, we obtain the sorted eigenvalues
$\lambda_{jk1}\geq\lambda_{jk2}\geq\cdots>0$ and the associated orthonormal
eigenfunctions $\varphi_{jk1}(\cdot,\cdot),\varphi_{jk2}(\cdot,\cdot),\dots$
satisfying

$\iint\Gamma_{jk}(u_{1},v_{1},u_{2},v_{2})\varphi_{jkr}(u_{2},v_{2})\,{\rm
d}u_{2}{\rm
d}v_{2}=\lambda_{jkr}\varphi_{jkr}(u_{1},v_{1})~{}~{}\text{for}~{}~{}r\in\mathbb{N}.$
(5)

We compute the sample estimator of $\Gamma_{jk}$ by

$\displaystyle\widehat{\Gamma}_{jk}(u_{1},v_{1},u_{2},v_{2})=\frac{1}{n}\sum_{i=1}^{n}\\{\widetilde{X}_{ij}(u_{1})\widetilde{X}_{ik}(v_{1})-\widehat{\Sigma}_{jk}(u_{1},v_{1})\\}\\{\widetilde{X}_{ij}(u_{2})\widetilde{X}_{ik}(v_{2})-\widehat{\Sigma}_{jk}(u_{2},v_{2})\\}$

with $\widetilde{X}_{ij}(\cdot)=X_{ij}(\cdot)-\widebar{X}_{j}(\cdot)$ and then
carry out an eigenanlysis of $\widehat{\Gamma}_{jk}$ in the form of (5), which
leads to the estimated eigenvalue/eigenfunction pairs
$\\{\hat{\lambda}_{jkr},\hat{\varphi}_{jkr}(\cdot,\cdot)\\}_{r\in\mathbb{N}}.$

Before investigating the theoretical properties, we need some regularity
conditions.

###### Condition 1

The covariance functions satisfy $\int\Sigma_{jj}(u,u)\,{\rm d}u=O(1)$ and
$\sup_{u\in{\cal U}}\Sigma_{jj}(u,u)=O(1)$ for any $j\in[p]$.

###### Condition 2

The subject-effect functions satisfy ${\mathbb{E}}(\|Z_{1j}\|^{4})=O(1)$ for
any $j\in[p]$ and $\sup_{(u,v)\in{\cal
U}^{2}}{\mathbb{E}}\\{Z_{1j}^{2}(u)Z_{1k}^{2}(v)\\}=O(1)$ for any $j,k\in[p].$

Conditions 1 and 2 are standard in functional data analysis literature, see,
e.g., Chapter 10 of Zhang (2013). We are now ready to present the theorem
regarding the limiting null distribution and the asymptotic power.

###### Theorem 1

Suppose that Conditions 1 and 2 hold, then the following two assertions hold.
(i) Under the null hypothesis $H_{0,jk}$, we have that
$T_{n,jk}\overset{d}{\rightarrow}T_{0,jk}$ as $n\rightarrow\infty,$ where
$T_{0,jk}\overset{d}{=}\sum_{r=1}^{\infty}\hat{\lambda}_{jkr}A_{r}$ with
$\\{A_{r}\\}_{r\in\mathbb{N}}$ being i.i.d. $\chi^{2}_{1}$ random variables.
(ii) Under the alternative hypothesis $H_{1,jk}$ with $\|\Sigma_{jk}\|_{{\cal
S}}>0$, we have that ${\mathbb{P}}\\{T_{n,jk}\geq
T_{0,jk}(\alpha)\\}\rightarrow 1$ as $n\rightarrow\infty$, where
$T_{0,jk}(\alpha)$ is the $100\alpha\%$ upper percentile of $T_{0,jk}$.

###### Remark 2

Theorem 1 indicates that $T_{n,jk}$ is asymptotically an infinite
$\chi^{2}$-type mixture under the null hypothesis $H_{0,jk}.$ Its distribution
can be approximated by a noncentral chi-squared distribution with parameters
determined by $\sum_{r=1}^{\infty}\hat{\lambda}_{jkr},$
$\sum_{r=1}^{\infty}\hat{\lambda}_{jkr}^{2},$
$\sum_{r=1}^{\infty}\hat{\lambda}_{jkr}^{3}$ and
$\sum_{r=1}^{\infty}\hat{\lambda}_{jkr}^{4}$, see, e.g., Liu et al. (2009). We
will use this approximation in our numerical studies. It is noteworthy that,
instead of calculating each individual eigenvalue $\hat{\lambda}_{jkr}$
explicitly, it is more convenient to compute the aforementioned infinite sums.
For example, $\sum_{r=1}^{\infty}\hat{\lambda}_{jkr}$ and
$\sum_{r=1}^{\infty}\hat{\lambda}_{jkr}^{2}$ can be obtained through the
discretization of functional data and the respective simple Riemann sum
approximations to the integrals $\iint\widehat{\Gamma}_{jk}(u,v,u,v)\,{\rm
d}u{\rm d}v$ and
$\iiiint\widehat{\Gamma}_{jk}(u_{1},v_{1},u_{2},v_{2})^{2}\,{\rm d}u_{1}{\rm
d}v_{1}{\rm d}u_{2}{\rm d}v_{2}.$ In this sense, our procedure is fully
functional without relying on any dimension reduction techniques that may
result in information loss.

###  2.2 Multiple testing of cross-covariance functions

To develop a general multiple testing procedure, we let ${\cal
H}=\\{(j,k):j,k\in[p],j<k\\}$ and reformulate (1) as simultaneously testing
$Q=p(p-1)/2$ hypotheses:

$\displaystyle H_{0}^{(q)}:\Sigma_{j_{q}k_{q}}(u,v)=0\mbox{ for any
}(u,v)\in{\cal U}^{2}\mbox{ and }(j_{q},k_{q})\in{\cal
H}~{}~{}\mbox{vs}~{}~{}H_{1}^{(q)}:H_{0}^{(q)}\mbox{ is not true}\,,$ (6)

for $q\in[Q]$. The test statistic of the $q$th hypothesis test $H_{0}^{(q)}$
is constructed as

$T_{n}^{(q)}:=T_{n,j_{q}k_{q}}=n\|\widehat{\Sigma}_{j_{q}k_{q}}\|_{{\cal
S}}^{2}~{},~{}~{}q\in[Q].$ (7)

Then we reject $H_{0}^{(q)}$ when $T_{n}^{(q)}$ takes some large values. Let
${\cal H}_{0}=\\{q\in[Q]:H_{0}^{(q)}~{}\textrm{is true}\\}$ be the set of all
true null hypotheses with $Q_{0}=|{\cal H}_{0}|.$ By Theorem 1(i), we obtain
that

${\mathbb{P}}\\{T_{n}^{(q)}\leq x\\}-{\mathbb{P}}\\{T_{0}^{(q)}\leq
x\\}=o(1)\text{ for any }q\in{\cal H}_{0},$ (8)

where
$T_{0}^{(q)}\overset{d}{=}\sum_{r=1}^{\infty}\hat{\lambda}_{r}^{(q)}A_{r}$
with $\hat{\lambda}_{r}^{(q)}:=\hat{\lambda}_{j_{q}k_{q}r}$ for notational
simplicity. Then the p-value of $H_{0}^{(q)}$ is given by ${\rm
pv}^{(q)}={\mathbb{P}}\\{T_{0}^{(q)}\geq T_{n}^{(q)}\\}$, and its normal
quantile transformation takes the form $V_{q}=\Phi^{-1}\\{1-{\rm
pv}^{(q)}\\},$ where $\Phi(\cdot)$ is the cumulative distribution function of
the standard normal distribution. For any threshold level $t\in{\mathbb{R}}$
such that $H_{0}^{(q)}$ is rejected if $V_{q}\geq t,$ denote the total number
of false positives by $R_{0}(t)=\sum_{q\in\mathcal{H}_{0}}I(V_{q}\geq t)$ and
the total number of rejections by $R(t)=\sum_{q=1}^{Q}I(V_{q}\geq t)$. Then,
the false discovery proportion (FDP) and the false discovery rate (FDR) are
respectively defined as

$\displaystyle{\rm FDP}(t)=\frac{R_{0}(t)}{R(t)\vee
1}~{}~{}~{}\textrm{and}~{}~{}~{}{\rm FDR}(t)=\mathbb{E}\\{{\rm FDP}(t)\\}\,.$

Given a prespecified level $\alpha\in(0,1)$, the main target is to find the
smallest $\hat{t}$ such that ${\rm FDP}(\hat{t})\leq\alpha.$ To achieve this,
we first consider ${\rm FDP}(t),$ the numerator of which needs to be estimated
due to the unknown set ${\cal H}_{0}.$ Using (8), we can approximate
${\mathbb{P}}(V_{q}\geq t)$ as $1-\Phi(t)$ for any $q\in\mathcal{H}_{0}.$ An
ideal estimate for ${\rm FDR}(t)$ can be obtained by
$[Q_{0}\\{1-\Phi(t)\\}]/\\{1\vee\sum_{q=1}^{Q}I(V_{q}\geq t)\\}$, where
$Q_{0}$ is unknown in practice. Thus we propose to estimate ${\rm FDP}(t)$ in
a more conservative way by replacing $Q_{0}$ with $Q:$

$\widehat{\rm FDP}(t)=\frac{Q\\{1-\Phi(t)\\}}{1\vee\sum_{q=1}^{Q}I(V_{q}\geq
t)}.$ (9)

We summarize the proposed multiple testing procedure in Algorithm 1.

Algorithm 1 Multiple testing procedure

Step 1. Calculate the test statistics $T_{n}^{(q)}$ in (7) and its
transformation $V_{q}$ for all $q\in[Q]$.

Step 2. Estimate the FDP by $\widehat{\rm FDP}(t)$ in (9).

Step 3. For a given $0<\alpha<1$, choose
$\displaystyle\hat{t}=\inf\big{\\{}0\leq t\leq(2\log Q-2\log\log
Q)^{1/2}:\,\widehat{\rm FDP}(t)\leq\alpha\big{\\}}\,.$

If $\hat{t}$ does not exist, set $\hat{t}=(2\log Q)^{1/2}$.

Step 4. For $q\in[Q]$, reject $H_{0}^{(q)}$ with $V_{q}\geq\hat{t}.$

Before presenting the theoretical results, we impose some regularity
conditions.

###### Condition 3

For each $j\in[p]$ and $i\in[n],$ $Z_{ij}(\cdot)$ is a mean-zero sub-Gaussian
process, i.e., there exists some constant $\check{c}_{j}>0$ such that for all
$x\in\mathbb{H},$ $\mathbb{E}\big{\\{}\exp(\langle
x,Z_{ij}\rangle)\big{\\}}\leq\exp\big{\\{}2^{-1}\check{c}_{j}^{2}\iint
x(u)\Sigma_{jj}(u,v)x(v)\rm{d}u{\rm d}v\big{\\}}.$

To accommodate the formulation in this section, for any $q\in\mathcal{H}_{0},$
the Karhunen-Loève expansion allows us to represent
$z_{iq}(u,v):=Z_{ij_{q}}(u)Z_{ik_{q}}(v)=\sum_{r=1}^{\infty}a_{ir}^{(q)}\varphi_{r}^{(q)}(u,v),$
where $\lambda_{r}^{(q)}:=\lambda_{j_{q}k_{q}r}$,
$\varphi_{r}^{(q)}(u,v):=\varphi_{j_{q}k_{q}r}(u,v)$, and $a_{ir}^{(q)}=\iint
z_{iq}(u,v)\varphi_{r}^{(q)}(u,v)\,{\rm d}u{\rm d}v.$

###### Condition 4

For any $q\in\mathcal{H}_{0},$ $\lambda_{1}^{(q)}>\lambda_{2}^{(q)}>\cdots>0,$
and there exist some constant $\eta_{1}>1$ such that
$\lambda_{r}^{(q)}\lesssim r^{-\eta_{1}}$ for $r\in\mathbb{N}_{+}.$

Condition 3 is an infinite-dimensional analog of the sub-Gaussian condition
within Hilbert space, which is satisfied by Gaussian processes. This condition
ensures that $z_{iq}(\cdot,\cdot)$ exhibits sub-exponential behavior, as
reflected in the associated basis coefficients
$\\{a_{ir}^{(q)}\\}_{r\in{\mathbb{N}}}.$ See Lemma A1 of the Supplementary
Material. These sub-exponentialities facilitate the Gaussian approximations of
the proposed test statistics. Condition 4 allows the polynomial decay of the
upper bounds on eigenvalues with parameter $\eta_{1}$ determining the decay
rate, and can be used to control the truncation error rate
$\sum_{r=M+1}^{\infty}\lambda_{r}^{(q)}\lesssim M^{-\eta_{1}+1}$ for
$M\rightarrow\infty.$ Larger values of $\eta_{1}$ correspond to faster decay
rate and smaller truncation errors. It can be easily verified that Condition 4
is fulfilled when the eigenvalues for the covariance function of
$Z_{ij_{q}}(\cdot)$ decay polynomially Kong et al. (2016); Qiao et al. (2019).

To investigate the theoretical properties of the proposed multiple testing
procedure, we need to give some regularity conditions on the dependency among
the test statistics $\\{T_{n}^{(q)}\\}_{q=1}^{Q}$. Notice that
$\zeta_{q}=\Phi^{-1}[F_{q}\\{T_{n}^{(q)}\\}]$ follows the standard normal
distribution, where $F_{q}(\cdot)$ denotes the distribution of $T_{n}^{(q)}$.
We use the correlation between two normal random variables $\zeta_{q}$ and
$\zeta_{q^{\prime}}$ to characterize the dependence between test statistics
$T_{n}^{(q)}$ and $T_{n}^{(q^{\prime})}$. For some constant $\gamma>0$, define
the set

$\mathcal{S}_{q}(\gamma)=\\{q^{\prime}\in[Q]:q^{\prime}\neq q,|{\rm
Corr}(\zeta_{q},\,\zeta_{q^{\prime}})|\geq(\log Q)^{-2-\gamma}\\}\,.$

For a given $q\in[Q]$, the remaining $Q-1$ test statistics
$\\{T_{n}^{(q^{\prime})}\\}_{q^{\prime}\in[Q]/\\{q\\}}$ can be categorized
into two scenarios: (i) for $q^{\prime}\in\mathcal{S}_{q}(\gamma)$,
$T_{n}^{(q^{\prime})}$ has a relatively strong dependence with $T_{n}^{(q)},$
and (ii) for $q^{\prime}\notin\mathcal{S}_{q}(\gamma)$, $T_{n}^{(q^{\prime})}$
has a quite weak dependence with $T_{n}^{(q)}$. We establish the theoretical
guarantee under these two scenarios separately using different technical
tools.

###### Theorem 2

Suppose that Conditions 1–4 hold, $\max_{1\leq q\neq q^{\prime}\leq Q}|{\rm
Corr}(\zeta_{q},\zeta_{q^{\prime}})|\leq c_{\zeta}$ for some constant
$c_{\zeta}\in(0,1),$ and $\max_{q\in[Q]}|\mathcal{S}_{q}(\gamma)|=o(p^{2\nu})$
for some constants $\gamma>0$ and $0<\nu<(1-c_{\zeta})/(1+c_{\zeta})$. If
$p\lesssim n^{\kappa}$ for some constant $\kappa>0$, then
$\limsup_{n,Q\rightarrow\infty}{\rm FDR}(\hat{t})\leq\alpha Q_{0}/Q$ and
$\lim_{n,Q\rightarrow\infty}{\mathbb{P}}\\{{\rm FDP}(\hat{t})\leq\alpha
Q_{0}/Q+\varkappa\\}=1$ for any $\varkappa>0$.

Theorem 2 shows that both FDP and FDR can be controlled below the level
$\alpha Q_{0}/Q$ asymptotically. The conditions on ${\rm
Corr}(\zeta_{q},\zeta_{q^{\prime}})$ and $\mathcal{S}_{q}(\gamma)$ are imposed
to bound the variance of $\sum_{q\in{\cal H}_{0}}I(V_{q}\geq t)$ in ${\rm
FDP}(t),$ see similar conditions in Liu (2013) and Chang et al. (2023).
Specifically, the condition $\max_{1\leq q\neq q^{\prime}\leq Q}|{\rm
Corr}(\zeta_{q},\zeta_{q^{\prime}})|\leq c_{\zeta}$ places a constraint on the
strength of dependence between different $\zeta_{q}$ and $\zeta_{q^{\prime}}$.
The condition $\max_{q\in[Q]}|\mathcal{S}_{q}(\gamma)|=o(p^{2\nu})$ controls
the number of pairs of $(\zeta_{q},\zeta_{q^{\prime}})$ with relatively strong
dependencies.

##  3 Test with error-contaminated functional data

In practical applications, ${\bf X}_{i}(\cdot)$ for $i\in[n]$ are rarely fully
observed. Instead, we estimate them using certain methods, resulting in
estimated curves $\widehat{{\bf
X}}_{i}(\cdot)=\\{\widehat{X}_{i1}(\cdot),\ldots,\widehat{X}_{ip}(\cdot)\\}^{{\mathrm{\scriptscriptstyle
T}}},$ denoted as error-contaminated functional data satisfying (2). This
leads to functional errors ${\mathbf{e}}_{i}(\cdot)=\widehat{{\bf
X}}_{i}(\cdot)-{\bf X}_{i}(\cdot)$ with
${\mathbf{e}}_{i}(\cdot)=\\{e_{i1}(u),\ldots,e_{ip}(\cdot)\\}^{{\mathrm{\scriptscriptstyle
T}}}.$ See below for two examples that are related to applications in Sections
4.1 and 4.2, respectively.

  * •

When $\\{{\bf X}_{i}(\cdot)\\}$ are partially observed with errors, we can
apply nonparametric smoothing methods to the observed data, thereby obtaining
functional estimates $\widehat{{\bf X}}_{i}(\cdot).$

  * •

When $\\{{\bf X}_{i}(\cdot)\\}$ represent the functional regression errors in
function-on-function linear regressions, we can estimate functional
coefficients and thus obtain functional fitted values, resulting in functional
residuals $\widehat{\bf X}_{i}(\cdot).$

In the above examples, it is apparent that each ${\mathbf{e}}_{i}(\cdot)$ may
not be independent, and the properties of ${\mathbf{e}}_{i}(\cdot)$ need to be
investigated for specific models. This section is devoted to the corresponding
tests for the cross-covariance functions $\\{\Sigma_{jk}(u,v)\\}_{1\leq j\neq
k\leq p}$ based on the estimated curves $\widehat{{\bf X}}_{i}(\cdot)$ under
some general conditions.

Analogous to Section 2.1, we first consider the marginal testing problem (1)
for error-contaminated functional data. To this end, we replace each ${\bf
X}_{i}(\cdot)$ by its error-contaminated version $\widehat{\bf X}_{i}(\cdot)$
in (3)–(4), and propose a new test statistic for $H_{0,jk}$ for a given $1\leq
j\neq k\leq p:$

$\displaystyle\widetilde{T}_{n,jk}=n\|\widetilde{\Sigma}_{jk}\|_{{\cal
S}}^{2}\,,$ (10)

where
$\widetilde{\Sigma}_{jk}(u,v)=(n-1)^{-1}\sum_{i=1}^{n}\\{\widehat{X}_{ij}(u)-{\widebar{\widehat{X}}}_{j}(u)\\}\\{\widehat{X}_{ik}(v)-\widebar{\widehat{X}}_{k}(v)\\}$
and $\widebar{\widehat{X}}_{j}(u)=n^{-1}\sum_{i=1}^{n}\widehat{X}_{ij}(u).$ To
control the difference between $\widetilde{T}_{n,jk}$ and $T_{n,jk}$ specified
in (4), we need the cross-covariance functions
$\boldsymbol{\Sigma}^{Xe}(u,v)={\text{Cov}}\\{{\bf
X}(u),{\mathbf{e}}(v)\\}=\\{\Sigma_{jk}^{Xe}(u,v)\\}_{j,k\in[p]}$ and
$\boldsymbol{\Sigma}^{ee}(u,v)={\text{Cov}}\\{{\mathbf{e}}(u),{\mathbf{e}}(v)\\}=\\{\Sigma_{jk}^{ee}(u,v)\\}_{j,k\in[p]},$
whose sample versions are

$\displaystyle\widehat{\boldsymbol{\Sigma}}^{Xe}(u,v)=\\{\widehat{\Sigma}_{jk}^{Xe}(u,v)\\}_{j,k\in[p]}=$
$\displaystyle~{}\frac{1}{n-1}\sum_{i=1}^{n}\\{{\bf X}(u)-\widebar{\bf
X}(u)\\}\\{{\mathbf{e}}(v)-\bar{\mathbf{e}}(v)\\}^{{\mathrm{\scriptscriptstyle
T}}}\,,$
$\displaystyle\widehat{\boldsymbol{\Sigma}}^{ee}(u,v)=\\{\widehat{\Sigma}_{jk}^{ee}(u,v)\\}_{j,k\in[p]}=$
$\displaystyle~{}\frac{1}{n-1}\sum_{i=1}^{n}\\{{\mathbf{e}}(u)-\bar{\mathbf{e}}(u)\\}\\{{\mathbf{e}}(v)-\bar{\mathbf{e}}(v)\\}^{{\mathrm{\scriptscriptstyle
T}}}\,,$

respectively. Before presenting the asymptotic properties of the new test
statistic $\widetilde{T}_{n,jk}$ under the null and alternative hypotheses, we
need the following condition, which will be verified by deriving the
corresponding convergence rates for two applications in Section 4.

###### Condition 5

There exist some constants $a_{1},a_{2}>1/2$ such that, for any $1\leq j\neq
k\leq p,$

$\displaystyle\|\widehat{\Sigma}^{Xe}_{jk}\|_{{\cal S}}=O_{\rm
p}(n^{-a_{1}})~{}~{}{\rm and}~{}~{}\|\widehat{\Sigma}^{ee}_{jk}\|_{{\cal
S}}=O_{\rm p}(n^{-a_{2}}).$ (11)

###### Theorem 3

Suppose that Conditions 1, 2, 5 hold, then the following two assertions hold.
(i) Under the null hypothesis $H_{0,jk}$, we have that
$\widetilde{T}_{n,jk}\overset{d}{\rightarrow}T_{0,jk}$ as
$n\rightarrow\infty,$ where
$T_{0,jk}\overset{d}{=}\sum_{r=1}^{\infty}\hat{\lambda}_{jkr}A_{r}$ with
$A_{r}\overset{i.i.d.}{\sim}\chi^{2}_{1}$.
(ii) Under the alternative hypothesis $H_{1,jk}$ with $\|\Sigma_{jk}\|_{{\cal
S}}>0$, we have ${\mathbb{P}}\\{\widetilde{T}_{n,jk}\geq
T_{0,jk}(\alpha)\\}\rightarrow 1$ as $n\rightarrow\infty.$

Analogous to Section 2.2, we next consider the multiple testing problem (6)
under the functional error-contamination setting. For each $q\in[Q]$, we
propose the test statistic for the $q$th hypothesis testing as follows:

$\displaystyle\widetilde{T}_{n}^{(q)}:=\widetilde{T}_{n,j_{q}k_{q}}=n\|\widetilde{\Sigma}_{j_{q}k_{q}}\|_{{\cal
S}}^{2}\,.$

By Theorem 3(i), we have ${\mathbb{P}}\\{\widetilde{T}_{n}^{(q)}\leq
x\\}-{\mathbb{P}}\\{T_{0}^{(q)}\leq x\\}=o(1)$ for any $q\in{\cal H}_{0}$. Let
the resulting p-value of $H_{0}^{(q)}$ be $\widetilde{\rm
pv}^{(q)}={\mathbb{P}}\\{T_{0}^{(q)}\geq\widetilde{T}_{n}^{(q)}\\}$, its
normal quantile transformation $\widetilde{V}_{q}=\Phi^{-1}\\{1-\widetilde{\rm
pv}^{(q)}\\},$ and ${t}$ be the threshold level such that $H_{0}^{(q)}$ is
rejected if $\widetilde{V}_{q}\geq{t}$. Given $\alpha\in(0,1)$, the ideal
choice of ${t}$ is the smallest ${t}$ such that we would reject as many true
positives as possible while controlling the false discovery proportion at the
prescribed level $\alpha$. To achieve this, similar to (9) and Step 3 in
Algorithm 1, we let $\widetilde{{\rm
FDP}}(t)=Q\\{1-\Phi(t)\\}/\\{1\vee\sum_{q=1}^{Q}I(\widetilde{V}_{q}\geq t)\\}$
and

$\tilde{t}=\inf\\{0\leq t\leq(2\log Q-2\log\log Q)^{1/2}:\,\widetilde{\rm
FDP}(t)\leq\alpha\\}.$

We then propose a new algorithm to implement the multiple testing procedure by
replacing $T_{n}^{(q)},$ $\widehat{{\rm FDP}}(t),$ $\hat{t}$ and $V_{q}$ in
Algorithm 1 with $\widetilde{T}_{n}^{(q)},$ $\widetilde{{\rm FDP}}(t),$
$\tilde{t}$ and $\widetilde{V}_{q},$ respectively.

To impose additional regularity conditions similar to those in Theorem 2, we
let $\widetilde{F}_{q}(\cdot)$ be the distribution of
$\widetilde{T}_{n}^{(q)},$ and write
$\widetilde{\zeta}_{q}=\Phi^{-1}[\widetilde{F}_{q}\\{\widetilde{T}_{n}^{(q)}\\}]$,
which follows the standard normal distribution. Accordingly, we define

$\displaystyle\widetilde{\mathcal{S}}_{q}(\tilde{\gamma})=\\{{q^{\prime}}\in[Q]:{q^{\prime}}\neq
q,|{\rm
Corr}(\widetilde{\zeta}_{q},\,\widetilde{\zeta}_{{q^{\prime}}})|\geq(\log
Q)^{-2-\tilde{\gamma}}\\}~{}\text{ for
}\tilde{\gamma}>0~{}\text{and}~{}q\in[Q].$

Similar to Theorem 2 established for fully observed functional data, the next
theorem indicates that our proposed multiple testing procedure for error-
contaminated functional data is guaranteed to control both FDR and FDP below
the nominal level asymptotically.

###### Theorem 4

Suppose that Conditions 1–5 hold, $\max_{q,q^{\prime}\in[Q]}|{\rm
Corr}(\tilde{\zeta}_{q},\tilde{\zeta}_{q^{\prime}})|\leq c_{\tilde{\zeta}}$
for some constant $c_{\tilde{\zeta}}\in(0,1)$, and
$\max_{q\in[Q]}|\widetilde{\mathcal{S}}_{q}(\tilde{\gamma})|=o(p^{2\tilde{\nu}})$
for some constants $\tilde{\gamma}>0$ and
$0<\tilde{\nu}<(1-c_{\tilde{\zeta}})/(1+c_{\tilde{\zeta}})$. If $p\lesssim
n^{\tilde{\kappa}}$ for some constant $\tilde{\kappa}>0$, then
$\limsup_{n,Q\rightarrow\infty}{\rm FDR}({\tilde{t}})\leq\alpha Q_{0}/Q$ and
$\lim_{n,Q\rightarrow\infty}{\mathbb{P}}\\{{\rm FDP}({\tilde{t}})\leq\alpha
Q_{0}/Q+\widetilde{\varkappa}\\}=1$ for any $\widetilde{\varkappa}>0$.

###### Remark 3

Note that our proposed multiple testing procedure considers $Q=p(p-1)/2$
hypotheses of cross-covariance functions induced from $p(p-1)/2$ pairs of
$\\{\widehat{X}_{ij}(\cdot),\widehat{X}_{ik}(\cdot)\\}$ based on $p$ estimated
curves $\widehat{X}_{i1}(\cdot),\dots,\widehat{X}_{ip}(\cdot).$ The
theoretical guarantee of this framework can be easily extended to a more
generalized multiple testing setting, where multiple hypotheses are induced
from $Q=p(p-1)/2$ pairs based on $p(p-1)/2$ instead of $p$ estimated curves.
For example, in Section 4.2, we transform the multiple testing problem of
conditional dependence structures to that of the error cross-covariance
patterns in (15), which is induced from $p(p-1)/2$ pairs of functional
residuals as specified in (16), i.e.
$\\{\hat{\varepsilon}_{i,jk}(\cdot),\hat{\varepsilon}_{i,kj}(\cdot)\\}$ for
$1\leq j<k\leq p.$ Not only does our proposal for the generalized multiple
testing problem remain unchanged, but also, following similar arguments as
those in the proofs of Theorems 2 and 4, it can still achieve effective
control of false discoveries asymptotically. Further details are available in
the application to the functional graphical model in Section 4.2.

##  4 Applications

###  4.1 Partially observed functional data

For any $i\in[n]$ and $j\in[p],$ each $X_{ij}(\cdot)$ is not directly
observable in practice. Instead, it is observed with errors at $T_{ij}$ random
time points, $U_{ij1},\dots,U_{ijT_{ij}}\in{\cal U}.$ Let $W_{ijt}$ be the
observed value of $X_{ij}(U_{ijt})$ satisfying

$W_{ijt}=X_{ij}(U_{ijt})+\varsigma_{ijt}\,,$ (12)

where the errors $\varsigma_{ijt}$’s, independent of $X_{ij}(\cdot)$’s, are
i.i.d. with $\mathbb{E}(\varsigma_{ijt})=0$ and ${\rm
Var}(\varsigma_{ijt})=\sigma_{j}^{2}<\infty.$ The sampling frequency $T_{ij}$
plays a vital role in the estimation as it may affect the choice of the
estimation procedure. For densely sampled functional data with $T_{ij}$’s
larger than some order of $n,$ the conventional approach is to implement
nonparametric smoothing to the observations from each subject, thus
reconstructing each individual curve before subsequent analysis Zhang and Chen
(2007). For sparsely sampled functional data with bounded $T_{ij}$’s, pre-
smoothing is no longer viable; it is standard practice to consider pooling
data from all subjects to borrow strength across all observations Zhang and
Wang (2016).

Since brain signals in neuroimaging data of our interest are typically
recorded at a dense set of points, we make use of the pre-smoothing to
facilitate the presentation of methodology and theory. In what follows, denote
$K_{h}(\cdot)=h^{-1}K(\cdot/h)$ for a univariate kernel with bandwidth $h>0.$
We first apply the frequently-adopted local linear smoothers to the observed
data $\\{(U_{ijt},W_{ijt}):t\in[T_{ij}]\\}$ for each $i\in[n]$ and $j\in[p],$
and thus can obtain the estimated curves $\widehat{X}_{ij}(u)=\hat{g}_{0ij},$
where

$(\hat{g}_{0ij},\hat{g}_{1ij})=\arg\min_{g_{0ij},g_{1ij}}\,\sum_{t=1}^{T_{ij}}\big{\\{}W_{ijt}-g_{0ij}-g_{1ij}(U_{ijt}-u)\big{\\}}^{2}K_{h_{ij}}(U_{ijt}-u).$

For any $j\in[p],$ individual functions across $i\in[n]$ often admit similar
smoothness properties and sometimes similar patterns, it is thus reasonable to
use a common bandwidth $h_{j}$ for all of them. We then employ the multiple
testing procedure in Section 3 to the estimated curves
$\\{\widehat{X}_{ij}(\cdot)\\},$ which leads to
$e_{ij}(\cdot)=\widehat{X}_{ij}(\cdot)-X_{ij}(\cdot).$

To provide the theoretical guarantee for such procedure based on Theorem 4, we
need to verify Condition 5, i.e. specifying the convergence rates of
$\|\widehat{\Sigma}_{jk}^{Xe}\|_{{\cal S}}$ and
$\|\widehat{\Sigma}_{jk}^{ee}\|_{{\cal S}}$ for any $1\leq j\neq k\leq p.$ To
this end, we impose the following regularity conditions.

###### Condition 6

For any $j\in[p],$ the average sampling frequency
$\widebar{T}_{j}=(n^{-1}\sum_{i=1}^{n}T_{ij}^{-1})^{-1}\rightarrow\infty,$
$h_{j}\rightarrow 0$ and $h_{j}\widebar{T}_{j}\rightarrow\infty$ as
$n\rightarrow\infty.$

###### Condition 7

(i) Let $\\{U_{ijt};i\in[n],j\in[p],t\in[T_{ij}]\\}$ be i.i.d. copies of a
random variable $U$ defined on ${\cal U}$; (ii) for any given $u$ within the
interior of ${\cal U}$, the first derivative of the density function
$f_{U}(u)$ of $U$ exists and is continuous over ${\cal U};$ (iii)
$X_{ij}(\cdot)$, $U_{ijt}$ and $\varsigma_{ij}$ are independent for each
$i,j,t.$

###### Condition 8

The kernel $K(\cdot)$ is bounded symmetric probability density function on
$[-1,1]$ with $\int u^{2}K(u)\,{\rm d}u<\infty$ and $\int K^{2}(u)\,{\rm
d}u<\infty$.

###### Condition 9

Let $u$ and $v$ be any two interior time points of ${\cal U}$. The individual
functions $X_{ij}(u)$ and their mean function $\mu_{j}(u)$ have up to 2-times
continuous derivatives for any $i\in[n]$ and $j\in[p]$. The covariance
function $\Sigma_{jk}(u,v)$ also has up to 2-times continuous derivatives for
any $j,k\in[p]$.

Conditions 6–9 are standard for the nonparametric estimation in functional
data analysis, see Zhang and Chen (2007). The following proposition indicates
that Condition 5 holds for partially observed functional data.

######  Proposition 1 (verify Condition 5)

Under Conditions 6–9, it holds that $\|\widehat{\Sigma}_{jk}^{Xe}\|_{{\cal
S}}=O_{\rm p}(h_{k}^{2})$ and $\|\widehat{\Sigma}_{jk}^{ee}\|_{{\cal
S}}=O_{\rm
p}\big{\\{}h_{j}^{2}h_{k}^{2}+(h_{j}\widebar{T}_{j})^{-1/2}(h_{k}\widebar{T}_{k})^{-1/2}\big{\\}}$.

###### Remark 4

To facilitate further discussion, we consider a simplified common scenario
where there exist two constants $\eta_{2}>\eta_{3}>0$ such that
$\widebar{T}_{j}\asymp n^{\eta_{2}}$ and $h_{j}\asymp n^{-\eta_{3}}$ for any
$j\in[p]$. Some calculations show that Condition 5 is satisfied by setting
$a_{1}=2\eta_{3}$ and $a_{2}=\min(4\eta_{3},\eta_{2}-\eta_{3}),$ provided that
$\eta_{3}>1/4$ and $\eta_{2}>\eta_{3}+1/2.$ With verifiable Condition 5 and
other conditions required in Theorem 4, an application of Theorem 4 yields
that our multiple testing procedure for partially observed functional data
ensures the effective control of both FDR and FDP asymptotically.

###  4.2 Functional graphical model

Suppose that we observe $p$-vector of functional data ${\bf
Y}_{i}(\cdot)=\\{Y_{i1}(\cdot),\ldots,Y_{ip}(\cdot)\\}^{{\mathrm{\scriptscriptstyle
T}}}$ for $i\in[n]$ defined on ${\cal U}$. Assume that ${\bf Y}_{i}(\cdot)$
follows a mean zero multivariate Gaussian process with covariance function
matrix $\boldsymbol{\Xi}(u,v)=\text{Cov}\\{{\bf Y}_{i}(u),{\bf Y}_{i}(v)\\}$
for $(u,v)\in{\cal U}^{2}.$ Denote ${\bf Y}_{i,-j,-k}(\cdot)$ be the subvector
of ${\bf Y}_{i}(\cdot)$ with $Y_{ij}(\cdot)$ and $Y_{ik}(\cdot)$ being
removed. In this section, we formulate the estimation of functional graphical
model depicting the conditional dependence structure among $p$ components of
${\bf Y}_{i}(\cdot)$ as simultaneous testing of $Q=p(p-1)/2$ hypotheses,

$\displaystyle\widetilde{H}_{0,jk}:Y_{ij}(\cdot)\perp\\!\\!\\!\perp
Y_{ik}(\cdot)\,|\,{\bf
Y}_{i,-j,-k}(\cdot)~{}~{}~{}\mbox{vs}~{}~{}~{}\widetilde{H}_{1,jk}:Y_{ij}(\cdot)\not\perp\\!\\!\\!\perp
Y_{ik}(\cdot)\,|\,{\bf Y}_{i,-j,-k}(\cdot)$ (13)

for $1\leq j<k\leq p$, and we identify a significant edge between nodes $j$
and $k$ if and only if $\widetilde{H}_{0,jk}$ is rejected.

For each $1\leq j<k\leq p$, consider regressing each pair $Y_{ij}(\cdot)$ and
$Y_{ik}(\cdot)$ on the remaining $p-2$ functional variables, i.e.,

$\displaystyle Y_{ij}(u)=\sum_{\ell\neq j,k}\int_{{\cal
U}}Y_{i\ell}(v)\beta_{jk,\ell}(u,v)\,{\rm
d}v+\varepsilon_{i,jk}(u)\,,~{}~{}u\in{\cal U}\,,$ (14) $\displaystyle
Y_{ik}(u)=\sum_{\ell\neq j,k}\int_{{\cal
U}}Y_{i\ell}(v)\beta_{kj,\ell}(u,v)\,{\rm
d}v+\varepsilon_{i,kj}(u)\,,~{}~{}u\in{\cal U}\,,$

where $\varepsilon_{i,jk}(\cdot)$ and $\varepsilon_{i,kj}(\cdot)$ represent
the functional regression errors, independent of ${\bf Y}_{i,-j,-k}(\cdot)$.
Under such functional regression setup, we have
$\text{Cov}\big{\\{}Y_{ij}(u),Y_{ik}(v)\big{|}{\bf
Y}_{i,-j-k}(\cdot)\big{\\}}$

$\displaystyle=$
$\displaystyle{\mathbb{E}}\Big{(}\big{[}Y_{ij}(u)-{\mathbb{E}}\\{Y_{ij}(u)|{\bf
Y}_{i,-j-k}(\cdot)\\}\big{]}\big{[}Y_{ik}(v)-{\mathbb{E}}\\{Y_{ik}(v)|{\bf
Y}_{i,-j-k}(\cdot)\\}\big{]}\Big{|}{\bf Y}_{i,-j-k}(\cdot)\Big{)}$
$\displaystyle=$
$\displaystyle{\mathbb{E}}\big{\\{}\varepsilon_{i,jk}(u)\varepsilon_{i,kj}(v)|{\bf
Y}_{i,-j,-k}(\cdot)\big{\\}},$

which possesses the property that
$\text{Cov}\big{\\{}Y_{ij}(u),Y_{ik}(v)\big{|}{\bf
Y}_{i,-j-k}(\cdot)\big{\\}}=0$ for any $(u,v)\in{\cal U}^{2},$ and hence
$\Sigma_{jk}^{\varepsilon}(u,v)={\text{Cov}}\\{\varepsilon_{i,jk}(u),\varepsilon_{i,kj}(v)\\}=0$
for any $(u,v)\in{\cal U}^{2}$ whenever $Y_{ij}(\cdot)\perp\\!\\!\\!\perp
Y_{ik}(\cdot)\,|\,{\bf Y}_{i,-j,-k}(\cdot).$ Thus, the hypotheses of
conditional independence in (13) reduces to

$\displaystyle\widetilde{H}_{0}^{(q)}:\Sigma_{j_{q}k_{q}}^{\varepsilon}(u,v)=0~{}\mbox{for
any }(u,v)\in{\cal U}^{2}\mbox{ and }(j_{q},k_{q})\in{\cal
H}~{}\mbox{vs}~{}\widetilde{H}_{1}^{(q)}:\widetilde{H}_{0}^{(q)}\mbox{ is not
true}\,,$ (15)

for $q\in[Q]=p(p-1)/2$.

Let
$\widehat{\boldsymbol{\beta}}_{jk}(\cdot,\cdot)=\\{\hat{\beta}_{jk,\ell}(\cdot,\cdot)\\}_{\ell\neq
j,k}^{{\mathrm{\scriptscriptstyle T}}}$ be some consistent estimator of
$\boldsymbol{\beta}_{jk}(\cdot,\cdot)=\\{\beta_{jk,\ell}(\cdot,\cdot)\\}_{\ell\neq
j,k}^{{\mathrm{\scriptscriptstyle T}}}$. Define the functional residuals by

$\hat{\varepsilon}_{i,jk}(u):=Y_{ij}(u)-\int_{{\cal U}}{\bf
Y}_{i,-j,-k}^{{\mathrm{\scriptscriptstyle
T}}}(v)\widehat{\boldsymbol{\beta}}_{jk}(u,v)\,{\rm
d}v={\varepsilon}_{i,jk}(u)+\text{deviation},$ (16)

which takes a similar form as (2) within the error-contamination framework.
Let the sample covariance functions between residuals be
$\widetilde{\Sigma}_{jk}^{\varepsilon}(u,v)=(n-1)^{-1}\sum_{i=1}^{n}\\{\hat{\varepsilon}_{i,jk}(u)-\bar{\hat{\varepsilon}}_{jk}(u)\\}\\{\hat{\varepsilon}_{i,kj}(v)-\bar{\hat{\varepsilon}}_{kj}(v)\\}$
with
$\bar{\hat{\varepsilon}}_{jk}(u)=n^{-1}\sum_{t=1}^{n}\hat{\varepsilon}_{i,jk}(u)$
for any $j,k\in[p]$. Building upon the methodology in Section 3, for each
$q\in[Q]$, we propose the test statistic for $\widetilde{H}_{0}^{(q)}$ as

$\displaystyle\widetilde{T}_{n}^{\varepsilon,(q)}=n\|\widetilde{\Sigma}_{j_{q}k_{q}}^{\varepsilon}\|_{{\cal
S}}^{2}\,.$

The above formulation then allows for the direct adaptation of the notation
and the multiple testing procedure from Section 3 to the functional graphical
model setting by substituting the true signal pairs
$\\{X_{ij}(\cdot),X_{ik}(\cdot)\\}$ with
$\\{\varepsilon_{i,jk}(\cdot),\varepsilon_{i,kj}(\cdot)\\}$, and the estimated
signal pairs $\\{\widehat{X}_{ij}(\cdot),\widehat{X}_{ik}(\cdot)\\}$ with
$\\{\hat{\varepsilon}_{i,jk}(\cdot),\hat{\varepsilon}_{i,kj}(\cdot)\\}$ for
any $1\leq j<k\leq p$. To save space, we annotate quantities related to the
functional graphical model with a superscript $\varepsilon$. For notational
consistency, we will continue to use FDP and FDR to denote the false discovery
proportion and rate, respectively, where ${\rm
FDP}(t)={R_{0}^{\varepsilon}(t)}/\\{R^{\varepsilon}(t)\vee 1\\}$ in this
example.

We then follow the proposed multiple testing procedure in Section 3 based on
$\widetilde{T}_{n}^{\varepsilon,(q)}$. Similar to Section 4.1, we need to
verify Condition 5 by establishing the convergence rates in (11) within the
functional graphical model framework. Denote
$\widehat{\boldsymbol{\Xi}}(u,v)=(n-1)^{-1}\sum_{i=1}^{n}\\{{\bf
Y}_{i}(u)-\widebar{\bf Y}(u)\\}\\{{\bf Y}_{i}(v)-\widebar{\bf
Y}(v)\\}^{{\mathrm{\scriptscriptstyle T}}}$ and $\widebar{\bf
Y}(u)=n^{-1}\sum_{i=1}^{n}{\bf Y}_{i}(u)$. Let $\boldsymbol{\Xi}_{-j,-k}(u,v)$
and $\widehat{\boldsymbol{\Xi}}_{-j,-k}(u,v)$ denote the submatrices of
$\boldsymbol{\Xi}(u,v)$ and $\widehat{\boldsymbol{\Xi}}(u,v)$ by removing the
$(j,k)$th rows and $(j,k)$th columns. Suppose that, for any $1\leq j\neq k\leq
p,$ $\widehat{\boldsymbol{\beta}}_{jk}(\cdot,\cdot)$ satisfies

$\displaystyle\|\widehat{\boldsymbol{\beta}}_{jk}-\boldsymbol{\beta}_{jk}\|_{{\cal
S},1}=O_{\rm p}$ $\displaystyle(\delta_{1n})\,,$ (17)
$\displaystyle\Big{\|}\int_{\cal
U}\widehat{\boldsymbol{\Xi}}_{-j,-k}(\cdot,s)(\widehat{\boldsymbol{\beta}}_{jk}-\boldsymbol{\beta}_{jk})(\cdot,s)\,$
$\displaystyle{\rm d}s\Big{\|}_{{\cal S},\max}=O_{\rm p}(\delta_{2n})\,.$ (18)

for some rates $\delta_{1n}$ and $\delta_{2n}$. Condition 5 then holds by
fulfilling

$\delta_{1n}\sqrt{\log p/n}\asymp n^{-a_{1}}\text{ and
}\delta_{1n}\delta_{2n}\asymp n^{-a_{2}}.$ (19)

See (S.15) and (S.16) in Section A.6 of the Supplementary Material for
technical details. Therefore, to verify Condition 5, we only need to specify
the rates in (17) and (18). To achieve this, we will introduce an estimation
procedure to obtain $\widehat{\boldsymbol{\beta}}_{jk}(\cdot,\cdot)$ and then
derive the corresponding rates in (17) and (18) respectively in Propositions 2
and 3 below.

Recall that our framework replaces $\\{X_{ij}(\cdot),X_{ik}(\cdot)\\}$ and
$\\{\widehat{X}_{ij}(\cdot),\widehat{X}_{ik}(\cdot)\\}$ in Section 3 by
$\\{\varepsilon_{i,jk}(\cdot),\varepsilon_{i,kj}(\cdot)\\}$ and
$\\{\hat{\varepsilon}_{i,jk}(\cdot),\hat{\varepsilon}_{i,kj}(\cdot)\\}$,
respectively, for any $1\leq j<k\leq p$. Applying a generalized version of
Theorem 4 as discussed in Remark 3, we establish the next theorem, which
guarantees that our proposed multiple testing procedure for functional
graphical model can effectively control false discoveries asymptotically.

###### Theorem 5

Suppose that Conditions 1–4 hold for
$\\{\varepsilon_{i,jk}(\cdot)\\}_{j,k\in[p]},$ $\max_{q,q^{\prime}\in[Q]}|{\rm
Corr}(\tilde{\zeta}_{q}^{\varepsilon},\tilde{\zeta}_{q^{\prime}}^{\varepsilon})|\leq
c_{\tilde{\zeta}^{\varepsilon}}$ for some constant
$c_{\tilde{\zeta}^{\varepsilon}}\in(0,1)$, and
$\max_{q\in[Q]}|\tilde{\mathcal{S}}_{q}^{\varepsilon}(\tilde{\gamma}^{\varepsilon})|=o(p^{2\tilde{\nu}^{\varepsilon}})$
for some constants $\tilde{\gamma}^{\varepsilon}>0$ and
$0<\tilde{\nu}^{\varepsilon}<(1-c_{\tilde{\zeta}^{\varepsilon}})/(1+c_{\tilde{\zeta}^{\varepsilon}})$.
If $p\lesssim n^{\tilde{\kappa}^{\varepsilon}}$ for some constant
${\tilde{\kappa}^{\varepsilon}}>0$ and (19) is satisfied with
$a_{1},a_{2}>1/2$, then it holds that $\limsup_{n,Q\rightarrow\infty}{\rm
FDR}({\tilde{t}}^{\varepsilon})\leq\alpha Q_{0}^{\varepsilon}/Q$ and
$\lim_{n,Q\rightarrow\infty}{\mathbb{P}}\\{{\rm
FDP}({\tilde{t}}^{\varepsilon})\leq\alpha
Q_{0}^{\varepsilon}/Q+\widetilde{\varkappa}^{\varepsilon}\\}=1$ for any
$\widetilde{\varkappa}^{\varepsilon}>0$.

We next develop a three-step procedure to estimate the functional coefficients
${\boldsymbol{\beta}}_{jk}(\cdot,\cdot)$ for $1\leq j\neq k\leq p$. To ensure
a feasible solution under a high-dimensional regime, we impose a sparsity
assumption on $\\{\boldsymbol{\beta}_{jk}(\cdot,\cdot)\\}$, i.e.,
$\boldsymbol{\beta}_{jk}(\cdot,\cdot)$ is functional sparse with support set
$S_{jk}=\big{\\{}\ell\in[p]\setminus\\{j,k\\}:\|\beta_{jk,\ell}\|_{\mathcal{S}}\neq
0\big{\\}}$ and its cardinality $|S_{jk}|=s_{jk}\ll p.$ Let
$s=\max_{jk}s_{jk}$.

In the first step, due to the infinite-dimensionality of functional data, we
approximate each $Y_{ij}(\cdot)$ for $j\in[p]$ under the Karhunen-Loève
expansion truncated at $d_{j}$ (to be specified in Section 5), serving as the
foundation of functional principal components analysis (FPCA),
$Y_{ij}(\cdot)\approx\sum_{m=1}^{d_{j}}\xi_{ijm}\phi_{jm}(\cdot)=\boldsymbol{\xi}_{ij}^{{\mathrm{\scriptscriptstyle
T}}}\boldsymbol{\phi}_{j}(\cdot),$ where
$\boldsymbol{\xi}_{ij}=(\xi_{ij1},\dots,\xi_{ijd_{j}})^{{\mathrm{\scriptscriptstyle
T}}},$
$\boldsymbol{\phi}_{j}(\cdot)=\\{\phi_{j1}(\cdot),\dots,\phi_{jd_{j}}(\cdot)\\}^{{\mathrm{\scriptscriptstyle
T}}}.$ Here $\xi_{ijm}=\langle Y_{ij},\phi_{jm}\rangle,$ namely FPC scores,
correspond to a sequence of random variables with
$\text{Cov}(\xi_{ijm},\xi_{ijm^{\prime}})=\omega_{jm}I(m=m^{\prime}),$ where
$\omega_{j1}\geq\omega_{j2}\geq\dots>0$ are the eigenvalues of
$\Xi_{jj}(\cdot,\cdot),$ the $(j,j)$th entry of
$\boldsymbol{\Xi}(\cdot,\cdot)$, and $\phi_{j1}(\cdot),\phi_{j2}(\cdot),\dots$
are the corresponding eigenfunctions. To implement FPCA based on $n$
observations, we perform eigenanlysis of $\widehat{\Xi}_{jj}(\cdot,\cdot),$
the $(j,j)$th entry of $\widehat{\boldsymbol{\Xi}}(\cdot,\cdot),$ and obtain
estimated eigenvalue/eigenvector pairs
$\\{\hat{\omega}_{jm},\widehat{\phi}_{jm}(\cdot)\\}_{m\in[d_{j}]}.$ The
estimated FPC scores are $\hat{\xi}_{ijm}=\langle
Y_{ij},\hat{\phi}_{jm}\rangle$ for $m\in[d_{j}].$

In the second step, some calculations lead to the representation of (14) as

$\displaystyle\boldsymbol{\xi}_{ij}^{{\mathrm{\scriptscriptstyle
T}}}=\sum_{\ell\neq j,k}\boldsymbol{\xi}_{i\ell}^{{\mathrm{\scriptscriptstyle
T}}}\boldsymbol{\Psi}_{jk,\ell}+{\mathbf{r}}_{i,jk}^{{\mathrm{\scriptscriptstyle
T}}}+\boldsymbol{\epsilon}_{i,jk}^{{\mathrm{\scriptscriptstyle T}}},$

where $\boldsymbol{\Psi}_{jk,\ell}=\int_{{\cal U}}\int_{{\cal
U}}\boldsymbol{\phi}_{\ell}(v)\beta_{jk,\ell}(u,v)\boldsymbol{\phi}_{j}(u)^{{\mathrm{\scriptscriptstyle
T}}}{\rm d}u{\rm d}v\in\mathbb{R}^{d_{\ell}\times d_{j}}$ for
$\ell\in[p]\setminus\\{j,k\\},$ ${\mathbf{r}}_{i,jk}$ and
$\boldsymbol{\epsilon}_{i,jk}$ are formed by truncation and random errors,
respectively. See (S.19) of the Supplementary Material for detailed
expressions of ${\mathbf{r}}_{i,jk}$ and $\boldsymbol{\epsilon}_{i,jk}$.
Notice that $\boldsymbol{\Psi}_{jk,\ell}=\boldsymbol{0}$ if
$\beta_{jk,\ell}(\cdot,\cdot)=0,$ we can rely on the block sparsity in
$\\{\boldsymbol{\Psi}_{jk,\ell}\\}$ to recover the functional sparsity in
$\\{\beta_{jk,\ell}(\cdot,\cdot)\\}.$ Towards this, we let
$\widehat{\boldsymbol{\xi}}_{ij}\\!=\\!(\hat{\xi}_{ij1},\dots,\hat{\xi}_{ijd_{j}})^{{\mathrm{\scriptscriptstyle
T}}}$ and
$\widehat{\boldsymbol{\phi}}_{j}(\cdot)\\!=\\!\\{\hat{\phi}_{j1}(\cdot),\dots,\hat{\phi}_{jd_{j}}(\cdot)\\}^{{\mathrm{\scriptscriptstyle
T}}},$ and then propose to solve a standardized group lasso Simon and
Tibshirani (2012) problem by minimizing the following penalized least squares
criterion over
$\\{\boldsymbol{\Psi}_{jk,\ell}\\}_{\ell\in[p]\setminus\\{j,k\\}}$:

$\frac{1}{2}\sum_{i=1}^{n}\|\widehat{\boldsymbol{\xi}}_{ij}^{{\mathrm{\scriptscriptstyle
T}}}-\sum_{\ell\neq
j,k}\widehat{\boldsymbol{\xi}}_{i\ell}^{{\mathrm{\scriptscriptstyle
T}}}\boldsymbol{\Psi}_{jk,\ell}\|^{2}+\tau_{n}\sum_{\ell\neq
j,k}\|\widehat{\bf V}_{\ell}\boldsymbol{\Psi}_{jk,\ell}\|_{\text{F}},$ (20)

where $\widehat{\bf
V}_{\ell}=(\widehat{\boldsymbol{\xi}}_{1\ell},\dots,\widehat{\boldsymbol{\xi}}_{n\ell})^{{\mathrm{\scriptscriptstyle
T}}}\in\mathbb{R}^{n\times d_{\ell}}$ for $\ell\in[p]\setminus\\{j,k\\}$ and
$\tau_{n}\geq 0$ is the regularization parameter. Let
$\\{\widehat{\boldsymbol{\Psi}}_{jk,\ell}\\}_{\ell\in[p]\setminus\\{j,k\\}}$
be the minimizer of (20).

In the third step, we estimate functional coefficients by

$\hat{\beta}_{jk,\ell}(u,v)=\widehat{\boldsymbol{\phi}}_{\ell}(v)^{{\mathrm{\scriptscriptstyle
T}}}\widehat{\boldsymbol{\Psi}}_{jk,\ell}\widehat{\boldsymbol{\phi}}_{j}(u)~{}~{}\text{for}~{}(u,v)\in{\cal
U}^{2}~{}\text{and}~{}\ell\in[p]\setminus\\{j,k\\}.$

Before presenting the theoretical results, we list some regularity conditions.
To simplify notation, we assume the same $d=d_{j}$ across $j\in[p]$, but our
theoretical results below can be extended naturally to the more general
setting where $d_{j}$’s are different.

###### Condition 10

For each $j\in[p]$, $\omega_{j1}>\omega_{j2}>\cdots>0$. There exists some
constant $\varpi>1$ such that $\omega_{jm}-\omega_{j(m+1)}\gtrsim
m^{-\varpi-1}$ and $\omega_{jm}\lesssim m^{-\varpi}$ for $m\in\mathbb{N}_{+}$.

###### Condition 11

Denote the diagonal function matrix by $\widetilde{\bf
D}_{0}=\text{diag}(\Xi_{11},\dots,\Xi_{pp}).$ The infimum
$\underline{\mu}=\inf\limits_{\boldsymbol{\Phi}\in\bar{\mathbb{H}}_{0}^{p}}\frac{\iint\boldsymbol{\Phi}^{{\mathrm{\scriptscriptstyle
T}}}(u)\boldsymbol{\Xi}(u,v)\boldsymbol{\Phi}(v)\,{\rm d}u{\rm
d}v}{\iint\boldsymbol{\Phi}^{{\mathrm{\scriptscriptstyle
T}}}(u)\widetilde{{\bf D}}_{0}(u,v)\boldsymbol{\Phi}(v)\,{\rm d}u{\rm
d}v}\gtrsim sd^{\varpi+3}\\{{\log(pd)}/{n}\\}^{1/2},$ where
$\boldsymbol{\Phi}\in\bar{\mathbb{H}}_{0}^{p}=\\{\boldsymbol{\Phi}\in\mathbb{H}^{p}:\iint\boldsymbol{\Phi}^{{\mathrm{\scriptscriptstyle
T}}}(u)\widetilde{{\bf D}}_{0}(u,v)\boldsymbol{\Phi}(v)\,{\rm d}u{\rm
d}v\in(0,\infty)\\}.$

###### Condition 12

For each $j,k\in[p]$ and $\ell\in S_{jk},$
$\beta_{jk,\ell}(u,v)=\sum_{m_{1},m_{2}=1}^{\infty}b_{jk,\ell
m_{1}m_{2}}\phi_{jm_{1}}(u)\phi_{\ell m_{2}}(v)$ and there exists some
constant $\upsilon>\varpi/2+1$ such that $|b_{jk,\ell
m_{1}m_{2}}|\lesssim(m_{1}+m_{2})^{-\upsilon-1/2}$ for $m_{1},m_{2}\geq 1.$

Conditions 10 and 12 are standard in functional linear regression literature
Kong et al. (2016); Guo and Qiao (2023) with parameter $\varpi$ capturing the
decay rate of eigenvalues and parameter $\upsilon$ controlling the level of
smoothness in nonzero coefficient functions. Condition 11 can be understood as
requiring a lower bound on the minimum eigenvalue of the correlation function
matrix of ${\bf Y}_{i}(\cdot).$ The following propositions verify the rates in
(17) and (18), respectively, for
$\widehat{\boldsymbol{\beta}}_{jk}(\cdot,\cdot)$ obtained through the proposed
estimation procedure. This implies that Condition 5 holds for the functional
graphical model.

######  Proposition 2 (verify the rate in (17))

Suppose that Conditions 10–12 hold. If $n\gtrsim\log(pd)d^{4\varpi+4},$ then
there exist some positive constants $c_{1}$ and $c_{2}$ such that, for any
regularization parameter, $\tau_{n}\gtrsim
s[d^{\varpi+2}\\{{\log(pd)}/{n}\\}^{1/2}+d^{-\upsilon+1}]$, the estimate
$\widehat{\boldsymbol{\beta}}_{jk}(\cdot,\cdot)$ satisfies

$\begin{split}\|\widehat{\boldsymbol{\beta}}_{jk}-\boldsymbol{\beta}_{jk}\|_{{\cal
S},1}\lesssim\underline{\mu}^{-1}{d^{\varpi/2}s\tau_{n}}\asymp\delta_{1n},\end{split}$

with probability greater than $1-c_{1}(pd)^{-c_{2}}.$

######  Proposition 3 (verify the rate in (18))

Suppose that the conditions in Proposition 2 hold, then with probability
greater than $1-c_{1}(pd)^{-c_{2}},$

$\Big{\|}\int_{\cal
U}\widehat{\boldsymbol{\Xi}}_{-j,-k}(\cdot,s)(\widehat{\boldsymbol{\beta}}_{jk}-\boldsymbol{\beta}_{jk})(\cdot,s)\,{\rm
d}s\Big{\|}_{{\cal
S},\max}\lesssim\tau_{n}+\underline{\mu}^{-1}d^{-\varpi/2+1}s\tau_{n}\asymp\delta_{2n}.$

###### Remark 5

Our procedure can be extended naturally to handle the functional graphical
model with partial observations by applying the methodology developed in
Section 4.2 to reconstructed curves obtained through the pre-smoothing step
introduced in Section 4.1. To provide theoretical support for our multiple
testing procedure, we still need to verify Condition 5 by establishing the
corresponding rates in (17) and (18). However, the initial pre-smoothing step
results in additional estimation errors, which makes the subsequent
verification much more involved. Empirically, we demonstrate the effectiveness
of such procedure in controlling false discoveries through extensive
simulations in Section 5.2.

##  5 Simulation studies

We conduct a number of simulations to illustrate the proposed multiple testing
procedures for the functional covariance model and functional graphical model
in Sections 5.1 and 5.2, respectively, for both fully and partially observed
functional data.

###  5.1 Functional covariance model

In each scenario, we generate functional variables by
$X_{ij}(u)={\mathbf{s}}(u)^{{\mathrm{\scriptscriptstyle
T}}}\boldsymbol{\theta}_{ij}$ for $i\in[n],j\in[p]$ and $u\in{\cal U}=[0,1],$
where ${\mathbf{s}}(u)$ is a $10$-dimensional Fourier basis function and
$\boldsymbol{\theta}_{i}=(\boldsymbol{\theta}_{i1}^{{\mathrm{\scriptscriptstyle
T}}},\dots,\boldsymbol{\theta}_{ip}^{{\mathrm{\scriptscriptstyle
T}}})^{{\mathrm{\scriptscriptstyle T}}}\in\mathbb{R}^{10p}$ is generated from
a mean zero multivariate Gaussian distribution with block covariance matrix
$\boldsymbol{\Lambda}\in\mathbb{R}^{10p\times 10p},$ whose $(j,k)$th block is
$\boldsymbol{\Lambda}_{jk}\in\mathbb{R}^{10\times 10}$ for $j,k\in[p].$ The
functional sparsity pattern in
$\boldsymbol{\Sigma}(\cdot,\cdot)=\\{\Sigma_{jk}(\cdot,\cdot)\\}_{j,k\in[p]}$
with its $(j,k)$th entry
$\Sigma_{jk}(u,v)={\mathbf{s}}(u)^{{\mathrm{\scriptscriptstyle
T}}}\boldsymbol{\Lambda}_{jk}{\mathbf{s}}(v)$ can thus be characterized by the
block sparsity structure in $\boldsymbol{\Lambda}.$ Define
$\boldsymbol{\Lambda}_{jk}=\Pi_{jk}\boldsymbol{\Delta}$ for $j,k\in[p]$ with
$\boldsymbol{\Delta}=\text{diag}(1^{-2},\dots,10^{-2}),$ which results in
polynomially decaying eigenvalues of $\Sigma_{jj}(\cdot,\cdot)$ for each $j.$

We start by comparing the performance of our proposed multiple testing
procedure for the functional covariance model (MFC) with Benjamini–Hochberg
procedure (BH) (Benjamini and Hochberg, 1995) and Bonferroni correction
procedure (BC) under sparse covariance settings. To this end, we generate
$\boldsymbol{\Pi}\equiv(\Pi_{jk})_{p\times p}={\bf A}\odot{\bf
B}+c^{\prime}{\bf I}_{p},$ where $\odot$ denotes the Hadamard product, ${\bf
A}=(A_{jk})_{p\times p}$ is a symmetric (0,1) matrix with 1 on the diagonal,
${\bf B}\in\mathbb{R}^{p\times p}$ is a symmetric matrix with entries being
i.i.d. sampled from $\text{Unif}[0.2,0.3]$, and
$c^{\prime}=\big{(}-\lambda_{\min}({\bf A}\odot{\bf B})\big{)}_{+}+0.01$ to
guarantee the positive definiteness of $\boldsymbol{\Pi},$ with
$\lambda_{\min}({\bf M})$ denoting the minimum eigenvalue of matrix ${\bf M}.$
Let the set $S_{A}=\\{(j,k):1\leq j<k\leq p,A_{jk}\neq 0\\}$ with its degree
$|S_{A}|=s_{A}$. We consider $n=200$, $p=30$ and examine the empirical FDR by

$\displaystyle\frac{1}{B}\sum_{b=1}^{B}\frac{\sum_{q\in{\cal
H}_{0}}I\\{V_{q}^{(b)}\geq\hat{t}^{(b)}\\}}{1\vee\sum_{q=1}^{Q}I\\{V_{q}^{(b)}\geq\hat{t}^{(b)}\\}}\,$
(21)

over $B=1000$ replications, where $V_{q}^{(b)}$ and $\hat{t}^{(b)}$ are
specified in Algorithm 1 for the $b$th replication. Figure 1 presents the
empirical FDRs of the three competing procedures and the frequency that
$\hat{t}$ exists for the proposed method at the $5\%$ nominal level when
$s_{A}\in\\{10,11,\dots,35\\}$. It is evident that BH tends to result in
inflated FDRs and thus cannot control the FDP if the number of true
alternatives is fixed (very sparse case). Similar patterns can also be found
in Liu and Shao (2014); Xia, Cai and Cai (2018).

Figure 1: Empirical FDRs of MFC, BH and BC (left) and the frequency that
$\hat{t}$ exists (right) over 1000 simulation runs.

We next assess the performance of the proposed MFC, BC and two commonly-used
thresholding strategies, i.e., hard and soft functional thresholding methods,
under both fully and partially observed functional scenarios. The $(j,k)$th
entries of hard and soft functional thresholding estimators are respectively
defined as $\widehat{\Sigma}_{jk}I(\|\widehat{\Sigma}_{jk}\|_{\cal
S}\geq\tilde{\tau})$ and
$\widehat{\Sigma}_{jk}(1-\tilde{\tau}/\|\widehat{\Sigma}_{jk}\|_{{\cal
S}})_{+},$ where $\widehat{\Sigma}_{jk}$ is specified in (3) and the
thresholding parameter $\tilde{\tau}>0$ is selected by the cross-validation
approach (Fang et al., 2023).

For each method, we first generate fully observed curves $X_{ij}(\cdot)$ under
different functional sparsity patterns in $\boldsymbol{\Sigma}(\cdot,\cdot)$,
i.e., block sparsity structures in $\boldsymbol{\Pi}$ as follows.

  * •

Model 1 (block banded). For $j,k\in[p],$ we generate
$\Pi_{jk}=(1-{|j-k|}/{3})_{+}$.

  * •

Model 2 (block sparse without any special structure). We generate
$\boldsymbol{\Pi}={\bf B}+c^{\prime\prime}{\bf I}_{p},$ where entries of ${\bf
B}$ are sampled independently from $\text{Unif}[0.3,0.8]$ with probability
$3/p$ or $0$ with probability $1-3/p,$ and
$c^{\prime\prime}=\big{(}-\lambda_{\min}({\bf B})\big{)}_{+}+0.01.$

We then generate the partially observed curves satisfying (12), where the
observational time points $U_{ijt}$ and errors $\varsigma_{ijt}$ are sampled
independently from $\text{Unif}[0,1]$ and ${\cal N}(0,1),$ respectively. We
consider settings of $T_{ij}\in\\{25,51\\}$ and ultilize the Gaussian kernel
with the optimal bandwidths proportional to $T_{ij}^{-1/5},$ as suggested in
Zhang and Chen (2007), when implementing the local linear smoothers to obtain
estimated curves $\widehat{X}_{ij}(\cdot)$. Let ${\cal H}_{1}={\cal
H}\setminus{\cal H}_{0}$ be the set of true alternatives. Define the empirical
power as

$\displaystyle\frac{1}{100}\sum_{b=1}^{100}\frac{1}{|{\cal
H}_{1}|}\sum_{q\in{\cal H}_{1}}I\\{V_{q}^{(b)}\geq\hat{t}^{(b)}\\}\,.$ (22)

We finally report the empirical FDRs and the empirical powers over $100$
replications under both fully and partially observed functional scenarios with
$n\in\\{100,200\\}$ and $p\in\\{30,60\\}$ at $\alpha=0.1$ in Table
LABEL:table.sim.12 and at $\alpha=0.05$ in Table S2 of the Supplementary
Material. Note that the empirical FDRs and powers for the partially observed
cases are calculated with $V_{q}^{(b)}$ and $\hat{t}^{(b)}$ in (21) and (22),
respectively, replaced by $\widetilde{V}_{q}^{(b)}$ and $\tilde{t}^{(b)},$ as
specified in Section 3. Several conclusions can be drawn from Tables
LABEL:table.sim.12 and S2. First, the empirical FDRs of the proposed MFC are
well maintained below the target FDR levels for both fully and partially
observed functional data. In contrast, BC tends to be conservative with
empirical FDRs substantially lower than the target levels. However, the soft
thresholding method deteriorates significantly with highly elevated FDRs.
Interestingly, the hard thresholding method seems to work well under these two
models, with lower FDRs and powers. Second, MFC consistently achieves the
highest empirical powers among four competitors across all settings. As
expected, we observe enhanced powers as the number of subjects $n$ increases
or as the dimension $p$ declines. Third, all methods exhibit higher powers
with increasing sampling frequency $T_{ij}$ and demonstrate comparable
performance when curves are densely observed, compared to the fully observed
functional case.

###  5.2 Functional graphical model

We now turn to the functional graphical model estimation. Different from the
data generating process in Section 5.1, we first generate the functional
components
$\tilde{\varepsilon}_{ij}(u)={\mathbf{s}}(u)^{{\mathrm{\scriptscriptstyle
T}}}\widetilde{\boldsymbol{\theta}}_{ij}$ for $i\in[n]$ and $j=[p],$ where
each $\widetilde{\boldsymbol{\theta}}_{ij}$ is sampled independently from
${\cal N}({\bf 0},\boldsymbol{\Delta}),$ and then sequentially generate fully
observed curves $Y_{i1}(\cdot),\dots,Y_{ip}(\cdot)$. To be specific, we begin
by establishing directed acyclic graphs according to the following functional
structural models and then moralize the directed graph to obtain the
undirected graph Cowell et al. (2007).

  * •

Model 3 (banded sparse). We generate $Y_{i1}(u)=\tilde{\varepsilon}_{i1}(u),$
$Y_{i2}(u)=\tilde{\varepsilon}_{i2}(u)$ and

$Y_{ij}(u)=\sum_{k=1}^{2}\int_{{\cal
U}}Y_{i(j-k)}(v)\beta_{j(j-k)}(u,v)dv+\tilde{\varepsilon}_{ij}(u)~{}~{}\text{for}~{}~{}j\in[p]\setminus[2].$

  * •

Model 4 (randomly sparse). For $j\in[p/3],$ we generate
$Y_{ij}(u)=\tilde{\varepsilon}_{ij}(u).$ To determine the directed edge set
$E_{D},$ we then randomly select one or two directed edges with equal
probability from the candidate directed edge set $\\{(k,j):1\leq k<j\\}$ for
each node $j\in\\{p/3+1,\dots,p\\}$ in a sequential way. Then

$Y_{ij}(u)=\sum_{(k,j)\in E_{D}}\int_{{\cal
U}}Y_{ik}(v)\beta_{jk}(u,v)dv+\tilde{\varepsilon}_{ij}(u)~{}~{}\text{for}~{}~{}j\in[p]\setminus[p/3],$

which leads to an overall edge connection probability of approximately $3/p$
in the undirected graph.

We generate functional coefficients
$\beta_{jk}(u,v)={\mathbf{s}}(u)^{{\mathrm{\scriptscriptstyle T}}}{\bf
B}_{jk}{\mathbf{s}}(v)$ for $(u,v)\in{\cal U}^{2}$ with ${\bf
B}_{jk}=(B_{jklm})_{l,m\in[10]}$ for $(k,j)\in E_{D}$ and
$B_{jklm}=(-1)^{l+m}s_{D,j}^{-1}c_{B}(l+m)^{-2},$ which results in
polynomially decaying basis coefficients. Here $c_{B}$’s are sampled
independently from $\text{Unif}[4,6],$ and the normalization term is
$s_{D,j}=|E_{D,j}|,$ where $E_{D,j}=\\{(k,j)\in E_{D},1\leq k<j\\}$ is defined
as the directed edge set for node $j.$ We finally generate partially observed
curves according to (12) and the same procedure adopted in Section 5.1 with
sampling frequency of $T_{ij}$ in $\\{25,51\\}$.

We consider $n=200,400$ subjects of $p=30,60$ functional variables, and repeat
each simulation 100 times. To select the truncated dimension $d_{j}$ for each
$j\in[p]$, we take the standard approach by selecting the largest $d_{j}$
eigenvalues of $\widehat{\Xi}_{jj}(\cdot,\cdot)$ such that the cumulative
percentage of selected eigenvalues exceeds $95\%$. Following the proposal of
Xia, Cai and Cai (2018), we select the optimal regularization parameter
$\hat{\tau}_{n}$ with the principle of making
$\sum_{q\in\mathcal{H}_{0}^{\varepsilon}}I(\widetilde{V}_{q}^{\varepsilon}\geq
t)$ and $Q\\{1-\Phi(t)\\}$ as close as possible. Specifically, for a sequence
of $\tau_{n}$ values, we obtain the estimated functional coefficients
$\\{\widehat{\boldsymbol{\beta}}_{jk,\tau_{n}}(\cdot,\cdot)\\}_{1\leq j\neq
k\leq p}$, construct the corresponding transformed statistics
$\widetilde{V}_{q,\tau_{n}}^{\varepsilon}$ for each $\tau_{n}$ and choose the
optimal $\hat{\tau}_{n}$ as the minimizer of

$\sum_{l=1}^{10}\bigg{\\{}\frac{\sum_{q\in\mathcal{H}}I(\widetilde{V}_{q,\tau_{n}}^{\varepsilon}\geq\Phi^{-1}[1-l\\{1-\Phi(\sqrt{\log
p})\\}/10])}{Q\cdot l\\{1-\Phi(\sqrt{\log p})\\}/10}-1\bigg{\\}}^{2}.$

We examine the performance of our proposed multiple testing procedure for the
functional graphical model (MFG) based on empirical FDR and power, and compare
it with the performance of BC, hard and soft functional thresholding
procedures. The numerical results at the $10\%$ and $5\%$ nominal levels are
summarized in Table LABEL:table.sim.34 and Table S2 of the Supplementary
Material, respectively. Similar conclusions can be drawn compared to the
numerical summaries obtained for the functional covariance model in Section
5.1. Furthermore, we observe that the hard thresholding approach fails
completely under both model settings. One possible explanation is that the
hard and soft thresholding methods adopt a uniform threshold level across all
pairs $(j,k)$ and thus are not able to handle the heteroscedastic problem of
the large-scale cross-covariance function estimation, i.e.
$\\{\Sigma_{jk}^{\varepsilon}(u,v)\\}_{1\leq j<k\leq p}$.

##  6 Real data analysis

###  6.1 HCP data

In this section, we apply our proposed MFC to functional magnetic resonance
imaging (fMRI) data from the Human Connectome Project (HCP). This dataset
includes resting-state fMRI scans of subjects, recorded every 0.72 seconds at
$L=1200$ measurement locations (14.4 minutes), along with their corresponding
fluid intelligence gF scores. These scores are measured by Raven’s Progressive
Matrices and serve as indicators of the subjects’ ability to process new
information, learn, and solve problems Cattell (1987). Specifically, we
consider $n_{\scriptstyle\textup{low}}=73$ subjects with $\textit{gF}\leq 8$
and $n_{\scriptstyle\textup{high}}=85$ subjects with $\textit{gF}\geq 23$, and
employ a $30$-second observation window ($T_{ij}=43$) starting from the 12th
minute to investigate the relationship between the short time scale functional
connectivity (Hutchison et al., 2013) and gF. The preprocessing of the raw
fMRI data involved the ICA-FIX preprocessed pipeline Glasser et al. (2013) and
a standard band-pass filter at $[0.01,0.08]$ Hz to filter out frequency bands
irrelevant to resting state functional connectivity Biswal et al. (1995).
Following the proposals of Zapata et al. (2022) and Miao et al. (2022), we
treat signals from different regions of interest (ROIs) as multivariate
functional data. We then implement MFC to construct the resting state
functional connectivity networks depicting the marginal correlation structures
among $p=83$ ROIs of three well-acknowledged modules in neuroscience study
Finn et al. (2015): the medial frontal (29 ROIs), frontoparietal (34 ROIs) and
default mode modules (20 ROIs), for two groups of subjects, respectively.

Figure 2 displays the functional connectivity networks based on MFC at the
target FDR level $\alpha=10\%$. Out of a total of $3403$ pairs, we identify
826 connections in subjects with $\textit{gF}\geq 23$ and $365$ in the group
with $\textit{gF}\leq 8$. We use the p-value, i.e. $\widetilde{\rm
pv}^{(l)}={\mathbb{P}}\\{T_{0}^{(l)}\geq\widetilde{T}_{n}^{(l)}\\}$ defined in
Section 3 as the connectivity measure for the $l$th identified pairs, with a
larger p-value indicating reduced connectivity strength. It is observable that
individuals with higher gF scores have significantly stronger connectivity
within the medial frontal and frontoparietal modules. Notably, despite having
a larger number of connections, the functional connectivity in the default
mode network appears to be weaker for the higher gF group compared to those
with lower gF scores. These patterns nicely align with the findings in Finn et
al. (2015).

(a) $\textit{gF}\leq 8$: the medial frontal module

(b) $\textit{gF}\leq 8$: the frontoparietal module

(c) $\textit{gF}\leq 8$: the default mode module

(d) $\textit{gF}\geq 23$: the medial frontal module

(e) $\textit{gF}\geq 23$: the frontoparietal module

(f) $\textit{gF}\geq 23$: the default mode module

Figure 2:  The connectivity strengths at fluid intelligence $\textit{gF}\leq
8$ and $\textit{gF}\geq 23$. Salmon, orange and yellow nodes represent the
ROIs in the medial frontal, frontoparietal and default mode modules,
respectively. The edge color from cyan to blue corresponds to the $p$-value
from large to small.

###  6.2 EEG data

In this section, we illustrate our proposed MFG using the electroencephalogram
(EEG) data from an alcoholism study Zhang et al. (1995). The study consists of
122 subjects, out of which $n_{a}=77$ are in the alcoholic group and
$n_{c}=45$ are in the control group. Each subject was exposed to a picture
stimulus while brain activities were measured at 256 time points over a one-
second time interval at $p=64$ electrodes/nodes. We follow the preprocessing
procedure in Zhu et al. (2016) and Qiao et al. (2019), which involved
averaging and filtering the EEG signals at $\alpha$ frequency bands across all
trials. In our analysis, we removed three anatomical landmarked electrodes, X,
Y and nd. Since the duration for the stimulus in each trial was 300 ms and the
subjects’ average response time for each stimulus was 653.45 ms, we consider
the time from 300–655 ms, which covers the time period of interest, with
$T_{ij}=90$ sampled time points in total.

We apply MFG to build brain networks across $p=61$ nodes. To approximate the
$\chi^{2}$-type mixture in Theorem 3, we adopted a fast dimension-reduction-
based approach Rice and Shum (2019). Figure 3 plots the brain networks
constructed at the target FDR level of $\alpha=2\%$, revealing that $6\%$ and
$5.3\%$ of node pairs are recognized as significant connections for the
alcohol and control groups, respectively. Unlike the symmetric connectivity
observed in the control group, the alcohol group exhibits enhanced
connectivity in the left frontal lobe area. This is consistent with existing
neuroscience literature documenting a compensatory elevation in brain
activation and functional connectivity, especially in the left frontal region,
among individuals with alcohol dependence to achieve comparable levels of task
performance to non-alcoholics (Desmond et al., 2003; Chanraud et al., 2011).
Additionally, the corpus callosum area (midline nodes in Figure 3)
demonstrates a significant number of connections for both groups, highlighting
its crucial role in facilitating information transfer between the two
hemispheres of the brain (Bloom and Hynd, 2005; van der Knaap and van der Ham,
2011). These observations again validate the scientific reliability of our
proposed method.

(a) alcohol group

(b) control group

Figure 3: Brain networks based on MFG for the alcoholic and control groups.

##  7 Discussion

Our test statistic $T_{n,jk}$ (or $\widetilde{T}_{n,jk}$) is constructed based
on the $L_{2}$-norm of the difference between the estimated cross-covariance
function $\widehat{\Sigma}_{jk}(u,v)$ (or $\widetilde{\Sigma}_{jk}(u,v)$) and
the true one $\Sigma_{jk}(u,v)$. For functional data with local spikes, one
may use the supermum-norm based test statistics,

${T}_{n,jk}’=n^{1/2}\sup_{(u,v)\in{\cal
U}^{2}}|\widehat{\Sigma}_{jk}(u,v)|\,,~{}~{}~{}\widetilde{T}_{n,jk}^{\prime}=n^{1/2}\sup_{(u,v)\in{\cal
U}^{2}}|\widetilde{\Sigma}_{jk}(u,v)|\,,$

upon which our multiple testing procedures in Sections 2.2 and 3 remain valid
under the associated technical analysis. We next provide a sketch for
establishing the theoretical properties of the supermum-norm based test
statistics. For fully observed functional data, under Conditions 1–2 and
regular continuity condition for $Z_{ij}(u)$, we can use the central limit
theorem for empirical process and continuous mapping theorem to obtain the
asymptotic properties of ${T}_{n,jk}^{\prime}$ under the null and alternative
hypotheses. It is not difficult to know the limiting null distribution of
${T}_{n,jk}^{\prime}$ is $\sup_{(u,v)\in{\cal U}^{2}}|G_{jk}(u,v)|$, where
$G_{jk}(u,v)$ follows a mean-zero Gaussian process with covariance function
$\Gamma_{jk}(u_{1},v_{1},u_{2},v_{2})$. As $\sup_{(u,v)\in{\cal
U}^{2}}|G_{jk}(u,v)|$ is usually intractable, a parametric bootstrap method
could be employed to determine critical values, which would also be used to
construct the multiple testing. To investigate the theoretical properties of
the multiple testing, we could apply the partition technique to reducing the
problem from supremum over ${\cal U}^{2}$ to the maximum over a grid of pairs,
and then adopt the idea of Gaussian approximation in Chernozhukov et al.
(2017) and the theory of multiple testing in Chang et al. (2023) to derive the
theoretical results. Furthermore, with the aid of initial pre-smoothing, this
procedure can be extended naturally to the densely observed functional
scenario. Although similar technical tools could be used, we need to address
the additional complexity that arises from the estimation errors in the pre-
smoothing step.

## References

  * (1)
  * Benjamini and Hochberg (1995) Benjamini, Y. and Hochberg, Y. (1995).  Controlling the false discovery rate: A practical and powerful approach to multiple testing, Journal of the Royal Statistical Society: Series B 57: 289–300.
  * Biswal et al. (1995) Biswal, B., Zerrin Yetkin, F., Haughton, V. M. and Hyde, J. S. (1995).  Functional connectivity in the motor cortex of resting human brain using echo-planar MRI, Magnetic resonance in medicine 34: 537–541.
  * Bloom and Hynd (2005) Bloom, J. S. and Hynd, G. W. (2005).  The role of the corpus callosum in interhemispheric transfer of information: Excitation or inhibition?, Neuropsychology review 15: 59–71.
  * Cai and Liu (2016) Cai, T. and Liu, W. (2016).  Large-scale multiple testing of correlations, Journal of the American Statistical Association 111: 229–240.
  * Cai (2017) Cai, T. T. (2017).  Global testing and large-scale multiple testing for high-dimensional covariance structures, Annual Review of Statistics and Its Application 4: 423–446.
  * Cattell (1987) Cattell, R. B. (1987).  Intelligence: Its structure, growth and action, Elsevier.
  * Chang et al. (2024) Chang, J., Chen, C., Qiao, X. and Yao, Q. (2024).  An autocovariance-based learning framework for high-dimensional functional time series, Journal of Econometrics 239: Paper No. 105385, 25.
  * Chang et al. (2023) Chang, J., He, J., Kang, J. and Wu, M. (2023).  Statistical inference for complex dependence of multimodal imaging data, Journal of the American Statistical Association 119: 1486–1499.
  * Chanraud et al. (2011) Chanraud, S., Pitel, A.-L., Pfefferbaum, A. and Sullivan, E. V. (2011).  Disruption of functional connectivity of the default-mode network in alcoholism, Cerebral cortex 21: 2272–2281.
  * Chernozhukov et al. (2017) Chernozhukov, V., Chetverikov, D. and Kato, K. (2017).  Central limit theorems and bootstrap in high dimensions, The Annals of Probability 45: 2309–2352.
  * Cowell et al. (2007) Cowell, R. G., Dawid, P., Lauritzen, S. L. and Spiegelhalter, D. J. (2007).  Probabilistic Networks and Expert Systems, Information Science and Statistics, Springer New York.
  * Desmond et al. (2003) Desmond, J. E., Chen, S. A., DeRosa, E., Pryor, M. R., Pfefferbaum, A. and Sullivan, E. V. (2003).  Increased frontocerebellar activation in alcoholics during verbal working memory: An fMRI study, Neuroimage 19: 1510–1520.
  * Fan et al. (2015) Fan, Y., James, G. M. and Radchenko, P. (2015).  Functional additive regression, The Annals of Statistics 43: 2296–2325.
  * Fang et al. (2023) Fang, Q., Guo, S. and Qiao, X. (2023).  Adaptive functional thresholding for sparse covariance function estimation in high dimensions, Journal of the American Statistical Association 119: 1473–1485.
  * Finn et al. (2015) Finn, E. S., Shen, X., Scheinost, D., Rosenberg, M. D., Huang, J., Chun, M. M., Papademetris, X. and Constable, R. T. (2015).  Functional connectome fingerprinting: Identifying individuals using patterns of brain connectivity, Nature Neuroscience 18: 1664–1671.
  * Glasser et al. (2013) Glasser, M. F., Sotiropoulos, S. N., Wilson, J. A., Coalson, T. S., Fischl, B., Andersson, J. L., Xu, J., Jbabdi, S., Webster, M., Polimeni, J. R. et al. (2013).  The minimal preprocessing pipelines for the human connectome project, Neuroimage 80: 105–124.
  * Guo and Qiao (2023) Guo, S. and Qiao, X. (2023).  On consistency and sparsity for high-dimensional functional time series with application to autoregressions, Bernoulli 29: 451–472.
  * Hu and Yao (2022) Hu, X. and Yao, F. (2022).  Sparse functional principal component analysis in high dimensions, Statistica Sinica 32: 1939–1960.
  * Hutchison et al. (2013) Hutchison, R. M., Womelsdorf, T., Allen, E. A., Bandettini, P. A., Calhoun, V. D., Corbetta, M., Della Penna, S., Duyn, J. H., Glover, G. H., Gonzalez-Castillo, J., Handwerker, D. A., Keilholz, S., Kiviniemi, V., Leopold, D. A., de Pasquale, F., Sporns, O., Walter, M. and Chang, C. (2013).  Dynamic functional connectivity: Promise, issues, and interpretations, NeuroImage 80: 360–378.  Mapping the Connectome.
  * Kong et al. (2016) Kong, D., Xue, K., Yao, F. and Zhang, H. H. (2016).  Partially functional linear regression in high dimensions, Biometrika 103: 147–159.
  * Lee et al. (2023) Lee, K.-Y., Ji, D., Li, L., Constable, T. and Zhao, H. (2023).  Conditional functional graphical models, Journal of the American Statistical Association 118: 257–271.
  * Li and Solea (2018) Li, B. and Solea, E. (2018).  A nonparametric graphical model for functional data with application to brain networks based on fMRI, Journal of the American Statistical Association 113: 1637–1655.
  * Liu et al. (2009) Liu, H., Tang, Y. and Zhang, H. H. (2009).  A new chi-square approximation to the distribution of non-negative definite quadratic forms in non-central normal variables, Computational Statistics & Data Analysis 53: 853–856.
  * Liu (2013) Liu, W. (2013).  Gaussian graphical model estimation with false discovery rate control, The Annals of Statistics 41: 2948––2978.
  * Liu and Shao (2014) Liu, W. and Shao, Q.-M. (2014).  Phase transition and regularized bootstrap in large-scale $t$-tests with false discovery rate control, The Annals of Statistics 42: 2003–2025.
  * Luo and Qi (2017) Luo, R. and Qi, X. (2017).  Function-on-function linear regression by signal compression, Journal of the American Statistical Association 112: 690–705.
  * Miao et al. (2022) Miao, R., Zhang, X. and Wong, R. K. (2022).  A wavelet-based independence test for functional data with an application to MEG functional connectivity, Journal of the American Statistical Association 118: 1876–1889.
  * Qiao et al. (2019) Qiao, X., Guo, S. and James, G. (2019).  Functional graphical models, Journal of the American Statistical Association 114: 211–222.
  * Rice and Shum (2019) Rice, G. and Shum, M. (2019).  Inference for the lagged cross-covariance operator between functional time series, Journal of Time Series Analysis 40: 665–692.
  * Simon and Tibshirani (2012) Simon, N. and Tibshirani, R. (2012).  Standardization and the group Lasso penalty, Statistica Sinica 22: 983–1001.
  * Solea and Li (2022) Solea, E. and Li, B. (2022).  Copula Gaussian graphical models for functional data, Journal of the American Statistical Association 117: 781–793.
  * Tsai et al. (2024) Tsai, K., Zhao, B., Koyejo, S. and Kolar, M. (2024).  Latent multimodal functional graphical model estimation, Journal of the American Statistical Association, in press .
  * van der Knaap and van der Ham (2011) van der Knaap, L. J. and van der Ham, I. J. (2011).  How does the corpus callosum mediate interhemispheric transfer? A review, Behavioural brain research 223: 211–221.
  * Xia, Cai and Cai (2018) Xia, Y., Cai, T. and Cai, T. T. (2018).  Multiple testing of submatrices of a precision matrix with applications to identification of between pathway interactions, Journal of the American Statistical Association 113: 328–339.
  * Xia, Cai and Li (2018) Xia, Y., Cai, T. T. and Li, H. (2018).  Joint testing the false discovery rate control in high-dimensional multivairate regression, Biometrika 105: 249–269.
  * Xia et al. (2020) Xia, Y., Lockhart, S. N. and Jagust, W. J. (2020).  Simultaneous covariance inference for multimodal integrative analysis, Journal of the American Statistical Association 115: 1279–1291.
  * Xue et al. (2023) Xue, K., Yang, J. and Yao, F. (2023).  Optimal linear discriminant analysis for high-dimensional functional data, Journal of the American Statistical Association 119: 1055–1064.
  * Zapata et al. (2022) Zapata, J., Oh, S. Y. and Petersen, A. (2022).  Partial separability and functional graphical models for multivariate Gaussian processes, Biometrika 109: 665–681.
  * Zhang (2013) Zhang, J.-T. (2013).  Analysis of Variance for Functional Data, CRC Press.
  * Zhang and Chen (2007) Zhang, J.-T. and Chen, J. (2007).  Statistical inferences for functional data, The Annals of Statistics 35: 1052–1079.
  * Zhang et al. (1995) Zhang, X. L., Begleiter, H., Porjesz, B., Wang, W. and Litke, A. (1995).  Event related potentials during object recognition tasks, Brain research bulletin 38: 531–538.
  * Zhang and Wang (2016) Zhang, X. and Wang, J.-L. (2016).  From sparse to dense functional data and beyond, The Annals of Statistics 44: 2281–2321.
  * Zhu et al. (2016) Zhu, H., Strawn, N. and Dunson, D. B. (2016).  Bayesian graphical models for multivariate functional data, Journal of Machine Learning Research 17: 1–27.

Supplementary material to “Large-scale multiple testing of cross-covariance
functions with applications to functional network models”

Qin Fang, Qing Jiang and Xinghao Qiao

This supplementary material contains proofs of main theoretical results in
Section A, additional technical proofs in Section B, and additional empirical
results in Section C.

##  Appendix A Proofs of main theoretical results

###  A.1 Proof of Theorem 1

Proof of Theorem 1(i). Notice that for any $j,k\in[p]$,

$\displaystyle\widehat{\Sigma}_{jk}(u,v)=$
$\displaystyle~{}\frac{1}{n-1}\sum_{i=1}^{n}\\{X_{ij}(u)-\widebar{X}_{j}(u)\\}\\{X_{ik}(v)-\widebar{X}_{k}(v)\\}$
$\displaystyle=$
$\displaystyle~{}\frac{1}{n-1}\sum_{i=1}^{n}\\{Z_{ij}(u)-\bar{Z}_{j}(u)\\}\\{Z_{ik}(v)-\bar{Z}_{k}(v)\\}$
$\displaystyle=$
$\displaystyle~{}\frac{1}{n-1}\sum_{i=1}^{n}Z_{ij}(u)Z_{ik}(v)-\frac{n}{n-1}\bar{Z}_{j}(u)\bar{Z}_{k}(v)$
$\displaystyle\overset{\triangle}{=}$
$\displaystyle~{}\frac{1}{n-1}\sum_{i=1}^{n}z_{ijk}(u,v)-\frac{n}{n-1}\bar{Z}_{j}(u)\bar{Z}_{k}(v)$

where $\widebar{X}_{j}(u)=n^{-1}\sum_{i=1}^{n}X_{ij}(u)$,
$\bar{Z}_{j}(u)=n^{-1}\sum_{i=1}^{n}Z_{ij}(u)$ and $z_{ijk}(u,v),i\in[n]$ are
i.i.d. with ${\mathbb{E}}\\{z_{ijk}(u,v)\\}=\Sigma_{jk}(u,v)$ and

$\displaystyle\Gamma_{jk}(u_{1},v_{1},u_{2},v_{2})=$
$\displaystyle~{}{\text{Cov}}\\{z_{1jk}(u_{1},v_{1}),z_{1jk}(u_{2},v_{2})\\}$
$\displaystyle=$
$\displaystyle~{}{\mathbb{E}}\\{Z_{1j}(u_{1})Z_{1k}(v_{1})Z_{1j}(u_{2})Z_{1k}(v_{2})\\}-\Sigma_{jk}(u_{1},v_{1})\Sigma_{jk}(u_{2},v_{2})\,.$

By Cauchy–Schwarz inequality and Condition 2, it holds that

$\displaystyle{\mathbb{E}}(\|z_{1jk}\|_{{\cal S}}^{2})=$
$\displaystyle~{}{\mathbb{E}}\bigg{\\{}\iint z_{1jk}^{2}(u,v)\,{\rm d}u{\rm
d}v\bigg{\\}}={\mathbb{E}}\bigg{\\{}\iint Z_{1j}^{2}(u)Z_{1k}^{2}(v)\,{\rm
d}u{\rm d}v\bigg{\\}}$ $\displaystyle\leq$
$\displaystyle~{}\bigg{(}{\mathbb{E}}\bigg{[}\bigg{\\{}\int
Z_{1j}^{2}(u)\,{\rm
d}u\bigg{\\}}^{2}\bigg{]}\bigg{)}^{1/2}\bigg{(}{\mathbb{E}}\bigg{[}\bigg{\\{}\int
Z_{1k}^{2}(u)\,{\rm d}u\bigg{\\}}^{2}\bigg{]}\bigg{)}^{1/2}$ $\displaystyle=$
$\displaystyle~{}\\{{\mathbb{E}}(\|Z_{1j}\|^{4})\\}^{1/2}\\{{\mathbb{E}}(\|Z_{1k}\|^{4})\\}^{1/2}<\infty\,.$

Let ${\rm GP}(\mu,\Sigma)$ denotes a Gaussian process with mean function
$\mu(\cdot)$ and covariance function $\Sigma(\cdot,\cdot)$. By the central
limit theorem of i.i.d. stochastic processes [Theorem 4.12 of Zhang (2013)],
we have

$\displaystyle\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\\{z_{ijk}(u,v)-\Sigma_{jk}(u,v)\\}\overset{d}{\rightarrow}{\rm
GP}(0,\Gamma_{jk})\,,$ (S.1)

as $n\rightarrow\infty$. Notice that ${\mathbb{E}}\\{\bar{Z}_{j}(u)\\}=0$ for
any $j\in[p]$. Cauchy–Schwarz inequality and Condition 1 imply that for any
$j\in[p]$,

$\displaystyle\sup_{u,v\in{\cal
U}}{\text{Cov}}\\{\bar{Z}_{j}(u),\bar{Z}_{j}(v)\\}=\frac{1}{n}\sup_{u,v\in{\cal
U}}\Sigma_{jj}(u,v)\leq\frac{\sup_{u\in{\cal
U}}\Sigma_{jj}(u,u)}{n}\rightarrow
0,~{}~{}~{}\mbox{as}~{}~{}n\rightarrow\infty\,.$

Then it holds that $\bar{Z}_{j}(u)$ converges to 0 in probability uniformly
over ${\cal U}$ for any $[p]$ and
$\sqrt{n}\,\bar{Z}_{j}(u)\overset{d}{\rightarrow}{\rm GP}(0,\Sigma_{jj})$.
Combining with (S.1), we have

$\displaystyle\sqrt{n}\\{\widehat{\Sigma}_{jk}(u,v)-\Sigma_{jk}(u,v)\\}\overset{d}{\rightarrow}{\rm
GP}(0,\Gamma_{jk})$

as $n\rightarrow\infty$. Following from Condition 2, it holds that

$\displaystyle{\rm tr}(\Gamma_{jk})=$
$\displaystyle~{}\iint\Gamma_{jk}\\{(u,v),(u,v)\\}\,{\rm d}u{\rm d}v$
$\displaystyle\leq$
$\displaystyle~{}\iint[{\mathbb{E}}\\{Z_{1j}^{2}(u)Z_{1k}^{2}(v)\\}-\Sigma_{jk}^{2}(u,v)]\,{\rm
d}u{\rm d}v$ $\displaystyle\leq$
$\displaystyle~{}2\iint{\mathbb{E}}\\{Z_{1j}^{2}(u)Z_{1k}^{2}(v)\\}\,{\rm
d}u{\rm d}v<\infty\,.$

Then applying the continuous mapping theorem for random elements taking values
in a Hilbert space (Billingsley, 1999) and Theorem 4.2 of Zhang (2013), it
holds that

$\displaystyle
n\iint\\{\widehat{\Sigma}_{jk}(u,v)-\Sigma_{jk}(u,v)\\}^{2}\,{\rm d}u{\rm
d}v\overset{d}{\rightarrow}\sum_{r=1}^{\infty}\lambda_{jkr}A_{r}~{}~{}~{}{\rm
as}~{}~{}~{}n\rightarrow\infty\,,$ (S.2)

where $A_{r}\overset{i.i.d.}{\sim}\chi^{2}_{1}$. The results of Dauxois et al.
(1982) imply that
$|\hat{\lambda}_{jkr}-\lambda_{jkr}|\lesssim\lambda_{jkr}O_{\rm p}(n^{-1/2})$
for any $r\in\mathbb{N}_{+}$ and $j,k\in[p]$. Together with Condition 2, it
holds that
$|\sum_{r=1}^{\infty}\hat{\lambda}_{jkr}A_{r}-\sum_{r=1}^{\infty}\lambda_{jkr}A_{r}|\leq\sum_{r=1}^{\infty}\lambda_{jkr}O_{\rm
p}(n^{-1/2})=O_{\rm p}(n^{-1/2})=o_{\rm p}(1)$. Notice that
$\Sigma_{jk}(u,v)=0$ for any $u,v\in{\cal U}$ under $H_{0,jk}$. Thus we
complete the proof of Theorem 1(i). $\hfill\Box$
Proof of Theorem 1(ii). Notice that for any $j,k\in[p]$,

$\displaystyle T_{n,jk}=$
$\displaystyle~{}\iint[\sqrt{n}\\{\widehat{\Sigma}_{jk}(u,v)-\Sigma_{jk}(u,v)\\}]^{2}\,{\rm
d}u{\rm d}v+n\iint\Sigma_{jk}^{2}(u,v)\,{\rm d}u{\rm d}v$
$\displaystyle~{}+2\sqrt{n}\iint\Sigma_{jk}(u,v)[\sqrt{n}\\{\widehat{\Sigma}_{jk}(u,v)-\Sigma_{jk}(u,v)\\}]\,{\rm
d}u{\rm d}v\,,$

Since $\iint\Sigma_{jk}^{2}(u,v)\,{\rm d}u{\rm d}v>0$, then following from
(S.2), Cauchy–Schwarz inequality and $H_{1,jk}$, it holds that
$\iint[\sqrt{n}\\{\widehat{\Sigma}_{jk}(u,v)-\Sigma_{jk}(u,v)\\}]^{2}\,{\rm
d}u{\rm d}v=O_{\rm p}(1)$, $n\iint\Sigma_{jk}^{2}(u,v)\,{\rm d}u{\rm d}v=O(n)$
and
$2\sqrt{n}\iint\Sigma_{jk}(u,v)[\sqrt{n}\\{\widehat{\Sigma}_{jk}(u,v)-\Sigma_{jk}(u,v)\\}]\,{\rm
d}u{\rm d}v=O_{\rm p}(n^{1/2})$, which implies that

$\displaystyle{\mathbb{P}}\\{T_{n,jk}\geq T_{0,jk}(\alpha)\\}\rightarrow 1$

as $n\rightarrow\infty$ under $H_{1,jk}$. $\hfill\Box$

###  A.2  Proof of Theorem 2

To prove Theorem 2, we need the following lemma whose proof is given in
Section B.1.

###### Lemma A1

Suppose Condition 3 hold. For any $q\in\mathcal{H}_{0},$ the Karhunen-Loève
expansion allows us to represent
$z_{iq}(u,v):=Z_{ij_{q}}(u)Z_{ik_{q}}(v)=\sum_{r=1}^{\infty}a_{ir}^{(q)}\varphi_{r}^{(q)}(u,v),$
where $\lambda_{r}^{(q)}:=\lambda_{j_{q}k_{q}r}$,
$\varphi_{r}^{(q)}(u,v):=\varphi_{j_{q}k_{q}r}(u,v)$, and the coefficients
$a_{ir}^{(q)}=\iint z_{iq}(u,v)\varphi_{r}^{(q)}(u,v)\,{\rm d}u{\rm d}v$ are
uncorrelated random variables with mean zero and ${\rm
Cov}\\{a_{ir}^{(q)},a_{ir^{\prime}}^{(q)}\\}=\lambda_{r}^{(q)}I(r=r^{\prime})$.
Then there exist some constant $M>0$ such that
$\max_{i\in[n]}\max_{q\in\mathcal{H}_{0}}\sup_{r\in\mathbb{N}_{+}}{\mathbb{E}}\big{(}\exp[M|a_{ir}^{(q)}|/\\{\lambda_{r}^{(q)}\\}^{1/2}]\big{)}=O(1)$.

Let
$\mathring{T}_{n}^{(q)}=\iint\\{n^{-1/2}\sum_{i=1}^{n}z_{iq}(u,v)\\}^{2}\,{\rm
d}u{\rm d}v$, where $z_{iq}(u,v)=Z_{ij_{q}}(u)Z_{ik_{q}}(v)$ for any $i\in[n]$
and $q\in[Q]$. Recall $Q=p(p-1)$ and $p\leq Cn^{\kappa}$ for some constants
$C>0$ and $\kappa>0$. By the proof of Theorem 1, we have
$\max_{q\in[Q]}|T_{n}^{(q)}-\mathring{T}_{n}^{(q)}|=O_{\rm p}\\{n^{-1/2}(\log
p)^{1/2}\\}=o_{\rm p}\\{(\log p)^{-2}\\}$. Notice that
$z_{iq}(u,v)=\sum_{r=1}^{\infty}a_{ir}^{(q)}\varphi_{r}^{(q)}(u,v)$, which
implies

$\displaystyle\mathring{T}_{n}^{(q)}=\iint\bigg{\\{}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\sum_{r=1}^{\infty}a_{ir}^{(q)}\varphi_{r}^{(q)}(u,v)\bigg{\\}}^{2}\,{\rm
d}u{\rm
d}v=\frac{1}{n}\sum_{i_{1}=1}^{n}\sum_{i_{2}=1}^{n}\sum_{r=1}^{\infty}a_{i_{1}r}^{(q)}a_{i_{2}r}^{(q)}=\sum_{r=1}^{\infty}\bigg{\\{}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}a_{ir}^{(q)}\bigg{\\}}^{2}\,.$

Denote
$\mathring{T}_{n}^{(q,K)}=\sum_{r=1}^{K}\\{n^{-1/2}\sum_{i=1}^{n}a_{ir}^{(q)}\\}^{2}$.
By Markov’s inequality and Condition 4, we have

$\displaystyle~{}{\mathbb{P}}\bigg{\\{}\max_{q\in{\cal
H}_{0}}|\mathring{T}_{n}^{(q)}-\mathring{T}_{n}^{(q,K)}|\geq(\log
p)^{-1}\bigg{\\}}={\mathbb{P}}\bigg{[}\max_{q\in{\cal
H}_{0}}\sum_{r=K+1}^{\infty}\bigg{\\{}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}a_{ir}^{(q)}\bigg{\\}}^{2}\geq(\log
p)^{-1}\bigg{]}$ $\displaystyle\lesssim$ $\displaystyle~{}p^{2}(\log
p)\max_{q\in{\cal
H}_{0}}\sum_{r=K+1}^{\infty}{\mathbb{E}}\bigg{[}\bigg{\\{}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}a_{ir}^{(q)}\bigg{\\}}^{2}\bigg{]}\lesssim\frac{p^{2}(\log
p)}{n}\max_{q\in{\cal
H}_{0}}\sum_{r=K+1}^{\infty}\sum_{i=1}^{n}\lambda_{r}^{(q)}$
$\displaystyle\lesssim$ $\displaystyle~{}\frac{p^{2}(\log
p)}{K^{\eta_{1}-1}}=O(p^{-1})$

provided that $K=p^{(3+\epsilon)/(\eta_{1}-1)}$ for some sufficiently small
constant $\epsilon>0$. Let $\lambda_{0}=\max_{q\in{\cal
H}_{0}}\sup_{r\in\mathbb{N}_{+}}\lambda_{r}^{(q)}\lesssim 1$. Define
$\hat{a}_{ir}^{(q)}=a_{ir}^{(q)}I\\{|a_{ir}^{(q)}|\leq\pi_{n}\\}-{\mathbb{E}}[a_{ir}^{(q)}I\\{|a_{ir}^{(q)}|\leq\pi_{n}\\}]$
with
$\pi_{n}=\lambda_{0}^{1/2}\log(p+n)(4+4\eta_{1}+2\epsilon^{\prime})/\\{(\eta_{1}-1)M\\}$
for some sufficiently small constant $\epsilon^{\prime}>\epsilon$. Let
$\hat{T}_{n}^{(q,K)}=\sum_{r=1}^{K}\\{n^{-1/2}\sum_{i=1}^{n}\hat{a}_{ir}^{(q)}\\}^{2}$.
By Cauchy–Schwarz inequality and Lemma A1,

$\displaystyle\max_{q\in{\cal
H}_{0}}\sup_{r\in\mathbb{N}_{+}}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}{\mathbb{E}}[|a_{ir}^{(q)}|I\\{|a_{ir}^{(q)}|\geq\pi_{n}\\}]$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}=\max_{q\in{\cal
H}_{0}}\sup_{r\in\mathbb{N}_{+}}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\\{\lambda_{r}^{(q)}\\}^{1/2}{\mathbb{E}}\bigg{(}\frac{|a_{ir}^{(q)}|}{\\{\lambda_{r}^{(q)}\\}^{1/2}}I\bigg{[}\frac{|a_{ir}^{(q)}|}{\\{\lambda_{r}^{(q)}\\}^{1/2}}\geq\frac{\lambda_{0}^{1/2}(4+4\eta_{1}+2\epsilon^{\prime})}{\\{\lambda_{r}^{(q)}\\}^{1/2}(\eta_{1}-1)M}\log(p+n)\bigg{]}\bigg{)}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\lesssim
n^{1/2}\max_{q\in{\cal
H}_{0}}\sup_{r\in\mathbb{N}_{+}}\max_{i\in[n]}{\mathbb{E}}\bigg{(}\frac{|a_{ir}^{(q)}|}{\\{\lambda_{r}^{(q)}\\}^{1/2}}I\bigg{[}\frac{|a_{ir}^{(q)}|}{\\{\lambda_{r}^{(q)}\\}^{1/2}}\geq\frac{4+4\eta_{1}+2\epsilon^{\prime}}{(\eta_{1}-1)M}\log(p+n)\bigg{]}\bigg{)}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\lesssim\frac{n^{1/2}}{(p+n)^{(2+2\eta_{1}+\epsilon^{\prime})/(\eta_{1}-1)}}\max_{q\in{\cal
H}_{0}}\sup_{r\in\mathbb{N}_{+}}\max_{i\in[n]}{\mathbb{E}}\bigg{(}\frac{|a_{ir}^{(q)}|}{\\{\lambda_{r}^{(q)}\\}^{1/2}}\exp\bigg{[}\frac{M|a_{ir}^{(q)}|}{2\\{\lambda_{r}^{(q)}\\}^{1/2}}\bigg{]}\bigg{)}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\lesssim\frac{n^{1/2}}{(p+n)^{(2+2\eta_{1}+\epsilon^{\prime})/(\eta_{1}-1)}}\,.$

Hence,

$\displaystyle{\mathbb{P}}\bigg{\\{}\max_{q\in{\cal
H}_{0}}\max_{r\in[K]}\bigg{|}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}a_{ir}^{(q)}-\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\hat{a}_{ir}^{(q)}\bigg{|}\geq
p^{-(3+\eta_{1}+\epsilon^{\prime})/(\eta_{1}-1)}\bigg{\\}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}\leq{\mathbb{P}}\bigg{\\{}\max_{i\in[n]}\max_{q\in{\cal
H}_{0}}\max_{r\in[K]}|a_{ir}^{(q)}|\geq\pi_{n}\bigg{\\}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}\lesssim
np^{2}K\max_{i\in[n]}\max_{q\in{\cal
H}_{0}}\max_{r\in[K]}{\mathbb{P}}\bigg{[}\frac{|a_{ir}^{(q)}|}{\\{\lambda_{r}^{(q)}\\}^{1/2}}\geq\frac{4+4\eta_{1}+2\epsilon^{\prime}}{(\eta_{1}-1)M}\log(p+n)\bigg{]}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}\lesssim\frac{np^{2}K}{{(p+n)^{(4+4\eta_{1}+2\epsilon^{\prime})/(\eta_{1}-1)}}}\max_{i\in[n]}\max_{q\in{\cal
H}_{0}}\max_{r\in[K]}{\mathbb{E}}\bigg{(}\exp\bigg{[}\frac{M|a_{ir}^{(q)}|}{\\{\lambda_{r}^{(q)}\\}^{1/2}}\bigg{]}\bigg{)}=O(p^{-1})\,.$

Triangle inequality yields

$\displaystyle\max_{q\in{\cal
H}_{0}}|\mathring{T}_{n}^{(q,K)}-\hat{T}_{n}^{(q,K)}|\lesssim$
$\displaystyle~{}K\max_{q\in{\cal
H}_{0}}\max_{r\in[K]}\bigg{|}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\hat{a}_{ir}^{(q)}\bigg{|}\cdot\max_{q\in{\cal
H}_{0}}\max_{r\in[K]}\bigg{|}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}a_{ir}^{(q)}-\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\hat{a}_{ir}^{(q)}\bigg{|}$
$\displaystyle~{}+K\max_{q\in{\cal
H}_{0}}\max_{r\in[K]}\bigg{|}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}a_{ir}^{(q)}-\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\hat{a}_{ir}^{(q)}\bigg{|}^{2}$
$\displaystyle=$ $\displaystyle~{}o_{\rm p}\\{(\log p)^{-2}\\}\,.$

Therefore, we have ${\mathbb{P}}[\max_{q\in{\cal
H}_{0}}|{T}_{n}^{(q)}-\hat{T}_{n}^{(q,K)}|=o\\{(\log p)^{-1}\\}]=1$. Recall
$T_{0}^{(q)}\overset{d}{=}\sum_{r=1}^{\infty}\hat{\lambda}_{r}^{(q)}A_{r}$
with $A_{r}\overset{i.i.d}{\sim}\chi^{2}_{1}$. Denote
$T_{0}^{(q,K)}\overset{d}{=}\sum_{r=1}^{K}\hat{\lambda}_{r}^{(q)}A_{r}$ and
$\mathring{T}_{0}^{(q,K)}\overset{d}{=}\sum_{r=1}^{K}\lambda_{r}^{(q)}A_{r}$.
Since $|\hat{\lambda}_{r}^{(q)}-\lambda_{r}^{(q)}|=\lambda_{r}^{(q)}O_{\rm
p}(n^{-1/2})$ for any $r\in\mathbb{N}_{+}$ and $q\in[Q]$ based on Dauxois et
al. (1982), by Markov’s inequality, triangle inequality and Lemma A1, it holds
that

$\displaystyle{\mathbb{P}}\bigg{\\{}\max_{q\in{\cal
H}_{0}}|T_{0}^{(q,K)}-\mathring{T}_{0}^{(q,K)}|\geq(\log
p)^{-1}\bigg{\\}}\lesssim(\log p)\mathbb{E}\bigg{[}\max_{q\in{\cal
H}_{0}}\bigg{|}\sum_{r=1}^{K}\\{\hat{\lambda}_{r}^{(q)}-\lambda_{r}^{(q)}\\}A_{r}\bigg{|}\bigg{]}\lesssim\frac{\log
p}{n^{1/2}}=o(1)$

provided that $p\lesssim n^{\kappa}$ for some positive constant $\kappa$.
Following the similar arguments, we can also obtain
${\mathbb{P}}[\max_{q\in{\cal H}_{0}}|{T}_{0}^{(q)}-{T}_{0}^{(q,K)}|=o\\{(\log
p)^{-1}\\}]=1$. By Theorem 1.1 in Zaitsev (1987), for $q\in{\cal H}_{0}$, we
have

$\displaystyle{\mathbb{P}}\\{\hat{T}_{n}^{(q,K)}\geq t\\}\leq$
$\displaystyle~{}{\mathbb{P}}\\{\mathring{T}_{0}^{(q,K)}\geq
t-\epsilon_{n}(\log
p)^{-1/2}\\}+\tilde{c}_{1}K^{5/2}\exp\bigg{\\{}-\frac{n^{1/2}\epsilon_{n}}{\tilde{c}_{2}K^{3}\pi_{n}(\log
p)^{1/2}}\bigg{\\}}\,,$ $\displaystyle{\mathbb{P}}\\{\hat{T}_{n}^{(q,K)}\geq
t\\}\geq$ $\displaystyle~{}{\mathbb{P}}\\{\mathring{T}_{0}^{(q,K)}\geq
t+\epsilon_{n}(\log
p)^{-1/2}\\}-\tilde{c}_{1}K^{5/2}\exp\bigg{\\{}-\frac{n^{1/2}\epsilon_{n}}{\tilde{c}_{2}K^{3}\pi_{n}(\log
p)^{1/2}}\bigg{\\}}\,,$

where $\tilde{c}_{1},\tilde{c}_{2}>0$ are constants and $\epsilon_{n}=o(1)$.
Recall $V_{q}=\Phi^{-1}\\{1-{\rm pv}^{(q)}\\}$ with ${\rm
pv}^{(q)}={\mathbb{P}}\\{T_{0}^{(q)}\geq T_{n}^{(q)}\\}$, and
$F_{q}(z)=\mathbb{P}\\{T_{n}^{(q)}\leq z\\}$ for any $z\in\mathbb{R}$. Then on
the event $\\{\max_{q\in{\cal
H}_{0}}|{T}_{n}^{(q)}-\hat{T}_{n}^{(q,K)}|=o\\{(\log
p)^{-1}\\}\\}\bigcap\\{\max_{q\in{\cal
H}_{0}}|{T}_{0}^{(q)}-{T}_{0}^{(q,K)}|=o\\{(\log
p)^{-1}\\}\\}\bigcap\\{\max_{q\in{\cal
H}_{0}}|{T}_{0}^{(q,K)}-\mathring{T}_{0}^{(q,K)}|=o\\{(\log p)^{-1}\\}\\}$, by
the fact that $G(t+o\\{(\log p)^{-1/2}\\})/G(t)=1+o(1)$ uniformly in $0\leq
t\leq 2(\log p)^{1/2}$ with $G(t)=1-\Phi(t)$, it holds that for $q\in{\cal
H}_{0}$,

$\displaystyle{\mathbb{P}}\\{V_{q}\geq
t\\}={\mathbb{P}}\big{[}{\mathbb{P}}\\{T_{0}^{(q)}\leq
T_{n}^{(q)}\\}\geq\Phi(t)\big{]}=G(t)\\{1+o(1)\\}\,.$ (S.3)

Let $t_{\max}=(2\log Q-2\log\log Q)^{1/2}$. We consider two cases: (i) there
exists $t\in[0,t_{\max}]$ such that $\widehat{{\rm FDP}}(t)\leq\alpha$, and
(ii) $\widehat{{\rm FDP}}(t)>\alpha$ for any $t\in[0,t_{\max}]$.

In Case (i), by the definition of $\hat{t}$, it holds that
$\widehat{\mathrm{FDP}}(t)>\alpha$ for any $t<\hat{t}$. Then we have

$\displaystyle\frac{QG(t)}{1\vee\sum_{q\in[Q]}I(V_{q}\geq\hat{t})}\geq\frac{QG(t)}{1\vee\sum_{q\in[Q]}I(V_{q}\geq
t)}=\widehat{\mathrm{FDP}}(t)>\alpha\,,$

based on the fact that $I(V_{q}>t)\geq I(V_{q}>\hat{t})$ for $t<\hat{t}$. By
letting $t\uparrow\hat{t}$ in the numerator of the first term in above
inequality, we have $\widehat{\mathrm{FDP}}(\hat{t})\geq\alpha$. On the other
hand, based on the definition of $\hat{t}$, there exists a sequence
$\\{t_{i}\\}$ with $t_{i}\geq\hat{t}$ and $t_{i}\downarrow\hat{t}$ such that
$\widehat{\mathrm{FDP}}(t_{i})\leq\alpha$. Thus we have
$I(V_{q}\geq\hat{t})\geq I(V_{q}\geq t_{i})$, which implies that

$\displaystyle\frac{QG(t_{i})}{1\vee\sum_{q\in[Q]}I(V_{q}\geq\hat{t})}\leq\frac{QG(t_{i})}{1\vee\sum_{q\in[Q]}I(V_{q}\geq
t_{i})}\leq\alpha\,.$

Letting $t_{i}\downarrow\hat{t}$ in the numerator of the first term in above
inequality, we have $\widehat{\mathrm{FDP}}(\hat{t})\leq\alpha$. Therefore, we
have $\widehat{\mathrm{FDP}}(\hat{t})=\alpha$ in Case (i).

In Case (ii), we first show that the threshold of $\hat{t}$ at $(2\log
Q)^{1/2}$ leads to no false rejection with probability tending to 1. Notice
that for any $t>0$, $G(t)<\phi(t)/t$, where $\phi(\cdot)$ is the density
function of the standard normal distribution $N(0,1)$. Then if $\hat{t}=(2\log
Q)^{1/2}$, following from (S.3), we have

$\displaystyle{\mathbb{P}}\Bigg{\\{}\sum_{q\in\mathcal{H}_{0}}I(V_{q}\geq\hat{t})\geq
1\Bigg{\\}}\leq
Q_{0}\max_{q\in\mathcal{H}_{0}}{\mathbb{P}}(V_{q}\geq\hat{t})\leq
Q_{0}G\\{(2\log Q)^{1/2}\\}\\{1+o(1)\\}=o(1)$

as $Q\rightarrow\infty$, which implies that ${\mathbb{P}}\\{{\rm
FDP}(\hat{t})=0\\}\rightarrow 1$ in Case (ii), that is, the probability of
false rejection is tending to 0 asymptotically. Then

$\displaystyle\mathbb{P}\bigg{\\{}{\rm FDP}(\hat{t})\leq\frac{\alpha
Q_{0}}{Q}+\varepsilon\bigg{\\}}=$ $\displaystyle~{}\mathbb{P}\bigg{\\{}{\rm
FDP}(\hat{t})\leq\frac{\alpha Q_{0}}{Q}+\varepsilon,~{}{\rm
Case~{}(i)~{}holds}\bigg{\\}}$ $\displaystyle+\mathbb{P}\bigg{\\{}{\rm
FDP}(\hat{t})\leq\frac{\alpha Q_{0}}{Q}+\varepsilon,~{}{\rm
Case~{}(ii)~{}holds}\bigg{\\}}$ $\displaystyle=$
$\displaystyle~{}\mathbb{P}\big{\\{}{\rm
Case~{}(i)~{}holds}\big{\\}}+\mathbb{P}\big{\\{}{\rm
Case~{}(ii)~{}holds}\big{\\}}$ $\displaystyle-\mathbb{P}\bigg{\\{}{\rm
FDP}(\hat{t})>\frac{\alpha Q_{0}}{Q}+\varepsilon,~{}{\rm
Case~{}(i)~{}holds}\bigg{\\}}$ $\displaystyle=$
$\displaystyle~{}1-\mathbb{P}\bigg{\\{}{\rm FDP}(\hat{t})>\frac{\alpha
Q_{0}}{Q}+\varepsilon,~{}{\rm Case~{}(i)~{}holds}\bigg{\\}}\,.$ (S.4)

Notice that $\widehat{{\rm FDP}}(\hat{t})=\alpha$ and $\hat{t}\in[0,t_{\max}]$
under Case (i). If we can show

$\displaystyle\sup_{t\in[0,t_{\max}]}\bigg{|}\frac{\mathrm{FDP}(t)}{\widehat{\mathrm{FDP}}(t)}-\frac{Q_{0}}{Q}\bigg{|}\to
0~{}~{}\mbox{in probability}$ (S.5)

as $n,Q\rightarrow\infty$, then we have
$\lim_{n,Q\rightarrow\infty}\mathbb{P}\\{{\rm FDP}(\hat{t})\leq\alpha
Q_{0}/Q+\varepsilon\\}=1$ by (A.2).

To prove (S.5), it suffices to show that

$\displaystyle\sup_{t\in[0,t_{\max}]}\Bigg{|}\frac{\sum_{q\in\mathcal{H}_{0}}\\{I(V_{q}\geq
t)-G(t)\\}}{QG(t)}\Bigg{|}\to 0~{}~{}\mbox{in probability}$

as $n,Q\rightarrow\infty$. Let $0=t_{0}<t_{1}<\cdots<t_{s}=t_{\max}$ such that
$t_{i}-t_{i-1}=\tilde{v}$ for $i\in[s-1]$ and $t_{s}-t_{s-1}\leq\tilde{v}$,
where $\tilde{v}=\\{(\log Q)(\log\log Q)^{1/2}\\}^{-1/2}$. Then we have
$s\asymp t_{\max}/\tilde{v}$. For any $t\in[t_{i-1},t_{i}]$, it holds that

$\displaystyle\frac{\sum_{q\in\mathcal{H}_{0}}I(V_{q}\geq
t_{i})}{QG(t_{i})}\frac{G(t_{i})}{G(t_{i-1})}\leq\frac{\sum_{q\in\mathcal{H}_{0}}I(V_{q}\geq
t)}{QG(t)}\leq\frac{\sum_{q\in\mathcal{H}_{0}}I(V_{q}\geq
t_{i-1})}{QG(t_{i-1})}\frac{G(t_{i-1})}{G(t_{i})}\,.$

Notice that there exists a universal constant $C_{2}>0$ such that
$e^{-t^{2}/2}\leq\max(C_{2},2t)\int_{t}^{\infty}e^{-x^{2}/2}\,{\rm d}x$ for
any $t>0$. Then

$\displaystyle
0<1-\frac{G(t_{i})}{G(t_{i-1})}=\frac{\int_{t_{i-1}}^{t_{i}}e^{-x^{2}/2}\,{\rm
d}x}{\int_{t_{i-1}}^{\infty}e^{-x^{2}/2}\,{\rm
d}x}\leq\frac{\tilde{v}e^{-t_{i-1}^{2}/2}}{\int_{t_{i-1}}^{\infty}e^{-x^{2}/2}\,{\rm
d}x}\leq\tilde{v}\max(C_{2},2t_{i-1})\,,$

which implies that $\max_{i\in[s]}|1-G(t_{i})/G(t_{i-1})|\leq
2\tilde{v}t_{\max}\rightarrow 0$ as $Q\rightarrow\infty$. Thus, to prove
(S.5), it suffices to show that

$\displaystyle\max_{0\leq i\leq
s}\Bigg{|}\frac{\sum_{q\in\mathcal{H}_{0}}\\{I(V_{q}\geq
t_{i})-G(t_{i})\\}}{QG(t_{i})}\Bigg{|}\to 0~{}~{}\mbox{in probability}$

as $n,Q\rightarrow\infty$. Following from Bonferroni inequality and Markov’s
inequality, we have

$\displaystyle{\mathbb{P}}\bigg{[}\max_{0\leq i\leq
s}\bigg{|}\frac{\sum_{q\in\mathcal{H}_{0}}\\{I(V_{q}\geq
t_{i})-G(t_{i})\\}}{QG(t_{i})}\bigg{|}\geq\varkappa\bigg{]}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\leq\sum_{i=0}^{s}{\mathbb{P}}\bigg{[}\bigg{|}\frac{\sum_{q\in\mathcal{H}_{0}}\\{I(V_{q}\geq
t_{i})-G(t_{i})\\}}{QG(t_{i})}\bigg{|}\geq\varkappa\bigg{]}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\leq\frac{C}{\tilde{v}}\int_{0}^{t_{\max}}{\mathbb{P}}\bigg{[}\bigg{|}\frac{\sum_{q\in\mathcal{H}_{0}}\\{I(V_{q}\geq
t)-G(t)\\}}{QG(t)}\bigg{|}\geq\varkappa\bigg{]}\,\mathrm{d}t$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\leq\frac{C}{\tilde{v}}\int_{0}^{t_{\max}}{\mathbb{P}}\bigg{[}\bigg{|}\frac{\sum_{q\in\mathcal{H}_{0}}\\{I(V_{q}\geq
t)-{\mathbb{P}}(V_{q}\geq
t)\\}}{QG(t)}\bigg{|}\geq\frac{\varkappa}{2}\bigg{]}\,\mathrm{d}t$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\frac{C}{\tilde{v}}\int_{0}^{t_{\max}}{\mathbb{P}}\bigg{[}\bigg{|}\frac{\sum_{q\in\mathcal{H}_{0}}\\{{\mathbb{P}}(V_{q}\geq
t)-G(t)\\}}{QG(t)}\bigg{|}\geq\frac{\varkappa}{2}\bigg{]}\,\mathrm{d}t$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\leq\frac{C}{\tilde{v}\varkappa^{2}}\int_{0}^{t_{\max}}\mathbb{E}\bigg{[}\bigg{|}\frac{\sum_{q\in\mathcal{H}_{0}}\\{I(V_{q}\geq
t)-{\mathbb{P}}(V_{q}\geq t)\\}}{QG(t)}\bigg{|}^{2}\bigg{]}\,{\rm d}t$ (S.6)

for any $\varkappa>0$ and sufficiently large $n$, where $C$ is some positive
constant. The second step is based on the relationship between integration and
its associated Riemann sum, and the forth step is due to (S.3). Thus it
suffices to prove the following statement is true,

$\displaystyle~{}\int_{0}^{t_{\max}}\mathbb{E}\bigg{[}\bigg{|}\frac{\sum_{q\in\mathcal{H}_{0}}\\{I(V_{q}\geq
t)-{\mathbb{P}}(V_{q}\geq t)\\}}{QG(t)}\bigg{|}^{2}\bigg{]}\,{\rm d}t$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}=\int_{0}^{t_{\max}}\sum_{q,q^{\prime}\in\mathcal{H}_{0}}\frac{{\mathbb{P}}(V_{q}\geq
t,V_{q^{\prime}}\geq t)-{\mathbb{P}}(V_{q}\geq
t){\mathbb{P}}(V_{q^{\prime}}\geq t)}{Q^{2}G^{2}(t)}\,{\rm
d}t=o(\tilde{v})\,.$ (S.7)

We define three subsets

$\displaystyle\mathcal{H}_{01}=$
$\displaystyle~{}\\{(q,q^{\prime}):q,q^{\prime}\in\mathcal{H}_{0},\,q=q^{\prime}\\}\,,$
$\displaystyle\mathcal{H}_{02}=$
$\displaystyle~{}\\{(q,q^{\prime}):q,q^{\prime}\in\mathcal{H}_{0},\,q\neq
q^{\prime},q\in\mathcal{S}_{q^{\prime}}(\gamma)\mbox{~{}or~{}}q^{\prime}\in\mathcal{S}_{q}(\gamma)\\}\,,$
$\displaystyle\mathcal{H}_{03}=$
$\displaystyle~{}\\{(q,q^{\prime}):q,q^{\prime}\in\mathcal{H}_{0}\\}\backslash(\mathcal{H}_{01}\cup\mathcal{H}_{02})\,.$

Then
$\sum_{q,q^{\prime}\in\mathcal{H}_{0}}=\sum_{(q,q^{\prime})\in\mathcal{H}_{01}}+\sum_{(q,q^{\prime})\in\mathcal{H}_{02}}+\sum_{(q,q^{\prime})\in\mathcal{H}_{03}}$.

When $(q,q^{\prime})\in\mathcal{H}_{01}$, ${\mathbb{P}}(V_{q}\geq
t,V_{q^{\prime}}\geq t)-{\mathbb{P}}(V_{q}\geq
t){\mathbb{P}}(V_{q^{\prime}}\geq t)={\mathbb{P}}(V_{q}\geq
t)\\{1-{\mathbb{P}}(V_{q}\geq t)\\}\leq{\mathbb{P}}(V_{q}\geq t)$. By (S.3),
it holds that

$\displaystyle\sum_{(q,q^{\prime})\in\mathcal{H}_{01}}\frac{{\mathbb{P}}(V_{q}\geq
t,V_{q^{\prime}}\geq t)-{\mathbb{P}}(V_{q}\geq
t){\mathbb{P}}(V_{q^{\prime}}\geq
t)}{Q^{2}G^{2}(t)}\lesssim\frac{1}{QG(t)}\,.$

Due to $\int_{0}^{t_{\max}}\\{G(t)\\}^{-1}\,{\rm
d}t\lesssim\exp(2^{-1}t_{\max}^{2})=Q(\log Q)^{-1}$, we have

$\displaystyle\int_{0}^{t_{\max}}\sum_{(q,q^{\prime})\in\mathcal{H}_{01}}\frac{{\mathbb{P}}(V_{q}\geq
t,V_{q^{\prime}}\geq t)-{\mathbb{P}}(V_{q}\geq
t){\mathbb{P}}(V_{q^{\prime}}\geq t)}{Q^{2}G^{2}(t)}\,{\rm d}t\lesssim$
$\displaystyle~{}\int_{0}^{t_{\max}}\frac{{\rm d}t}{QG(t)}$
$\displaystyle\lesssim$ $\displaystyle~{}\frac{1}{\log Q}=o(\tilde{v})\,.$
(S.8)

Since $F_{q}\\{T_{n}^{(q)}\\}\sim U[0,1]$ and
$F_{q^{\prime}}\\{T_{n}^{(q^{\prime})}\\}\sim U[0,1]$, we have
${\mathbb{P}}(V_{q}\geq t,V_{q^{\prime}}\geq
t)\leq{\mathbb{P}}\big{[}F_{q}\\{T_{n}^{(q)}\\}\geq\Phi(t),F_{q^{\prime}}\\{T_{n}^{(q^{\prime})}\\}\geq\Phi(t)\big{]}\\{1+o(1)\\}\leq{\mathbb{P}}(\zeta_{q}\geq
t,\zeta_{q^{\prime}}\geq t)\\{1+o(1)\\}$, where
$\zeta_{q}=\Phi^{-1}[F_{q}\\{T_{n}^{(q)}\\}]\sim N(0,1)$ and
$\zeta_{q^{\prime}}=\Phi^{-1}[F_{q^{\prime}}\\{T_{n}^{(q^{\prime})}\\}]\sim
N(0,1)$. When $(q,q^{\prime})\in\mathcal{H}_{02}$, since $\max_{1\leq q\neq
q^{\prime}\leq Q}|\mathrm{Corr}(\zeta_{q},\zeta_{q^{\prime}})|\leq
c_{\zeta}<1$, Lemma 2 in Berman (1962) implies that
${\mathbb{P}}(\zeta_{q}\geq t,\zeta_{q^{\prime}}\geq t)\lesssim
t^{-2}\exp\\{-t^{2}/(1+c_{\zeta})\\}$ for any $t>C_{3}$, where $C_{3}>0$ is a
universal constant. Notice that
$e^{-t^{2}/2}\leq\max(C_{2},2t)\int_{t}^{\infty}e^{-x^{2}/2}\,{\rm d}x$ for
any $t>0$. Then ${\mathbb{P}}(V_{q}\geq t,V_{q^{\prime}}\geq t)\lesssim
t^{-2}\exp\\{-t^{2}/(1+c_{\zeta})\\}\lesssim
t^{-2}\\{\max(C_{2},2t)\\}^{2/(1+c_{\zeta})}\\{G(t)\\}^{2/(1+c_{\zeta})}$ for
any $t>\max\\{C_{3},C_{2}/2\\}$. Since
$\max_{q\in[Q]}|\mathcal{S}_{q}(\gamma)|=o(Q^{\nu})$, we have
$|\mathcal{H}_{02}|=O(Q^{1+\nu})$. Note that
$\nu<(1-c_{\zeta})/(1+c_{\zeta})<1$ and ${\mathbb{P}}(V_{q}\geq
t,V_{q^{\prime}}\geq t)\leq 1$ for any $0<t<\max\\{C_{3},C_{2}/2\\}$. We have

$\displaystyle\int_{0}^{t_{\max}}\sum_{(q,q^{\prime})\in\mathcal{H}_{02}}\frac{{\mathbb{P}}(V_{q}\geq
t,V_{q^{\prime}}\geq t)-{\mathbb{P}}(V_{q}\geq
t){\mathbb{P}}(V_{q^{\prime}}\geq t)}{Q^{2}G^{2}(t)}\,{\rm
d}t\leq\int_{0}^{t_{\max}}\sum_{(q,q^{\prime})\in\mathcal{H}_{02}}\frac{{\mathbb{P}}(V_{q}\geq
t,V_{q^{\prime}}\geq t)}{Q^{2}G^{2}(t)}\,{\rm d}t$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}\lesssim\int_{0}^{\max\\{1,C_{5},C_{3}/2\\}}\frac{{\rm
d}t}{Q^{1-\nu}G^{2}(t)}+\int_{\max\\{1,C_{3},C_{2}/2\\}}^{t_{\max}}\frac{{\rm
d}t}{Q^{1-\nu}\\{G(t)\\}^{2c_{\zeta}/(1+c_{\zeta})}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}\lesssim\frac{1}{Q^{1-\nu}}+\frac{Q^{(c_{\zeta}-1)/(1+c_{\zeta})+\nu}}{(\log
Q)^{(1+3c_{\zeta})/(2+2c_{\zeta})}}=o(\tilde{v})\,.$ (S.9)

When $(q,q^{\prime})\in\mathcal{H}_{03}$, denote $\rho_{q,q^{\prime}}={\rm
Corr}(\zeta_{q},\zeta_{q^{\prime}})$. By Theorem 2.1.e of Lin and Bai (2010),

${\mathbb{P}}(\zeta_{q}\geq t,\zeta_{q^{\prime}}\geq
t)\leq\left\\{\begin{aligned}
G(t)G\bigg{\\{}\frac{(1-\rho_{q,q^{\prime}})t}{(1-\rho_{q,q^{\prime}}^{2})^{1/2}}\bigg{\\}}\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}&\mbox{if}~{}-1<\rho_{q,q^{\prime}}\leq
0\,;\\\
(1+\rho_{q,q^{\prime}})G(t)G\bigg{\\{}\frac{(1-\rho_{q,q^{\prime}})t}{(1-\rho_{q,q^{\prime}}^{2})^{1/2}}\bigg{\\}}\,,~{}~{}~{}~{}&\mbox{if}~{}0\leq\rho_{q,q^{\prime}}<1\,.\end{aligned}\right.$
(S.10)

Note that $|\rho_{q,q^{\prime}}|\leq(\log Q)^{-2-\gamma}$ for any
$(q,q^{\prime})\in\mathcal{H}_{03}$. When $-(\log
Q)^{-2-\gamma}\leq\rho_{q,q^{\prime}}\leq 0$, due to
$(1-\rho_{q,q^{\prime}})/(1-\rho_{q,q^{\prime}}^{2})^{1/2}\geq 1$, we have
$G\\{(1-\rho_{q,q^{\prime}})t/(1-\rho_{q,q^{\prime}}^{2})^{1/2}\\}\leq G(t)$,
which implies that ${\mathbb{P}}(\zeta_{q}\geq t,\zeta_{q^{\prime}}\geq t)\leq
G^{2}(t)\leq\\{1+(\log Q)^{-1-\gamma}\\}G^{2}(t)$. When
$0<\rho_{q,q^{\prime}}\leq(\log Q)^{-2-\gamma}$, by the mean-value theorem,
there exists $t^{\prime}$ satisfying
$(1-\rho_{q,q^{\prime}})t/(1-\rho_{q,q^{\prime}}^{2})^{1/2}<t^{\prime}<t$ such
that

$\displaystyle
G\bigg{\\{}\frac{(1-\rho_{q,q^{\prime}})t}{(1-\rho_{q,q^{\prime}}^{2})^{1/2}}\bigg{\\}}=G(t)+\phi(t^{\prime})\bigg{\\{}t-\frac{(1-\rho_{q,q^{\prime}})t}{(1-\rho_{q,q^{\prime}}^{2})^{1/2}}\bigg{\\}}\,.$

Then we have

$\displaystyle
G\bigg{\\{}\frac{(1-\rho_{q,q^{\prime}})t}{(1-\rho_{q,q^{\prime}}^{2})^{1/2}}\bigg{\\}}\\{G(t)\\}^{-1}=$
$\displaystyle~{}1+\frac{\phi(t^{\prime})}{G(t)}\bigg{\\{}t-\frac{(1-\rho_{q,q^{\prime}})t}{(1-\rho_{q,q^{\prime}}^{2})^{1/2}}\bigg{\\}}$
$\displaystyle\leq$
$\displaystyle~{}1+\frac{t\phi(t^{\prime})}{t\phi(t)/(1+t^{2})}\bigg{\\{}1-\frac{1-\rho_{q,q^{\prime}}}{(1-\rho_{q,q^{\prime}}^{2})^{1/2}}\bigg{\\}}$
$\displaystyle\leq$ $\displaystyle~{}1+\frac{\phi(t^{\prime})}{\phi(t)}\cdot
2\rho_{q,q^{\prime}}(1+t^{2})$

for any $t>0$. For $0<t<t_{\max}$, it holds that
$\phi(t^{\prime})/\phi(t)=\exp\\{(t-{t}^{\prime})(t+{t}^{\prime})/2\\}\leq\exp[t_{\max}^{2}\\{1-(1-\rho_{q,q^{\prime}})/(1+\rho_{q,q^{\prime}}^{2})^{1/2}\\}]\leq\exp(2\rho_{q,q^{\prime}}t_{\max}^{2})\leq\exp\\{4(\log
Q)^{-1-\gamma}\\}\lesssim 1$, which implies that
$G\\{(1-\rho_{q,q^{\prime}})t/(1-\rho_{q,q^{\prime}}^{2})^{1/2}\\}\leq
G(t)[1+O\\{(\log Q)^{-1-\gamma}\\}]$ for any $t\in[0,t_{\max}]$ if
$0<\rho_{q,q^{\prime}}\leq(\log Q)^{-2-\gamma}$, where the term $O\\{(\log
Q)^{-1-\gamma}\\}$ holds uniformly over $t\in[0,t_{\max}]$. Then
$(1+\rho_{q,q^{\prime}})G(t)G\\{(1-\rho_{q,q^{\prime}})t/(1-\rho_{q,q^{\prime}}^{2})^{1/2}\\}\leq
G^{2}(t)[1+O\\{(\log Q)^{-1-\gamma}\\}]$ for any $t\in[0,t_{\max}]$ if
$0<\rho_{q,q^{\prime}}\leq(\log Q)^{-2-\gamma}$. By (S.10),
${\mathbb{P}}(\zeta_{q}\geq t,\zeta_{q^{\prime}}\geq t)\leq[1+O\\{(\log
Q)^{-1-\gamma}\\}]G^{2}(t)$ for $t\in[0,t_{\max}]$. Thus we have
$\max_{(q,q^{\prime})\in\mathcal{H}_{03}}{\mathbb{P}}(V_{q}\geq
t,V_{q^{\prime}}\geq t)\leq[1+O\\{(\log Q)^{-1-\gamma}\\}]G^{2}(t)$ for
$t\in[0,t_{\max}]$. It then holds that

$\displaystyle\int_{0}^{t_{\max}}\sum_{(q,q^{\prime})\in\mathcal{H}_{03}}\frac{{\mathbb{P}}(V_{q}\geq
t,V_{q^{\prime}}\geq t)-{\mathbb{P}}(V_{q}\geq
t){\mathbb{P}}(V_{q^{\prime}}\geq t)}{Q^{2}G^{2}(t)}\,{\rm d}t$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\leq O\\{(\log
Q)^{-1-\gamma}\\}\cdot\int_{0}^{t_{\max}}1\,{\rm d}t=o(\tilde{v})\,.$

Together with (S.8) and (S.9), we know (A.2) holds. Then
$\lim_{n,Q\rightarrow\infty}\mathbb{P}\\{{\rm FDP}(\hat{t})\leq\alpha
Q_{0}/Q+\varkappa\\}=1$ for any $\varkappa>0$. Since ${\rm
FDR}(t)={{\mathbb{E}}}\\{{\rm FDP}(t)\\}$, it holds that
$\limsup_{n,Q\to\infty}{\rm FDR}(\hat{t})\leq\alpha Q_{0}/Q$. We complete the
proof of Theorem 2. $\hfill\Box$

###  A.3 Proof of Theorem 3

Proof of Theorem 3(i). Notice that
$\widetilde{T}_{n,jk}=n\|\widetilde{\Sigma}_{jk}\|_{{\cal S}}^{2}$ and
$\widetilde{\Sigma}_{jk}(u,v)=(n-1)^{-1}\sum_{i=1}^{n}\\{\widehat{X}_{ij}-\widebar{\widehat{X}}_{j}(u)\\}\\{\widehat{X}_{ik}(v)-\widebar{\widehat{X}}_{k}(v)\\}$
with $\widehat{X}_{ij}(u)=X_{ij}(u)+e_{ij}(u)$. Then it holds that

$\displaystyle\widetilde{\Sigma}_{jk}(u,v)=$
$\displaystyle~{}\frac{1}{n-1}\sum_{i=1}^{n}\\{X_{ij}(u)+e_{ij}(u)-\widebar{X}_{j}(u)-\bar{e}_{j}(u)\\}\\{X_{ik}(v)+e_{ik}(v)-\widebar{X}_{k}(v)-\bar{e}_{k}(v)\\}$
$\displaystyle=$
$\displaystyle~{}\frac{1}{n-1}\sum_{i=1}^{n}\\{X_{ij}(u)-\widebar{X}_{j}(u)\\}\\{X_{ik}(v)-\widebar{X}_{k}(v)\\}$
$\displaystyle~{}+\frac{1}{n-1}\sum_{i=1}^{n}\\{X_{ij}(u)-\widebar{X}_{j}(u)\\}\\{e_{ik}(v)-\bar{e}_{k}(v)\\}$
$\displaystyle~{}+\frac{1}{n-1}\sum_{i=1}^{n}\\{e_{ij}(u)-\bar{e}_{j}(u)\\}\\{X_{ik}(v)-\widebar{X}_{k}(v)\\}$
$\displaystyle~{}+\frac{1}{n-1}\sum_{i=1}^{n}\\{e_{ij}(u)-\bar{e}_{j}(u)\\}\\{e_{ik}(v)-\bar{e}_{k}(v)\\}$
$\displaystyle\overset{\triangle}{=}$
$\displaystyle~{}\widehat{\Sigma}_{jk}(u,v)+\widehat{\Sigma}^{Xe}_{jk}(u,v)+\widehat{\Sigma}^{Xe}_{kj}(u,v)+\widehat{\Sigma}^{ee}_{jk}(u,v)$

for any $u,v\in{\cal U}$. By (11),

$\displaystyle n\|\widehat{\Sigma}^{Xe}_{jk}\|_{{\cal S}}^{2}=O_{\rm
p}(n^{1-2a_{1}})=n\|\widehat{\Sigma}^{Xe}_{kj}\|_{{\cal S}}^{2}\,,$

and $n\|\widehat{\Sigma}^{ee}_{jk}\|_{{\cal S}}^{2}=O_{\rm p}(n^{1-2a_{2}})$,
where $a_{1},a_{2}>1/2$. Following from Cauchy–Schwarz inequality, it holds
that $\widetilde{T}_{n,jk}=T_{n,jk}+O_{\rm p}(n^{1/2-a_{1}})+O_{\rm
p}(n^{1/2-a_{2}})=T_{n,jk}+o_{\rm p}(1)$. Under Conditions 1–4 and $H_{0,jk}$,
following the same arguments as the proof of Theorem 1(i), we can derive the
results of Theorem 3(i). $\hfill\Box$
Proof of Theorem 3(ii). Following the similar arguments as the proof of
Theorem 3(i), based on (11) and Cauchy–Schwarz inequality, we have
$\widetilde{T}_{n,jk}=T_{n,jk}+o_{\rm p}(1)$. Then it holds that

$\displaystyle\widetilde{T}_{n,jk}=$
$\displaystyle~{}\iint[\sqrt{n}\\{\widehat{\Sigma}_{jk}(u,v)-\Sigma_{jk}(u,v)\\}]^{2}\,{\rm
d}u{\rm d}v+n\iint\Sigma_{jk}^{2}(u,v)\,{\rm d}u{\rm d}v$
$\displaystyle~{}+2\sqrt{n}\iint\Sigma_{jk}(u,v)[\sqrt{n}\\{\widehat{\Sigma}_{jk}(u,v)-\Sigma_{jk}(u,v)\\}]\,{\rm
d}u{\rm d}v+o_{\rm p}(1)\,,$

Since $\iint\Sigma_{jk}^{2}(u,v)\,{\rm d}u{\rm d}v>0$, then following from
(S.2), Cauchy–Schwarz inequality and $H_{1,jk}$, it holds that
$\iint[\sqrt{n}\\{\widehat{\Sigma}_{jk}(u,v)-\Sigma_{jk}(u,v)\\}]^{2}\,{\rm
d}u{\rm d}v=O_{\rm p}(1)$, $n\iint\Sigma_{jk}^{2}(u,v)\,{\rm d}u{\rm d}v=O(n)$
and
$2\sqrt{n}\iint\Sigma_{jk}(u,v)[\sqrt{n}\\{\widehat{\Sigma}_{jk}(u,v)-\Sigma_{jk}(u,v)\\}]\,{\rm
d}u{\rm d}v=O_{\rm p}(n^{1/2})$, which implies that
${\mathbb{P}}\\{\widetilde{T}_{n,jk}\geq T_{0,jk}(\alpha)\\}\rightarrow 1$ as
$n\rightarrow\infty$ under $H_{1,jk}$. $\hfill\Box$

###  A.4 Proof of Theorem 4

Based on (11) and triangle inequality, we have
$\max_{q\in[Q]}|\widetilde{T}_{n}^{(q)}-\mathring{T}_{n}^{(q)}|=O_{\rm
p}\\{n^{-c}(\log p)^{1/2}\\}=o_{\rm p}\\{(\log p)^{-2}\\}$ for some
sufficiently small constant $c>0$ provided that $p\lesssim n^{\kappa}$ for
some constant $\kappa>0$, where
$\mathring{T}_{n}^{(q)}=\iint\\{n^{-1/2}\sum_{i=1}^{n}Z_{ij_{q}}(u)Z_{ik_{q}}(v)\\}^{2}\,{\rm
d}u{\rm d}v$ for each $q\in[Q]$. Following the same arguments as the proof of
Theorem 2, it is easy to derive the results of Theorem 4. $\hfill\Box$

###  A.5 Proof of Proposition 1

By local linear smoothing method, for fixed $i,j$ and any $u\in{\cal U}$, we
have

$\displaystyle\widehat{X}_{ij}(u)=\frac{R_{ij0}S_{ij2}-R_{ij1}S_{ij1}}{S_{ij0}S_{ij2}-S_{ij1}^{2}}\,,$

where
$R_{ijr}=T_{ij}^{-1}\sum_{t=1}^{T_{ij}}K_{h}(U_{ijt}-u)\\{(U_{ijt}-u)/h_{j}\\}^{r}W_{ijt}$
and
$S_{ijr}=T_{ij}^{-1}\sum_{t=1}^{T_{ij}}K_{h}(U_{ijt}-u)\\{(U_{ijt}-u)/h_{j}\\}^{r}$
for $r=0,1,2$. Here, $K(\cdot)$ is the kernel function and $h_{j}$ is the
bandwidth. For any fixed $i\in[n]$ and $j\in[p]$, we can rewrite
$\widehat{X}_{ij}(u)$ as

$\displaystyle\widehat{X}_{ij}(u)=\sum_{t=1}^{T_{ij}}K_{h_{j}}^{T_{ij}}(U_{ijt}-u)W_{ijt}\,,$
(S.11)

where $K^{T_{ij}}(u)$ is the empirical equivalent kernels for the local
smoothing method, see Fan and Gijbels (1996). For convenience, we define the
following widely-used functionals of a kernel $K$:

$\displaystyle B_{r}(K)=\int K(u)u^{r}\,{\rm d}u\,,~{}~{}~{}~{}~{}~{}V(K)=\int
K^{2}(u)\,{\rm d}u\,,~{}~{}~{}~{}~{}K^{(1)}(u)=\int K(s)K(s+t)\,{\rm d}s\,.$

Denote the set of all the time points as ${\cal
T}=\\{U_{ijt},i\in[n],j\in[p],t\in[T_{ij}]\\}$, and
$\Sigma_{jk}^{l_{1},l_{2}}(u,v)$ as the $(l_{1},l_{2})$-times partial
derivative of $\Sigma_{jk}(u,v)$, that is,
$\Sigma_{jk}^{l_{1},l_{2}}(u,v)=\frac{\partial^{l_{1}+l_{2}}\Sigma_{jk}(u,v)}{\partial^{l_{1}}u\,\partial^{l_{2}}v}$.
Following the similar arguments as the proof of Theorem 1 in Zhang and Chen
(2007), we also have the average MSE (mean squared errors) of the local
smoothing estimator $\widehat{X}_{ij}(u)$, $i\in[n]$ and $j\in[p]$, satisfying

$\displaystyle{\mathbb{E}}[\\{\widehat{X}_{ij}(u)-X_{ij}(u)\\}^{2}\,|\,{\cal
T}]$
$\displaystyle~{}~{}~{}~{}~{}~{}=\bigg{\\{}\frac{B_{2}^{2}(K^{*})[\\{\mu_{j}^{(2)}(u)\\}^{2}+\Sigma_{jj}^{(2,2)}(u,u)]}{4}h_{j}^{4}+\frac{V(K^{*})\sigma_{j}^{2}}{f_{U}(u)}(T_{ij}h_{j})^{-1}\bigg{\\}}\\{1+o_{\rm
p}(1)\\}\,,$

where $K^{*}$ is the equivalent kernel of the local smoothing method (Fan and
Gijbels, 1996, page 64). Since the coefficients of $h_{j}^{4}$ and
$(T_{ij}h_{j})^{-1}$ are uniformly bounded over the finite interval ${\cal
U}$. Then, ${\mathbb{E}}[\\{\widehat{X}_{ij}(u)-X_{ij}(u)\\}^{2}\,|\,{\cal
T}]=O_{\rm UP}\\{h_{j}^{4}+(T_{ij}h_{j})^{-1}\\}$, where $O_{\rm UP}(1)$ means
bounded in probability uniformly for any $t$ within the interior of ${\cal U}$
and all $i\in[n]$ and $j\in[p]$. Therefore, we have

$\displaystyle\widehat{X}_{ij}(u)=X_{ij}(u)+\\{h_{j}^{2}+(T_{ij}h_{j})^{-1/2}\\}O_{\rm
UP}(1)\,,$ (S.12)

which implies that for any $j,k\in[p]$,

$\displaystyle\widehat{\Sigma}_{jk}^{ee}(u,v)=$
$\displaystyle~{}\frac{1}{n-1}\sum_{i=1}^{n}[\\{\widehat{X}_{ij}(u)-\widebar{\widehat{X}}_{j}(u)\\}-\\{X_{ij}(u)-\widebar{{X}}_{j}(u)\\}][\\{\widehat{X}_{ik}(v)-\widebar{\widehat{X}}_{k}(v)\\}-\\{X_{ik}(v)-\widebar{X}_{k}(v)\\}]$
$\displaystyle=$
$\displaystyle~{}\frac{1}{n-1}\sum_{i=1}^{n}\\{\widehat{X}_{ij}(u)-X_{ij}(u)\\}\\{\widehat{X}_{ik}(v)-X_{ik}(v)\\}$
$\displaystyle~{}-\frac{n}{n-1}\bigg{[}\frac{1}{n}\sum_{i=1}^{n}\\{\widehat{X}_{ij}(u)-X_{ij}(u)\\}\bigg{]}\bigg{[}\frac{1}{n}\sum_{i=1}^{n}\\{\widehat{X}_{ik}(v)-X_{ik}(v)\\}\bigg{]}$
$\displaystyle=$
$\displaystyle~{}\frac{1}{n}\sum_{i=1}^{n}\\{h_{j}^{2}+(T_{ij}h_{j})^{-1/2}\\}\\{h_{k}^{2}+(T_{ik}h_{k})^{-1/2}\\}O_{\rm
UP}(1)$
$\displaystyle~{}+\bigg{[}\frac{1}{n}\sum_{i=1}^{n}\\{h_{j}^{2}+(T_{ij}h_{j})^{-1/2}\\}\bigg{]}\bigg{[}\frac{1}{n}\sum_{i=1}^{n}\\{h_{k}^{2}+(T_{ik}h_{k})^{-1/2}\\}\bigg{]}O_{\rm
UP}(1)$ $\displaystyle=$
$\displaystyle~{}\bigg{[}h_{j}^{2}h_{k}^{2}+h_{j}^{-1/2}h_{k}^{-1/2}\bigg{\\{}\frac{1}{n}\sum_{i=1}^{n}T_{ij}^{-1/2}T_{ik}^{-1/2}+\bigg{(}\frac{1}{n}\sum_{i=1}^{n}T_{ij}^{-1/2}\bigg{)}\bigg{(}\frac{1}{n}\sum_{i=1}^{n}T_{ik}^{-1/2}\bigg{)}\bigg{]}O_{\rm
UP}(1)\,.$

Therefore, it holds that

$\displaystyle\|\widehat{\Sigma}_{jk}^{ee}\|_{{\cal S}}^{2}=$
$\displaystyle~{}O_{\rm
p}\bigg{[}h_{j}^{4}h_{k}^{4}+h_{j}^{-1}h_{k}^{-1}\bigg{(}\frac{1}{n}\sum_{i=1}^{n}T_{ij}^{-1/2}T_{ik}^{-1/2}\bigg{)}^{2}+h_{j}^{-1}h_{k}^{-1}\bigg{(}\frac{1}{n}\sum_{i=1}^{n}T_{ij}^{-1/2}\bigg{)}^{2}\bigg{(}\frac{1}{n}\sum_{i=1}^{n}T_{ik}^{-1/2}\bigg{)}^{2}\bigg{]}$
$\displaystyle=$ $\displaystyle~{}O_{\rm
p}\bigg{\\{}h_{j}^{4}h_{k}^{4}+h_{j}^{-1}h_{k}^{-1}\bigg{(}\frac{1}{n}\sum_{i=1}^{n}T_{ij}^{-1}\bigg{)}\bigg{(}\frac{1}{n}\sum_{i=1}^{n}T_{ik}^{-1}\bigg{)}\bigg{\\}}=O_{\rm
p}\big{\\{}h_{j}^{4}h_{k}^{4}+(h_{j}\widebar{T}_{j})^{-1}(h_{k}\widebar{T}_{k})^{-1}\big{\\}}\,,$

where $\widebar{T}_{j}=(n^{-1}\sum_{i=1}^{n}T_{ij}^{-1})^{-1}$ for any
$j\in[p]$.

Recall $W_{ijt}=X_{ij}(U_{ijt})+\varsigma_{ijt}$. By (S.11) and the fact that
$\sum_{t=1}^{T_{ij}}K_{h_{j}}^{T_{ij}}(U_{ijt}-u)=1$, then

$\displaystyle\widehat{X}_{ij}(u)-X_{ij}(u)=\sum_{t=1}^{T_{ij}}K_{h_{j}}^{T_{ij}}(U_{ijt}-u)\varsigma_{ijt}+\sum_{t=1}^{T_{ij}}K_{h_{j}}^{T_{ij}}(U_{ijt}-u)\\{X_{ij}(U_{ijt})-X_{ij}(u)\\}$

for any $i\in[n]$, $j\in[p]$ and $u\in{\cal U}$. Analogous to Lemma 3 in Zhang
and Chen (2007), we also have

$\displaystyle{\text{Cov}}\\{\widehat{X}_{ij}(u)-X_{ij}(u),X_{ik}(v)\,|\,{\cal
T}\\}=$
$\displaystyle~{}\sum_{t=1}^{T_{ij}}K_{h_{j}}^{T_{ij}}(U_{ijt}-u)\\{\Sigma_{jk}(U_{ijt},v)-\Sigma_{jk}(u,v)\\}$
$\displaystyle=$
$\displaystyle~{}\frac{B_{2}(K^{*})\Sigma_{jk}^{(2,0)}(u,v)}{2}h_{j}^{2}\\{1+o_{\rm
p}(1)\\}=O_{\rm UP}(h_{j}^{2})\,.$ (S.13)

Together with the law of large numbers,

$\displaystyle\widehat{\Sigma}_{jk}^{Xe}(u,v)=$
$\displaystyle~{}\frac{1}{{n-1}}\sum_{i=1}^{n}\\{X_{ij}(u)-\widebar{X}_{j}(u)\\}[\\{\widehat{X}_{ik}(v)-\widebar{\widehat{X}}_{k}(v)\\}-\\{X_{ik}(v)-\widebar{X}_{k}(v)\\}]$
$\displaystyle=$
$\displaystyle~{}{\mathbb{E}}\bigg{\\{}\frac{1}{n}\sum_{i=1}^{n}{\mathbb{E}}\big{(}\\{X_{ij}(u)-\widebar{X}_{j}(u)\\}[\\{\widehat{X}_{ik}(v)-\widebar{\widehat{X}}_{k}(v)\\}-\\{X_{ik}(v)-\widebar{X}_{k}(v)\\}]\,|\,{\cal
T}\big{)}\bigg{\\}}O_{\rm p}(1)$ $\displaystyle=$
$\displaystyle~{}{\mathbb{E}}[{\text{Cov}}\\{X_{ij}(u),\widehat{X}_{ik}(v)-X_{ik}(v)\,|\,{\cal
T}\\}]O_{\rm p}(1)$ $\displaystyle=$ $\displaystyle~{}h_{k}^{2}O_{\rm
UP}(1)\,,$

which implies

$\displaystyle\|\widehat{\Sigma}_{jk}^{Xe}\|_{{\cal S}}^{2}=O_{\rm
p}(h_{k}^{4})\,.$

Hence, we complete the proof of Proposition 1. $\hfill\Box$

###  A.6  Proof of Theorem 5

Notice that $\hat{\varepsilon}_{i,jk}(u)=Y_{ij}(u)-\int_{{\cal U}}{\bf
Y}_{i,-j,-k}^{{\mathrm{\scriptscriptstyle
T}}}(v)\widehat{\boldsymbol{\beta}}_{jk}(u,v)\,{\rm d}v$ and
$\bar{\hat{\varepsilon}}_{jk}(u)=n^{-1}\sum_{i=1}^{n}\hat{\varepsilon}_{i,jk}(u)$.
Then

$\displaystyle\hat{\varepsilon}_{i,jk}(u)-\bar{\hat{\varepsilon}}_{jk}(u)=$
$\displaystyle~{}Y_{ij}(u)-\widebar{Y}_{j}(u)-\int_{{\cal U}}\\{{\bf
Y}_{i,-j,-k}(v)-\widebar{{\bf Y}}_{-j,-k}(v)\\}^{{\mathrm{\scriptscriptstyle
T}}}\widehat{\boldsymbol{\beta}}_{j}(u,v)\,{\rm d}v$ $\displaystyle=$
$\displaystyle~{}\int_{{\cal U}}\\{{\bf Y}_{i,-j,-k}(v)-\widebar{{\bf
Y}}_{-j,-k}(v)\\}^{{\mathrm{\scriptscriptstyle
T}}}\\{\boldsymbol{\beta}_{jk}(u,v)-\widehat{\boldsymbol{\beta}}_{jk}(u,v)\\}\,{\rm
d}v+\varepsilon_{i,jk}(u)-\bar{\varepsilon}_{jk}(u)$ $\displaystyle=$
$\displaystyle~{}\int_{{\cal U}}\\{{\bf Y}_{i,-j,-k}(v)-\widebar{{\bf
Y}}_{-j,-k}(v)\\}^{{\mathrm{\scriptscriptstyle
T}}}\\{\boldsymbol{\beta}_{jk}(u,v)-\widehat{\boldsymbol{\beta}}_{jk}(u,v)\\}\,{\rm
d}v+\tilde{\varepsilon}_{i,jk}(u)\,,$

where $\widebar{Y}_{j}(u)=n^{-1}\sum_{i=1}^{n}Y_{ij}(u)$, $\widebar{\bf
Y}_{-j,-k}=n^{-1}\sum_{i=1}^{n}{\bf Y}_{i,-j,-k}$,
$\bar{\varepsilon}_{jk}(u)=n^{-1}\sum_{i=1}^{n}\varepsilon_{i,jk}(u)$ and
$\tilde{\varepsilon}_{i,jk}(u)=\varepsilon_{i,jk}(u)-\bar{\varepsilon}_{jk}(u)$