# Doubly ranked tests for grouped functional data

Mark J. Meyer
Department of Mathematics and Statistics, Georgetown University
Washington, DC USA

###### Abstract

Nonparametric tests for functional data are a challenging class of tests to
work with because of the potentially high dimensional nature of functional
data. One of the main challenges for considering rank-based tests, like the
Mann-Whitney or Wilcoxon Rank Sum tests (MWW), is that the unit of observation
is a curve. Thus any rank-based test must consider ways of ranking curves.
While several procedures, including depth-based methods, have recently been
used to create scores for rank-based tests, these scores are not constructed
under the null and often introduce additional, uncontrolled for variability.
We therefore reconsider the problem of rank-based tests for functional data
and develop an alternative approach that incorporates the null hypothesis
throughout. Our approach first ranks realizations from the curves at each time
point, then summarizes the ranks for each subject using a sufficient statistic
we derive, and finally re-ranks the sufficient statistics in a procedure we
refer to as a doubly ranked test. As we demonstrate, doubly rank tests are
more powerful while maintaining ideal type I error in the two sample, MWW
setting. We also extend our framework to more than two samples, developing a
Kruskal-Wallis test for functional data which exhibits good test
characteristics as well. Finally, we illustrate the use of doubly ranked tests
in functional data contexts from material science, climatology, and public
health policy.

##  1 Introduction

Expanded research in functional data analysis has produced a variety of tests
couched in the framework of classic nonparametric tests including work by Hall
and Van Keilegom (2007); López-Pintado and Romo (2009); López-Pintado et al.
(2010); Chakraborty and Chaudhuri (2015); Pomann et al. (2016); López-Pintado
and Wrobel (2017); Abramowicz et al. (2018); Berrett et al. (2021) and
Meléndez et al. (2021). These various authors have built out hypothesis
testing procedures for functional versions of Mann-Whitney or Wilcoxon Rank
Sum, Anderson-Darling tests, and permutation tests, among others. The primary
feature of functional data that makes the development of such nonparametric
tests challenging is that the unit of observation is a curve. For example,
$X_{g,1}(s),\ldots,X_{g,n_{g}}(s)$ are a sample of $n_{g}$ curves from groups
$g=1,\ldots,G$. Each curve is typically measured on a fine grid,
$\mathcal{S}=\\{s:s=s_{1},\ldots,s_{S}\\}$ for a total of $S$ measurements.
Any rank-based procedure must employ a way to rank curves directly or at least
rank summary metrics of each subject’s curve.

Depth is a popular and useful concept for ranking functions and many authors
have utilized depth for Mann-Whitney-Wilcoxon (MWW) type tests (López-Pintado
and Romo, 2009; López-Pintado et al., 2010; López-Pintado and Wrobel, 2017),
functional boxplots (Sun and Genton, 2011; Genton et al., 2014; Wrobel et al.,
2016), outlier detection (Sun and Genton, 2012; Dai and Genton, 2018; Huang
and Sun, 2019; Dai et al., 2020; Alemán-Gómez et al., 2022), and ranking
curves (López-Pintado and Romo, 2009, 2011; Sun et al., 2012; Sguera and
López-Pintado, 2021). The basic idea behind depth is to characterize the
proportion of time a curve is bounded above and below by another curve in the
sample. The resulting scores, when ordered, rank the curves from the inside
out with the smallest values representing the extremes, the maximum or
minimum, and the largest value representing the median (López-Pintado and
Romo, 2009; Sun et al., 2012). These scores are effectively a measure of the
distance a curve is from the middle. MWW-type tests using depth to test for a
difference between two groups must construct a third, artificial reference
group (López-Pintado and Romo, 2009; López-Pintado et al., 2010; López-Pintado
and Wrobel, 2017). Depth is first determined by group, then the test procedure
counts the number of depth values in each group that are larger than the depth
in the reference group. Code is publicly available for the depth-based MWW
test, see López-Pintado and Wrobel (2017).

Alternative approaches include a spatial sign-based statistic by Chakraborty
and Chaudhuri (2015) and a random projections based approach by Meléndez et
al. (2021). In the former, Chakraborty and Chaudhuri (2015) develop an MWW-
type statistic using the spatial sign for infinite-dimensional data that
weakly convergences to a Gaussian distribution under limited assumptions. The
spatial sign is used as a distance metric to construct the test, taking
pairwise distances between groups with data assumed to arise from a random
variable defined on a separable Hilbert space. The latter alternative
approach, a random projection-based test by Meléndez et al. (2021), maps the
points from the high dimensional functional space to a randomly chosen low-
dimensional space defined using Brownian motion. This approach has application
in parametric tests for grouped functional as well, see for example Cuesta-
Albertos and Febrero-Bande (2010). The random projections-based approach
effectively generates scores for each subject via a random basis function and
an integral approximation of the subject-specific curve. It then treats these
scores as the data in a traditional MWW test, thus the scores are ranked and
then summed by group. Neither of these alternative methods have publicly
available code.

Both depth- and random projection-based approaches perform well in larger
samples (López-Pintado and Romo, 2009; López-Pintado et al., 2010; Meléndez et
al., 2021). However, each method also introduces an additional and
uncontrolled for source of variation: the random artificial reference group
for the depth-based approach and the brownian motion-based basis functions for
the random projections. While depth is useful for ranking curves particularly
when identifying medians and outliers, it does not produce a score that ranks
in a traditional fashion, from smallest to largest. When using random
projections, the integration is performed before constructing ranks which
ignores the fact that rank may be dynamic over time. To handle the curve
ranking problem, both procedures rely on subject-specific summary scores,
although neither constructs their scores under the null. Thus, the tests are
only conducted under the null when comparing the scores between groups.

The spatial sign-based approach differs from the depth- and random projection-
based methods in that, for finite-dimensional approximations of functional
data, it effectively creates a time-point specific score instead of subject-
specific scores over time. Chakraborty and Chaudhuri (2015) show that the
spatial sign-based test performs well empirically when the underlying data is
Gaussian with a relatively small sample. The test does consider the MWW null
throughout its construction, but appears to be underpowered when the data is
non-Gaussian. While all of these methods are formulated for the MWW or two
group case, to our knowledge no prior work considers the Kruskal-Wallis (KW)
or three or more group case. Thus, with respect to these authors, we believe
that a further examination of this problem is warranted.

In this manuscript, we re-examine rank-based test settings for comparing
groups of functional data. We propose a general testing procedure working
within the MWW or KW assumptions, developing our tests under the relevant null
hypotheses at all stages—a general overview of MWW and KW tests as well as
functional data is given in Section 2. Our approach first constructs ranks of
either the raw data or preprocessed functional data at each discretely
measured time point giving every subject a curve of ranks. Under the MWW and
KW null hypotheses, we then derive the exact and approximate distribution of
each subject’s $r$th order statistic at a given time point (Section 3). The
resulting approximate distribution can be shown to be an exponential family
member from which we obtain sufficient statistics for each subject’s value of
$r$. The summaries resulting from the sufficient statistic become our new
data, or rank-based scores, which are then ranked and analyzed in the usual
way by the relevant tests. Because the MWW and KW tests ultimately re-rank our
rank-based scores, we refer to our testing procedures as doubly ranked tests.

We empirically demonstrate the power of doubly ranked tests in Section 4. Our
approach can accommodate curves with a potentially large number of
measurements and groups of size $G\geq 2$. Further, we show that doubly ranked
tests perform well in both small and larger sample sizes. We illustrate their
use with an analysis of data from several different studies across different
substantive fields (Section 5). The illustrations include data on resin
viscosity, Canadian weather patterns, and changes in driving requests in the
United States prior to and just after the initiation of COVID-19 policy
measures in various states. While we formulate doubly ranked tests in the
functional data context with the preprocessed data using functional principal
components, the tests are applicable to any high dimensional data testing
setting. We discuss this, and our work in general, in Section 6.

##  2 Statistical Framework

Because our method blends several statistical frameworks, we now briefly
describe each beginning with the MWW and KW tests in their scalar forms. We
then discuss some general concepts for functional data. Functional data
analysis is a large area of research. Review articles by Morris (2015), Wang
et al. (2016), and Greven and Scheipl (2017), for example, or the text by
Ramsay and Silverman (2005) will provide more in-depth discussions of this
topic.

###  2.1 Mann-Whitney-Wilcoxon

Let $X_{1},\ldots,X_{n_{1}}$ be observed data from group 1,
$Y_{1},\ldots,Y_{n_{2}}$ be observed data from group 2, and $n=n_{1}+n_{2}$.
Assuming the groups are independent, the null and alternative hypotheses can
have one of two forms. The first tests for stochastic ordering where the null
is $H_{0}:F_{X}(c)=F_{Y}(c)$ and the alternative is $H_{A}:F_{X}(c)\leq
F_{Y}(c)$ for the CDFs $F_{X}(c)$ and $F_{Y}(c)$ in each group. Assuming these
distributions are same up to location shift $\Delta$, i.e.
$F_{X}(c)=F_{Y}(c-\Delta)$, gives a second form of the hypotheses which is
$H_{0}:\Delta=0$ and $H_{A}:\Delta\neq 0$.

The test statistic to test either form of $H_{0}$ is
$T^{+}=T-\frac{n_{2}(n_{2}+1)}{2}$ where $T=\sum_{i=1}^{n_{2}}R[Y_{i}]$ for
$R[Y_{i}]$ equal to the rank of $Y_{i}$ among the combined samples, i.e. among
$X_{1},\ldots,X_{n_{1}},$ $Y_{1},\ldots,Y_{n_{2}}$. The sum of the ranks, $T$,
is the test statistic for the Wilcoxon rank sum test (Wilcoxon, 1945; Kloke
and McKean, 2015). $T^{+}$ is also equivalent to the test statistic described
by Mann and Whitney (1947). Consequently, we refer to this test as the MWW
test.

The exact distribution of $T$ is non-standard but can be easily computed by
standard statistical software for samples up to $n=50$, provided there are no
tied ranks. For larger samples and in the presence of ties, a normal
approximation to the distribution of the ranks is used instead, see the
wilcox.test() documentation in the stats package in R (R Core Team, 2022). The
test statistic based off of the normal approximation is

$\displaystyle
T_{z}=\left[T-\frac{n_{2}(n+1)}{2}\right]\bigg{/}\sqrt{\frac{n_{1}n_{2}(n+1)}{12}},$

which has a standard normal distribution.

###  2.2 Kruskal-Wallis

The KW test generalizes the MWW to three or more groups (Kruskal and Wallis,
1952; Kloke and McKean, 2015). Let $G$ denote the total number of independent
groups with data $X_{g,1},\ldots,X_{g,n_{g}}$ for $g=1,\ldots,G$. The KW test
assumes the distributions of the data are the same except possibly for some
location parameter or parameter(s). The null hypothesis for the KW test
generalizes the second MWW null where groups differ only up to some set of
location shifts, $\Delta_{g}$. The null then has the form
$H_{0}:\Delta_{1}=\ldots=\Delta_{G}$ with the alternative hypothesis,
$H_{A}:\Delta_{g}\neq\Delta_{g^{\prime}}$ for some $g\neq g^{\prime}$.

For $n=\sum_{g=1}^{G}n_{g}$, the test statistic to test $H_{0}$ is

$\displaystyle
H=\frac{12}{n(n+1)}\sum_{g=1}^{G}n_{g}\left(\bar{R}_{g}-\frac{n+1}{2}\right)^{2},$

where $\bar{R}_{g}$ is the average rank in group $g$, i.e.
$\bar{R}_{g}=n_{g}^{-1}\sum_{i=1}^{n_{g}}R[X_{g,i}]$.Under $H_{0}$, $H$ is
distribution-free. Tables for its exact distribution do exist, see for example
Hollander and Wolfe (1999), Chapter 6. $H$ is asymptotically $\chi_{G-1}^{2}$
under $H_{0}$ as well, see Kruskal and Wallis (1952) or Kloke and McKean
(2015) which is what the the stats package in R uses when kruskal.test() is
called (R Core Team, 2022).

###  2.3 Functional Data

Figure 1: Blue curves denote the stations in the Arctic region, tan curves
indicate stations the Atlantic region, red curves for stations in the
Continental region, and green curves are the Pacific region.

Functional data are data observations whose units are curves, often measured
over time but they can be measured over location as well. These curves are
typically sampled at a high frequency over a pre-defined interval. The data
comes sampled on a grid, $\\{s:s=s_{1},\ldots,s_{S}\\}$ for example as we note
in Section 1. The data used for analysis may be the raw values themselves,
$X_{1}(s),\ldots,X_{n_{1}}(s)$, or they may be preprocessed and projected into
a well-behaved space using a known basis function. The idea behind
preprocessing is to model the underlying function from which the raw data
values were observed and then analyze the functions themselves. We denote the
preprocessed data with $\tilde{X}_{1}(s),\ldots,\tilde{X}_{n_{1}}(s)$. There
are many ways to preprocess functional data including functional principal
components, wavelet-based techniques, and Fourier transforms (Wang et al.,
2016). Our method is compatible with either raw or pre-processed data, but for
this work we will preprocess functional data using a sandwich smoother-based
approach referred to as Fast Covariance Estimation or FACE (Xiao et al., 2013,
2016). FACE is a nonparametric approach for constructing functional principal
components based on a set of curves with a potentially large sampling density.
It can be implemented using the refund package in R (Goldsmith et al., 2022).
As an example, Figure 1 contains the preprocessed (using FACE) daily
temperature curves from the Canadian Weather dataset (Ramsay and Silverman,
2005; Ramsay et al., 2009). The dataset contains daily temperature and
precipitation measurements, averaged over the years 1960 to 1994 from 35
different sites across Canada. We analyze this data in Section 5.

##  3 Doubly Ranked Tests

Suppose we obtain a sample of curves from two or more groups:
$X_{g,1}(s),\ldots,X_{g,n_{g}}(s)$ for $g=1,\ldots,G$, $s=s_{1},\ldots,s_{S}$,
and $n_{1}+\ldots+n_{g}+\ldots+n_{G}=n$. Measurements are discretely sampled
on the same grid, $\mathcal{S}=\\{s:s=s_{1},\ldots,s_{S}\\}$, and come from
continuous process. Under the null, we assume no difference in the curves
between the $G$ groups across all $s$. To avoid having to make a
distributional assumption on the process that generates $X_{g,i}(s)$, we
evaluate this null by analyzing the ranks of $X_{g,i}(s)$ at each measurement
occurrence $s$, i.e. time point or location depending on what the function is
measured over. Thus, we first rank either the raw or preprocessed data at each
measurement occurrence, ignoring group assignment. This results in a curve of
ranks for each subject which we denote as $z_{g,i}(s)$. If there is no
difference in the curves across $s$ under the null, then we should also
observe no difference in the ranks across $s$. Assuming each subject’s curve
has a true order in the sample, we can consider
$\\{z_{g,i}(s_{1}),\ldots,z_{g,i}(s_{S})\\}$ as a set of draws from the
distribution describing that true order. We then construct a summary of each
curve and re-rank the summaries for analysis using either the MWW or KW tests,
depending on the number of groups. Assuming the summaries are reasonable
representations of the subject’s order in the sample, no difference in the
curves and subsequent ranks over $s$ between the groups implies no difference
in the summarized ranks between the groups.

More formally, the doubly ranked test inherits the null and alternatives from
the MWW and KW tests. Since the summaries are based on a sample of observed
ranks, it is reasonable to assume the distributions are the same up to some
location shift. Let $t_{1},t_{2},\ldots t_{g},\ldots,t_{G}$ denote the
observed summarized ranks from their corresponding random variables
$T_{1},T_{2},\ldots T_{g},\ldots,T_{G}$. The doubly ranked MWW test then has
the null and alternative of $H_{0}:\theta=0$ and $H_{0}:\theta\neq 0$ where
$\theta$ is a parameter such that $F_{T_{1}}(t_{1})=F_{T_{2}}(t_{2}-\theta)$.
When $G\geq 3$, the null and alternative hypotheses for the doubly ranked KW
test are $H_{0}:\theta_{1}=\ldots=\theta_{G}$ with the alternative hypothesis,
$H_{A}:\theta_{g}\neq\theta_{g^{\prime}}$ for some $g\neq g^{\prime}$. In
other words, under $H_{0}$ we assume there is no difference between the groups
in the location parameters of the distributions of summarized ranks. The
alternative is that the there is a difference, when $G=2$, and at least one
difference, when $G\geq 3$.

###  3.1 Test Derivation

As in the univariate tests, the ranks are randomly assigned under $H_{0}$. We
apply this both at the measurement occurrence level and again at the summary
level when the summaries are re-ranked. Thus, for each subject and at each
measurement occurrence, $z_{g,i}(s)\sim\mathcal{U}\\{1,n\\}$ where $n$ is the
total number of curves and $\mathcal{U}$ denotes the discrete uniform
distribution. To obtain the summaries, consider the sampling distribution of
the $r$th order statistic for any given subject at time $s$ which we denote
$Z_{(r)}$. We let $z$ generically denote an observed rank, i.e. $z_{g,i}(s)$,
for ease of derivation. The PMF for $Z_{(r)}$ is given by

$\displaystyle P\bigg{[}Z_{(r)}=z\bigg{]}=1-\sum_{k=0}^{r-1}{n\choose
k}\left(\frac{z}{n}\right)^{k}\left(1-\frac{z}{n}\right)^{n-k}$ (1)
$\displaystyle-\left[1-\sum_{k=0}^{r-1}{n\choose
k}\left(\frac{z-1}{n}\right)^{k}\left(1-\frac{z-1}{n}\right)^{n-k}\right].$

The summations in Equation (1) are the CDFs of binomial random variables with
probabilities of success $\frac{z}{n}$ and $\frac{z-1}{n}$, respectively. Let
$I_{q}(a,b)$ denote the regularized incomplete beta function evaluated at $q$
with parameters $a$ and $b$. Equation (1) can be re-expressed as

$\displaystyle P\bigg{[}$ $\displaystyle Z_{(r)}=z\bigg{]}$ (2)
$\displaystyle=I_{z/n}(r,n-r+1)-I_{(z-1)/n}(r,n-r+1).$

The distribution in Equation (2) is non-standard and difficult to work with
analytically, although we can compute values from it using standard
statistical software. The relationship between the beta and binomial
distributions can, however, be exploited to find an approximation of
$P\left[Z_{(r)}=z\right]$ with a more workable form. Specifically, the
regularized incomplete beta functions in Equation (2) are also equivalent to
the CDFs of beta random variables evaluated at $\frac{z}{n}$ and
$\frac{z-1}{n}$, respectively, with parameters $r$ and $n-r+1$. After
expressing the beta CDFs in their integral form, we re-write
$P\left[Z_{(r)}=z\right]$ as

$\displaystyle P\bigg{[}$ $\displaystyle Z_{(r)}=z\bigg{]}$
$\displaystyle=\frac{\Gamma(n+1)}{\Gamma(r)\Gamma(n-r+1)}\left[\int_{\frac{z}{n}-\frac{1}{n}}^{\frac{z}{n}}t^{r-1}\left(1-t\right)^{n-r}dt\right].$

A full derivation of $P\left[Z_{(r)}=z\right]$ is in Appendix A.

We approximate the integral using the midpoint rule and select one interval,
$[\frac{z}{n}-\frac{1}{n},\frac{z}{n}]$, with a midpoint at
$\frac{z}{n}-\frac{1}{2n}$. This value is also approximately equal to the
optimal numerical solution to the mean value theorem problem for this
integral, see Section 1 of Supplement A (Meyer, 2024a) for a numerical
comparison. The mass function in (2) can then be approximated by

$\displaystyle P\bigg{[}Z_{(r)}=z$
$\displaystyle\bigg{]}\approx\frac{\Gamma(n+1)}{\Gamma(r)\Gamma(n-r+1)}\frac{1}{n}$
(3)
$\displaystyle\times\left(\frac{z}{n}-\frac{1}{2n}\right)^{r-1}\left(1-\frac{z}{n}+\frac{1}{2n}\right)^{n-r},$

with the approximation improving as $n$ increases. We observe this graphically
in Section 1 of Supplement A (Meyer, 2024a). When $r=\frac{n+1}{2}$, the
expected rank under $H_{0}$, the approximation is quite good.

######  Claim 1.

The distribution described by the mass function in Equation 3 is an
exponential family member.

###### Proof.

The proof of Claim 1 is in Appendix B. ∎

Since Equation (3) is an exponential family,
$t(z)=\log\left[\left(\frac{z}{n}-\frac{1}{2n}\right)\bigg{/}\left(1-\frac{z}{n}+\frac{1}{2n}\right)\right]$
is a sufficient statistic for $r$. Noting that, under $H_{0}$,
$z\sim\mathcal{U}\\{1,n\\}$, the statistic has a useful property when applied
to the sample of ranks.

######  Claim 2.

Under $H_{0}$, $E[(t(z)]=0$.

###### Proof.

The proof of Claim 2 is in Appendix C. ∎

It is straightforward to show that $t(z)=0$ only when $z=\frac{n+1}{2}$, which
is the expected value of a discrete uniform distribution on the interval
$\\{1,n\\}$. Thus, the sufficient statistic preserves this property in its
expectation under the null. In the context of the sample of ranks for subject
$i$ in group $g$, $r$ is the true order of subject $i$’s curve within the
sample, under $H_{0}$. Subsequently, we summarize each $z_{g,i}(s)$ over $s$
using the sufficient statistic

$\displaystyle
t\left[z_{g,i}(\mathcal{S})\right]=\frac{1}{S}\sum_{s=1}^{S}\log\left[\frac{\frac{z_{g,i}(s)}{n}-\frac{1}{2n}}{1-\frac{z_{g,i}(s)}{n}+\frac{1}{2n}}\right].$
(4)

Our testing procedure for grouped functional data proceeds as follows: 1)
preprocess the data using FACE; 2) generate ranks at each time point $s$
ignoring group assignment; 3) calculate $t[z_{g,i}(\mathcal{S})]$; 4) perform
the relevant test (MWW or KW), depending on $G$. The doubly ranked MWW test
statistic is

$\displaystyle
T_{DR}^{+}=\sum_{j=1}^{n_{2}}R\left(t\left[z_{2,i}(\mathcal{S})\right]\right)-\frac{n_{2}(n_{2}+1)}{2}.$

For $G\geq 3$, the doubly ranked KW test statistic is

$\displaystyle
H_{DR}=\frac{12}{n(n+1)}\sum_{g=1}^{G}n_{g}\left[\overline{R\left(t\left[z_{g,i}(\mathcal{S})\right]\right)}-\frac{n+1}{2}\right]^{2},$

where
$\overline{R\left(t\left[z_{g,i}(\mathcal{S})\right]\right)}=\frac{1}{n_{g}}\sum_{i=1}^{n_{g}}R\left(t\left[z_{g,i}(\mathcal{S})\right]\right)$.
Because the sufficient statistic is a one-to-one transformation of the ranks,
when $S=1$, the doubly ranked tests will be equivalent to their univariate
analogue, i.e. the standard MWW or KW test. The statistic
$t\left[z_{g,i}(\mathcal{S})\right]$ is a univariate score for each subject
which use as the “new” data in each backend test. Thus, the distributional
results from Section 2 hold for both the doubly ranked MWW and KW tests. That
is, $T_{DR}^{+}$ is approximately normal while $H_{DR}$ is asymptotically
$\chi_{G-1}^{2}$. Implementation of our method is available for download as an
R package, titled runDRT, at https://github.com/markjmeyer/runDRT.

##  4 Simulation Study

In our empirical study, we assume a balanced design with groups of size
$n_{1}=n_{2}=10,25$ and 50 and, in the KW case, $n_{3}$ defined similarly. The
total sample sizes are $n=20,50,$ and 100 (MWW) and $n=30,75,$ and 150 (KW).
The doubly ranked MWW test is for two groups and we consider the three group
case for the doubly ranked KW test. To simulate the functional data, we use
the approach employed by Chakraborty and Chaudhuri (2015) who generate data in
the MWW setting using a Karhunen-Loève expansion of the form

$\displaystyle
X_{1,i}(s)=\sum_{k=1}^{K}\sqrt{2}\left[(k-0.5)\pi\right]^{-1}Z_{ik}\sin\left[(k-0.5)\pi
s\right],$

for those in group $X$. The value $Z_{ik}$ is a random variable which we
generate from one of three distributions: Gaussian, $t_{2}$ (a $t$
distribution with 2 degrees of freedom), and log-normal. Formally, the value
of $K$ could be infinite. In simulation, we take it to be large, $K=1000$
basis functions. To generate the data in the other group(s), we use:

$\displaystyle X_{g,j}(s)=$
$\displaystyle\mu(s)+\sum_{k=1}^{K}\sqrt{2}\left[(k-0.5)\pi\right]^{-1}Z_{jk}\sin\left[(k-0.5)\pi
s\right],$

where $\mu(s)$ is a time-dependent function for $s\in[0,1]$, $g=2$ for the MWW
case, and $g=2,3$ for the KW setting. These generative models produce
functional data with high levels of autocorrelation, see Section 2 of
Supplement A (Meyer, 2024a).

We consider three different functional forms for $\mu(s)$:

$\displaystyle\mu_{1}(s)=\xi s,\ $ $\displaystyle\mu_{2}(s)=\xi 4s(1-s),\text{
and}$ $\displaystyle\mu_{3}(s)=\xi m^{-1}$ $\displaystyle
B(2,6)^{-1}s^{2-1}(1-s)^{6-1},$

for $m$ equal to $\max_{s}B(2,6)^{-1}s^{2-1}(1-s)^{6-1}$ and $B(2,6)$ equal to
the Beta function evaluated at 2 and 6. The parameter $\xi$ is set to zero to
evaluate Type I Error and incremented up to 3, by steps of size 0.12, to
evaluate power. The length of the sampling grid, $s$, varies from sparser to
denser: $S=40,80,160,$ and 320 and is taken to be equally spaced over the unit
interval. We limit our direct comparison to the depth-based method, since it
alone has publicly available code. All curves are preprocessed using FACE,
retaining 99% of the variability. Code to generate the empirical studies is in
Supplement B (Meyer, 2024b).

Table 1: Type I Error for functional Mann-Whitney-Wilcoxon tests by distribution of $Z_{k}$ and sample size when $S=40$ and 320. Each value is based on 10000 simulated datasets. G stands for Gaussian, T for $t_{2}$, and L for log-normal.  $S$ | $Z_{k}$ | $n$ | Doubly | Depth-based
---|---|---|---|---
Ranked
40 | G | 20 | 0.042 | 0.031
|  | 50 | 0.052 | 0.041
|  | 100 | 0.050 | 0.049
| T | 20 | 0.042 | 0.032
|  | 50 | 0.048 | 0.044
|  | 100 | 0.047 | 0.050
| L | 20 | 0.041 | 0.032
|  | 50 | 0.051 | 0.042
|  | 100 | 0.052 | 0.049
320 | G | 20 | 0.041 | 0.031
|  | 50 | 0.052 | 0.041
|  | 100 | 0.049 | 0.049
| T | 20 | 0.041 | 0.031
|  | 50 | 0.048 | 0.045
|  | 100 | 0.049 | 0.051
| L | 20 | 0.042 | 0.033
|  | 50 | 0.050 | 0.045
|  | 100 | 0.050 | 0.048
Table 2: Type I Error for doubly ranked Kruskal-Wallis tests by distribution of $Z_{k}$ and sample size when $S=40$ and 320. Each value is based on 10000 simulated datasets. G stands for Gaussian, T for $t_{2}$, and L for log-normal.  $S$ | $n$ | $Z_{k}$
---|---|---
G | T | L
40 | 30 | 0.051 | 0.046 | 0.050
| 75 | 0.049 | 0.050 | 0.049
| 150 | 0.051 | 0.048 | 0.051
320 | 30 | 0.052 | 0.046 | 0.049
| 75 | 0.048 | 0.051 | 0.049
| 150 | 0.051 | 0.048 | 0.051

Type I error estimates for the both the doubly ranked MWW test and the depth-
based test are presented in Table 1 while Table 2 contains type I error for
the doubly ranked KW test. Because type I error did not vary significantly by
$S$, both tables contain values under the sparsest sampling density ($S=40$)
and are broken down by distribution of $Z_{k}$ and total sample size, $n$.
Each empirical type I error rate is based on 10000 simulated datasets. Testing
was performed at the $\alpha=0.05$ level. The doubly ranked MWW test achieves
or is within Monte Carlo error of the nominal $\alpha$-level for the larger
sample sizes ($n=50,100$), regardless of the distribution of $Z_{k}$ (Table
1). Under smaller sample sizes, the empirical type I error rate is below
nominal, near $0.041$ for each $Z_{k}$ (Table 1). Type I error is more stable
in the doubly ranked KW test simulation with all but one setting either
achieving nominal or within Monte Carlo error of it; when $Z_{k}\sim t_{2}$,
the error is under nominal, 0.046 (Table 2). Notably, from Table 1, the error
rates are similar between doubly ranked and depth-based MWW tests with depth-
based test error rates tending to be under nominal for $n=20$ ($0.031$ to
$0.032$) and nearing or achieving nominal as $n$ increases. Additional results
for type I error under other sampling densities ($S=80,160$ and 320) are in
Section 2 of Supplement A (Meyer, 2024a).

Figure 2: Power curves for functional Mann-Whitney-Wilcoxon tests under the
sparsest sampling density, $S=40$. Rows index varying distributions of
$Z_{k}$, columns index values of $\mu(s)$. Curves for the doubly ranked test
are in shades of blue, curves for depth-based tests are in shades of tan.
Solid curves are for when $n=100$, dashed curves for when $n=50$, and dotted-
dashed curves for when $n=20$.

Power curves are based on 5000 simulated datasets for each combination of
$\xi$, $S$, $n$, $\mu(s)$, and distribution of $Z_{k}$. Figure 2 presents the
curves under the sparsest sampling density ($S=40$) since increasing the
sampling density does not significantly impact the power of the doubly ranked
MWW test. Solid curves denote the settings where $n=100$, dashed were $n=50$,
and dotted-dashed for when $n=20$. Curves for the doubly ranked MWW test are
in blue while curves for the depth-based MWW test are in tan. The rows of
Figure 2 index the distribution of $Z_{k}$ and the columns index the varying
functional forms of $\mu(s)$. Overall, power increases for both methods as $n$
increases, no matter the distribution of $Z_{k}$ and the form of $\mu(s)$.
Comparing the two and at a given $n$, the doubly ranked MWW test has
consistently higher power than the depth-based MWW test. Similar power curves
to those in Figure 2 when $S=80,160$ or 320 are in Section 2 of Supplement A
(Meyer, 2024a).

Figure 3: Power curves for doubly ranked Kruskal-Wallis tests under the
sparsest sampling density, $S=40$. Rows index varying distributions of
$Z_{k}$, columns index values of $\mu(s)$. Solid curves are for when $n=100$,
dashed curves for when $n=50$, and dotted-dashed curves for when $n=20$.

Figure 3 contains the power curves for the doubly ranked KW test when $S=40$.
As with the MWW tests, the curves are based on 5000 simulated datasets per
combination of $\xi$, $n$, $\mu(s)$, and distribution of $Z_{k}$. Solid curves
now denote the settings where $n=150$, dashed where $n=75$, and dotted-dashed
for when $n=30$. Overall, power increases as $n$ increases across $Z_{k}$ and
forms of $\mu(s)$. The form of $\mu(s)$ does have a slight impact on power
when $n$ is small but the distribution of $Z_{k}$ does not appear to have much
of an impact on power. Additional power curves when $S=80,160,$ and 320 are in
Section 3 of Supplement A (Meyer, 2024a).

##  5 Data Illustrations

To illustrate both doubly ranked MWW and KW tests, we examine functional data
from three different substantive fields with varying numbers of observations
and differing lengths of functions. Two of the datasets are available in R
packages and the third will be made available upon request. Each dataset has
at least one factor of interest that may differentiate the outcome curves.
While we will draw light conclusions based on the findings, the purpose of
these illustrations is to demonstrate the use of doubly ranked tests and not
to elaborate on the scientific findings (or lack-thereof) that result from our
analysis. Results from all three data sources for various combinations of
outcomes and factors are in Table 3. All outcomes were first pre-processed
using FACE, retaining 99% of the variation within the curves. Graphs of the
functional curves from each data source are in Section 4 of Supplement A
(Meyer, 2024a).

Table 3:  Test statistics and $p$-values for doubly ranked tests applied to the three data illustrations: resin viscosity, Canadian weather, and COVID mobility. Temp. is short for temperature, Cur. for curing, Rot. for rotational, and Precip. for precipitation. Doubly ranked KW tests were performed for both Canadian weather outcomes and the COVID mobility factor MD, VA, & WV. The remaining tests were doubly ranked MWW tests.  Data Set | Outcome | Factor | Test | $p$-value
---|---|---|---|---
Statistic
Resin | Viscosity | Resin Temp. | 194 | $<0.001$
|  | Cur. Agent Temp. | 440 | 0.337
|  | Tool Temp. | 207 | $<0.001$
|  | Rot. Speed | 419.5 | 0.217
|  | Mass Flow | 507.5 | 0.957
Weather | Temp. | Region | 21.44 | $<0.001$
| Precip. | Region | 22.46 | $<0.001$
COVID | Requests | CO & UT | 132 | 0.001
|  | IA & MN | 1871 | $<0.001$
|  | MD, VA, & WV | 2.214 | 0.331

###  5.1 Resin Viscosity

The resin viscosity data comes from an experimental setting conducted at the
Technical University of Munich’s Institute for Carbon Composites and is freely
available in the R package FDBoost (Brockhaus et al., 2020). The data set
contains measurements of the viscosity of resin over the course of 838 seconds
with the goal assessing factors that influence the curing process in a resin
mold. The grid is not equally spaced due to technical reasons whereby
viscosity can be measured every two seconds at first but, due to hardening,
can only be measured every ten seconds after the 129th second. In total, 64
different molds were poured under five different experimental conditions:
temperature of resin, temperature of the curing agent, temperature of the
tools, rotational speed, and mass flow. Each condition is a binary factor,
with a “low” and “high” levels.

Figure 4: Blue curves were measured under the “low” temperature condition
while tan curves were measured under the “high” temperature condition.

From Table 3, we observe that in two of five factors the location parameters
in the distributions of the ranks differ significantly. Specifically, the
resin temperature ($W=194$, $p=<0.001$) and tool temperature ($W=207$,
$p<0.001$) have $p$-values below nominal. This suggests that the ranks of the
curves that describe the viscosity of the resin differ significantly between
the groups defined by these factors. To determine which group is higher or
lower, we turn to a graphical assessment. When considering the tool
temperature, Figure 4 displays clear clustering in the curves with low
temperature tools appearing to induce slower curing than high temperature
tools. This separation lessens around 400 seconds where the sets of curves
begin to overlap. The analysis by temperature of resin captures greater
separation after 400 seconds which can be seen in the additional figures in
Section 4 of Supplement A (Meyer, 2024a). The remaining factors did not
significantly distinguish the curves as can also be seen in the graphs in
Section 4 of Supplement A (Meyer, 2024a).

###  5.2 Canadian Weather

The Canadian weather data is a classic dataset for illustrating functional
data analysis techniques having appeared in texts by Ramsay and Silverman
(2005) and Ramsay et al. (2009) as well as in many articles that are too
numerous to list here. It contains daily average temperature (C) and
precipitation (mm) measurements, averaged over the years 1960 to 1994, taken
at 35 different sites across Canada. The factor of interest is the region from
which the measurement was taken: Arctic, Atlantic, Continental, or Pacific. It
is available in the R package fda (Ramsay et al., 2022). For the purpose of
this illustration, the goal of the analysis is to see if the doubly ranked KW
test can detect differences in temperature and precipitation curves by region
as might be expected given their geographic differences.

Figure 5: Blue curves denote the stations in the Arctic region, tan curves
indicate stations the Atlantic region, red curves for stations in the
Continental region, and green curves are the Pacific region.

Table 3 contains the results of the analysis for both outcomes. We see that
the doubly ranked KW returns significant tests at the nominal level for both
temperature curves ($X^{2}=21.44$, $p=<0.001$) and precipitation curves
($X^{2}=22.46$, $p<0.001$) differing by region. Formally, we would say the
location parameters of the distributions of the summarized ranks for
temperature and precipitation over time differ by region. To observe the
nature of the difference by group, we examine Figures 1 and 5. We see that the
temperature curves from the Arctic region are the coldest while those from the
Pacific region tend to be warmest. The precipitation curves show inverted
patterns based on region with the Pacific region tending to see more
precipitation during the winter months while the Continental and Arctic
regions experience more during the summer months.

###  5.3 COVID Mobility

The COVID mobility data was extracted from publicly available data provided by
Apple Inc. regarding daily mapping requests for driving, walking, and transit
directions during the onset of the COVID-19 pandemic and for some time
thereafter. A cached version of the data is available from Gassen (2022). The
data itself represents the percent change in requests for directions from a
baseline date of January 13, 2020. Requests are categorized by whether they
were for driving, walking, or transit directions. In our analysis, we focus on
the percent change in daily requests for driving directions. For the United
States, the data is available at the county level which is our unit of
observation.

Rather than aligning the data by calendar time, we align it by the first date
of a major COVID policy implementation at the state level, e.g. shelter-in-
place or safer-at-home orders. Policy dates were obtained from the publicly
available COVID AMP data (Katz et al., 2023). The timeline is centered at the
date of the first COVID-policy implementation in each state and we take 30
days prior and 30 days post. The main objective of this analysis is to compare
different states in similar regions to each other, focusing on whether or not
there is a difference post-policy implementation. Specifically, we compare
county level percent changes in driving requests between counties in Iowa (IA)
and Minnesota (MN), Colorado (CO) and Utah (UT), and finally between Maryland
(MD), Virginia (VA), and West Virginia (WV).

Figure 6: Blue curves denote counties in Utah while tan curves denote counties
in Colorado.

When comparing the percent change in county-level requests for driving
directions between Colorado and Utah, we find a significant difference between
the groups of curves by state ($W=132$, $p=0.001$). Figure 6 shows that the
county-level curves have substantial overlap pre-initial policy
implementation. After implementation, however, Colorado counties tend to
experience a lower percent change in requests. Graphs for the remaining
comparisons are in Section 4 of Supplement A (Meyer, 2024a). We see a similar
trend when comparing Minnesota and Iowa with Minnesota counties tending to
experience lower percent change in driving requests which continues until
three to four weeks post-implementation. This is reflected in the doubly
ranked MWW test which suggests a significant difference in the location
parameters of distributions of the summarized ranks between Minnesota and Iowa
($W=1871$, $p<0.001$). In contrast to the first two comparisons, we lack
evidence to suggest a difference in the location parameters of the
distributions the ranks between Maryland, Virginia, and West Virginia
($X^{2}=2.214$, $p=0.331$).

##  6 Discussion

In this manuscript, we investigate two novel nonparametric tests for comparing
groups of functional curves: the doubly ranked MWW test and the doubly ranked
KW test. The tests arise from the same sufficient statistic derived using the
distribution of the ranks under $H_{0}$ over time. Treating the ranks from
individual curves as a sample over time and using the sufficient statistic to
summarize them allows us to reduce the dimensionality of the testing problem
and improve efficiency. While other approaches reduce dimensionality using a
summary score for the data, only ours uses a score that is constructed under
$H_{0}$. Doubly ranked tests are distribution free since they rely only on the
ranks of the sufficient statistics derived in Equation 4. We demonstrate
empirically that both the doubly ranked MWW and KW tests perform well in terms
of type I error and power for a variety of true error distributions under
varying curve types.

The test procedures we propose here are global tests of differences between
curves, as are the other methods by Hall and Van Keilegom (2007); López-
Pintado and Romo (2009); López-Pintado et al. (2010); Chakraborty and
Chaudhuri (2015) and Meléndez et al. (2021). With any global test, one must be
careful to examine the underlying data for where the differences lie. These
tests can instruct us that there is difference but cannot be used to identify
exactly where the difference is nor can they tell us the magnitude of
difference. For that, we would need to model the group location curves
themselves and rely on point-wise testing, ideally adjusted for multiple
testing. See, for example, some of the work on functional regression reviewed
in Morris (2015) and Greven and Scheipl (2017).

Point-wise testing is not goal of this work, however our method can be useful
as a global test in a variety of data contexts. In our data illustrations, we
demonstrate the use of the doubly ranked MWW and KW tests to study problems in
material science, climatology, and public health policy. These tests can be
performed as a main analysis or a first pass analysis in conjunction with
other modeling steps. Our tests are also easy to implement and easy for
practitioners to interpret—in particular, for those who are already familiar
with univariate MWW and KW tests. Doubly ranked testing is not limited to
functional data applications nor is it even limited to FACE-based preprocessed
data. The assumptions for the tests apply broadly to both functional and
multivariate data.

Our method assumes the sampling grid is the same for each subject. Subjects
may, however, have mistimed measurements and therefore asynchronous grids. A
full investigation of this data structure is left to future work. However, in
the current context, an FPCA that accounts for asynchronous grids can be
employed should such data arise. For example, the FPCA by smoothed covariance,
see Yao et al. (2005) or Di et al. (2009) among others, can be used in the
context of sparse longitudinal data of which asynchronous functional data is a
more general case. Such an approach could be used for asynchronously sampled
profiles.

##  Appendix A Derivation of the PMF

The full derivation of the exact distribution in Equations (1) and (2) begins
with the discrete PMF:

$\displaystyle P\bigg{[}Z_{(r)}=z\bigg{]}$
$\displaystyle=\sum_{k=r}^{n}{n\choose
k}\bigg{[}\left(\frac{z}{n}\right)^{k}\left(1-\frac{z}{n}\right)^{n-k}$
$\displaystyle-\left(\frac{z-1}{n}\right)^{k}\left(1-\frac{z-1}{n}\right)^{n-k}\bigg{]}$
$\displaystyle=\sum_{k=r}^{n}{n\choose
k}\left(\frac{z}{n}\right)^{k}\left(1-\frac{z}{n}\right)^{n-k}$
$\displaystyle-\sum_{k=r}^{n}{n\choose
k}\left(\frac{z-1}{n}\right)^{k}\left(1-\frac{z-1}{n}\right)^{n-k},$

from which we obtain Equation (1). To get to Equation (2), we note that the
summations in Equation (1) are equivalent to the CDF of the binomial
distribution and can be rewritten as regularized incomplete beta functions:

$\displaystyle P\bigg{[}Z_{(r)}=z\bigg{]}$
$\displaystyle=1-I_{1-z/n}[n-(r-1),r-1+1]$
$\displaystyle-\left[1-I_{1-(z-1)/n}\\{n-(r-1),r-1+1\\}\right],$

where $I_{q}(a,b)$ is as defined in Section 3. Regularized incomplete beta
functions have the property that $I_{q}(a,b)=1-I_{1-q}(b,a)$. After applying
this, we arrive at Equation (2). To obtain the approximate PMF in Equation
(3), we note that the CDF of the beta distribution is also a regularized
incomplete beta function giving the expression

$\displaystyle P\bigg{[}$ $\displaystyle Z_{(r)}=z\bigg{]}$
$\displaystyle=\int_{0}^{\frac{z}{n}}\frac{1}{B(r,n-r+1)}\left(t\right)^{r-1}\left(1-t\right)^{n-r+1-1}dt$
$\displaystyle-\int_{0}^{\frac{z-1}{n}}\frac{1}{B(r,n-r+1)}\left(t\right)^{r-1}\left(1-t\right)^{n-r+1-1}dt$
$\displaystyle=\frac{1}{B(r,n-r+1)}\bigg{[}\int_{0}^{\frac{z}{n}}t^{r-1}\left(1-t\right)^{n-r+1-1}dt$
$\displaystyle-\int_{0}^{\frac{z-1}{n}}t^{r-1}\left(1-t\right)^{n-r+1-1}dt\bigg{]}.$

Writing the beta function in terms of gamma functions and combining the
integrals, we arrive at Equation (3).

##  Appendix B Proof of Claim 1

###### Proof.

Consider the mass function in Equation (3) and rewrite it as

$\displaystyle
P\bigg{[}Z_{(r)}=z\bigg{]}\approx\frac{\Gamma(n+1)}{\Gamma(r)\Gamma(n-r+1)}\frac{1}{n}$
$\displaystyle\hskip
15.0pt\times\left(\frac{z}{n}-\frac{1}{2n}\right)^{r-1}\left(1-\frac{z}{n}+\frac{1}{2n}\right)^{n-r}$
$\displaystyle=\frac{\Gamma(n+1)}{\Gamma(r)\Gamma(n-r+1)}\frac{1}{n}$
$\displaystyle\hskip
15.0pt\times\exp\bigg{[}n\log\left(1-\frac{z}{n}+\frac{1}{2n}\right)-\log\left(\frac{z}{n}-\frac{1}{2n}\right)\bigg{]}$
$\displaystyle\hskip
15.0pt\times\exp\bigg{[}r\log\left\\{\left(\frac{z}{n}-\frac{1}{2n}\right)\bigg{/}\left(1-\frac{z}{n}+\frac{1}{2n}\right)\right\\}\bigg{]}$

Noting that $n$ is fixed, define the following functions:

$\displaystyle c(r)$
$\displaystyle=\frac{1}{n}\frac{\Gamma(n+1)}{\Gamma(r)\Gamma(n-r+1)},\
w(r)=r,$ $\displaystyle h(z)$
$\displaystyle=\exp\bigg{[}n\log\left(1-\frac{z}{n}+\frac{1}{2n}\right)-\log\left(\frac{z}{n}-\frac{1}{2n}\right)\bigg{]},$
and $\displaystyle t(z)$
$\displaystyle=\log\left\\{\left(\frac{z}{n}-\frac{1}{2n}\right)\bigg{/}\left(1-\frac{z}{n}+\frac{1}{2n}\right)\right\\}.$

The mass function can then be written as

$\displaystyle P\bigg{[}Z_{(r)}=z\bigg{]}$
$\displaystyle=h(z)c(r)\exp\left[w(r)t(z)\right],$

which is the form of an exponential family. ∎

##  Appendix C Proof of Claim 2

###### Proof.

First note the mass function of a discrete uniform random variable over
$\\{1,n\\}$ is $P(Z=z)=\frac{1}{n}\ \forall\ z\in\\{1,n\\}$,
$n\in\mathbb{Z}^{+}$. The expectation of $t(z)$ is

$\displaystyle E\left[t(z)\right]$
$\displaystyle=E\left[\log\left\\{\left(\frac{z}{n}-\frac{1}{2n}\right)\bigg{/}\left(1-\frac{z}{n}+\frac{1}{2n}\right)\right\\}\right]$
$\displaystyle=\sum_{z=1}^{n}\log\left\\{\left(\frac{z}{n}-\frac{1}{2n}\right)\bigg{/}\left(1-\frac{z}{n}+\frac{1}{2n}\right)\right\\}\frac{1}{n}$
$\displaystyle=\frac{1}{n}\log\left\\{\prod_{z=1}^{n}\left(\frac{z}{n}-\frac{1}{2n}\right)\bigg{/}\left(1-\frac{z}{n}+\frac{1}{2n}\right)\right\\}$

Expanding out the first few terms and last few terms of the product, we see
that it is a telescoping product:

$\displaystyle\prod_{z=1}^{n}\left(\frac{z}{n}-\frac{1}{2n}\right)\bigg{/}\left(1-\frac{z}{n}+\frac{1}{2n}\right)=$
$\displaystyle\hskip
10.0pt\frac{\frac{1}{n}-\frac{1}{2n}}{1-\frac{1}{n}+\frac{1}{2n}}\times\frac{\frac{2}{n}-\frac{1}{2n}}{1-\frac{2}{n}+\frac{1}{2n}}\times\frac{\frac{3}{n}-\frac{1}{2n}}{1-\frac{3}{n}+\frac{1}{2n}}\times\cdots$
$\displaystyle\hskip
10.0pt\cdots\times\frac{\frac{n-2}{n}-\frac{1}{2n}}{1-\frac{n-2}{n}+\frac{1}{2n}}\times\frac{\frac{n-1}{n}-\frac{1}{2n}}{1-\frac{n-1}{n}+\frac{1}{2n}}\times\frac{\frac{n}{n}-\frac{1}{2n}}{1-\frac{n}{n}+\frac{1}{2n}}$
$\displaystyle=\frac{\frac{1}{2n}}{1-\frac{1}{2n}}\times\frac{\frac{1}{n}+\frac{1}{n}-\frac{1}{2n}}{1-\frac{1}{n}-\frac{1}{n}+\frac{1}{2n}}\times\frac{\frac{1}{n}+\frac{2}{n}-\frac{1}{2n}}{1-\frac{1}{n}-\frac{2}{n}+\frac{1}{2n}}\times\cdots$
$\displaystyle\hskip
10.0pt\cdots\times\frac{1-\frac{2}{n}-\frac{1}{2n}}{1-1+\frac{2}{n}+\frac{1}{2n}}\times\frac{1-\frac{1}{n}-\frac{1}{2n}}{\frac{1}{n}+\frac{1}{2n}}\times\frac{1-\frac{1}{2n}}{\frac{1}{2n}}$
$\displaystyle=\frac{\frac{1}{2n}}{1-\frac{1}{2n}}\times\frac{\frac{1}{n}+\frac{1}{2n}}{1-\frac{1}{n}-\frac{1}{2n}}\times\frac{\frac{1}{n}+\frac{2}{n}-\frac{1}{2n}}{1-\frac{2}{n}-\frac{1}{2n}}\times\cdots$
$\displaystyle\hskip
10.0pt\cdots\times\frac{1-\frac{2}{n}-\frac{1}{2n}}{\frac{1}{n}+\frac{2}{n}-\frac{1}{2n}}\times\frac{1-\frac{1}{n}-\frac{1}{2n}}{\frac{1}{n}+\frac{1}{2n}}\times\frac{1-\frac{1}{2n}}{\frac{1}{2n}}$

Examining the last two lines, we see that the first fraction is the reciprocal
of the last fraction, the second fraction is the reciprocal of the second to
last fraction, the third fraction is the reciprocal of the third to last
fraction, and so forth. When $n$ is even, the first $\frac{n}{2}$ fractions
will cancel with the last $\frac{n}{2}$ fractions.

When $n$ is odd, the first $\left\lfloor\frac{n}{2}-1\right\rfloor$ fractions
will cancel with the last $\left\lfloor\frac{n}{2}-1\right\rfloor$ fractions.
The remaining fraction will be the $z=\frac{n+1}{2}$ case which equals 1.
Thus, for any $n\in\mathbb{Z}^{+}$, the product will equal 1. Substituting, we
have

$\displaystyle E\left[t(z)\right]$
$\displaystyle=\frac{1}{n}\log\left\\{\prod_{z=1}^{n}\left(\frac{z}{n}-\frac{1}{2n}\right)\bigg{/}\left(1-\frac{z}{n}+\frac{1}{2n}\right)\right\\}$
$\displaystyle=\frac{1}{n}\log(1)=0$

Thus, the expected value of the sufficient statistic is zero. ∎

## References

  * Abramowicz et al. [2018] K. Abramowicz, C. K. Häger, A. Pini, L. Schelin, S. Sjöstedt de Luna, and S. Vantini.  Nonparametric inference for functional-on-scalar linear models applied to knee kinematic hop data after injury of the anterior cruciate ligament.  _Scandinavian Journal of Statistics_ , 45:1036–1061, 2018\.  doi: 10.1111/sjos.12333.
  * Alemán-Gómez et al. [2022] Y. Alemán-Gómez, A. Arribas-Gil, M. Desco, A. Elías, and J. Romo.  Depthgram: Visualizing outliers in high-dimensional functional data with application to fmri data exploration.  _Statistics in Medicine_ , 41:2005–2024, 2022.  doi: https://doi.org/10.1002/sim.9342.
  * Berrett et al. [2021] T. B. Berrett, I. Kontoyiannis, and R. J. Samworth.  Optimal rates for independence testing via U-statistic permutation tests.  _The Annals of Statistics_ , 49:2457–2490, 2021.  doi: 10.1214/20-AOS2041.
  * Brockhaus et al. [2020] S. Brockhaus, D. Rügamer, and S. Greven.  Boosting functional regression models with FDboost.  _Journal of Statistical Software_ , 94(10):1–50, 2020.  doi: 10.18637/jss.v094.i10.
  * Chakraborty and Chaudhuri [2015] A. Chakraborty and P. Chaudhuri.  A wilcoxon–mann–whitney-type test for infinite-dimensional data.  _Biometrika_ , 102:239–246, 2015.  doi: 10.1093/biomet/asu072.
  * Cuesta-Albertos and Febrero-Bande [2010] J. A. Cuesta-Albertos and M. Febrero-Bande.  A simple multiway ANOVA for functional data.  _Test_ , 19:537–557, 2010.  doi: 10.1007/s11749-010-0185-3.
  * Dai and Genton [2018] W. Dai and M. G. Genton.  Multivariate functional data visualization and outlier detection.  _Journal of Computational and Graphical Statistics_ , 27:923–934, 2018.  doi: 10.1080/10618600.2018.1473781.
  * Dai et al. [2020] W. Dai, T. Mrkvi$\check{\text{c}}$ka, Y. Sun, and M. G. Genton.  Functional outlier detection and taxonomy by sequential transformations.  _Computational Statistics & Data Analysis_, 149:106960, 2020.  doi: https://doi.org/10.1016/j.csda.2020.106960.
  * Di et al. [2009] C.-Z. Di, D. M. Crainiceanu, B. S. Caffo, and N. M. Punjabi.  Multilevel functional principal component analysis.  _The Annals of Applied Statistics_ , 3:458–488, 2009.  doi: 10.1214/08-AOAS206.
  * Gassen [2022] J. Gassen.  Download Apple mobility trend reports data.  _tidycovid19_ https://joachim-gassen.github.io/tidycovid19/reference/download_apple_mtr_data.html, 2022\.  Accessed June, 2023.
  * Genton et al. [2014] M. G. Genton, C. Johnson, K. Potter, G. Stenchikov, and Y. Sun.  Surface boxplots.  _Stat_ , 3:1–11, 2014.  doi: https://doi.org/10.1002/sta4.39.
  * Goldsmith et al. [2022] J. Goldsmith, F. Scheipl, L. Huang, J. Wrobel, C. Di, J. Gellar, J. Harezlak, M. W. McLean, B. Swihart, L. Xiao, C. Crainiceanu, and P. T. Reiss.  _refund: Regression with Functional Data_ , 2022.  URL https://CRAN.R-project.org/package=refund.  R package version 0.1-28.
  * Greven and Scheipl [2017] S. Greven and F. Scheipl.  A general framework for functional regression modelling.  _Statistical Modelling_ , 17:1–35, 2017.  doi: 10.1177/1471082X16681317.
  * Hall and Van Keilegom [2007] P. Hall and I. Van Keilegom.  Two sample tests in functional data analysis starting from discrete data.  _Statistica Sinica_ , 17:1511–1531, 2007.  URL https://www.jstor.org/stable/24307686.
  * Hollander and Wolfe [1999] M. Hollander and D. A. Wolfe.  _Nonparametric statistical methods_.  John Wiley & Sons, Inc., New York, 1999.
  * Huang and Sun [2019] H. Huang and Y. Sun.  A decomposition of total variation depth for understanding functional outliers.  _Technometrics_ , 61:445–458, 2019.  doi: 10.1080/00401706.2019.1574241.
  * Katz et al. [2023] R. Katz, M. Boyce, and E. Graedon.  Visualizing the impact of policies on COVID response.  _COVID AMP_ https://covidamp.org/, 2023.  Accessed June, 2023.
  * Kloke and McKean [2015] J. Kloke and J. W. McKean.  _Nonparametric Statistical Methods using R_.  CRC Press, Taylor & Francis Group, Boca Raton, Florida, 2015.
  * Kruskal and Wallis [1952] W. H. Kruskal and W. A. Wallis.  Use of ranks in one-criterion variance analysis.  _Journal of the American Statistical Association_ , 47(260):583–621, 1952.  doi: 10.2307/2280779.
  * López-Pintado and Romo [2009] S. López-Pintado and J. Romo.  On the concept of depth for functional data.  _Journal of the American Statistical Association_ , 104:718–734, 2009.  doi: 10.1198/jasa.2009.0108.
  * López-Pintado and Romo [2011] S. López-Pintado and J. Romo.  A half-region depth for functional data.  _Computational Statistics & Data Analysis_, 55:1679–1695, 2011.  doi: 10.1016/j.csda.2010.10.024.
  * López-Pintado and Wrobel [2017] S. López-Pintado and J. Wrobel.  Robust non-parametric tests for imaging data based on data depth.  _Stat_ , 6:405–419, 2017.  doi: 10.1002/sta4.168.
  * López-Pintado et al. [2010] S. López-Pintado, J. Romo, and A. Torrente.  Robust depth-based tools for the analysis of gene expression data.  _Biostatistics_ , 11:254–264, 2010.  doi: 10.1093/biostatistics/kxp056.
  * Mann and Whitney [1947] H. B. Mann and D. R. Whitney.  On a test of whether one of two random variables is stochastically larger than the other.  _The Annals of Mathematical Statistics_ , 18:50–60, 1947\.  doi: 10.1214/aoms/1177730491.
  * Meléndez et al. [2021] R. Meléndez, R. Giraldo, and V. Leiva.  Sign, Wilcoxon and Mann-Whitney tests for functional data: An approach based on random projections.  _Mathematics_ , 9:44, 2021.  doi: 10.3390/math9010044.
  * Meyer [2024a] M. J. Meyer.  Supplement A to “Doubly ranked tests for grouped functional data”.  2024a.  doi: 10.1214/.
  * Meyer [2024b] M. J. Meyer.  Supplement B to “Doubly ranked tests for grouped functional data”.  2024b.  doi: 10.1214/.
  * Morris [2015] J. S. Morris.  Functional regression.  _Annual Review of Statistics and Its Application_ , 2:321–359, 2015.  doi: 10.1146/annurev-statistics-010814-020413.
  * Pomann et al. [2016] G.-M. Pomann, A.-M. Staicu, and S. Ghosh.  A two-sample distribution-free test for functional data with application to a diffusion tensor imaging study of multiple sclerosis.  _Journal of the Royal Statistical Society, Series C_ , 65:395–414, 2016.  doi: 10.1111/rssc.12130.
  * R Core Team [2022] R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2022.  URL https://www.R-project.org/.
  * Ramsay and Silverman [2005] J. O. Ramsay and B. W. Silverman.  _Functional Data Analysis_.  Springer, New York, 2${}^{\text{nd}}$ edition, 2005.
  * Ramsay et al. [2009] J. O. Ramsay, G. Hooker, and S. Graves.  _Functional Data Analysis in R and Matlab_.  Springer, New York, 2009.
  * Ramsay et al. [2022] J. O. Ramsay, S. Graves, and G. Hooker.  _fda: Functional Data Analysis_ , 2022.  URL https://CRAN.R-project.org/package=fda.  R package version 6.0.5.
  * Sguera and López-Pintado [2021] C. Sguera and S. López-Pintado.  A notion of depth for sparse functional data.  _Test_ , 30:630–649, 2021.  doi: 10.1007/s11749-020-00734-y.
  * Sun and Genton [2011] Y. Sun and M. G. Genton.  Functional boxplots.  _Journal of Computational and Graphical Statistics_ , 20:316–334, 2011.  doi: 10.1198/jcgs.2011.09224.
  * Sun and Genton [2012] Y. Sun and M. G. Genton.  Adjusted functional boxplots for spatio-temporal data visualization and outlier detection.  _Environmetrics_ , 23:54–64, 2012.  doi: 10.1002/env.1136.
  * Sun et al. [2012] Y. Sun, M. G. Genton, and D. W. Nychka.  Exact fast computation of band depth for large functional datasets: How quickly can one million curves be ranked?  _Stat_ , 1:68–74, 2012.  doi: https://doi.org/10.1002/sta4.8.
  * Wang et al. [2016] J.-L. Wang, J.-M. Chiou, and H.-G. Müller.  Functional data analysis.  _Annual Review of Statistics and Its Application_ , 3:257–295, 2016.  doi: 10.1146/annurev-statistics-041715-033624.
  * Wilcoxon [1945] F. Wilcoxon.  Individual comparisons by ranking methods.  _Biometrics Bulletin_ , 1:80–83, 1945.  doi: 10.2307/3001968.
  * Wrobel et al. [2016] J. Wrobel, S. Y. Park, A.-M. Staicu, and J. Goldsmith.  Interactive graphics for functional data analyses.  _Stat_ , 5:108–118, 2016.  doi: https://doi.org/10.1002/sta4.109.
  * Xiao et al. [2013] L. Xiao, Y. Li, and D. Ruppert.  Fast bivariate p-splines: the sandwich smoother.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 75:577–599, 2013.  doi: https://doi.org/10.1111/rssb.12007.
  * Xiao et al. [2016] L. Xiao, V. Zipunnikov, D. Ruppert, and C. Crainiceanu.  Fast covariance estimation for high-dimensional functional data.  _Statistics and Computing_ , 26:409–421, 2016.  doi: 10.1007/s11222-014-9485-x.
  * Yao et al. [2005] F. Yao, H.-G. Müller, and J.-L. Wang.  Functional linear regression analysis for longitudinal data.  _Annals of Statistics_ , 33:2873–2903, 2005.  doi: 10.1214/009053605000000660.