Testing homogeneity: the trouble with sparse functional data Open Access

Abstract

1 Introduction

Two-sample testing of equality of distributions, which is the homogeneity hypothesis, is fundamental in statistics and has a long history that dates back to Cramér (1928), Von Mises (1928), Kolmogorov (1933), and Smirnov (1939). The literature on this topic can be categorized by the data type considered. Classical tests designed for low-dimensional data include Bickel (1969), Friedman and Rafsky (1979), Bickel and Breiman (1983), Schilling (1986), Henze (1988) among others. For recent developments that are applicable to data of arbitrary dimension, we refer to energy distance (ED) (Székely & Rizzo, 2004) and maximum mean discrepancy (Gretton et al., 2012; Sejdinovic et al., 2013). To suit high-dimensional regimes (Aoshima et al., 2018; Hall et al., 2005; Zhong et al., 2021), extensions of ED and maximum mean discrepancy were studied in Chakraborty and Zhang (2021), Gao and Shao (2021), and Zhu and Shao (2021). Some other interesting developments of ED for data residing in a metric space include Klebanov (2006) and Lyons (2013). In this paper, we focus on functional data that are random samples of functions on a real interval, e.g., $[0,1]$ (Davidian et al., 2004; Hsing & Eubank, 2015; Ramsay & Silverman, 2005; Wang et al., 2016).

Two-sample inference for functional data is gaining more attention due to the explosion of data that can be represented as functions. A substantial literature has been devoted to comparing the mean and covariance functions between two groups of curves or testing the nullity of model coefficients, see Fan and Lin (1998), Cuevas et al. (2004), Cox and Lee (2008), Panaretos et al. (2010), Zhang et al. (2010), Horváth and Kokoszka (2012), Zhang and Liang (2014), Staicu et al. (2015), Pini and Vantini (2016), Paparoditis and Sapatinas (2016), Wang et al. (2018), Guo et al. (2019), Zhang et al. (2019), He et al. (2020), Yuan et al. (2020), and Wang (2021). However, the underlying distributions of two random functions can have the same mean and covariance function, but differ in other aspects. Testing homogeneity, which refers to a hypothesis testing procedure to determine the equality of the underlying distributions of any two random objects, is thus of particular interest and practical importance. The literature on testing the homogeneity for functional data is much smaller and restricted to the two-sample case based on either fully observed (Benko et al., 2009; Cabaña et al., 2017; Krzyśko & Smaga, 2021; Wynne & Duncan, 2022) or intensively measured functional data (Hall & Keilegom, 2007; Jiang et al., 2019).

For intensively measured functional data, a presmoothing step is often adopted to each individual curve in order to construct a smooth curve before carrying out subsequent analysis, which may reduce mean square error (Ferraty et al., 2012) or, with luck, remove the noise, a.k.a. measurement error, in the observed discrete data. This results in a two-stage procedure, smoothing first and then testing homogeneity based on the presmoothed curves. For example, Hall and Keilegom (2007) and Jiang et al. (2019) adopted such an approach. In addition, the ED in Székely and Rizzo (2004) could be extended to the space of $L2$ functions to characterize the distribution differences of random functions, provided the functional data are fully observed without errors (Klebanov, 2006).

Specifically, the ED between two random functions X and Y is defined as

where $X′,Y′$ are i.i.d copies of $X,Y$⁠, respectively, and $‖X−Y‖L2=(∫01(X(t)−Y(t))2dt)1/2$⁠. According to the results of Klebanov (2006) and Lyons (2013), $ED(X,Y)$ fully characterizes the distributions of X and Y in the sense that $ED(X,Y)≥0$ and $ED(X,Y)=0⇔X=dY$⁠, where we use $X=dY$ to indicate that $X,Y$ are identically distributed. Given two samples of functional data ${Xi}i=1n$ and ${Yi}i=n+1n+m,$ which are intensively measured at some discrete time points, the reconstructed functions, denoted as ${X^i}i=1n$ and ${Y^i}i=n+1n+m$⁠, can be obtained using the aforementioned presmoothing procedure. Then, $ED(X,Y)$ can be estimated by the U-type statistic

While the above presmoothing procedure to reconstruct the original curves may be promising for intensely measured functional data, it has not yet been utilized to our knowledge, perhaps due to the technical and practical challenges to implement this approach as the level of intensity in the measurement schedule and the proper amount of smoothing are both critical. First, the distance between the reconstructed functional data and its target (the true curve) needs to be tracked and reflected in the subsequent calculations. Second, such a distance would depend on the intensity of the measurements and the bandwidth used in the presmoothing stage. Neither is easy to nail down in practice.

Furthermore, in real world applications, such as in longitudinal studies, each subject often may only have a few measurements, leading to sparse functional data (Yao et al., 2005a). Here, presmoothing individual data no longer works and one must borrow information from all subjects to reconstruct the trajectory of an individual subject. The principal analysis through conditional estimation (PACE) approach in Yao et al. (2005a) offers such an imputation method, yet it does not lead to consistent estimates of the true curve for sparsely observed functional data as there are not enough data available for each individual. Consequently, the quantities $E[‖X−Y‖L2],E[‖X−X′‖L2],andE[‖Y−Y′‖L2]$ in equation (1) are not consistently estimable as these expectations are outside the corresponding $L2$ norms. Pomann et al. (2016) reduced the problem to testing the homogeneity of the scores of the two processes by assuming that the random functions are finite dimensional. Such an approach would not be consistent either as the scores still cannot be consistently estimated for sparse functional data. To our knowledge, there exists no consistent test of homogeneity for sparse functional data. In fact, it is not feasible to test full homogeneity based on sparsely observed functional data as there are simply not enough data for such an ambitious goal.

This seems disappointing, since although it has long been recognized that sparse functional data are much more challenging to handle than intensively measured functional data, much progress has been made to resolve this challenge. For instance, both the mean and covariance function can be estimated consistently at a certain rate (Li & Hsing, 2010; Yao et al., 2005a; Zhang & Wang, 2016) for sparsely observed functional data. Moreover, the regression coefficient function in a functional linear model can also be estimated consistently with rates (Yao et al., 2005b). This motivated us to explore a less stringent concept of homogeneity that can be tested consistently for sparse functional data. In this paper, we provide the answer by proposing a test of marginal homogeneity for two independent samples of functional data. For ease of presentation, we assume that the random functions are defined on the unit interval $[0,1]$⁠.

Definition of marginal homogeneity. Two random functions X and Y defined on [0, 1] are marginal homogeneous if

From this definition we can see that unlike testing homogeneity that involves testing the entire distribution of functional data, testing marginal homogeneity only involves simultaneously testing the marginal distributions at all time points. This is a much more manageable task that works for all sampling designs, be it intensively or sparsely observed functional data, and it is often adequate in many applications. Testing marginal homogeneity is not new in the literature and has been investigated by Chakraborty and Zhang (2021) and Zhu and Shao (2021) for high-dimensional data. They show that the marginal tests can be more powerful than their joint counterparts under the high-dimensional regime. In a larger context, the idea of aggregating marginal information originates from Zhu et al. (2020), where they consider a related problem of testing the independence between two high-dimensional random vectors.

For real applications of testing marginal homogeneity, taking the analysis of biomarkers over time in clinical research as an example, comparing differences between marginal distributions of the treatment and control groups may be sufficient to establish the treatment effect. To contrast stocks in two different sectors, the differences between marginal distributions might be more important than the differences between joint distributions. In addition, differences between marginal distributions can be seen as the main effect of differences between distributions. Thus, it makes good sense to test marginal homogeneity, especially in situations where joint distribution testing is not feasible or inefficient.

Let λ be the Lebesgue measure on $R$⁠. The focus of this paper is to test

This can be accomplished through the marginal energy distance (MED) defined as

where $X′$ and $Y′$ are independent copies of X and Y, respectively. Indeed, MED is a metric for marginal distributions in the sense that $MED(X,Y)≥0$ and $MED(X,Y)=0⇔X(t)=dY(t)$ for almost all $t∈[0,1]$⁠. A key feature of our approach is that it can consistently estimate $MED(X,Y)$ for all types of sampling plans. Moreover, $E[|X(t)−Y(t)|]$⁠, $E[|X(t)−X′(t)|]$⁠, and $E[|Y(t)−Y′(t)|]$ can all be reconstructed consistently for both intensively and sparsely observed functional data. Such a unified procedure for all kinds of sampling schemes may be more practical as the separation between intensively and sparsely observed functional data is usually unclear in practical applications. Moreover, it could happen that while some of the subjects are intensively observed, others are sparsely observed. In the extremely sparse case, our method can still work if each subject only has one measurement.

Measurement errors (or noise) are common for functional data, so it is important to accommodate them. If noise is left unattended, there will be bias in the estimates of MED as the observed distributions are no longer the true distributions of X and Y. One might hope that the measurement errors can be averaged out during the estimation of the function $E[|X(t)−Y(t)|]$ in $MED(X,Y)$ in analogy to estimating the mean function $μ(t)=E[X(t)]$ or covariance function $C(s,t)=E[(X(t)−μ(t))(X(s)−μ(s))]$ (Yao et al., 2005a). However, this is not the case. To see why, let $e1$⁠, $e1′$⁠, $e2$⁠, and $e2′$ be independent white noise. When estimating mean or covariance function at any fixed time $t,s∈[0,1]$⁠, it holds that $μ(t)=E[X(t)+e1]$ and $C(s,t)=E[(X(t)−μ(t)+e1)(X(s)−μ(s)+e1′)]$ for $t≠s$⁠. But $E[|X(t)−Y(s)+e1−e2|]≠E[|X(t)−Y(s)|]$⁠. Likewise, we can see that the ED in (1) would have the same challenge to handle measurement errors unless these errors were removed in a presmoothing step before carrying out the test. So the challenges with measurement errors is not triggered by the use of the $L1$ norm in MED. The $L2$ norm in ED will face the same challenge.

For intensely measured functional data, a presmoothing step is often used to handle measurement errors in the observed data in the hope that smoothing will remove the error. However, this is a delicate issue, as it is difficult to know the amount of smoothing needed in order for the subsequent analysis to retain the same convergence rate as if the true functional data were fully observed without errors. For instance, Zhang and Chen (2007) study the effects of smoothing to obtained reconstructed curves and show that in order to retain the same convergence rate of mean estimation for fully observed functional data, the number of measurements per subject that generates the curves must be of higher order than the number of independent subjects. This requires functional data that are intensively sampled well beyond ultra-dense (or dense) functional data that have been studied in the literature (Zhang & Wang, 2016).

In this paper, we propose a new way of handling measurement errors so that the $MED$-based testing procedure is still consistent in the presence of measurement errors. The key idea is to show that when the measurement errors $e1$ and $e2$ of X and Y, respectively, are identically distributed, i.e., $e1=de2$⁠, the $MED(X,Y)$-based approach can still be applied to the contaminated data with consistency guaranteed under mild assumptions (cf. Corollary 3). When $e1≠de2$⁠, we propose an error-augmentation approach, which can be applied jointly with our unified estimation procedure.

To conduct hypothesis testing, $MED(X,Y)$-based approaches can be implemented as permutation tests. We refer the book by Lehmann and Romano (2005) for an extensive introduction to permutation test. The use of permutation test for distance/kernel-based tests is not new in the literature. The consistency of permutation test for distance covariance (Székely et al., 2007) or Hilbert–Schmidt independence criterion (Gretton et al., 2007) has been investigated by Pfister et al. (2018), Rindt et al. (2021), and Kim et al. (2022) for low-dimensional data, where distance covariance and Hilbert–Schmidt independence criterion are distance-based and kernel-based independence measures, respectively. Under the high dimension, low sample size setting (Hall et al., 2005; Zhu & Shao, 2021) studied the size and power behaviours of permutation test for ED. The difference in the proof of consistency of permutation test between longitudinal data and vector-valued data is substantial. For instance, MED estimated from longitudinal data is an integration of local weighted averages, where the measurement time points in the weights are also random and hence requires special handling in the proof.

The rest of the paper is organized as follows. Section 2 contains the main methodology and supporting theory about testing marginal homogeneity. Numerical studies are presented in Section 3. The conclusion is in Section 4. All technical details are postponed to Section 5.

2 Testing marginal homogeneity

We first consider the case where there are no measurement errors and postpone the discussion of measurement errors to the end of this section. Let ${Xi}i=1n$ and ${Yi}i=n+1n+m$ be i.i.d copies of X and Y, respectively. In practice, the functions are only observed at some discrete points:

This sampling plan allows the two samples to be measured at different schedules and additionally each subject within the sample has its own measurement schedule. This is a realistic assumption but the consequence is that two-dimensional smoothers will be needed to estimate the targets. Fortunately, the convergence rate of our estimator attains the same convergence rate as that of a one-dimensional smoothing method. This intriguing phenomenon will be explained later.

For notational convenience, denote $Zi=Xi,ifi=1,2,…,n$ and $Zi=Yi,ifi=n+1,…,n+m$⁠. Let $Z=(z1;…;zn+m)$ be the combined observations, where for $i=1,2,…,n+m$⁠, $zi$ is a vector of length $Ni$⁠,

The observations corresponding to X and Y are defined as $X=(z1;…;zn)$ and $Y=(zn+1;…;zn+m)$⁠, respectively.

To estimate $MED(X,Y)$ in (4), note that we actually have no observations for the one-dimensional functions $E[|X(t)−Y(t)|]$⁠, $E[|X(t)−X′(t)|]$⁠, and $E[|Y(t)−Y′(t)|]$⁠, due to the longitudinal design where X and Y are observed at different time points. Thus, the sampling schedule for X and Y are not synchronized. A consequence of such asynchronized functional/longitudinal data is that a one-dimensional smoothing method that has typically been employed to estimate a one-dimensional target function, e.g., $E[|X(t)−Y(t)|],$ does not work here because $X(t)−Y(t)$ cannot be observed at any point t. However, a workaround is to estimate the following two-dimensional functions first:

then set $t1=t2=t$ in all three estimators. Since $G1,G2,G3$ can all be recovered by some local linear smoother, $MED(X,Y)$ admits consistent estimates for both intensively and sparsely observed functional data. For instance, $G1(t1,t2)$ can be estimated by $G^1(t1,t2)=β^0$⁠, where

and $G2(t1,t2)$ can be estimated by $G^2(t1,t2)=α^0$⁠, where

$G3(t1,t2)$ can be estimated similarly as $G2(t1,t2)$ by an estimator $G^3(t1,t2)$⁠. In the above, $Ni1$ and $Ni2$ should be understood as the respective length of the vector $zi1$ and $zi2$ and $Kh(⋅)=K(⋅/h)/h$ is a one-dimensional kernel with bandwidth h. The sample estimate of $MED(X,Y)$ can then be constructed as

For hypothesis testing, the critical value or p-value can be determined by permutations (Lehmann & Romano, 2005). To be more specific, let $π:{1,2,…,n+m}→{1,2,…,n+m}$ be a permutation. There are $(n+m)!$ number of permutations in total and we denote the set of permutations as $Pn+m={πl:l=1,2,…,(n+m)!}$⁠. For $l=1,2,…,(n+m)!$⁠, define the permutation of $πl$ on $Z$ as

Write the statistic that is based on the permuted sample $πl⋅Z$ as $MEDn(πl⋅Z)$ and let $Π1,…,ΠS−1$ be i.i.d and uniformly sampled from $Pn+m$⁠, we define the permutation-based p-value as

Then, the level-α permutation test w.r.t. $MEDn(Z)$ can be defined as

2.1 Convergence theory

In this subsection, we show that $MEDn(Z)$ is a consistent estimator and develop its convergence rate.

The next assumption specifies the relationship of the number of observations per subject and the decay rate of the bandwidth parameters $hx,hy$⁠.

The following theorem states that we can consistently estimate $MED(X,Y)$ with sparse observations.

Given two sequences of positive real numbers $an$ and $bn$⁠, we say that $an$ and $bn$ are of the same order as $n→∞$ (denoted as $an≍bn$⁠) if there exists constants $0<c1<c2<∞$ such that $c1≤limn→∞an/bn≤c2$ and $c1≤limn→∞bn/an≤c2$⁠. The convergence rates of $MEDn(Z)$ for different sampling plans are provided in the following corollary.

2.2 Validity of the permutation test and power analysis

We now justify the permutation-based test for sparsely observed functional data. Under the null hypothesis and the mild assumption that ${Ni}$ are i.i.d across subjects, the size of the test can be guaranteed by the fact that the distribution of the sample is invariant under permutation. For a rigorous argument, see Theorem 15.2.1 in Lehmann and Romano (2005). Thus, the permutation test based on the test statistic (7) produces a legitimate size of the test. The power analysis is much more challenging and will be presented below.

Let $π∈Pn+m$ be a fixed permutation and $G^π,I(t1,t2),I=1,2,3$⁠, be the estimated functions from algorithms (5) and (6) using permuted samples $π⋅Z$⁠. The conditions on the decay rate of bandwidth parameters $hx,hy$ that ensure the convergence of $G^π,I$ for any fixed permutation π are summarized below.

Let Π be a random permutation uniformly sampled from $Pn+m$⁠. If the sample is randomly shuffled, it holds that $ZΠ(i)=dZΠ(j)$ and $MED(ZΠ(i),ZΠ(j))=0$ for any $i,j=1,2,…,n+m$⁠. For the sample estimate $MEDn(Π⋅Z)$ based on the permuted sparse observations, we show that $MEDn(Π⋅Z)$ converges to 0 in probability.

For the original data that have not been permuted, $MEDn(Z)→pMED(X,Y)$⁠, which is strictly positive under the alternative hypothesis. On the other hand, we know from Theorem 2 that the permuted statistics converges to 0 in probability. This suggests that under mild assumptions, the probability of rejecting the null approaches 1 as $n,m→∞$⁠. We make this idea rigorous in the following theorem.

2.3 Handling of measurement errors

With the presence of measurement errors, the actual observed data are

where ${eij:i=1,2,…,n,j=1,2,…,Ni}∼i.i.de1$⁠, ${eij:i=n+1,…,n+m,j=1,2,…,Ni}∼i.i.de2$⁠, and $e1,e2$ are mean 0 independent univariate random variables. Denote the combined noisy observations by $Z~=(z~1;…;z~n+m)$⁠, where

The local linear smoothers described in equations (5) and (6) are then applied with the input data ${zi1j1}$ and ${zi2j2}$ replaced, respectively, by ${z~i1j1}$ and ${z~i2j2}$⁠. The resulting outputs are denoted as $H^1,H^2,H^3$⁠, leading to the estimator

Correspondingly, the proposed test with contaminated data $Z~$ is

where $p~=(1/S){1+∑l=1S−1I{MEDn(Πl⋅Z~)≥MEDn(Z~)}}.$ To study the convergence of $MEDn(Z~)$⁠, define the two-dimensional functions $H1(s,t)$⁠, $H2(s,t)$⁠, and $H3(s,t)$ as

where $e1′$ and $e2′$ are independent and identical copies of $e1$ and $e2$⁠, respectively. The target of $MEDn(Z~)$ is shown to be

To show the approximation error of $MEDn(Z~)$⁠, we need the following assumptions.

By the property of ED, it holds that $MED~(X,Y)≥0$ and $MED~(X,Y)=0⇔X(t)+e1=dY(t)+e2$ for almost all $t∈[0,1]$⁠. Then, we show that under the following assumptions, the condition that $X(t)+e1=dY(t)+e2$ would imply the homogeneity of $X(t)$ and $Y(t)$⁠.

For common distributions, such as Gaussian and Cauchy, their characteristic functions are of exponential form and have no real zeros. Some other random variables, such as exponential, Chi-square, and gamma, have characteristic functions of the form $ψ(t)=(1−itθ)−k$ with only a finite number of real zeros. Since it is common to assume Gaussian measurement errors, the restriction on the real zeros of the characteristic function in Assumption 5(E.2) is very mild.

Based on the above theorem, we have the important property that $MED~(X,Y)=0⇔X(t)=dY(t)$ for almost all $t∈[0,1]$⁠. As discussed before, $MED~(X,Y)$ can be consistently estimated by $MEDn(Z~)$ and the test can be conducted via permutations. Consequently, the data contaminated with measurement errors can still be used to test marginal homogeneity as long as $e1=de2$⁠. We make this statement rigorous below.