Continuous testing for Poisson process intensities: A new perspective on scanning statistics

TL;DR

This paper introduces a continuous testing framework for Poisson process intensities, enabling rigorous, non-asymptotic control of error rates and extending traditional scanning window methods with new test statistics and Monte Carlo estimation.

Contribution

It presents a novel continuous testing approach for Poisson processes that controls error rates non-asymptotically and introduces kernel-based test statistics.

Findings

Framework effectively controls family-wise error rate and false discovery rate.

Method successfully applied to Neurosciences and Genomics data.

Provides a Monte Carlo-based estimation for p-value processes.

Abstract

We propose a novel continuous testing framework to test the intensities of Poisson Processes. This framework allows a rigorous definition of the complete testing procedure, from an infinite number of hypothesis to joint error rates. Our work extends traditional procedures based on scanning windows, by controlling the family-wise error rate and the false discovery rate in a non-asymptotic manner and in a continuous way. The decision rule is based on a \pvalue process that can be estimated by a Monte-Carlo procedure. We also propose new test statistics based on kernels. Our method is applied in Neurosciences and Genomics through the standard test of homogeneity, and the two-sample test.

Tables

Table 1. Table 1: Summary of the required computations to implement the continuous tests in practice for the homogeneity and two-sample test.

Procedure	test statistic	computations	homogeneity test	two-sample test
$\min$-$p$	count	$p_{\eta }$	exact
		full procedure	1-step MC
	kernel	$p_{\eta }$	1-step MC
		full procedure	2-step MC
wBH	count	$p_{\eta }$	exact
		full procedure	exact
	kernel	$p_{\eta }$	1-step MC
		full procedure	1-step MC

Equations

$$\lambda ={\int }_0^1\nu (s)ds.$$

$$\forall t\in [0,1],\theta (t)=\frac{\nu (t)}{{\int }_0^1\nu (s)ds}-1.$$

$${\epsilon }_T|N\sim 2\mathcal{B}\left(\frac{\theta (T)+1}{2}\right)-1,$$

$$\forall t\in [0,1],\theta (t)=\frac{{\nu }_A(t)-{\nu }_B(t)}{{\nu }_A(t)+{\nu }_B(t)}.$$

$$\displaystyle\forall t\in[0,1],\,\,H_{0,t}:\big{\{}\theta(t)=0\big{\}}$$

$$\displaystyle\text{or }\forall t\in[0,1],\,\,H_{0,t}:\big{\{}\theta(t)=0\big{\}}$$

$$J_{0}^{\eta}(\theta)=\big{\{}x\in\mathcal{X}_{\eta}\>:\>\forall t\in I_{\eta}(x),\theta(t)=0\big{\}}.$$

$$\forall x\in {\mathcal{X}}_{\eta },p_{\eta }(x)=F_{\mathcal{B}(N([0,1]),\frac{1}{2})}(S_{\eta }(x)),$$

$$\forall x\in {\mathcal{X}}_{\eta },p_{\eta }(x)=F_{\mathcal{P}(\eta \lambda )}(S_{\eta }(x))$$

$$\forall x\in {\mathcal{X}}_{\eta },p_{\eta }(x)=F_{\mathcal{B}(N(I_{\eta }(x)),\frac{1}{2})}(S_{\eta }(x)).$$

$$x\in {\mathcal{X}}_{\eta }\mapsto Y\cap I_{\eta }(x),$$

$$(\mathcal{P})$$

$$\forall x\in J_0^{\eta },F_{\mathfrak{N},x}(s)=P_{{\theta }_0}(S_{\eta }(x)\ge s|\mathfrak{N})=P_{\theta }(S_{\eta }(x)\ge s|\mathfrak{N}),$$

$$p_{\eta }(x)=F_{\mathfrak{N},x}(S_{\eta }(x)),x\in {\mathcal{X}}_{\eta }.$$

$$\mbox{if $\mathcal{H}_{0}\left\{I_{\eta}(x)\right\}$ is true, }\>\>\forall\alpha\in[0,1],P_{\theta,\lambda}({p}_{\eta}(x)\leq\alpha)\leq\alpha.$$

$${\mathcal{R}}_{\eta }(u):=\left\{x\in {\mathcal{X}}_{\eta }:p_{\eta }(x)\le u\right\},$$

$$\text{mbox}{FWER}_{\theta ,\lambda }^{\eta }({\mathcal{R}}_{\eta })=P_{\theta ,\lambda }\left(J_0^{\eta }\cap {\mathcal{R}}_{\eta }\ne \varnothing \right).$$

$$\text{mbox}{FDR}_{\theta ,\lambda }^{\eta }({\mathcal{R}}_{\eta })=E_{\theta ,\lambda }\left(\frac{\Lambda \left(J_0^{\eta }\cap {\mathcal{R}}_{\eta }\right)}{\Lambda \left({\mathcal{R}}_{\eta }\right)}\right),$$

$$\forall (\theta ,\lambda ),\text{mbox}{FWER}_{\theta ,\lambda }^{\eta }({\mathcal{R}}_{\eta }(u^{\alpha }))\le \alpha .$$

$$\begin{gathered} \left\{J_0^{\eta }\cap {\mathcal{R}}_{\eta }(u)\ne \varnothing \right\}=\left\{\exists x\in J_0^{\eta }:p_{\eta }(x)\le u\right\} \\ =\left\{\underset{x\in J_0^{\eta }}{\inf }\left\{p_{\eta }(x)\right\}\le u\right\}, \end{gathered}$$

$$\text{mbox}{FWER}_{\theta ,\lambda }^{\eta }({\mathcal{R}}_{\eta }(u))$$

$$\mbox{FWER}_{\theta,\lambda}^{\eta}(\mathcal{R}_{\eta}(u))\leq E_{\lambda}\left[P_{\theta_{0}}\left(\underset{x\in\mathcal{X}_{\eta}}{\inf}\left\{p_{\eta}(x)\right\}\leq u\>\big{|}\>\mathfrak{N}\right)\right].$$

$$q_{\eta }(x)=F_{\mathfrak{N}}^{\text{m}}\left(p_{\eta }(x)\right),x\in {\mathcal{X}}_{\eta },$$

$${\mathcal{R}}_{\eta }^{\text{m}}:=\{x\in {\mathcal{X}}_{\eta }:q_{\eta }(x)\le \alpha \}$$

$$\underset{x\in {\mathcal{X}}_{\eta }}{\inf }{p}_{\eta }(x)=F_{\mathcal{B}(N([0,1]),\frac{1}{2})}\left(\underset{x\in {\mathcal{X}}_{\eta }}{\sup }{S}_{\eta }(x)\right).$$

$$\forall (\theta ,\lambda ),\text{mbox}{FDR}_{\theta ,\lambda }^{\eta }(v^{\alpha })\le \alpha .$$

$$\begin{gathered} \text{mbox}{FDR}_{\theta ,\lambda }^{\eta }(v)=E_{\theta ,\lambda }\left(\frac{\Lambda \left(J_0^{\eta }\cap {\mathcal{R}}_{\eta }(v)\right)}{\Lambda \left({\mathcal{R}}_{\eta }(v)\right)}\right) \\ ={\int }_{J_0^{\eta }}{E}_{\theta ,\lambda }\left(\frac{1_{\left\{p_{\eta }(x)\le v\right\}}}{\Lambda \left({\mathcal{R}}_{\eta }(v)\right)}\right)d\Lambda (x). \end{gathered}$$

$$\frac{\Lambda \left({\mathcal{R}}_{\eta }(v)\right)}{\Lambda \left({\mathcal{X}}_{\eta }\right)}\ge \frac{v}{\alpha }$$

$$\text{mbox}{FDR}_{\theta ,\lambda }^{\eta }(v)\le \alpha {\int }_{J_0^{\eta }}{v}^{-1}{E}_{\theta ,\lambda }\left(\frac{1_{\left\{p_{\eta }(x)\le v\right\}}}{\Lambda \left({\mathcal{X}}_{\eta }\right)}\right)d\Lambda (x),$$

$$\text{mbox}{FDR}_{\theta ,\lambda }^{\eta }(v)$$

$$V^{\alpha }=\max \left\{v\ge 0:\frac{\Lambda \left({\mathcal{R}}_{\eta }(v)\right)}{\Lambda \left({\mathcal{X}}_{\eta }\right)}\ge v/\alpha \right\}.$$

Full text

1. 1 Introduction

2. 2 The continuous testing framework for Poisson Processes

3. 2.1 Homogeneity and two-sample conditional tests

4. 2.2 An infinite set of local null hypothesis

5. 2.3 A continuum of scanning windows

6. 2.4 A first simple construction of the $p$-value process

7. 2.5 Making continuous testing computationally tractable in general

8. 2.6 General definition of the $p$-value process

9. 3 False positive control in continuous time

10. 3.1 Definition of the multiplicity error rates

11. 3.2 A $\min$-$p\,$procedure to control the FWER in continuous time

12. 3.3 Link with the scan statistic framework.

13. 3.4 A weighted BH procedure to control the FDR in continuous time

14. 4 Kernel-based statistics for sliding windows

15. 4.1 On the limitations of the count-based statistic

16. 4.2 The kernel vision in the two-sample case

17. 5 Monte-carlo procedures for $p$-values estimation and adjustment

18. 5.1 The distribution of the single test statistic being explicit.

19. 5.2 The distribution of the single test statistic being not explicit

20. 5.3 Estimation of the $p$-value process by Monte-Carlo

21. 5.4 Two-step Monte-Carlo method for the $\min$-$p\,$procedure.

22. 6 Simulations

23. 6.1 Homogeneity test

24. 6.2 Two sample test

25. 7 Applications

26. 7.1 Do neurons respond to some stimulus ?

27. 7.2 Differential analysis of replication origins between cell lines

28. 8 Conclusion

29. A Supplementary material

30. A.1 Step-down improvement

31. A.2 Other kernel-based single test statistics

32. A.2.1 Homogeneity test statistics

33. A.2.2 Two-sample test statistics (one-sided)

34. B Proofs

35. B.1 Proof of Theorem 3.1

36. B.2 Proof of Theorem 3.5

37. B.3 Proof of Theorem 5.1

1 Introduction

Continuous testing has recently emerged as a suitable theoretical framework to test a possibly infinite number of hypothesis. It is especially adapted to situations where observed data come from an underlying continuously-indexed random process, like in the classical white noise model (Dümbgen and Spokoiny, 2001), or when data can be modeled by a point process (Blanchard et al., 2014). Our work being motivated by data coming from Neurosciences and Genomics, we focus on point processes that have been very effective to model occurrences of events through time or space in many fields. More than a broad range of applications, point processes have the advantage to be discrete by essence, which allows a drastic reduction of the continuum number of hypotheses, into to a random finite number of tests. This property was exploited by Blanchard et al. (2014) in a particular context, and we elaborate upon this work to propose a continuous testing framework that allows the rigorous definition of the complete testing procedure, from multiple hypothesis to error rates. Our procedure is very general and fully operational thanks to the use of sliding windows. This allows us to we revisit the old-tale of scanning statistics, and we extend their traditional possibilities and performance by providing the theoretical controls for the Family Wise Error Rate (FWER) as well as for the False Discovery Rate (FDR). We focus on Poisson processes for their remarkable properties, but the main definitions, methodological developments or computational tricks are valid for general point processes.

Scanning windows probably constitute the most common method to perform tests on Poisson processes for the detection of unexpected clusters of observations (Kulldorf, 1997; Chan and Zhang, 2007; Perone Pacifico et al., 2007; Siegmund et al., 2011) or to compare samples (Walther, 2010). It consists in computing a test statistic on overlapping windows of a given length that are shifted by all possible delays, which has the advantage of avoiding empirical discretization of the data (Grün, 1996; Tuleau-Malot et al., 2014). Then the decision rule is based on the so-called scan statistic, that refers to the window with most signal. Therefore these procedures rely on the distribution of the maximum of the test statistics under the full null, that can be either approximated (Naus, 1982; Robin et al., 2005; Fu et al., 2012; Rivera and Walther, 2013; Walther, 2010; Chan and Zhang, 2007) or learned by a randomization procedure (Kulldorff et al., 2005; Kulldorff et al., 2009; Arias-Castro et al., 2015), which is common in standard multiple testing procedures (Westfall and Young, 1993; Romano and Wolf, 2005a). In this work we adopt the conditional testing approach (Loader, 1991; Rivera and Walther, 2013) that has the advantage of getting rid of nuisance parameters. It consists in deriving tests that are conditional to the number of occurrence (for the homogeneity test) or to the position of occurrences (for the two-sample test). Also, this strategy is consistent with the strong connections between continuous testing and non-parametric statistics (Blanchard et al., 2014), since the conditional test may depend on the unknown intensity functions of the observed processes. Therefore, our strategy is also to go beyond standard parametric approaches, that are mainly based on observed counts or on a likelihood ratio statistic. As we will show, these parametric statistics do not fully exploit the spatial/temporal information that lie within the data, but only use integrated information. Inspired by recent results on non-parametric testing (Gretton et al., 2012; Fromont et al., 2013), we enrich the scan-statistics toolbox with kernel-based statistics to increase the power of detection of sharp alternatives.

Another limitation of current window-based strategies is that they correspond to global procedures that focus on a global null hypothesis. This resumes to ask whether the process is globally homogeneous, or if the samples are globally identical, which inevitably results in a global yes/no answer. However, with the increasing availability of more and more complex data, there is a need for flexible and local approaches to improve sensitivity for local structures. As we will see, continuous testing allows us to revisit the scanning windows framework, by defining a (possibly) uncountably infinite set of local hypothesis, and most importantly by providing a rigorous framework for defining joint error rates.

Indeed, the well known statistical challenge that arises when dealing with sliding windows is the appropriate control of Type-I errors, which is made difficult by the induced dependency structure. Current strategies mostly rely on the asymptotic approximation of the scan distribution (Kulldorf, 1997; Chan and Zhang, 2007), which has two advantages : $i)$ it somewhat avoids the dependency matter by focusing on the scan statistic only, $ii)$ it is computationally efficient. We rather adopt a non-asymptotic point of view, by considering the complete $p$-value process that becomes here a stochastic process . For some particular statistics (the count statistic for instance), this process is easily computed. When the statistic distribution is unknown we propose a Monte-Carlo procedure to estimate the $p$-value process. Then we propose a $\min$-$p\,$procedure in continuous time to control the FWER in a non-asymptotic manner, based on a new computationally efficient double Monte-Carlo scheme. Simulations show that this non asymptotic procedure outperforms the one of Chan and Zhang (2007) in terms of FWER control.

The control of the False Discovery Rate (FDR), that is the expected proportion of errors among rejections, is more intricate, and tentatives are more recent. Perone Pacifico et al. (2004) propose a continuous analogue to the FDR in the context of random fields, where the amount of null hypotheses is assessed by the Lebesgue measure. However, their (continuous) FDR controlling procedure is based upon a FWER guarantee. Their strategy has also been considered in a setting close to our homogeneity framework (Perone Pacifico et al., 2007). Siegmund et al. (2011) also investigated the FDR for scanning statistics, but only under two strong conditions (Poisson distribution for the number of false discoveries and independence between null and alternative statistics), none of which being satisfied in standard models. Moreover, their procedure eventually resumes to discretized $p$-values , and is based on the tail approximation of the scan statistic, which is not valid for other statistics other than the scan statistic. Thanks to the continuous testing framework, we propose a weighted Benjamini-Hochberg procedure in continuous time for Poisson processes testing, which outperforms on simulations Siegmund et al. (2011). This procedure corresponds to a continuous analogue of the well-known step-up procedure of Benjamini and Hochberg (1995), and extends the results of Blanchard et al. (2014) for continuous FDR in the context of local hypothesis testing.

The present work is inspired by two applications of Poisson processes: the analysis of spike trains data in Neurosciences and analysis of Genomic data coming from Next Generation Sequencing (NGS) technologies. In Neurosciences, Poisson processes are used to model spike trains, that is the succession of action potentials of a given neuron (Pipa et al., 2013). In particular inhomogeneous Poisson processes constitute a powerful model to understand rapid changes in firing rates induced by some stimulii (Kass et al., 2003). While the estimation of the Poisson process intensity has been already deeply studied (Shimazaki and Shinomoto, 2010), the questions asked in practice are often somehow different from a pure curve reconstruction matter: does the stimulus affect the spike apparition process (global hypothesis)? When (local) ? Does it depend on the type of stimulus and if so when exactly are the differences? In Genomics, next generation sequencing technologies (NGS) now allow the fine mapping of genomic features along chromosomes, and this spatial organization has become central to better understand genomes architecture and regulation. Here we consider the spatial organization of human replication origins that constitute the starting points of chromosomes duplication. To maintain the stability and integrity of genomes, replication origins are under a very strict spatio-temporal control, and part of their positioning has been shown to be associated with cell differentiation (Picard et al., 2014). Their recent mapping on the human genome has allowed the comparative analysis of replication origins maps between cell lines. Unfortunately, there is a lack of statistical framework to compare genomic maps of spatially organized genomic features, and we show how our procedure can be applied to the identification of replication starting sites that are linked to cell differentiation.

This article is organized as follows. First, we define the continuous testing framework, and we show how it can be easily handled in the Poisson process set-up (Section 2), with a particular focus on the test of homogeneity and the two-sample test. We also define the infinite set of local test hypotheses, the continuum of scanning windows to test them and the associated $p$-value stochastic process. To proceed, and for the sake of simplicity, we first focus on the standard count statistic that directly gives access to the $p$-values . Then two types of control are discussed in Section 3: the classical FWER and the more involved control of the FDR in continuous time. Both procedures result in an adjusted $p$-value process, which makes our method operational and easily implemented. To increase the power of our procedure we propose in Section 4 local kernel-based statistics. Those kernel statistics require non explicit conditional distributions to be transformed into $p$-values and we discuss Monte-Carlo procedures to estimate the $p$-values , with a new computationally efficient double Monte-Carlo step to control the FWER in Section 5. In Section 6 we show through simulations that the resulting procedures are more powerful than existing ones (in particular those based on unconditioned scanning statistics or on discretization). Section 7.2 shows that the procedure works extremely well on experimental data, issued from Neurosciences or Genomics. Our method is available in the form of the contest R-package, and offers many perspectives in the field of continuous testing.

2 The continuous testing framework for Poisson Processes

In the sequel a point process is a random countable set of points usually denoted by $N$. The corresponding point measure is then denoted by $dN$, that is $\sum_{T\in N}\delta_{T}$ with $\delta$ denoting the Dirac mass. For any interval $I$, $N(I)$ is a random variable that corresponds to the number of points of $N$ in $I$. We focus on processes in $\mathbb{R}$ (which can be thought as the DNA strand for genomic applications and time line for neuroscience applications) and we consider processes in [0,1] for the sake of simplicity.

2.1 Homogeneity and two-sample conditional tests

In order to propose a unified framework, we use the following notations: we observe a random set of points $Y$ with a distribution of the form $P_{\theta,\lambda}$, where $\lambda$ is a nuisance parameter, while $\theta$ represents the signal of interest. In both cases, the aim is to test whether $\theta$ is locally null.

Homogeneity test.

We observe $Y=N$, a point process modeled by a heterogeneous Poisson Process on $[0,1]$ with unknown intensity $\nu$ (non negative and in $L^{2}[0,1]$) with respect to the Lebesgue measure. We would like to detect the regions that are particularly rich or poor in occurrences of the process, with respect to its mean/stationary behavior. To do so, we compare the intensity $\nu$ of the observed process to the mean behavior of the intensity

[TABLE]

Then we rewrite the distribution of $Y$ in terms of two different parameters $\lambda$ and $\theta$, with $\theta$, defined by

[TABLE]

The signal of interest is $\theta$, whereas $\lambda$ may be seen as a nuisance parameter because its value is unknown under the null hypothesis. The homogeneity null hypothesis is therefore reduced to ”$\theta=0$”: it is a composite hypothesis since under the null hypothesis, the distribution can still vary according to $\lambda$. As a particular easier problem, one can also assume that $\lambda$ is known. The problem is then a classical goodness-of-fit test and the null hypothesis is not composite anymore.

Two-sample test.

We observe here $Y=(N_{A},N_{B})$, where $N_{A}$ and $N_{B}$ are two heterogeneous independent Poisson Processes on $[0,1]$ with intensities $\nu_{A}$ and $\nu_{B}$ (non negative and in $L^{2}[0,1]$) with respect to the Lebesgue measure. In the following, for any interval $I$, we denote $Y(I)=(N_{A}(I),N_{B}(I))$, and by a slight abuse of notation, (Y\cap I)=\big{(}N_{A}\cap I,N_{B}\cap I\big{)}. Our aim here is to detect the regions where $\nu_{A}\not=\nu_{B}$ (or, alternatively, $\nu_{A}>\nu_{B}$).

To re-parametrize the distribution of $Y$, one can provide a one-to-one correspondence between $Y$ and the couple $(N,\varepsilon)$, where $N=N_{A}\cup N_{B}$ is the joint process and where $\varepsilon=(\varepsilon_{T})_{T\in N}$ is a set of marks with values in $\{-1,1\}$ such that $\varepsilon_{T}=1$ if $T$ actually belongs to $N_{A}$ and $\varepsilon_{T}=-1$ if $T$ actually belongs to $N_{B}$ (Kingman, 1993). An important point here is to note that $N$ is actually a Poisson process with intensity $\lambda=\nu_{A}+\nu_{B}$ and that conditionally to $N$, the $\varepsilon_{T}$’s are independent variables with distribution

[TABLE]

with

[TABLE]

Once again, $\theta$ stands for the amount of signal in the data and $\lambda$ is a nuisance parameter. The two-sample null hypothesis is then once again summarized by ”$\theta=0$”. Since $\lambda$ is an unknown function, this hypothesis is composite. This case is thus much more complex than the homogeneity case, and calls for a non-parametric procedure.

Conditional distributions.

To avoid in both cases the use of specific unknown values of $\lambda$ to derive the distribution under the null hypothesis, our strategy is based on conditional testing (Loader, 1991; Rivera and Walther, 2013). Indeed, for the homogeneity test, conditionally on the total number of points, $N([0,1])=n$, the observed process is actually an $n$ i.i.d. sample of density $1+\theta$ that does not depend on $\lambda$ (Loader, 1991; Rivera and Walther, 2013). Similarly, in the two-sample case, given the positions of the joint process $N$, the remaining variability under the null hypothesis lies in the distribution of the marks $(\varepsilon_{T})_{T\in N}$ that only depends on $\theta$ (and not on $\lambda$ anymore). Thus we introduce the variable $\mathfrak{N}$, such that $\mathfrak{N}=N([0,1])$ in the homogeneity case, and $\mathfrak{N}=N$ in the two-sample case and our testing procedure relies on the distribution of $Y$ conditionally to $\mathfrak{N}$. Since the conditional distribution of $Y$ given $\mathfrak{N}$ does not depend on $\lambda$, it is denoted by $P_{\theta}$ in the sequel, whereas $E_{\lambda}$ refers to the expectation w.r.t. $\mathfrak{N}$.

2.2 An infinite set of local null hypothesis

The continuous testing procedure developed here aims at determining where the signal $\theta$ is non zero. The classical strategy consists in testing the full null hypothesis, which corresponds to the single question “Is the signal $\theta$ null over [0,1]?”. Therefore a unique hypothesis is tested and the test produces a binary answer. Here, we rather consider local null hypothesis such as,

[TABLE]

In practice, one-sided tests of homogeneity are usually considered to detect regions that are too rich in occurrences, whereas two-sided tests of the two-sample problem seem more relevant when there is a complete equivalence between both processes $N_{A}$ and $N_{B}$ from a modeling point of view.

Since information at the single point level seems difficult to catch in general, we focus on finding intervals or unions of intervals on which the null hypothesis holds. We denote $\mathcal{H}_{0}\left\{I\right\}$, the null hypothesis corresponding to $``\forall t\in I,\theta(t)=0"$. As a particular case, $\mathcal{H}_{0}\left\{[0,1]\right\}$ is the full null hypothesis and corresponds to $"\theta=\theta_{0}"$, where $\theta_{0}$ is the null function on $[0,1]$.

2.3 A continuum of scanning windows

Our procedure is based on scanning windows with a given resolution $\eta\in]0,1[$ that is fixed. We envision here all the possible continuum of windows with such given length, meaning that we are performing a whole continuum of tests, one for each window. To avoid any confusion in the sequel, $x$ always denotes a window center whereas $t$ always denotes a point, that is a possible value for the observations $Y$. Thus we have $x\in\mathcal{X}_{\eta}=[\eta/2,1-\eta/2]$ (while the set of possible values for $t$ is $[0,1]$). The scanning windows are denoted by $I_{\eta}(x)=(x-\eta/2,x+\eta/2]$ for any $x\in\mathcal{X}_{\eta}$. Therefore, a typical relationship between $x$ and $t$ is $"t\in I_{\eta}(x)"$. In the sequel, all the proposed multiple testing procedures are actually based on single tests of the null hypothesis $\mathcal{H}_{0}\left\{I_{\eta}(x)\right\}$, for all the possible window centers $x\in\mathcal{X}_{\eta}$, which is different than the null hypotheses $H_{0,t}$ given by $\theta(t)=0$ for $t$ in $[0,1]$. The continuous testing procedure aims at accepting the set of ”true windows”, which is indexed by their centers:

[TABLE]

$J_{0}^{\eta}(\theta)$ formally depends on $\theta$, but is denoted by $J_{0}^{\eta}$ in the sequel, for short.

2.4 A first simple construction of the $p$-value process

Homogeneity test.

Consider as an individual test statistic $S_{\eta}(x)=N(I_{\eta}(x))$,that is the number of points in the window of center $x$. Then classically, a one-sided test based on this statistic rejects $\mathcal{H}_{0}\left\{I_{\eta}(x)\right\}$ if $S_{\eta}(x)$ is large enough. Under $\mathcal{H}_{0}\left\{I_{\eta}(x)\right\}$, it is easy to see that the conditional distribution of $S_{\eta}(x)$ given $\mathfrak{N}=n$ (or $N([0,1])=n$ here) is a Binomial variable with parameter $n$ and $1/2$. Hence, if we denote by $F_{\mathcal{B}(n,\frac{1}{2})}(.)$, the function such that $F_{\mathcal{B}(n,\frac{1}{2})}(z)=P(Z\geq z)$ with $Z\sim\mathcal{B}(n,\frac{1}{2})$, then

[TABLE]

can be seen as a $p$-value associated to the classical one-sided single test of $\mathcal{H}_{0}\left\{I_{\eta}(x)\right\}$. Since $F_{\mathcal{B}(n,\frac{1}{2})}$ is explicit and since the number of points, as a function of $x$, can be very efficiently computed, one can have access to a $p$-value process $x\mapsto p_{\eta}(x)$, with a very low computational cost.

If $\lambda$ is known, one does not need to use a conditional distribution. The unconditional distribution of $S_{\eta}(x)$ is a Poisson distribution with mean $\eta\lambda$, under $\mathcal{H}_{0}\left\{I_{\eta}(x)\right\}$. If $F_{\mathcal{P}(\eta\lambda)}(z)=P(Z\geq z)$ with $Z\sim\mathcal{P}(\eta\lambda)$, then

[TABLE]

also defines a $p$-value process. Note that this set-up is quite close to the one of Chan and Zhang (2007).

Two-sample test.

Consider now as individual test statistic $S_{\eta}(x)=N_{A}(I_{\eta}(x)).$ A classical one-sided test based on this statistic also rejects $\mathcal{H}_{0}\left\{I_{\eta}(x)\right\}$, if $S_{\eta}(x)$ is large enough. Under $\mathcal{H}_{0}\left\{I_{\eta}(x)\right\}$, it is easy to see that the conditional distribution of $S_{\eta}(x)$ given $\mathfrak{N}=N$ (the joint process here), is the same as the conditional distribution of $S_{\eta}(x)$ given $N(I_{\eta}(x))=N_{A}(I_{\eta}(x))+N_{B}(I_{\eta}(x))$. It is a Binomial variable with parameter $N(I_{\eta}(x))$ and $1/2$. The corresponding $p$-value associated to the classical one-sided single test of $\mathcal{H}_{0}\left\{I_{\eta}(x)\right\}$ can then be defined by:

[TABLE]

This process depends on $N$ and $x$, and once again can be computed very easily.

2.5 Making continuous testing computationally tractable in general

Even if the method is based on a continuum of tests, the two previous examples are computationally tractable because they are based on $Y(I_{\eta}(x))$, the number of points within windows $I_{\eta}(x)$, which can be computed very efficiently and which is piece-wise constant as a function of $x$. In full generality, the number of points may not be informative enough, and to gain power, one can consider more general single test statistics (see Section 4) that only depend on the composition of each window, that is on the function

[TABLE]

which records the random positions of points in $Y$ that lie in $I_{\eta}(x)$. Since the composition is piece-wise constant with changes only each time a point leaves or enters the scanning window, the single test statistic that only depends on the composition is a piece-wise constant function on a random finite partition $\tau=(\tau_{m})_{0\leq m\leq M}$. Therefore, the continuum of tests can be reduced to a random finite number of tests, one for each segment defined by $\tau$. An example of such a partition if provided in Fig. 1.

2.6 General definition of the $p$-value process

In the sequel, we assume that a test statistic $S_{\eta}(x)$ is given for each window $I_{\eta}(x)$, that only depends on the composition of the window. Because Poisson processes have nice properties of independence between disjoint sets, it is easy to show the following property, which is usually referred to as subset pivotality (Westfall and Young, 1993; Romano and Wolf, 2005a).

[TABLE]

This roughly says the the conditional distribution of $S_{\eta}(x)$ given $\mathfrak{N}$ is not modified by the truth or the falsehood of the remaining hypotheses. In particular, the conditional distribution of $S_{\eta}(x)$ given $\mathfrak{N}$ under $\mathcal{H}_{0}\left\{I_{\eta}(x)\right\}$ is the same as the one under the full null $\mathcal{H}_{0}\left\{[0,1]\right\}$.

Note that in full generality the conditional distribution of $S_{\eta}(x)$ given $\mathfrak{N}$ depends on $\mathfrak{N}$ but also on the composition of the window centered on $x$. Therefore we define:

[TABLE]

and the associated $p$-value process by

[TABLE]

It satisfies for all $x\in\mathcal{X}_{\eta}$, for all $(\theta,\lambda)$,

[TABLE]

Moreover since $p_{\eta}(x)$ is a deterministic function of $S_{\eta}(x)$ (conditionally to $\mathfrak{N}$), property $(\mathcal{P})$ also holds for $p_{\eta}(x)$, and not only for $S_{\eta}(x)$.

3 False positive control in continuous time

Except when specifically mentioned, the results presented here hold for any test statistic $S_{\eta}(x)$ that only depends on the composition of the window (as defined in Section 2.5). Their associated $p$-value process are as described previously.

A rejection set $\mathcal{R}_{\eta}$ is a subset of $\mathcal{X}_{\eta}$, which is a function of the $p$-value process $(p_{\eta}(x))_{x\in\mathcal{X}_{\eta}}$. Classically, it corresponds to small values of this process and a typical example is given by

[TABLE]

for a given threshold $u$ potentially depending on the data. Once the rejection set $\mathcal{R}_{\eta}$ given, the set of accepted windows is defined by its complement $\mathcal{A}_{\eta}:=(\mathcal{R}_{\eta})^{c}.$

3.1 Definition of the multiplicity error rates

If the procedure behaves properly, the set $\mathcal{A}_{\eta}$ should be a good approximation of $J_{0}^{\eta}$, or similarly $\mathcal{R}_{\eta}$ should be a good approximation of $(J_{0}^{\eta})^{c}$. To evaluate the quality of a procedure $\mathcal{R}_{\eta}$, we focus on false positive windows that correspond to elements of $J_{0}^{\eta}\cap\mathcal{R}_{\eta}$. The size of $J_{0}^{\eta}\cap\mathcal{R}_{\eta}$ can be gauged in many ways. In this work, we focus on the most famous criteria, namely the Family-Wise Error Rate (FWER) and the False Discovery Rate (FDR).

From a family-wise point of view, one wants to avoid any false positive with a large probability, which corresponds to a control of the quantity

[TABLE]

While a FWER control provides a strong assessment on the amount of false positives, this approach generally leads to conservative procedures. A more balanced view is to allow some false positives, in a pre-specified fraction of the total amount of positives. This can be achieved by controlling

[TABLE]

with $0/0=0$. Importantly, since the setting is continuous here, the “number” of (false) positives is quantified above by the Lebesgue measure, $\Lambda$, on $[0,1]$.

3.2 A $\min$-$p\,$procedure to control the FWER in continuous time

We would like to find a threshold $u^{\alpha}\in[0,1]$ such that:

[TABLE]

For all $u\in[0,1]$,

[TABLE]

with $p_{\eta}(x)$ defined in (3). The set $J_{0}^{\eta}$ being unknown, one can use some rough upper-bound. We have, thanks to subset pivotality $\left(\mathcal{P}\right)$, that for any $u>0$ (possibly depending on the conditioning variable $\mathfrak{N}$),

[TABLE]

Since $J_{0}^{\eta}\subset\mathcal{X}_{\eta}$, we obtain

[TABLE]

The distribution of $\underset{x\in\mathcal{X}_{\eta}}{\inf}\left\{p_{\eta}(x)\right\}$ under $\theta_{0}$ and given $\mathfrak{N}$ are in both cases (homogeneity and two sample cases) known/observable once the variable $\mathfrak{N}$ is fixed and it is possible to calibrate an adequate threshold $u^{\alpha}$. However, due to the classical subtleties of inversion of non continuous c.d.f., it is in fact more powerful to use adjusted $p$-values that are defined by

[TABLE]

with $F^{\text{m}}_{\mathfrak{N}}$ being the conditional c.d.f. of $\underset{x\in\mathcal{X}_{\eta}}{\inf}\left\{p_{\eta}(x)\right\}$ under the full null hypothesis, that is $\theta=̀\theta_{0}$.

Theorem 3.1**.**

Let $\alpha\in(0,1)$ and $q_{\eta}$ given by (7). Then the procedure defined by

[TABLE]

has a FWER controlled by $\alpha$.

Note that $\mathcal{R}_{\eta}^{\text{m}}$, which is a continuous version of the classical $\min$-$p\,$procedure, is coherent with the case where there is only one value in the $p$-value process because the resolution $\eta$ is equal to 1. In this extreme case, $\mathcal{X}_{\eta}$ is reduced to one point and the $p$-value is exactly the adjusted $p$-value , the control of $\mbox{FWER}(\mathcal{R}_{\eta}^{\text{m}})$ being then an obvious consequence of Section (4).

Remark 3.2**.**

The previous procedure is usually referred to as a single-step procedure. We also propose refinements by using a step-down procedure in the Supplementary Material.

Remark 3.3**.**

Unfortunately, even for the classical homogeneity or two-sample tests the c.d.f. $F^{\text{m}}_{\mathfrak{N}}$ is not explicit and needs to be estimated. This is the purpose of Section 5 in which we also derive meaningful theoretical results when this distribution is approximated by Monte-Carlo that are necessary to provide a completely operational method with strong theoretical guarantee.

3.3 Link with the scan statistic framework.

In the case of the homogeneity test, the $p$-value process does not depend on $x$ and because function $F_{\mathcal{B}(N([0,1]),\frac{1}{2})}$ is a non increasing function, we have

[TABLE]

Hence, rejecting for small values of the infinimum of the $p$-value process (similar to standard minP procedures), or rejecting for large of the supremum of the test statistic (i.e. the scan statistic, often called the maxT procedure) is strictly equivalent (Dudoit et al., 2003; Daley and Vere-Jones, 2008). Compared with the work of Chan and Zhang (2007), we mainly added two ingredients: ($i$) the distribution we are considering is conditional and exact, whereas theirs is asymptotic and still depends on the unknown nuisance parameter $\lambda$ (even if they provide a plug-in procedure in practice to avoid this issue), ($ii$) thanks to Property $(\mathcal{P})$, we are able to prove the FWER control, whereas their control stands under the full null only.

This minP-maxT equivalence is valid in the case of the homogeneity test only, but it does not hold as soon as the c.d.f. of the infinimum of the $p$-value under the full null depends on the composition of the window. This is typically the case for the two-sample test, where $F_{\mathcal{B}(N(I_{\eta}(x)),\frac{1}{2})}$ depends on $x$ so that Eq. (9) cannot hold anymore. In this case, the minP framework offers the advantage to scale all tests on the same level (Dudoit et al., 2003; Daley and Vere-Jones, 2008).

3.4 A weighted BH procedure to control the FDR in continuous time

To increase the detection capability, we propose a second procedure based on the control of the FDR criterion (6).

Heuristic.

We would like to find a threshold $v^{\alpha}\in[0,1]$ (see (3)) such that:

[TABLE]

Following Blanchard et al. (2014), let us first explain heuristically how one can reach this goal in a continuous framework. By Fubini’s theorem, $\forall v\in[0,1]$,

[TABLE]

Then, choosing $v$ such that

[TABLE]

and doing as if $v$ was deterministic leads to

[TABLE]

which leads to

[TABLE]

by using the property of $p$-values (4). The above heuristic suggests to choose $V^{\alpha}$ by using

[TABLE]

Weighted-BH procedure.

Quantile $V^{\alpha}$ is random (even conditionally on $\mathfrak{N}$) because it depends on the observed $p$-value process. Nevertheless, since the $p$-value process is piece-wise constant in our framework, finding threshold $V^{\alpha}$ is very simple. We propose a step-up algorithm that is inspired from the famous weighted BH threshold (Benjamini and Hochberg, 1997), with the essential difference that our algorithm is based on random weights. Recall that the processes $(S_{\eta}(x))_{x}$ and therefore $(p_{\eta}(x))_{x}$ are piece-wise constant on the random partition $\tau=(\tau_{m})_{0\leq m\leq M}$ ($\tau_{0}=0$, $\tau_{M}=1$). Therefore

[TABLE]

This implies that the threshold $V^{\alpha}$ defined in (3.4) can be derived as follows:

• •

Compute the weights,

[TABLE]

• •

Denote $p_{m}=p_{\eta}(\tau_{m-1})$, $1\leq m\leq M$,

and order these $p$-values in increasing order $p_{\sigma(1)}\leq\dots\leq p_{\sigma(M)}$ for an appropriate permutation $\sigma$ of $\{1,\dots,M\}$;

• •

Consider $\widehat{k}=\max\{k\in\{1,\dots,M\}\>:\>p_{\sigma(k)}\leq\alpha\sum_{\ell=1}^{k}w_{\sigma(\ell)}\}$ (with $\widehat{k}=0$ if the set is empty);

• •

Compute $V^{\alpha}$ as $\alpha\sum_{\ell=1}^{\widehat{k}}w_{\sigma(\ell)}.$

• •

Define $\mathcal{R}_{\eta}^{wBH}=\{x\in\mathcal{X}_{\eta}:p_{\eta}(x)\leq V^{\alpha}\}.$

By contrast with the classical BH step-up procedure, each $p$-value has its own weight $w_{m}$ which is not uniform but adaptive to the data. This weight means that the quantity of interest is not the number of tests, i.e. the cardinality of the partition $\tau$, but the proportion of ”time” that the $p$-value process spends at a particular value.

Remark 3.4**.**

It is worth to note that $\mathcal{R}_{\eta}^{wBH}=\{x\in\mathcal{X}_{\eta}:q_{\eta}(x)\leq\alpha\}$ where ${q}_{\eta}(x)$ stands for the adjusted $p$-value

[TABLE]

Theoretical result.

Our procedure is based on a heuristic that has been made fully theoretically valid by Blanchard et al. (2014) under suitable assumptions of measurability and positive dependence (namely, finite dimensional positively regressively dependent $p$-value process). The difficulty comes from the interplay between the numerator and the denominator within the expectation in (3.4) via $v$, because $v=V^{\alpha}$ is a random variable. While the measurability issue is not critical here (by the piece-wise constant property), the positive dependence condition seems difficult to check in general. For instance, this condition is known to exclude the two-sided test statistic (Yekutieli, 2008). We prove here that this positive dependence condition holds in the one-sided testing setting (both for the homogeneity test with a known $\lambda$ and the two-sample test), by using the peculiar properties of Poisson processes.

Theorem 3.5**.**

For all $\alpha\in(0,1)$, for a $p$-value process given either by Eq. 1 or Eq. 2, for one-sided null hypotheses, the corresponding procedure $\mathcal{R}_{\eta}^{wBH}$, satisfies

[TABLE]

This result is proved in the Supplementary File. Theorem 3.5 hence theoretically justifies the use of a weighted step-up procedure to control the continuous FDR. However, it is only valid for simple versions of our test statistic, which ensure the positive dependence structure. For more general cases, Blanchard et al. (2014) proposed a modification of the step-up algorithm. This modification consists in introducing a function $\beta(v)\leq v$ of a specific form, such that the quantity $\sum_{\ell=1}^{k}w_{\sigma(\ell)}$ is replaced by the smaller quantity $\beta\left(\sum_{\ell=1}^{k}w_{\sigma(\ell)}\right)$. However, this modification is reported to be quite conservative (Blanchard and Roquain, 2008). Furthermore, our simulation results show that the FDR control (3.4) is achieved in most of the situations we studied, even for more complex $p$-value processes that the ones we considered here (homogeneity test with known intensity, and two-sample test)

Note that the homogeneity test with known intensity in Theorem 3.5 is similar to Corollary 4.5 in Blanchard et al. (2014). However, a subtle distinction concerns the way the null hypotheses are defined (windows versus points). A consequence is that our approach does not rely on any knowledge regarding the regularity of the underlying process intensity.

4 Kernel-based statistics for sliding windows

4.1 On the limitations of the count-based statistic

Let us focus on the two-sample test to illustrate our motivations. In this context, our objective is to detect local differences between unknown functions $\nu_{A}$ and $\nu_{B}$ based on individual tests defined on a set of windows. In a first approximation, this can be achieved by comparing $N_{A}(I)$ with $N_{B}(I)$, the count of occurrences on window $I$, that are distributed as a couple of independent Poisson variables. However, given that $N_{A}(I)$ (resp. $N_{B}(I)$) is an approximation of $\int_{I}\nu_{A}(t)dt$ (resp. $\int_{I}\nu_{B}(t)dt$), the count based-statistic strategy can be considered as an hypothesis that tests the equality of function integrals rather than the equality of the functions themselves. Hence the induced testing procedure may lack of accuracy.

Therefore we developed a kernel-based strategy to propose single test statistics that do not depend on the sole number of occurrences in a window, but that consider its whole composition as defined in Section 2.5. The interest of kernel-based statistics is that they account the physical distances between occurrences, which explains the gain of performance.

4.2 The kernel vision in the two-sample case

The kernel framework for testing proposed by Gretton et al. (2012) and Fromont et al. (2013) has been especially designed to separate functions, and not just their integrals. It is based on the estimation of a distance between the functions $\nu_{A}$ and $\nu_{B}$ (in a RKHS, or in $L^{2}$). In particular, Fromont et al. (2013) show that these test statistics are minimax with respect to various smoothness classes, and in this sense, powerful against some alternatives that cannot be distinguished by count-based tests.

In the following we consider translation invariant kernels denoted by $K$, and in practice we will restrict ourselves (see Sections 6 and 7) to the Gaussian kernel that has been shown to be the most powerful on simulations (Fromont et al., 2013). The shape of the statistic is then given by

[TABLE]

where $N$ is the joint process and the $\varepsilon_{T}$’s are the marks (see Section 2.1). In particular, it only depends on the composition of the window. To understand, why this statistic is meaningful, let us compute its expectation. By classical computations for Poisson processes,

[TABLE]

If $K\star f$ denote the convolution product, we end up with

[TABLE]

In particular if $K$ is a convolution kernel with bandwidth $h$ we obtain, under classical assumptions, that

[TABLE]

This short derivation illustrates the motivation to use such a kernel-based statistic, as a good estimate of a distance between intensity functions. The one-sided version of this statistic, as well as the homogeneity case are presented in the Supplementary Material. Note that from a computational point of view, the complexity for computing $S_{\eta}(x)$ is in $O(N^{2}(I_{\eta}(x))$. However, thanks to the piece-wise constant property of $S_{\eta}(x)$ on partition $\tau$, the continuum $\{S_{\eta}(x),x\in\mathcal{X}_{\eta}\}$, can be computed using the in/out points of each scanning window, which drastically reduces the global computational burden.

5 Monte-carlo procedures for $p$-values estimation and adjustment

5.1 The distribution of the single test statistic being explicit.

For tests based on the count statistic the $p$-value process is explicitly computable. Hence, no need for Monte-Carlo in this case. Moreover, since our weighted BH procedure (Section 3.2) is completely explicit, $\mathcal{R}^{wBH}_{\eta}$ can be easily computed in practice. In this situation we have already proved the FDR control for some particular cases and we will show that it holds in more general cases by using simulations (Section 6).

However, the $\min$-$p\,$procedure is not explicit since it requires the conditional distribution of $\inf_{x\in\mathcal{X}_{\eta}}p_{\eta}(x)$ given $\mathfrak{N}$ under the full null hypothesis (Section 3.2). Therefore we need to use a Monte-Carlo procedure to approach this distribution. This can be achieved while maintaining a valid FWER control by following Romano and Wolf (2005a), as detailed in Section 5.4 below.

5.2 The distribution of the single test statistic being not explicit

This concerns kernel-based tests (Section 4) for which the $p$-value process is not even computable. One can use again the estimated $p$-values developed by Romano and Wolf (2005b) since their result shows that they satisfy (4) and $\mathcal{R}^{wBH}_{\eta}$ can be easily computed in practice. However, there is no theoretical FDR guarantee and we will only show the FDR control by simulations in Section 6. Consequently, to compute $\widehat{p}_{\eta}(x)$, the weighted-BH procedure only requires one Monte-Carlo step in this situation.

As for the $\min$-$p\,$procedure (Section 3.2), it requires an additional Monte-Carlo step. To this end, a classical and sufficient method is to perform an approximation using two independent Monte-Carlo samples. We detail this approach for the Homogeneity case in the supplementary material. For the two-sample case, the Monte-Carlo procedure is quite demanding and is close in spirit to bootstrap methods. In this case, we developed a double Monte-Carlo procedure, which only requires one Monte-Carlo sample and still guarantees a controlled FWER (Section 5.4).

5.3 Estimation of the $p$-value process by Monte-Carlo

Recall that this estimation is required when using test statistics different from the statistics based on counts. The full procedures are summarized in Table 1.

Homogeneity test.

Here, the conditioning variable $\mathfrak{N}=N([0,1])$ is the same for all scanning windows and the distribution of $S_{\eta}(x)$ under $\mathcal{H}_{0}\{I_{\eta}(x)\}$ does not depend on $x$. Let us focus on the first window $I_{\eta}(\eta/2)=(0,\eta]$ and simulate points inside this window under $\mathcal{H}_{0}\{I_{\eta}(\eta/2)\}$. Hence we consider $B$ random samples of the number of points $N^{b}(I_{\eta}(\eta/2))$ such that

[TABLE]

and then we draw

[TABLE]

to sample positions of a homogeneous Poisson process $N^{b}$ inside window $(0,\eta]$. For each $b$, one can compute thanks to the positions of the points inside $[0,\eta]$ the value of the test statistic for the first window, that we denote $S_{\eta}^{b}$. But since all the statistics are distribution invariant by translation under the null, one can also say that this sample has the same distribution as $S_{\eta}(x)$ under $\mathcal{H}_{0}\{I_{\eta}(x)\}$ for every window center $x\in\mathcal{X}_{\eta}$.

Therefore, whatever $B$, one has easily access to \big{(}S^{1}_{\eta}\ldots,S^{B}_{\eta}\big{)}, $B$ i.i.d. variables with the same distribution as $S_{\eta}(x)$ under $\mathcal{H}_{0}\{I_{\eta}(x)\}$ for every $x$. Then, following Romano and Wolf (2005b), we define:

[TABLE]

which amounts to use the empirical function $F_{\mathfrak{N},x}$ defined in Eq. 3 over a $B+1$ sample where $S^{0}_{\eta}(x)$ is the observed statistic.

Two sample test.

We label the observed marks $\varepsilon^{0}:=\varepsilon=(\varepsilon_{T})_{T\in N},$ such that it constitutes the first term of a $(B+1)$-sample of marks, filled by $B$ independent draws of i.i.d. Rademacher sets, $\varepsilon^{b}:=(\varepsilon^{b}_{T})_{T\in N}$, for $b=1,...,B$. The conditional distribution given $N=\mathfrak{N}$ of these Rademacher processes are:

[TABLE]

that is the distribution that the observed marks $\varepsilon_{T}$’s should have under the full null hypothesis. As previously, the Monte-Carlo approximated conditional $p$-value of window $I_{\eta}(x)$ is defined by:

[TABLE]

where $S^{0}_{\eta}(x)$ is the observed statistic and where $S_{\eta}^{b}(x)$ is the value of the test statistic on $I_{\eta}(x)$ that is computed with the resampled marks $\varepsilon^{b}$ and fixed joint process $N$. Let us underline that because in the two-sample case the partition $\tau$ only depend on $N$ and do not vary with the resampling scheme $\varepsilon^{b}$, the processes $(S^{b}_{\eta}(x))_{x\in\mathcal{X}_{\eta}}$, for $b=0,...,B$, and therefore the process $(\widehat{p}_{\eta}(x))_{x\in\mathcal{X}_{\eta}}$ are piece-wise constant on $\tau$.

5.4 Two-step Monte-Carlo method for the $\min$-$p\,$procedure.

Recall that when using the $\min$-$p\,$procedure with a $p$-value process that is not explicit, we need two Monte-Carlo approximations.

Homogeneity test.

In that case, the distribution of the test statistics under the null are independent of the window centers, while the distribution of the minimum of the $p$-values relies on the observations over all windows. Hence, it is quite natural to perform two separated Monte-Carlo schemes for these two steps.

Two sample test.

In the two-sample case, the trick is that one can use the same sample of $(N,\varepsilon^{b})$’s as the one used in (10) for both approximations. More precisely, the $p$-values $\widehat{p}_{\eta}^{b}(x)$, corresponding to the simulated process $(N,\varepsilon^{b})$, can be computed as

[TABLE]

that is we do the same computation as in Eq. 10 except that we replace the observed statistics $S^{0}_{\eta}(x)$ by the simulated one $S^{b}_{\eta}(x)$, hereby computing a ”resampled” p-value process corresponding to $\varepsilon^{b}$. Then, let us define, for all $b=0,...,B$

[TABLE]

The corresponding adjusted Monte-Carlo $p$-values are then given by

[TABLE]

Finally the rejection set is given by

[TABLE]

The previous double Monte-Carlo procedure is based on the use of the same random signs $\left(\varepsilon^{b}_{T}\right)_{T\in N}$, $0\leq b\leq B$ and not on two sets of different simulations. Indeed the $\varepsilon^{b}$’s are used twice: a first time to estimate the $p$-values of the observed process and a second time to give the $p$-values even on the simulated realizations leading to the adjusted $p$-values . This procedure has the double advantage of sparing computational resource, while enjoying the same control property as shown hereafter.

Theorem 5.1**.**

Let $\alpha\in(0,1)$. For the two-sample case, the continuous testing procedure defined by (5.4) satisfies

[TABLE]

This result holds whatever the choice of the single test statistic that only depends on the composition of the window.

6 Simulations

6.1 Homogeneity test

For a fixed parameter $r$, we simulate a homogeneous Poisson process corresponding to the following piecewise constant signal $\theta$ (illustrated in Figure 9, Supp. Material):

[TABLE]

with

[TABLE]

and $\mathcal{I}_{1}=\mathcal{I}_{1}^{+}\cup\mathcal{I}_{1}^{-}$, $\mathcal{I}_{0}=[0,1]\backslash\mathcal{I}_{1}$, the positions on which $\theta(.)$ is not null (or null, respectively). We consider this particular shape of $\theta$ to ensure that $\int_{0}^{1}\theta(t)dt=0$.

For a fixed background intensity parameter $\nu^{*}$, we therefore simulate 1000 times a Poisson process $N$ of intensity $\nu(.)=\nu^{*}(1+\theta(.))$. Parameter $\theta^{*}$ is used to tune the distance to homogeneity (which reflects the difficulty of the test), with $\theta^{*}=0.01$ standing for a situation where almost every position is under the homogeneous regime, $\theta^{*}=0.99$ standing for a situation with high separability. Parameter $\nu^{*}$ controls the number of occurrence and varies in $\{500,1000\}$. To compare with other scanning procedures, we opt for one-sided tests (the alternative corresponding to too large $\theta$), with fixed window size $\eta=2r$. More precisely, we compare our count-based $\min$-$p\,$and count-based weighted BH procedures in continuous time with asymptotic approaches based on the scan statistic that controls the FWER (Chan and Zhang, 2007) and the FDR (Siegmund et al., 2011).

We consider this homogeneity test to compare asymptotic approaches based on the scan statistic to control the FWER (Chan and Zhang, 2007) and the FDR (Siegmund et al., 2011), with our $\min$-$p\,$and weighted BH procedures in continuous time, based on the count statistic. Figure 2 clearly shows that our $\min$-$p\,$procedure in continuous time controls the FWER even for moderate sample sizes, whereas the scan statistic procedure requires a strong signal ($\theta^{*}\sim 1)$, or a higher number of occurrences ($\nu^{*}=5000$, not shown). This trend is even more clear for the FDR control, that is not ensured by method of Siegmund et al. (2011). Note that contrary to continuous testing, this method requires a discretization step on which the performance depend. In this case we choose fixed-length non-overlapping windows to compute their procedure.

6.2 Two sample test

To test the differences between intensities $\nu_{A}$ and $\nu_{B}$, we sample a joint point process $N$ over $[0,1]$ with constant intensity $\nu(.)=\nu^{*}\in\{500,1000\}$. Then we sample the labels $\varepsilon^{0}$ given the occurrences of the joint process $T\in N$, in an independent manner, such that

[TABLE]

Thus positions $T\in N$ with $\varepsilon_{T}^{0}=+1$ (resp. $-1$) stand for points of process $N_{A}$ (resp. $N_{B}$) with intensity $\nu_{A}$ (resp. $\nu_{B}$). Parameter $\theta^{*}$ represents the distance between $\nu_{A}$ and $\nu_{B}$ (which reflects the difficulty of the test), with $(\theta^{*}+1)/2=0.51$ standing for a situation where the two intensities are nearly equal, and $(\theta^{*}+1)/2=0.9$ standing for a situation with high separability. $\mathcal{I}_{0}$ and $\mathcal{I}_{1}$ are defined as in the homogeneity case, and each configuration is repeated 1000 times. For each value of $\theta^{*}$, bootstrap $p$-values are computed using $B=10^{5}$ Monte Carlo samples, along with the Gaussian kernel $K$ of bandwidth $h=\eta$. Decision rules are based either on the $\min$-$p\,$procedure, or on the weighted BH procedure, and we are comparing the performance of the count statistic (one-sided) with the statistic based on the Gaussian Kernel (one-sided). The size of the windows is fixed $\eta=2r$.

Even if the FWER and the FDR are well controlled by our two procedures, the Gaussian kernel statistic seems to be less conservative when controlling the FWER and seems more robust changes in the intensity $\nu^{*}$ of the process (Fig. 3). The control of the FDR are comparable between the two statistics. However, when the kernel and the count statistics are compared using sensitivity and specificity (Fig. 4), for the same level of signal (same $\theta^{*}$), the kernel statistics is globally more accurate than the count statistics (higher sensitivity and specificity). This can be explained by the fact that kernel-based statistics are in fact $U$-statistics, hereby using information from pairs of points through the kernel. The difference of these statistics will also be illustrated on the analysis of experimental data.

7 Applications

Description of the full procedure in practice.

In the following applications, we compute the $p$-value process as well as its $\min$-$p\,$and wBH adjusted versions. In order to propose a completely operational procedure, we need to clarify the distinction between the points space and the windows space on which the decision is made. Indeed, if window $I_{\eta}(x)$ is accepted by our procedure, and if there is no mistake, then one is sure that $H_{0,t}$ holds for all point t in $I_{\eta}(x)$, whereas if a window $I_{\eta}(x)$ is rejected without mistake, it is yet not clear whether there is not a smaller interval included in $I_{\eta}(x)$ on which the null hypothesis holds. Consequently we start by defining the set of accepted windows

[TABLE]

with $\widehat{q}_{\eta}(x)$ denoting the estimated adjusted $p$-value process. Finally, the graphical representation of our procedure is based on the estimated set of points that lie within rejected windows, namely: $\widehat{\mathcal{I}}_{1}^{\eta}=[0,1]\backslash\widehat{\mathcal{I}}_{0}^{\eta}$, where

[TABLE]

7.1 Do neurons respond to some stimulus ?

In Neurosciences, the succession of action potentials of a given neuron, also called spike trains, constitutes the key biological signal that carries information regarding the functioning of neurons. In particular, understanding the rapid changes in firing rates due to a certain stimulus (Kass et al., 2003) is a challenging task, and inhomogeneous Poisson processes are classically used to model such data (Pipa et al., 2013). Asking whether there is a change of firing rate after exposure typically refers to our homogeneity test, and we are seeking for time intervals on which the rate is modified with respect to the standard behaviour. To illustrate our procedure we analyzed public experimental data (Pouzat and Chaffiol, 2009) that consist in spike trains of a cockroach neuron stimulated by different odors like citronellal or terpineol. A first question is to determine if this neuron responds to the stimulus (one-sided homogeneity test), and if the response is different from one stimulus to another (two-sided two-sample test).

The neuronal activity starts to be recorded at time 0, and the citronellal is presented between 0.42 and 0.45 (normalized time, Fig. 5). Using the count statistic and the $\min$-$p\,$procedure, our method identifies a neuronal response only when the neuron is stimulated (Fig. 5, dark-grey rectangle), but not after the stimulus, whereas the weighted BH procedure identifies two firing rate changes posterior to the stimulus. Interestingly, the weighted BH procedure also identifies a deviation before the stimulus, which may correspond to habituation and therefore anticipation of the stimulus by the animal (see for a similar phenomenon with synchrony instead of firing rates modifications (Tuleau-Malot et al., 2014)). What is remarkable is that the Gauss Kernel statistic globally identifies the same intervals of changes when compared with the count statistic (Fig. 6), with the notable difference that the instantaneous response to the stimulus is located more precisely, which illustrates the higher specificity of the kernel-based statistics. Moreover, the kernel detects an additional interval of change in the firing rate. When studying the differential response of the same neuron to different stimuli (citronellal, terpineol) the differential change of fire rate occurs just after the stimulus (Fig. 7). This difference is identified by the weighted BH procedure, but not by the $\min$-$p\,$, which lacks of power, as already documented.

7.2 Differential analysis of replication origins between cell lines

Next Generation Sequencing technologies (NGS) are based on the massive parallel sequencing of short DNA fragments, which have given rise to a new area in the field of high throughput biology. The applications of such technologies are many, and we focus on the particular cases where NGS allow the fine mapping of genomic features along the genome, like transcription factors (using the ChIP-Seq technology), or copy number variations (CNV-Seq) for instance. When applied to the field of replication, NGS technologies have allowed the fine-scale mapping of replication origins along the human genome. Briefly, replication is the mechanism by which genomes are duplicated into two exact copies. Genomic stability is under the control of a spatio-temporal program that orchestrates both the positioning and the timing of firing 100,000 replication starting points, also called replication origins (Picard et al., 2014). This fine mapping has allowed us to identify potential regulators of the replication dynamics, and the comparative analysis of replication origins maps can help to identify a core of origins that are common to all human cells, compared with cell-type specific origins more related to cell differentiation. Moreover, pathologies like cancer being linked to replication instabilities, comparing replication origins in cancer cells vs. non-cancer cells could be promising.

The chromosomal organization of replication origins can easily be modeled by a heterogeneous Poisson process of unknown intensity, and we consider the comparison of origins between cancer cells (Hela) and embryonic stem cells (iPS). Since the total number of origins differ from one condition to another and that globally there is no real difference in shape between both intensities except at some particular positions of interests, we need to account for the difference in the total number of occurrences per condition. To proceed, we modify the conditional distribution of the marks under the null, which becomes conditional to the total (random) number of points coming from condition A and B respectively:

[TABLE]

with $\lambda_{A}=\int_{0}^{1}\nu_{A}(t)dt$ and $\lambda_{B}=\int_{0}^{1}\nu_{B}(t)dt$. Since in practice $\lambda_{A}$ and $\lambda_{B}$ are unknown, they are replaced by the total number of points in $N_{A}$ and $N_{B}$. This procedure resumes to testing the equality of density functions instead of intensities.

As expected, origins densities show strong variations along chromosomes, but those variations are globally comparable between cell lines (Fig. 8). One region shows significantly more origins in Hela cells (cancer cells), that is detected by both the count and the Gaussian kernel statistic after the continuous wBH procedure. The $\min$-$p\,$adjustment does not allow the detection of this region at the same nominal level, but the variations of the adjusted $p$-value process indicate a strong difference of densities between cell lines. This region corresponds to a gene-rich region that needs investigation for further biological characterization. More and more cell lines have been investigated to determine the core of replication origins that are common to any human cells, and the set of origins that are cell-type specific (Besnard et al., 2012). This context offers a nice perspective of the present work, by generalized our two-sample test to the more general framework of a $k$-sample (or ANOVA-like) test.

8 Conclusion

We proposed a new framework for testing Poisson processes intensities based on continuous testing. This framework appears very powerful, from both the theoretical and the practical points of views. It allows the complete definition of the testing procedure, based on local null hypothesis, which provides a clear definition of false positive windows, that was lacking in existing works. We introduced the $p$-value process that becomes a stochastic process that needs to be adjusted to account for joint error rates. Contrary to existing approaches, continuous testing avoid the arbitrary discretization of this process, which is coherent with the continuous nature of the underlying model. The continuous framework especially provides a rigorous definition of error rates in continuous time, and we proposed two procedures to control the FWER and the FDR. In addition to theoretical controls, we show the empirical properties of our procedure, that can be easily implemented, and is available in the form of the contest R-package. The application on experimental data inspires extensions of our framework, that could be adapted to the testing of more than 2 conditions, or to other tests like local independence tests that are essential in Neurosciences (Albert et al., 2015). Calibrating the window size is certainly one of the major perspective of this work. While the subject has been extensively investigated in the context of model estimation, only a few works considered the question in the testing framework (Gao and Gijbels, 2008). Consequently, one challenge ahead will be to determine if the continuous testing framework would be appropriate to derive a calibration method to automatically select the resolution of the test from the data. Other future extensions involve either spatial aspects or more general point processes from a theoretical point of view.

Acknowledgment

We are grateful to Anne-Laure Fougères for fruitful discussions on the subject, to Matthieu Lerasle for sharing with us the trick of the double Monte-Carlo step, and to Ivan Bardet for his work on the Chan and Zhang (2007) paper as an intern. This research was partly supported by the french Agence Nationale de la Recherche (ANR 2011-BS01-010-01 Calibration project, ANR-11-BINF-0001-06 ABS4NGS project, ANR-15-CE12-0004-02 OriMolMech project and ANR-16-CE40-0019 SansSouci project). This work was performed using the computing facilities of the computing center LBBE/PRABI.

Appendix A Supplementary material

A.1 Step-down improvement

Let us first introduce some notation. For any fixed subset $\mathcal{C}$ of $\mathcal{X}_{\eta}$, $F_{\theta_{0},\mathfrak{N}}^{\mathcal{C}}$ is the conditional c.d.f. given $\mathfrak{N}$ of $\inf_{x\in\mathcal{C}}p_{\eta}(x)$. Then we define

[TABLE]

The theoretical step-down algorithm can then be derived as follows:

• •

Step 0 : compute $\mathcal{R}^{0}=\emptyset$ given by (8).

• •

Step $j\geq 1$: compute $\mathcal{R}^{j}=\mathcal{N}\left((\mathcal{R}^{j-1})^{c}\right)$

The step-down rejection set $\mathcal{R}^{step-down}_{\eta}$ is then defined as the limit of $\mathcal{R}^{j}$ when $j$ tends to infinity. Note that $\mathcal{R}^{1}=\mathcal{R}_{\eta}^{\text{m}}$ is in fact given by (8) and that the sequence of $\mathcal{R}^{j}$ is increasing.

The fact that this algorithm ends almost surely in a finite number of steps is just an easy consequence of the fact that for a given realization, there is only a finite number of possible values for the $p$-value process, which is in fact upper bounded by the size of the partition $\tau$. Indeed, at each round, $\mathcal{R}^{j}$ cannot decrease and if it increases, it is by absorbing sets of the type $\{x\in\mathcal{X}_{\eta},p_{\eta}(x)=a\}$ for some $a$.

One can prove that such an algorithm guarantees a controlled FWER following the ideas of (Goeman and Solari, 2010). More precisely if one defines the event

[TABLE]

one has by $\left(\mathcal{P}\right)$ (see Section 2.6) that

[TABLE]

The last point comes from the fact that $F_{\theta_{0},\mathfrak{N}}^{J_{0}^{\eta}}(\inf_{x\in J_{0}^{\eta}}p_{\eta}(x))$ can be seen as the $p$-value of the test based on the statistics $-\inf_{x\in J_{0}^{\eta}}p_{\eta}(x)$ and described for instance in Lemma 1 of (Fromont et al., 2016), which implies that (4) holds even for the $p$-value $F_{\theta_{0},\mathfrak{N}}^{J_{0}^{\eta}}(\inf_{x\in J_{0}^{\eta}}p_{\eta}(x))$. This is the single-step property of (Goeman and Solari, 2010).

For the monotony property of (Goeman and Solari, 2010), remark that if $\mathcal{C}\subset\mathcal{C^{\prime}}$ then $F_{\theta_{0},\mathfrak{N}}^{\mathcal{C}}\leq F_{\theta_{0},\mathfrak{N}}^{\mathcal{C^{\prime}}}$ and $\mathcal{N}(\mathcal{C}^{\prime})\subset\mathcal{N}(\mathcal{C})$. Hence for all $j$, if $J_{0}^{\eta}\subset(\mathcal{R}^{j})^{c}$, $R^{j+1}=\mathcal{N}([\mathcal{R}^{j}]^{c})\subset\mathcal{N}(J_{0}^{\eta})$.

Now, on $\mathcal{E}$, $\mathcal{N}(J_{0}^{\eta})\cap J_{0}^{\eta}=\emptyset$ and therefore if $J_{0}^{\eta}\subset(\mathcal{R}^{j})^{c}$, $J_{0}^{\eta}\subset(\mathcal{R}^{j+1})^{c}$ on the same event. Since $J_{0}^{\eta}\subset(\mathcal{R}^{0})^{c}=\mathcal{X}_{\eta}$, one has with probability larger than $1-\alpha$ that $J_{0}^{\eta}\subset\lim_{j}\mathcal{R}^{j}=\mathcal{R}_{\eta}^{step-down},$ which means that $\mbox{FWER}_{\theta,\lambda}^{\eta}(\mathcal{R}_{\eta}^{step-down})\leq\alpha.$

A.2 Other kernel-based single test statistics

A.2.1 Homogeneity test statistics

Two sided test statistics.

The starting point is that, conditionally on the event $\{N([0,1])=n\}$, the points in $N$ form an $n$ i.i.d. sample with density $f=1+\theta$. Hence, for any translation invariant kernel $K_{h}:[0,1]\to\mathbb{R}$ of bandwith $h$, we have

[TABLE]

is an unbiased consistent estimate (when $n$ tends to infinity) of

[TABLE]

Consequently, the corresponding $U$-statistics given by

[TABLE]

is an unbiased estimate of:

[TABLE]

under standard conditions on the kernel $K_{h}$. Actually, in a testing framework, we do not want to estimate $\|f\|_{I_{\eta}(x)}^{2}$ stricto sensu: here, we do not need to consider small values for $h$, we even recommend to take $h=\eta$ when it is unknown. This choice would not be appropriate for estimation, but seems sufficient for testing. Moreover if $\mathcal{H}_{0}\left\{I_{\eta}(x)\right\}$ holds, $f=1$ on $I_{\eta}(x)$ and $\|f\|^{2}_{I_{\eta}(x)}=\eta$. So we propose the following statistics

[TABLE]

and the local null hypothesis is rejected for high values of $S_{\eta}(x)$.

One-sided test statistic.

In this case we focus on what is happening when $f>1$. Hence the target quantity to estimate is

[TABLE]

Here we consider the following density estimator:

[TABLE]

which leads to the following statistics:

[TABLE]

If there exists a $t$ in $I_{\eta}(x)$ such that $f(t)>1$ and if $f$ is smooth enough such that this holds on a neighborhood, $E[\max(f(T),1)\mathbf{1}_{\left\{T\in I_{\eta}(x)\right\}}]$ should be larger than $\eta$, whereas its value should be smaller than $\eta$ if there are no such $t$.

A.2.2 Two-sample test statistics (one-sided)

Applying a reasoning similar to above, the quantity

[TABLE]

is a good estimate of $\theta(T)\lambda(T)$, for each $T\in N$. If one wants to reject only when there exists a $t\in I_{\eta}(x)$ such that $\theta(t)>0$, then one needs to take the positive part of the previous quantity and integrate it with respect to $T\in N\cap I_{\eta}(x)$. We do not want to multiply again by $\varepsilon(T)$ because we do not want to multiply by $\theta$ twice. This leads to

[TABLE]

which estimates E_{\lambda}\big{[}\max(\theta(T),0)\mathbf{1}_{\left\{T\in I_{\eta}(x)\right\}}\big{]} for $T$ of density proportional to $\lambda$.

Appendix B Proofs

B.1 Proof of Theorem 3.1

Computations and arguments similar to Section A.1 gives that

[TABLE]

which concludes the proof.

B.2 Proof of Theorem 3.5

We follow the methodology introduced in Blanchard et al. (2014) by applying their Theorem 4.1 in the ”Poisson case”, as discussed in Example 2.5 and Section 4.2.4 therein. While the measurability conditions can be easily checked because the $p$-value process is càdlàg, the result will be proved if we show that the $p$-value process is finite-dimensional strong PRDS on $J_{0}^{\eta}$, as defined in Blanchard et al. (2014). For the homogeneity test, Lemma A.2 in Blanchard et al. (2014) shows that the $p$-value process is finite-dimensional strong PRDS on any subset, and thus also on $J_{0}^{\eta}$. For the two-sample test, this condition takes the following form: for any $q\geq 1$, for any $(x_{j})_{1\leq j\leq q}\in(\mathcal{X}_{\eta})^{q}$ such that $x_{1}\in J_{0}^{\eta}$, for any nonincreasing set $D\subset\mathbb{N}^{q}$, for any $n\geq 0$,

[TABLE]

Above, a nonincreasing set is defined as a subset $D\subset\mathbb{N}^{q}$ such that for all $(m,m^{\prime})\in\mathbb{N}^{q}$ with $m_{j}\leq m^{\prime}_{j}$ for $1\leq j\leq q$, we have $m^{\prime}\in D\Rightarrow m\in D$.

Now, to establish (12), we work conditionally on the joint Poisson process $N$. In that case, the random part of $S_{\eta}(x)=\sum_{T\in N\cap I_{\eta}(x)}(2\varepsilon_{T}-1)$ is only carried by the set of marks $(\varepsilon_{T})_{T\in N}$. Let $b_{T}=2\varepsilon_{T}-1$ for $T\in N$ and

[TABLE]

Since $D$ is a nonincreasing set, so is $D^{\prime}$ and thus there exists some subset $\mathcal{T}\subset N$ such that $D^{\prime}=\{(b_{T})_{T\in N}\in\{0,1\}^{N}\,:\,\forall T\in\mathcal{T},b_{T}=0\}$. Furthermore, since $(b_{T})_{T\in N\cap I_{\eta}(x_{1})}$ is independent of $(b_{T})_{T\in N\backslash I_{\eta}(x_{1})}$, we have

[TABLE]

by denoting $U=\#N\cap I_{\eta}(x_{1})$, $V=\#\mathcal{T}\cap I_{\eta}(x_{1})$ and by an elementary combinatorial argument using that $b_{T}$, $T\in N\cap I_{\eta}(x_{1})$, are i.i.d. $\mathcal{B}(1/2)$ variables. Finally, since $\binom{U-V}{n}/\binom{U}{n}=\binom{U-n}{V}/\binom{U}{V}$, the right-hand-side of (14) is nonincreasing w.r.t. $n$ and the relation (12) is proved.

B.3 Proof of Theorem 5.1

Thanks to $\left(\mathcal{P}\right)$ (see Section 2.6) that is used for the observed sample and the fact that the distribution of the resampled $\varepsilon^{b}$ for $b=1,...,B$, is fixed whatever $\theta$, we have that

[TABLE]

As a consequence, by using the notation $m^{0}$, we have

[TABLE]

Since, under $P_{\theta_{0}}$, the random variables $\varepsilon^{0},...,\varepsilon^{B}$ are i.i.d. Rademacher, the process vector $([\widehat{p}^{b}_{\eta}(x)]_{x\in\mathcal{X}_{\eta}})_{b=0,...,B}$ is exchangeable (under $P_{\theta_{0}}$ and cond. on $\mathfrak{N}$), which in turn implies that the vector $(m^{b})_{b=0,...,B}$ is exchangeable (under $P_{\theta_{0}}$ and cond. on $\mathfrak{N}$). Therefore, applying Lemma 1 of (Romano and Wolf, 2005b), we obtain

[TABLE]

and it remains to integrate in $\mathfrak{N}$ to conclude.