Profile and Globe Tests of Mean Surfaces for Two-Sample Bivariate Functional Data
Jin Yang Tao Zhang Chunling Liu Kam Chuen Yuen Aiyi Liu




#1 1




0
  

  
  Profile  and Globe Tests of Mean Surfaces for Two-Sample Bivariate Functional Data
    Jin Yang
    The authors gratefully acknowledge support of Post-doctoral Fellowship of Nankai University;

    School of Mathematics and Statistics, Nankai University, Tianjin, 300071, P. R. China
 (<michael.jin.yang@nankai.edu.cn>) 

    and 

    Tao Zhang 
    	The authors gratefully acknowledge support of National Natural Science Foundation of China (Grant No. 11561006,61462008), research projects of colleges and universities in Guangxi (Grant No. KY2015YB171), Scientific Research and Technology Development Project of Liuzhou(Grant No. 2016C05020;

School of Sciences, Guangxi University of Science and Technology, Liuzhou, 545006, P. R. China
 (<tzhangmaths@126.com>) 

        and 

    Chunling Liu 
	The authors gratefully acknowledge support of The Hong Kong Polytechnic University funding 4-BCC0;

Department of Applied Mathematics, The Hong Kong Polytechnic University, Hong Kong
 (<macliu@polyu.edu.hk>) 

and 

    Kam Chuen Yuen 
    	The authors gratefully acknowledge support of The University of Hong Kong Grant No. HKU17329216;

Department of Statistics and Actuarial Science, The University of Hong Kong, Hong Kong
 (<kcyuen@hku.hk>) 

            and 

    Aiyi Liu 
    	The authors gratefully acknowledge support of Eunice Kennedy Shriver National Institute of Child Health and Human Development, National Institutes of Health, USA.

Biostatistics and Bioinformatics Branch, NICHD, NIH, Bethesda, Maryland, USA.
 (<liua@mail.nih.gov>) 


    December 30, 2023
=============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================

 

1
  

    Profile  and Globe Tests of Mean Surfaces for Two-Sample Bivariate Functional Data


   
Multivariate functional data has received considerable attention but testing for equality of mean surfaces and its profile has limited progress. The existing literature has tested equality of either mean curves of univariate functional samples directly, or mean surfaces  of bivariate functional data samples but turn into functional curves comparison again. In this paper, we aim to develop both the profile and globe tests of mean  surfaces  for two-sample bivariate functional data.
We present valid approaches of tests by employing the idea of pooled projection and by developing a novel profile functional principal component analysis tool.
The proposed methodology enjoys the merit of readily interpretability and implementation. Under mild conditions, we derive the asymptotic behaviors of test statistics under null and alternative hypotheses.
Simulations show that the proposed tests have a good control of the type I error by the size and can detect difference in mean surfaces  and its profile effectively in terms of power in finite samples.
Finally, we apply the testing procedures to two real data sets associated with the precipitation change affected jointly by time and locations in the Midwest of USA, and the trends in human mortality from European period life tables.



Keywords:  
Asymptotic Chi-square; Bivariate functional data; Globe test; Mean surface; Profile test.


1.45 


§ INTRODUCTION

In multivariate functional stochastic process X(u), there has increasing research interest in data type that is both functional and multidimentional. That is, u=(s, t) has two arguments where s ∈𝒮⊂ℝ^d_1 and t ∈𝒯⊂ℝ^d_2 with d_1 and d_2 being positive integers. Here s and t inherently belong to distinct domains 𝒮 and 𝒯 in terms of scientific meaning or research design. For example, X(s,t) may be the mortality rate  of age s during year t in a given country.  A typical example of such data comes from neuroimaging studies using functional magnetic resonance imaging (fMRI), in which the so-called voxels data, i.e. brain activity like blood flow changes are discrepantly recorded at a large number of locations at irregular time units <cit.>. Spatiotemporal study is no doubt another important application of this kind of data where  t is defined on a temporal domain and  s is defined on a spatial domain. Although functional data of afore structure are encountered in many applications,  there is rare progress in inferential aspect for such data <cit.>.  In the present work, we plan to investigate the profile and globe   tests  of  mean surfaces for  two  bivariate functional samples.

A practical motivation for this research comes from precipitation  data in Midwest  of the United States,  where  the daily data  of precipitation  from 1941 to 2000 are  collected at 59 spatial locations scattered over 12 states in the Midwest of USA.  For ease of reference, we provide a map of Midwest states with the locations of the climate monitoring stations in Figure <ref>.
The Midwest  is a breadbasket of the United States and its agriculture has continued to play a major role in the economy of the region  <cit.>. The agriculture in the Midwest is vulnerably affected by the climate, of which precipitation is a vital component. To monitoring the future agricultural activities, it therefore has long been  recognized as an important problem to reveal how the change of precipitation takes place for different locations, different regions, or different years in the same region.







The study of the  precipitation  data has led to several interesting findings.  For instance, <cit.> detected no changes during the period 1941-2000 for only  individual station. However, it is difficult to implement if we sequently tested for every station when the number of stations were large.
<cit.> used cumulative sum paradigm to expose the fact that, the mean precipitation curves before and after 1966 were different over the whole region. Nevertheless their method was particularly designed to detect the temporal change
but not applicable to detect the difference in spatio domain, not to mention
the joint spatiotemporal effect on the precipitation.
Looking into analysis of heatmaps of yearly sample mean surfaces where X_i(s,t), i=1941,⋯, 1967,  corresponds to the  precipitation of the tth  day in the ith year
at the sth station,
intuitively we have observed that the yearly sample mean surface of precipitation in the Great Plains is different from that in the Great Lakes, refer to Figure <ref>. Also, we can recognize from Figure <ref> that some  profiles of mean surface are same but others are different.
These motivate us to develop more powerful inferential procedures to detect if mean surfaces or its profiles have significant difference for either different regions or different individual stations.






Tracking back testing procedures for the equality of mean functions in the functional data setting, existing works mainly focus on detecting the curve equality for univariate functional data. 
In the two-sample testing scenario,
<cit.> presented bootstrap procedures for testing the equality of mean curves through the eigenelements for two independent functional samples.
Under the Gaussian assumption,
<cit.> considered the two-sample test based on L^2-norm.
<cit.> derived mean functions comparison through a normal approximation method but only applicable to dense functional data samples.  
<cit.> still solved testing the curve equality problem though in
bivariate (two-dimensional by their words) functional data setting and for distribution function testing.
Regarding the k-sample testing or the  one-way ANOVA for functional data, works include HANOVA <cit.>, Cramér-von Mises type test <cit.>, F-type test <cit.>,  B-spline test <cit.>, and Mahalanobis distance <cit.>, among others.
In the case of within-curve dependence in each sample,
<cit.> detected the mean curve variation using  L^2-norm criterion. <cit.> and its  multiple group extension  <cit.> worked on parametric testing relying on quite strong assumptions. Notice that, throughout our literature review, since our awareness concentrates on testing the equality of mean functions, we leave out other inferential topics such as testing the equality of coefficient operators or testing independency within a sample, and etc.

It has series of work in functional time series literature on testing the equality of mean functions, where weak dependence between or within two samples are accommodated in reality. Testing mean function difference in such functional time series study  had still been on comparison of mean curve functions (, among others).

Aforementioned literature in both functional curve samples and functional time series have all inclined to testing the equality of mean curve functions, i.e. the inferential target is on univariate functional data.
However, for comparison between samples of multivariate functional data, there have been few works by far. Only <cit.> raised testing the equality of the mean surfaces of bivariate functional data, but eventually the equality of mean curves indexed at all locations were tested. Also to the best of our knowledge, the profile test of mean surfaces has not been considered for two bivariate functional data samples.
Although the profile test of mean surfaces may belong to the curve test scope, it attributes to two different topics due to the different subjects. 
Above dire need in real-world data analysis and literature review motivates us to develop valid tests for equality of means surfaces and the corresponding profile test for  bivariate functional data samples.

To address the problem in demand, 
firstly, we obtain the marginal eigen-function of the pooled sample by marginal functional principal component analysis (FPCA) and project the profiles of mean surfaces on marginal eigenfunctions.  
The profile testing statistic measures the distance of the profile of mean surfaces for two bivariate functional samples. Once the marginal eigenfunctions are obtained,  
the eigensurfaces of the pooled sample can be constructed by further FPCA. 
The distance between mean surfaces for two samples can be measured by the globe test statistic using the analogous projection ideas.       
Consequently, our proposed profile testing procedures can be implemented for every profile of the mean surface, which corresponds to simultaneously test whether mean precipitation curves have significant difference for every station.  
The globe test performs well in terms of both the size and the power in that it includes the information of two domains effectively.

The major contribution of this paper is threefold.
Firstly, the presented methodology may be the first one to detect difference of mean surfaces and its profile for two-sample bivariate functional data.
In contrast to the literature that we can search out by far, of which the focus has almost all been on testing the equality of mean curves as a matter of fact. When one argument is fixed, our profile test methodology can also simultaneously detect the mean difference in the other domain. 
Secondly, our testing procedures are interpretable and easily implemented. This will help fill out some theoretical gaps in functional inference and facilitate the real application and interpretation in statistical perspective.
Finally, asymptotic distributions of the test statistics under null hypotheses has been derived. The consistency of  test procedure has been proved. 
In addition, simulation studies show that the proposed tests have a good control of the type I error by the size and can detect difference in mean surfaces and its profile effectively in terms of power in finite samples.




The rest of the paper is organized as follows.
In Section <ref>, we describe the model and data structure.
The profile test procedure of mean surfaces for two bivariate functional data samples is presented in Section <ref>,
while globe test procedure is proposed in Section <ref>.
The finite sample performance for several representative scenarios is investigated in Section <ref>. In Section <ref>, we demonstrate two applications associated with the precipitation changes affected jointly by time and locations in the
Midwest of USA, and the trends in human mortality from European period life tables. The paper concludes with
a brief discussion in Section <ref>.   
Theory proofs are included in AppSupplementary. 



§ MODEL AND DATA STRUCTURE

Let L^2(𝒮×𝒯) be the separable
Hilbert space.
{X^(m)(s,t):(s,t)∈𝒮×𝒯} is a square integrable  stochastic  process    on L^2(𝒮×𝒯 )  with
mean function μ_m(s,t)=E{X^(m)(s,t)} and covariance function

    C^(m){(s,v),(u,t)}= E{X^(m)c(s,v)X^(m)c(u,t)},





where X^(m)c(s,t)=X^(m)(s,t)-μ_m(s,t), for m=1,2, respectively.
With this notation, we
can decompose X^(m)(s,t) into

    [ X^(m)(s, t)=μ_m(s, t)+ε^(m)(s,t),
    m=1,2, ]

where ε^(m)(s,t)  is the stochastic part of X^(m)(s, t) with
E{ε^(m)(s,t)}=0   and covariance function  C^(m){(s,v),(u,t)}.


Functional samples {X_i^(m)(s, t),  m=1,2;  i=1,⋯,n_m } may usually be modeled as independent realizations of the underlying  stochastic process X^(m)(s,t).
In practice,  {X_i^(m)(s, t),  m=1,2;  i=1,⋯,n_m }  can not be observed, but rather,  measurements are taken at discrete  time points.   In this paper,   we assume  {X_i^(m)(s, t),  m=1,2;  i=1,⋯,n_m }
are recorded on a regular and dense grid of  time points as follows,

    [ X_i^(m)(s_il_1, t_il_2)=μ_m(s_il_1, t_il_2)+ε_i^(m)(s_il_1,t_il_2);;                             m=1,2; i=1,⋯,n_m; l_1=1,⋯,N; l_2=1,⋯,M. ]


In this paper, we are  firstly interested in profile test of bivariate functional data samples,  i.e.  for every fixed t^*∈𝒯,

    H_0^𝒮: μ_1(s,t^*)=μ_2(s,t^*) vs.
    		H_1^𝒮: μ_1(s,t^*)≠μ_2(s,t^*), s ∈𝒮,

or for every fixed s^*∈𝒮,

    H_0^𝒯: μ_1(s^*,t)=μ_2(s^*,t) vs.
    		H_1^𝒯: μ_1(s^*,t)≠μ_2(s^*,t), t ∈𝒯.



Then we go to the second target to present a globe test procedure for bivariate functional data samples  with hypothesis  below,

    H_0: μ_1(s,t)=μ_2(s,t) vs.
    		H_1: μ_1(s,t)≠μ_2(s,t), s ∈𝒮,  t ∈𝒯.


The equality in hypothesis (<ref>) means that ∫_𝒮{μ_1(s,t^*)-μ_2(s,t^*)}^2 ds =0   for every fixed t^*∈𝒯, and the alternative means that ∫_𝒮{μ_1(s,t^*)-μ_2(s,t^*)}^2 ds >0.   Analogously meaning can be interpreted for (<ref>). However,  null hypothesis of (<ref>) implies ∫_𝒮∫_𝒯{μ_1(s,t)-μ_2(s,t)}^2dtds=0 while the alternative means that ∫_𝒮∫_𝒯{μ_1(s,t)-μ_2(s,t)}^2dtds>0.
For statistical inference of bivariate functional data,  marginal FPCA is a  widely used tool, which often assumes that  bivariate functional data can project onto  finite-dimensional eigensurfaces <cit.>. It is our start point for the proposed profile and globe test procedures.



§ PROFILE TEST OF BIVARIATE FUNCTIONAL DATA

Profile test of bivariate functional data is an important problem, as it allows to
provide multiple insight from multiple angles, and also is of interest in many
applications. For example, in analysis of precipitation, the testing problem (<ref>) corresponds to test whether mean precipitation curves have significant difference  before and after 1966 for every station, while the testing problem (<ref>) means to test  whether different stations have significant difference for every day. <cit.> considered detection the difference only on an individual station. However, it is difficult to implement when the number of stations is large if we sequentially test for every station by their method. 
So, we propose the profile test of mean functions which is easy to implement and can  simultaneously detect difference of all stations. In this section, we address the test problem (<ref>) only as (<ref>) can be analogously implemented.

As a first step, the marginal covariance function is denoted to be G_𝒮^(m)(s,u)=∫_𝒯C^(m){(s,t),(u,t)}dt, as the form of (5) in <cit.>, and may be estimated by

    G_𝒮^(m)(s_h,s_l)=1/n_mM∑_i=1^n_m∑_k_2=1^M   X_i^(m)c(s_h,t_ik_2)
       ×X_i^(m)c(s_l,t_ik_2),

where X_i^(m)c(s, t)=X_i^(m)(s,t)- X^(m)(s,t)
with 
 X^(m)(s,t)=1/n_m∑_i=1^n_mX_i^(m)(s,t).

Denote

    [ G_𝒮(s,u)=n_2/n_1+n_2G_𝒮^(1)(s,u)
    +n_1/n_1+n_2G_𝒮^(2)(s,u), s,u∈𝒮. ]


It is easy to see  G_𝒮(s,u) p⟶ (1-θ)G_𝒮^(1)(s,u)+θ G_𝒮^(2)(s,u)≡ G_𝒮(s,u), where  θ is defined in Assumption 6 stated in next section and  G_𝒮(s,u) is the pooled covariance function.
Consequently, it has orthogonal eigenfunctions {ψ_j}_j≥ 1 and
non-negative eigenvalues {ν_j}_j≥ 1 satisfying

    [                                                   ∫_𝒮^ G_𝒮(s,u)ψ_j(u)du=ν_jψ_j(s),  s,u∈𝒮, j=1,2,…. ]


Such eigencomponents can be numerically estimated by suitably discretized eigenequations,

    ∫_𝒮G_𝒮(s,u)ψ_j(u)du=
    	ν_jψ_j(s), j=1,2,…,

with orthogonal constraints on {ψ_j}_j≥ 1.

Once the estimators of marginal eigen-functions ψ_j(s), j=1,2,…,  are obtained, we project the observations onto the marginal eigenfunctions and obtain the
profile estimators of mean functions as follows:  for every fixed t^*∈𝒯,


    μ_m(s,t^*)=∑_j=1^Jη_j^(m)(t^*)
    	ψ_j(s), m=1,2,

with

    η_j^(m)(t^*)=1/n_m∑_i=1^n_mη_ij^(m)(t^*),
       η_ij^(m)(t^*)=1/N∑_l_1=1^NX_i^(m)(s_il_1,t^*)
    		ψ_j(s_il_1).

For practical implementation, one has to decide the magnitude of J.
A practical strategy is J=min{j:ν_1+ν_2+⋯
	+ν_k/ν_l+ν_2+⋯>q},  where  ν_l, l=1,2,⋯  are defined in (<ref>).
We find that  q=90% threshold works well for our numerical examples.

Based on  above discussion, we propose the following profile test statistic

    TP(t^*)=n_1n_2/n_1+n_2∑_j=1^J(η_j^(1)(t^*)- η_j^(2)(t^*))^2/λ_j(t^*),

where
λ_j(t^*)=n_2/n_1+n_2λ_j^(1)(t^*)
+n_1/n_1+n_2λ_j^(2)(t^*) with

λ_j^(m)(t^*)=n_m^-1∑_i=1^n_m{η_ij^(m)(t^*)-η_j^(m)(t^*)}^2, m=1,2.


	It is easy to see that n_1n_2/n_1+n_2∫[μ_1(s,t^*) 
 - μ_2(s,t^*)]^2dt p⟶ U_n_1, n_2= n_1n_2/n_1+n_2∑_l=1^ K(η_l(t^*)- η_2(t^*))^2. However,  the variance of U_n_1, n_2 may be unnecessarily inflated by the
	presence of, possibly many, very small estimates μ_1(s,t^*)- μ_2(s,t^*).  This drawback can be remedied by giving a divisor to λ_j(t^*).



We then establish asymptotic behaviors of the test statistic TP(t^*) under the null hypothesis H_0^𝒮 and the alternative one H_1^𝒮. To derive the asymptotic properties of profile test statistic, we make the following assumptions.


	ν_1>ν_2>⋯
	where {ν_j}_j=1,2,… are the eigenvalues of covariance operates  G_S(s,u).



	For every fixed t^*,   μ_m(s,t^*), m=1,2 may be written as μ_m(s,t^*)=∑_j=1^∞η_j^(m)(t^*)ψ_j(s),  where η_j^(m)(t^*)=∫_0^1μ_m(s,t^*)ψ_j(s)ds.



	Assume sup _(s,t)∈𝒮×𝒯μ_m^2(s,t), m=1,2   are bounded
	and E(sup|ε^m(s,t)|^4), m=1,2 are bounded.



	The grid point {t_il_1:l_1=1,…,N}  and
	{s_il_2: l_2=1,…,M} are equidistant. We assume
	n_1/N^2=o(1),  n_1/M^2=o(1),  n_2/N^2=o(1)   and n_2/M^2=o(1).



	min{n_1,n_2}→∞, n_1/(n_1+n_2)→θ  for a fixed constant θ∈ (0,1).


Assumptions <ref> and <ref> are regular conditions. One needs these conditions to uniquely (up to signs) choose ψ_j(s) and obtain  the bound of ψ_j(s)-ψ_j(s).  Assumption <ref> means that the profiles of mean surface are projected onto a space that is generated by a large set of basis functions.  
Assumption <ref> requires that functional data are recorded on dense grid.
Assumption <ref> is of standard for two-sample asymptotic inference.


	Under Assumptions <ref>-<ref> and H_0^𝒮, we have TP(t^*)d⟶χ_J^2,
	where χ_J^2 stands for a χ^2 -distributed random variable with J degrees of freedom. Under H_1^𝒮 and 0<θ<1, we have TP(t^*)p⟶∞.


From the expression of TP(t^*) and remark <ref>, we can see that TP(t^*) depends on sample sizes n_1, n_2, and η_j^(1)(t^*)- η_j^(2)(t^*), j=1,⋯,J,  which
reflects  the difference of profile mean functions μ_1(s,t^*) and μ_2(s,t^*).
Intuitively, √(n_1n_2/n_1+n_2)(η_j^(1)(t^*)- η_j^(2)(t^*))λ_j^-1/2  has a limiting standard  normal distribution under  H_0^𝒮.
Theorem <ref> shows that TP(t^*) asymptotically follows the chi-square distribution with J degrees of freedom if H_0^𝒮 holds. Furthermore, TP(t^*) is consistent under H_1^𝒮. The proof of this theorem is provided in AppSupplementary.



§ GLOBE TEST OF BIVARIATE FUNCTIONAL DATA

Compared with the profile test, the globe test of bivariate functional data attempts to detect the joint effects impacted by both domains. In this section, we develop a globe test method for bivariate functional data which aims to detect whether mean surfaces of precipitation have significant difference over a specific time window and/or a specific area, or whether two regions exist significant difference during different time windows.

Based on the estimated  marginal eigenfunctions ψ_j(s)  in Section <ref>, we next estimate the
marginal  functional  principal component scores
ξ^(m)_j,i(t).
The traditional integral estimates of
ξ^(m)_j,i(t)  based on the definition

    ξ^(m)_j,i(t)=∫_𝒮 X_i^(m)c(s,t)ψ_j(s)ds,
       i=1,…,n_m; j=1,2,….

are

    ξ^(m)_j,i(t)=∑_l=2^N X_i^(m)c(s_l,t)ψ_j(s_l)(s_l-s_l-1),
       i=1,…,n_m; j=1,2,….

where  N  is the number of measurements for  X_i^(m)c(s,t)  in the direction 𝒮.

Notice that each score function ξ^(m)_j,i(t) is a centered new random curve. Denote the covariance function of ξ_j,i^(m)(t) by G_𝒯,j^(m)(v,t)=E{ξ_j,i^(m)(v)ξ_j,i^(m)(t)}. Then,  the  estimator of G_𝒯,j^(m) is denoted as,

    G^(m)_𝒯,j(t_h,t_l)=1/n_m∑_i=1^n_mξ_j,i^(m)(t_h)ξ_j,i^(m)(t_l),
       t_h, t_l∈𝒯;  j=1,2,….

Let

    G_𝒯,j(v,t)=n_2/n_1+n_2G^(1)_𝒯,j(v,t)
    		+n_1/n_1+n_2G^(2)_𝒯,j(v,t),
       v,t∈𝒯; j=1,2,….

It is easy to see G_𝒯,j(v,t)p⟶
(1-θ)G^(1)_𝒯,j(v,t)+θ G^(2)_𝒯,j(v,t)
 ≡ G_𝒯,j(v,t)
where
G_𝒯,j(v,t) is the covariance function   and  has orthogonal eigenfunctions {ϕ_jk}_k≥ 1 and
non-negative eigenvalues {ν_jk}_k≥ 1 satisfying

    ∫_𝒯  G_𝒯,j(v,t)ϕ_jk(v)dv=ν_jkϕ_jk(t),
       v,t∈𝒯; k, j=1,2,….


Then estimators of eigenvalues and eigenfunctions {(ν_jk, ϕ_jk(t)): j,k≥ 1} are obtained by the following equations,

    ∫_𝒯G_𝒯,j(v,t)ϕ_jk(v)dv=
    	ν_jkϕ_jk(t), k, j=1,2,…,

with orthogonal constraints on {ϕ_jk}_k≥ 1.

Denote φ_jk(s,t)≡ϕ_jk(t)ψ_j(s) and its consistent estimator by φ_jk(s,t)=ϕ_jk(t)ψ_j(s). We propose estimators of the mean surfaces which
are projection of observations onto a hyperspace spanned from the pooled eigensurfaces {φ_jk(s,t): j,k≥ 1}, written as

    μ_m(s,t)=∑_j=1^J∑_k=1^K_jη_jk^(m)φ_jk(s,t), m=1,2,

with

    η_jk^(m)=1/n_m∑_i=1^n_mη_ijk^(m),
       η_ijk^(m)=1/MN∑_l_2=1^M∑_l_1=1^NX_i^(m)(s_il_1,t_il_2)φ_jk(s_il_1,t_il_2),

where selection of J is the same to in Section <ref> and K_j can be decided by analogous procedure. In details, we select  K_j=min{k:ν_j1+ν_j2+⋯
	+ν_jk/ν_j1+ν_j2+⋯>0.9}, where ν_jl, l=1,2,⋯ are defined in (<ref>).
























It is natural to take into consideration the term TC≡∫_𝒮∫_𝒯{μ_1(s,t)-μ_2(s,t)}^2dtds to measure the distance between two estimated mean surfaces.

It is readily seen that  TCp⟶∑_j=1^J∑_k=1^K_j(η_jk^(1)- η_jk^(2))^2.
Therefore, H_0  will be rejected if  TC is large.  Similarly,  the variance of TC may be unnecessarily inflated by the
presence of, possibly many, very small estimates η_jk^(1)- η_jk^(2). This drawback can also be remedied by giving a divisor to their variance.







Based on the above steps, we propose the following test statistic

    [ TM=n_1n_2/n_1+n_2∑_j=1^J∑_k=1^K_j(η_jk^(1)- η_jk^(2))^2/λ_jk, ]

where
λ_jk=n_2(n_1+n_2)^-1λ_jk^(1)
+n_1(n_1+n_2)^-1λ_jk^(2)  with
λ_jk^(m)=(n_m-1)^-1∑_i=1^n_m(η_ijk^(m)-η_jk^(m))^2, m=1,2.


From (<ref>), we can see that X_i^(1)(·,·) and  X_i^(2)(·,·)
are directly projected on  the common basis surface and obtain η_ijk^(1) and η_ijk^(2).  η_jk^(1) and η_jk^(2), which are the average of such projection, and hence can be viewed as the scores of projection that two  mean surfaces μ_1(s,t)   and μ_2(s,t)  project on the same basis function space, respectively.
The representation of TM measures the total such deviation between two samples. Therefore, the proposed method has a nice explanation and easy to implement.


Next we establish asymptotic behavior of the test statistic TM under hypotheses (<ref>). Additionally, we need the following assumptions.


	ν_j1>ν_j2>⋯ where {ν_jk}_k=1,2,…; j=1,2,… are the eigenvalues of the covariance function G_T(v,t).



	Assume μ_m(s,t), m=1,2  may be written as μ_m(s,t)=∑_j=1^∞∑_k=1^∞η_jk^(m)φ_jk(s,t),  where η_jk^(m)=∫_0^1∫_0^1μ_m(s,t)φ_jk(s,t)dsdt.


Assumption <ref> along with Assumption <ref> in Section <ref> ensures the bound of ϕ_jk(t)-ϕ_jk(t). The interpretation of Assumption <ref> is similar to Assumption <ref> in Section <ref>. 


	Under Assumptions <ref>-<ref> and H_0,  we have
	
    [ TMd⟶χ_∑_j=1^JK_j^2, ]

	where χ_∑_j=1^JK_j^2 stands for a χ^2-distributed random variable with ∑_j=1^JK_j degrees of freedom. Under H_1 and 0<θ<1, we have TMp⟶∞.


Intuitively √(n_1n_2/n_1+n_2)(η_jk^(1)- η_jk^(2))λ_jk^-1/2  has a limiting standard  normal distribution under H_0. Theorem <ref> shows that TM asymptotically follows the chi-square distribution with ∑_j=1^JK_j degrees of freedom under H_0.
The consistency of  TM  is also illustrated under H_1, which
together provides clear theoretical justification of the empirical properties of the proposed test.
The proof of this theorem is provided in AppSupplementary.


§ SIMULATION STUDIES

We conduct extensive simulation studies and report two representative examples here. Examples <ref> and <ref> evaluate two proposed testing procedures in terms of empirical size and power when covariance functions of two samples are identical or distinct, separately. The data grid for argument s consists of 100 equispaced points on [0, 1], and the grids for argument s consists of 50 equispaced points on [0,1]. Each pair of data-generated processes was replicated 1000 times.

	Identical covariance functions.

In this example, we consider the following model

    X_i^(1)(s, t)=ε_i^(1)(s, t), i=1,…,n_1,
       X_i^(2)(s, t)=δ (s+t)+ ε_i^(2)(s, t), i=1,…,n_2,

where ε_i^(1)(s, t)  and  ε_i^(2)(s, t) are independently generated from

    ε(s, t)=∑_j=1^2ξ_j(t)ψ_j(s), s∈ [0,1],t∈ [0,1],

with ψ_1(s)=s^2 and ψ_2(s)=s^3, s∈ [0,1].
ξ_j(t) is generated from

    [ ξ_j(t)=∑_k=1^2χ_jkϕ_jk(t), j=1,2, ]

with
ϕ_11(t)=ϕ_21(t)=-√2cos(2π t), ϕ_12(t)=ϕ_22(t)=√2sin(2π t), t∈ [0,1];
χ_11∼ N(0,3), χ_12∼ N(0,1.5), χ_21∼ N(0,2),  and  χ_22∼ N(0,1).



	Distinct covariance functions.

To compare with Example <ref>, we consider the following model

    X_i^(1)(s, t)=ε_i^(1)(s, t), i=1,…,n_1,
       X_i^(2)(s, t)=δ (s+t)+ ε_i^(2)(s, t), i=1,…,n_2,

where ε_i^(1)(s, t) is generated from

    [ ε^(1)(s, t)=∑_j=1^2ξ_j(t)ψ_j(s), s∈ [0,1],t∈ [0,1], ]

and  ε_i^(2)(s, t) from

    [ ε^(2)(s, t)=ξ_1(t)ψ_1(s), s∈ [0,1],t∈ [0,1], ]

with ψ_1(s)=s^2 and ψ_2(s)=s^3, s∈ [0,1].
ξ_j(t) is generated from

    [ ξ_j(t)=∑_k=1^2χ_jkϕ_jk(t), j=1,2, ]

with
ϕ_11(t)=√2cos(2π t),
ϕ_21(t)=√2sin(2π t),
ϕ_12(t)=2cos(4π t), ϕ_22(t)=2sin(4π t), t∈ [0,1];
χ_11∼ N(0,3), χ_12∼ N(0,1.5), χ_21∼ N(0,2),  and  χ_22∼ N(0,1).

Example <ref> can be seen as two-sample tests where covariance functions are identical, while covariance functions of Example <ref> are distinct. 
The sample size pair is taken to be (n_1,n_2)=(25,75), (50,150), (100,300), (50,50), (100,100), and (200,200), respectively.
The empirical sizes of profile test are computed for different s   and t.  To save space, we here only present the results of different s  for (n_1,n_2)= (100,100) in Figure <ref>. Next, we can also compute the empirical sizes of the globe test. The results are reported in Table <ref>.
The empirical power can be evaluated  when δ≠ 0. The empirical
power at δ= 0.4, 0.6, 0.8  of profile tests are displayed in Figure <ref> while the results of globe tests at δ=0.2, 0.4, 0.6, 0.8, 1.0, 1.2 are scatter plotted in Figure <ref>.



















Several observations can be concluded from Figures <ref> and <ref>. Firstly, the profile tests have a good control of the type I error. The empirical sizes of identical covariance scenarios are better than that of distinct covariance cases.
Secondly, the empirical power of the test becomes larger  when δ increases from 0.4 to 0.8,  which is expected. Lastly, the empirical power for the same covariance case is slightly
larger than that of  the  different covariance function cases.





















We may observe from Table <ref> and Figure <ref> that the globe test approach can keep steady empirical size even at pairs of small sample sizes (n_1, n_2)= (25,75) or (50, 50).
The empirical power of two test methods increases as the sample size increases.
When δ increases from 0.2 to 1.2, the empirical power of the test becomes more and more large, which is evidence of the consistency of the testing procedures. Also the empirical power of equal sample size scenario is slightly better than that of  unequal sample size one.


































§ REAL DATA EXAMPLES

To illustrate profile and globe tests methods, we analyse the historical precipitation data in the Midwest of USA  and the period lifetables in Europe for human mortality trend analysis.




 §.§ Precipitation data

The first example is used to analyze the changes of precipitation during 1941-2000 or in different regions
in the Midwest of USA. <cit.> detected no changes during the period 1941-2000 for only one station  while <cit.> detected the change of precipitation during 1941-2000 over the whole region.

The precipitation data is available from the global historical climatological network database. The comprehensive U.S. Climate Normals dataset includes various derived products including daily air temperature normals, precipitation normals and hourly normals.
The dataset that we analyzed in this paper can be downloaded directly from GHCN (Global Historical Climatology Network)-Daily, an integrated public database of NOAA <(https://www.ncdc.noaa.gov/oa/climate/ghcn-daily/)> by an R interface. Our interest is daily precipitation records from Midwestern states  including  Illinois, Indiana, Iowa, Kansas, Michigan, Minnesota, Missouri, Nebraska, North Dakota, Ohio, South Dakota, and Wisconsin.
In Figure <ref>, totally 59 locations of the climate monitoring stations are indicated with blue circles ∘ in 4 states from the Great Plains (light green region), and  with red triangles  in 5 states from the Great Lakes (yellow region). Notice that there is no climate monitoring stations in Iowa, Michigan, and Missouri. 
We target to detect whether the changes of average precipitation took place for different time phases or regions.

Let Y_i(s,t) be the precipitation of the tth day in the ith year
of the sth station.
Before we apply the proposed method, we need to do registration with the data.  To remove the effects due to the heavy tail distribution, we apply
the transformation

    Z_i(s,t)=log_10{Y_i(s,t)+1},

where {Y_i(s,t)} are original records.
After the transformation, we
pre-smooth data by using the cubic splines function. It is noted that  the data of every climate
monitoring stations from 1941 to 2000 can be constituted into a time series with length
21900(365 day by 60 year). Then, the data of the 59 climate monitoring stations can be
seen as a sample with sample size being 21900 and variables being 59. According to the empirical
Pearson correlation of 59 variables, the 59 climate monitoring stations is stringed into a
function
by the stringing method in <cit.>. Consequently the spatiotemporal data {Y_i(s,t)}  are converted into the bivariate functional data  {X_i(s,t)}. Notice that the difference between the spatiotemporal data Y_i(s,t) and the bivariate functional data X_i(s,t)  is that the argument s in the former expression has no order but it is ranked in the latter.

<cit.> studied the data Y_i(s,t) and detected out the change of  the average precipitation at about  1967.   In this subsection,
we firstly apply the profile test to check if the profile of mean surfaces  are equal during the periods 1941-1967 and 1968-2000. It   corresponds to  test  whether  the average precipitation of every
station has changes during these periods. The p-values  of  the profile  tests  are computed  and  results are displayed
in Figure <ref>.  As can be seen from Figure <ref>, most of the  p-values are less than 0.05 or significant except 11  stations.
For ease of reference,  we list  the latitude and longitude
in Table <ref> for 11 stations. This displays that  the average precipitation  of most locations had changed during  the periods 1941-1967 and 1968-2000.














Next,  we implement the following globe test

    H_0^Midwest: μ_1(s,t)=μ_2(s,t)
       vs. H_1^Midwest: μ_1(s,t)≠μ_2(s,t), s∈ℝ^59, t ∈ℝ^365.

From the globe  test procedure presented in Section <ref> together with the asymptotic distribution
of the test statistic   TM, we calculate the corresponding p-value to be  0.001.  This result is consistent  with the conclusion of Gromenko, Kokoszka and Reimherr. That is, the patterns of mean surfaces
are different over the whole Midwest region between  before 1967 and after 1967.  Intuitively,
according to the results of the profile test, the precipitation had changed  in  most of locations  which lead to  the variations  of whole region.









The heatmaps in Figure <ref> leak the information that sample mean values of annual precipitation in the Great Lakes (GL) based on 28 stations are more than that in the Great Plains (GP). This motivates us to further explore how the mean functions of bivariate functional data {X_i(s,t)} was affected  by temporal and spatial effects from both domains. It is natural to test the equality of two mean surfaces of the precipitation for the 31 stations located in the GP and the 28 stations located in the GL during the periods 1941-1967 and 1968-2000, respectively by

    H_0^1967-: μ^GP=μ^GLvs.
    H_1^1967-: μ^GP≠μ^GL,

and

    H_0^1967+: μ^GP=μ^GLvs.
    H_1^1967+: μ^GP≠μ^GL.

All the p-values by globe test procedures for above two hypotheses are tiny approaching to zero indicating rejecting the null hypotheses but in favor of the alternative one. It is consistent with the intuition that the mean patterns of precipitation at Great Plains and at Great Lakes are different.

Furthermore,
for the 28 stations located in the GL, we test the mean surfaces of precipitation before and after 1967,  denoted by

    H_0^GL: μ^1967-=μ^1967+vs.
    H_1^GL: μ^1967-≠μ^1967+.

The p-value is 0.0163. The null hypothesis would be rejected at 0.05 significance level.
Testing equality of the mean surfaces of precipitation before and after   1967  is also implemented for the 31 stations located in the GP,
denoted  by

    H_0 ^GP: μ^1967-=μ^1967+vs.
    H_1^GP: μ^1967-≠μ^1967+.

The p-values by globe testing method are 0.5677. The null hypothesis would not be rejected at 0.05 significance level.
That is, averagely speaking, the precipitation in the Great Lakes changed before 1967 and after 1967, whereas the mean pattern of precipitation in the Great Plains had no change before 1967 and after 1967. Therefore, our analysis provides evidence that change in the mean function of precipitation was mainly due to the Great Lakes but the Great Plains may be affected little.  By looking up the map,  we find that all the stations in Table <ref>  are located in the Great Plains. It further verify the reliability of the proposed methods. All testing results are presented in Table <ref>.











 §.§ European human mortality rate data

In the second example, we will analyse the trends in human mortality based on the records in the period life tables during the calendar years 1960-2006 for Europe countries. A period life table represents the mortality conditions at a specific moment in time. It is approachable from the Human Mortality Database via the website linkage <www.mortality.org> (). The analysis of trends in human mortality is important to recover  the demographic impacts.  Results of such research will benefit the prediction and forecasting of future cohort mortality <cit.>. We focus on  comparison of different countries or genders,   specifically   on the  older ages over 50 years old.

There are 32 countries included in the European period life tables. It contains five Eastern European countries, Belarus, Bulgaria, Russia, Ukraine and Lithuania, and the remaining 27 Western European countries.
Following the notation introduced in Section <ref>,
X_i^(1)(s, t), i=1,…,5,  denotes the mortality rate of
the  five  Eastern European countries for subjects at age s and calendar year t,  where 50 ≤ s ≤ 90, focusing on
the death rates of older individuals, and on a recent block of 47 years, 1960 ≤ t ≤ 2006.  Similarly,  X_i^(2)(s, t), i=1,…,27,  denotes the mortality rate for other countries.
The sample  mean function  μ_1(s, t)=∑_i=1^5X_i^(1)(s, t)  and  μ_2(s, t)=∑_i=1^27X_i^(2)(s, t)  for two clusters of countries are visualized in Figure <ref>.
The heatmaps and sample mean surfaces show obvious opposite trend of mortality rates particularly for very aged people in Eastern and Western European countries as the calendar year passed 1980 or so. We apply the profile and globe test procedures to test  if the two underlying mean surfaces and its profile are different.

















According to profile test method introduced in Section <ref>,  we implement the tests (<ref>) and (<ref>).
The p-values
for fixed s^* or t^* are calculated, respectively. The results are presented in Figure <ref>.
For every fixed age s^*, we find that all of p-values are approaching to zero. This indicates that the mean mortality rates of the Eastern and Western European
is different for every age s^*= 50,⋯,90. For every fixed year t^*, almost all p-values are less than 0.05
except for years t^*=1978 and 1986.

Sequentially, we implement the globe test for the mean mortality rates of the Eastern and Western European.
The numbers of included components is J = 2, K_1=2, K_2=2  are  chosen by the fraction of variance explained (FVE) criterion with the threshold 0.90.
Based on the asymptotic distribution of the test statistic   TM, the p-value is calculated to be 0. It coincides with the intuition on images in Figure <ref> and is evidence that the mean surfaces of the mortality rates  are different between  the Eastern and Western European countries. Also, it is consistent with  the conclusion of the profile test because almost of
H_0^𝒮  and H_0^𝒯  are rejected for fixed  s^* and t^*.











Next we examine the equality of mean surfaces and its profile  between female and male clusters in West Europe.
The heatmaps and sample mean surfaces for male and female clusters    are displayed in Figure <ref>. Intuitively it does not show obvious difference. However,
all the p-values of profile tests are zero for fixed s^* and t^*.
Furthermore, we also implement globe test and obtain the p-value that is 0. Therefore, the mean surface and its profile are different in Western Europe for aged people in different gender type.






















§ DISCUSSION

Bivariate functional data have been definitely presented  in <cit.>, <cit.> and
 <cit.>.  However, inferential in testing procedures of such data type has not been adequately dealt with in the literature.
This paper develops profile and globe  tests to detect if    mean surfaces and their profile are different for  two bivariate functional samples.

In this paper, for simplicity, we assume X^(m)(s,t) are recorded on a regular and dense grid of pair time point. It is noted that the mean function and its profile estimation in (<ref>) and (<ref>) can always be obtained in the sparse and irregular setting  by additional smoothing steps.
Our proposed profile and globe tests can also be implemented if marginal eigenfunctions  and mixed  eigenfunctions can be effectively estimated. Hence, the proposed testing methodology will have wider application and much
more flexible framework.

Additionally, for  bivariate functional data X^(m)(s,t),   product  functional principal component analysis (PFPCA) due to <cit.> and double functional principal component analysis (DFPCA) due to <cit.> have been developed. Using the methodology similar to Section <ref> and <ref>,  profile and globe tests  based on PFPCA and DFPCA can also be developed.

Higher dimensional functional data also occur in practice. For example, in a recent environment and conservation project, the raw water sample has been collected weekly from different branch streams of Dongjiang at Pearl River Delta, China. Quite a few  water quality indices were measured for each sample.  The  measurements database  form
naturally a trivariate functional data.  To save cost and to monitor water quality more effectively, we are interested in detecting whether water quality has significantly changed across time, locations and water quality indices. The corresponding null hypothesis is thus then

    H_0:  μ_1(b,s,t)=μ_2(b,s,t)  vs. H_1: μ_1(b,s,t)
    	≠μ_2(b,s,t),

where μ_i(b,s,t) is the mean function for the ith sample {X_i(b,s,t)}_i=1^n_i, i=1, 2 with sample size n_i  collected on time t at location s with water quality index b measured. Although as mentioned earlier that the methods developed can be similarly extended to more than two samples, technical derivations become tedious.
We are currently working on methods for trivariate or higher order multivariate functional data.



SUPPLEMENTARY MATERIAL






  Proof of Theorem

This file is to present the detail of the proof procedure of the corresponding theorems in the article. (file type: pdf)








ECA_jasa

