A Sequential Scheme for Large Scale Bayesian Multiple Testing

TL;DR

This paper introduces a scalable Sequential Monte Carlo method for large-scale Bayesian multiple testing, effectively controlling false positives while incorporating covariate information, demonstrated on neural data.

Contribution

It presents a computationally efficient alternative to MCMC for Bayesian multiple testing that scales well with large datasets and uses sequential data processing.

Findings

The method matches MCMC accuracy with single data pass.

It effectively incorporates covariate information.

Demonstrated on neural recording data.

Abstract

The problem of large scale multiple testing arises in many contexts, including testing for pairwise interaction among large numbers of neurons. With advances in technologies, it has become common to record from hundreds of neurons simultaneously, and this number is growing quickly, so that the number of pairwise tests can be very large. It is important to control the rate at which false positives occur. In addition, there is sometimes information that affects the probability of a positive result for any given pair. In the case of neurons, they are more likely to have correlated activity when they are close together, and when they respond similarly to various stimuli. Recently a method was developed to control false positives when covariate information, such as distances between pairs of neurons, is available. This method, however, relies on computationally-intensive Markov Chain Monte Carlo (MCMC). Here we develop an alternative, based on Sequential Monte Carlo, which scales well with the size of the dataset. This scheme considers data items sequentially, with relevant probabilities being updated at each step. Simulation experiments demonstrate that the proposed algorithm delivers results as accurately as the previous MCMC method with only a single pass through the data. We illustrate the method by using it to analyze neural recordings from extrastriate cortex in a macaque monkey.

Tables

Table 1. TABLE I: Parameter setting for the SMC sampler for both the simulated and the real data experiments

$N_{0,0,m}$	$N_{1,0,m}$	$M$	$K_m$	${\mu }_{1,m}$	${\sigma }_{1,m}$	${\sigma }_{0,m}$
9	1	$1e4$	1	3	$\sqrt{20}$	1.5

Table 2. TABLE II: The estimation result for parameters of the null and the alternative distributions

	${\sigma }_0$	${\mu }_1$	${\sigma }_1$
estimated by SMC	0.9945	3.3470	0.7930
estimated by BFDR	1.0152	3.1070	0.1012
the true answer	1	3	0.5

Table 3. TABLE III: Number of detections and errors reported by the SMC sampler algorithm

	Declared null	Declared alternative	Total
True null	9539	38	9577
True alternative	130	293	423
Total	9669	331	10000

Table 4. TABLE IV: Number of detections and errors reported by the BFDR method

	Declared null	Declared alternative	Total
True null	9541	36	9577
True alternative	138	285	423
Total	9679	321	10000

Equations

$$z_i$$

$$f_0(z)$$

$$f_1(z)$$

$$c(x_i)$$

$$s(x_i)$$

$$(z_i|{\theta }_i)$$

$$({\theta }_i|h_i)$$

$$\text{mbox}P(h_i=1)=c(x_i)$$

$$E_{{\eta }_{t+1}}\{q(\varphi )\}$$

$$=$$

$$M^{\frac{1}{2}}\left\{\sum _{m=1}^M{\omega }_{m}q({\varphi }_m)-E_{{\eta }_{t+1}}\{q(\varphi )\}\right\}\to \mathcal{N}(0,V(q))$$

$$\varphi =\{\beta ,K,{\mu }_0,{\sigma }_0,{\mu }_1,\dots ,{\mu }_K,{\sigma }_1,\dots ,{\sigma }_K,w_1,\dots ,w_{K-1}\}.$$

$${\mu }_{0,m}$$

$${\sigma }_{0,m}$$

$$w_{k,m}=(1-{\alpha }_1)\cdot w_{k,m}+{\alpha }_1\cdot M_{k,m},k=1,\dots ,K_m,$$

$${\mu }_{k,m}$$

$${\sigma }_{k,m}^2$$

$$\varrho =\frac{{\alpha }_1}{{\alpha }_1+w_{k,m}}.$$

$$p(\beta |z_1,\dots ,z_{t+1})\simeq \sum _{i=1}^M\mathcal{N}(\beta |{\hat{\beta }}_i,b^{2}Q)$$

Full text

\IEEEkeywords

Bayesian, massive data, multiple hypothesis testing, neuroscience, Sequential Monte Carlo

I Introduction

I-A Multiple Testing

Statistical hypothesis testing is an extremely powerful approach to knowledge discovery from data [1, 2, 3, 4]. Classically, a single proposition $h$ is tested and two opposing hypotheses, namely the null hypothesis ($h=0$) and the alternative hypothesis ($h=1$), are examined. It is often desirable to employ this approach to test not just a single proposition but multiple propositions $h_{1},\ldots,h_{n}$ simultaneously, especially in many burgeoning data-intensive research fields such as genomics [5, 6], neurophysiology [7], electronic commerce [8], machine learning [3], medical and epidemiological research [9, 10], etc.

In the simplest form of multiple testing, each test yields a summary statistic $z_{i}$ and the goal is to decide which of the $z_{i}$s are signals $(h_{i}=1)$ and which are null $(h_{i}=0)$ while limiting the number of false positives. The most widely applied approach in many fields is to control the false discovery rate (FDR), i.e., the expected proportion of false positives among those null hypotheses that are rejected [11]. It has been shown that the FDR solution is closely related to the Bayesian posterior $p(h_{i}=0|\mbox{Data})$, and this alternative formulation can offer advantages in some settings [11, 12, 13]. To improve testing performance, it is sometimes possible to introduce covariate information that can appropriately modify the prior probability of each particular proposition, and the recently developed Bayesian FDR regression (BFDR) provides a formalism and implementation for this [14]. However, it relies on a Markov Chain Monte Carlo (MCMC) based iterative sampling procedure for parameter inference. The MCMC mechanism requires a complete scan of the dataset for each iteration. Hence it does not scale well to large datasets.

I-B Neural Interaction Detection: A Large Scale Multiple Testing Problem

There is great interest in examining the correlated activity among neurons recorded from various networks in the brain while an animal performs a task, as this can furnish clues about the functional operation of these networks [15, 16, 17, 18, 19, 20, 21]. However, as the number of neurons being recorded simultaneously increases, statistical methods that scale well with dataset size are needed. In experiments that will take place in the next few years, tens of millions of pairs of neurons will be recorded simultaneously. Determining the ways in which patterns of correlation change with experimental conditions thus poses a large-scale multiple testing problem.

I-C Contribution of this paper

In this paper we first translate the problem of large scale multiple testing into a sequential parameter updating process.

We then design an efficient algorithm to implement the proposed sequential scheme based on a specific parametric model that incorporates covariates information in modeling the test statistics. This algorithm features novel particle rejuvenation operations designed for updating particles defined on a joint parameter space with both discrete and continuous parameters. Each particle corresponds to a specific hypothesis on the model parameter. We conduct simulated experiment to show proof-of-concept. Results demonstrate the desirable estimation and testing performance of our algorithm and its scalability to large dataset.

Finally, we apply the proposed algorithm to analyze neural recordings from extrastriate visual cortex in a macaque monkey. The pattern of neural functional connectivity is inferred from the data. The network structure is found to vary with the experimental condition. It is also shown that the sequential feature of our algorithm allows the analyst to inspect this variation in real time. The results we provide are illustrative, as we anticipate datasets that are orders of magnitude larger in the near future.

II Model

We adopt a parametric model to describe the relationship between the test statistics and the covariates. This model is designed as follows [14]

[TABLE]

where $i$ is the index of the observations, $c(\cdot)\in(0,1)$ denotes the prior probability that $z_{i}$ is a signal, $x_{i}$ denotes the covariates, which affects $c$ in a way as specified by Eqn.(4), $J$ is the dimensionality of $x$, $\beta=[\beta_{0},\beta_{1},\ldots,\beta_{J}]$ is the regression coefficient, $f_{0}$ and $f_{1}$ respectively denote the null ($h_{i}=0$) and alternative ($h_{i}=1$) distributions, and $\pi(\theta)$ is a $K-$component mixture of Gaussians defined to be $\pi(\theta)\triangleq\sum_{j=1}^{K}w_{j}\mathcal{N}(\theta|\mu_{j},\sigma_{j}^{2})$. The above model can be rewritten in a hierarchical form

[TABLE]

where $\delta_{a}$ denotes the Dirac-delta function located at $a$. Such a hierarchical model together with priors for $\beta$ and the unknown distribution $\pi(\theta)$, defines a joint posterior distribution over all model parameters. A full Bayesian approach, namely BFDR, has been developed to sample from the posterior [14]. As a Bayesian approach, the BFDR method provides a natural way to quantify uncertainty about the regression function $s(x)$ and $\pi(\theta)$ jointly, while this is counterbalanced by the additional computational complexity of the fully Bayesian method. As mentioned above in the Introduction section, the main building block of the BFDR method is MCMC, which consists of numerous iterations and each iteration requires the whole dataset to be processed. The nature of the MCMC methods also determines that, if new data items are added into the dataset, the whole sampling and computation process should be re-started from scratch to update the testing result. This cumbersome computing feature makes the BFDR method unsuitable for large-scale multiple testing problems. To this end, we propose a sequential computation scheme in the next section, which is also full Bayesian, while scaling well to large dataset.

III The Proposed Sequential Scheme for Multiple Testing

Given a parametric model of the test statistics, here we show how to perform multiple testing in a sequential manner. Denote $\phi$ as the unknown model parameter that must be estimated. Denote $t$, $t\in\mathbb{Z},t>0$, as a time variable and let $\eta_{t}\triangleq p(\phi|z_{1},\ldots,z_{t})$ represent the target probabilistic density function (pdf) at time step $t$. As the proposed scheme processes data items one by one in a sequential manner, $z_{t}$ essentially represents the $t$-th data item that is processed. If the sample size of the dataset to be analyzed is $n$, then the target pdf of our final concern is $\eta_{n}$.

We now consider a series of target densities, $\eta_{1},\ldots,\eta_{n}$, and focus on the question how to sample from them sequentially. Suppose that, at time $t+1$, a sampling mechanism is not readily available for the target density $\eta_{t+1}$, but one is available for another sampling density $\eta_{t}$, then we can use importance sampling and write [22]

[TABLE]

where $q$ is any measurable function such that the left hand side of Eqn.(9) exists, $\phi_{m}$ is a draw from $\eta_{t}(\phi)$ and $\hat{\omega}_{m}\propto\eta_{t+1}(\phi_{m})/\eta_{t}(\phi_{m})$. The weights $\hat{\omega}_{m}$ can be known only up to a multiplicative constant, which is cancelled by the denominator $\sum_{m=1}^{M}\hat{\omega}_{m}$. According to importance sampling theory, we have

[TABLE]

in distribution, where $V(q)=\mbox{var}_{\eta_{t}}[\{\eta_{t+1}(\phi)/\eta_{t}(\phi)\}\{q(\phi)-E_{\eta_{t+1}}q(\phi)\}]$ [22], $\omega_{m}=\hat{\omega}_{m}/\sum_{i=1}^{M}\hat{\omega}_{i}$. The weighting operation thus shifts the target density of the particle set from $\eta_{t}$ to $\eta_{t+1}$. Given an initial target density $\eta_{0}$, we perform the above sampling and weighting operations sequentially, shifting the target of interest from $\eta_{0}$ to $\eta_{1}$, from $\eta_{1}$ to $\eta_{2}$ and so on, until from $\eta_{n-1}$ to $\eta_{n}$. Such a sequential importance sampling (SIS) procedure constitutes the basic building block of a Sequential Monte Carlo (SMC) method, which has gained a lot of attentions in the fields of statistics and signal processing [23, 24].

In the basic SIS procedure, as long as the target density evolves one time step forward, say from $t$ to $t+1$, a reweighting operation $\omega_{m}^{\prime}=\hat{\omega}_{m}\eta_{t+1}(\phi_{m})/\eta_{t}(\phi_{m})$ will be performed. As this reweighting operation iterates, fewer and fewer particles will retain significant weights. This phenomenon is termed particle degeneracy or impoverishment. A common strategy to avoid that is to perform resampling followed by a rejuvenation step [25]. The resamping operation duplicates particles with higher weights and discards those with lower weights [26, 27]. The rejuvenation step is operated on the resampled particles, moving them according to a Markov transitional kernel with stationary distribution set to be the target density at current time step. In this way, the target density of the particles is not changed, but more refreshed particles are generated. Therefore the phenomenon of particle degeneracy may be reduced a lot. Under mild regularity conditions, rigorous convergence results for general SMC methods have been proved [28].

Here we propose to employ SMC, instead of MCMC, to perform multiple testing. One iteration of the proposed scheme is presented in Algorithm 1 as follows.

IV Algorithm Design

Here we design an efficient algorithm to implement the proposed scheme presented in Section 3 based on the model described in Section 2. Recall that the model is defined by Eqns. (1)-(5). The unknown parameter to be estimated is

[TABLE]

As $w_{K}=1-\sum_{j=1}^{K-1}w_{j}$, it is not included in $\phi$. In contrast with standard Bayesian sequential filtering problems, one notable feature of our problem lies in that the dimensionality of the model parameter is a variable, whose value depends on the number of mixing components $K$. Therefore, to implement the sequential scheme presented in Algorithm 1, we should take into account of the cases, wherein the particles are assigned with different $K$ values and thus have different dimensions. To our knowledge, these is seldom algorithm that can operate on a set of particles with different dimensions in the SMC literature. Additionally, when $K$ is a large number, the dimensionality of $\phi$ becomes large, then a traditional SMC method may collapse due to the curse of dimensionality [29]. Hence it is far from trivial to develop an efficient algorithm to implement the proposed scheme based on the model we are concerned here. We describe our algorithm design in what follows.

IV-A The weighting and resampling steps

The weighting step involves the calculation of $\omega_{m}$. Based on Bayesian theorem, we have $\eta_{t+1}(\phi)=\eta_{t}(\phi)\cdot p(z_{t+1}|\phi)/C$, where $C$ is the normalizing constant. Hence we just set $\hat{\omega}_{m}=p(z_{t+1}|\phi_{m})$, which can be calculated directly based on Eqns.(1)-(5), and then we get $\omega_{m}=\hat{\omega}_{m}/\sum_{j=1}^{M}\hat{\omega}_{j}$. To implement the resampling step, we select the unbiased residual resampling method [30]. Other widely used resampling methods include the multinomial selection method [31] and the stratified resampling method [32].

IV-B Particle rejuvenation step

We develop an efficient approach to implement the particle rejuvenation step of Algorithm 1. It is the ad hoc design of the particle rejuvenation algorithm that discriminates our SMC sampler from all the other related SMC methods in the literature. First we separate the components of $\phi$ into two parts. Let $\phi=\{\phi^{1},\phi^{2}\}$, where $\phi^{1}=\beta$ and $\phi^{2}$ is composed of the rest of parameters in $\phi$. So $\phi^{1}$ is of fixed dimensionality and the parameters included in $\phi^{2}$ correspond to a mixture model that consists of the null distribution and the alternative distribution. In the model defined with Eqns. (1)-(5), the null distribution is normal and the alternative distribution is a $K$-component normal mixture. Hence the dimensionality of $\phi^{2}$ is variable depending on the value of $K$. We design rejuvenation operations for $\phi^{1}$ and $\phi^{2}$, respectively. Suppose that $\phi_{m}=\{\phi_{m}^{1},\phi_{m}^{2}\}$, a random draw from $\eta_{t}$, is available, we now take this particle as an example to describe the designed particle rejuvenation operations that are performed at time step $t+1$.

IV-B1 Particle rejuvenation operations for $\phi_{m}^{2}$

Now we describe how to update $\phi_{m}^{2}$ based on the new data item $z_{t+1}$ at current time step. First we calculate the prior probability that $z_{t+1}$ is a signal, through Eqn. (4), based on the current hypothesis on the model parameter $\phi_{m}$. If this probability is less than 0.5, we allocate $z_{t+1}$ into the null distribution associated with the $m$th particle, represented by $f_{0,m}$, and update its parameters $\mu_{0,m}$ and $\sigma_{0,m}$ as follows

[TABLE]

where $\alpha_{0}=1/(1+N_{0,t,m})$ denotes the weight assigned to $z_{t+1}$. $N_{0,t,m}$ denotes the number of $z_{i}$s that were allocated into $f_{0,m}$ in the previous $t$ time steps.

If the aforementioned prior probability is greater than or equal to 0.5, we allocate $z_{t+1}$ into $f_{1,m}$. Now we present how to update the parameter value of $f_{1,m}$ if $z_{t+1}$ is allocated into it. We develop an online $K$-means approximation method to update the parameter value of $f_{1,m}$. First the test statistic $z_{t+1}$ is checked against the existing $K_{m}$ Gaussian distributions of $f_{1,m}$, until a match is found. A match is defined as the value of $z_{t+1}$ falling within 2.5 standard deviations of a distribution. The threshold 2.5 is also used in an online method for background modeling in visual tracking problems [33]. Then we adjust the prior weights of the $K_{m}$ distributions, $w_{1,m},\ldots,w_{K,m}$, as follows

[TABLE]

where $M_{k,m}$ is 1 for the model which matched and 0 for the remaining models, $\alpha_{1}=1/(1+N_{1,t,m})$ is the weight assigned to $z_{t+1}$. $N_{1,t,m}$ denotes the number of $z_{i}$s that were allocated into $f_{1,m}$ in the previous $t$ time steps. If none of the $K_{m}$ distributions match $z_{t+1}$, a new probable distribution is added to $f_{1,m}$, with $z_{t+1}$ as its mean value, an initially high variance, and a prior weight set at $\alpha_{1}$. As a new component is added, $K_{m}$ is updated with $K_{m}=K_{m}+1$. After the above updating operation, the weights of mixture components are renormalized. The $\mu$ and $\sigma$ parameters of unmatched mixture components remain the same. The parameters of the component that matches $z_{t+1}$ are updated as follows

[TABLE]

where

[TABLE]

The above design of the mixture updating operations is inspired by a $K$-means approximation method used for background modeling in visual object tracking problems [33], while our method is different from that in [33] in several ways. First, our method allows the number $K$ to be adaptable, while, the value of $K$ is fixed in [33]. Second, in our method, the learning rates, $\alpha_{0}$, $\alpha_{1}$ and $\varrho$, are controlled through $N_{0,t,m}$ and $N_{1,t,m}$ adaptively. In contrast, for the $K$-means approximation method in [33], these learning rates are set as free parameters, whose values are specified empirically without clear explanations. Our method uses $N_{0,t,m}$ and $N_{1,t,m}$ to record the numbers of $z_{i}$s that have been involved in the previous $t$ time steps for constructing $f_{0,m}$ and $f_{1,m}$, respectively. When $z_{t+1}$ arrives, the amount of the relative contribution of $z_{t+1}$ in updating $f_{0,m}$ (resp. $f_{1,m}$) is encoded by $\alpha_{0}$ (resp. $\alpha_{1}$), which is a function of $N_{0,t,m}$ (resp. $N_{1,t,m}$). As the values of $N_{0,t,m}$ and $N_{1,t,m}$ increase over time, the relative contribution of $z_{t+1}$ in updating the parameters of $f_{0,m}$ and $f_{1,m}$ is weakened automatically over time.

IV-B2 Particle rejuvenation operations for $\phi_{m}^{1}$

After updating $\phi_{m}^{2}$, we rejuvenate $\phi_{m}^{1}$, i.e., $\beta_{m}$. First, based on the set of resampled $\beta_{i}$s, $i=1,\ldots,M$, a kernel smoothing method [34] is adopted here to approximate the posterior $p(\beta|z_{1},\ldots,z_{t+1})$ by a series of Gaussian kernel functions:

[TABLE]

where $Q$ is the sample Monte Carlo variance, $b=\left(\frac{4}{(d+2)M}\right)^{\frac{1}{d+4}}$ is the kernel bandwidth, $d$ is the dimensionality of $\beta$, $\hat{\beta}_{i}=a\cdot\beta_{i}+(1-a)\cdot\bar{\beta}$, $a=\sqrt{(1-b^{2})}$ and $\bar{\beta}$ is the current Monte Carlo mean $\beta$ value. This choice of bandwidth has been proved to be asymptotically optimal for multivariate-normal distribution cases [35]. Now we set $\beta_{m}$ as a random draw from $\mathcal{N}(\beta|\hat{\beta}_{m},b^{2}V)$. All particle weights of these rejuvenated $\beta_{m}$s are equally set at $1/M$ because the resample step preceded this rejuvenation step. These rejuvenated $\beta_{m}$s now represent a more diverse set of $\beta$ values, and, under mild conditions, these samples converge to those drawn directly from the target density with theoretic guarantees [36, 34].

IV-C Initialization

To run Algorithm 1, $\eta_{0}(\phi)$ along with a weighted set of particles $\{\phi_{m},\omega_{m}\}_{m=1}^{M}$ drawn from it, $N_{0,0,m}$ and $N_{1,0,m}$ must be initialized. We consider two cases for initialization. In the first case, we assume that a historical dataset, which is collected under the same condition as that of the current dataset, has been processed and a corresponding posterior is available, then we just set $\eta_{0}(\phi)$, $N_{0,0,m}$ and $N_{1,0,m}$ based on this posterior. In another word, the posterior corresponding to the historical dataset is taken as the prior for analyzing the current dataset. In the other case, there is no information from a historical dataset available for use. Then we resort to domain knowledge to do initialization. If the domain knowledge contains little information, we can select a noninformative prior, e.g., a diffusing prior over the parameter space, for use. If there is rich domain knowledge, then prior elicitation techniques can be used for distilling the prior from the domain knowledge [37, 38]. Specifying the prior is an inevitable step to apply Bayesian methods, while, it is noteworthy to mention that, as more and more observations have been processed, the impact of the prior on the final posterior inference will become weaker and weaker, as long as the prior is properly formulated [39].

V Simulated Data

The proposed SMC sampler algorithm is first tested on synthetic dataset in order to verify its efficiency and merits, and confirm its ability to yield meaningful results when applied to experimental spike train data. A number $n=10000$ of test statistics $z_{i},i=1,\ldots,n$, are drawn according to the covariate-dependent mixture model defined with Eqns.(1)-(5). The dataset has two covariates $x=(x^{1},x^{2})$ and the regression coefficient $\beta=(-3.5,\frac{\sqrt{2}}{2},\frac{\sqrt{2}}{2})$. The null and the alternative distributions are specified as $f_{0}(z)=\mathcal{N}(z|0,1)$ and $f_{1}(z)=\mathcal{N}(z|3,0.5^{2})$, respectively.

We select the BFDR method as a baseline for performance comparison, as it represents the state-of-the-art multiple testing method for the neural interaction detection problem. We implement the BFDR method in full accordance with its companion R package FDRreg [14], where the number of mixing components $K$ in $f_{1}$ is chosen based on the Akaike information criterion via a preliminary run of an Expectation-Maximization algorithm; the other components of $\phi$ are estimated by MCMC. The MCMC procedure consists of 2200 iterations of sampling, and the first 200 iterations are taken as the burn-in period. For the proposed SMC sampler, its parameter setting is presented in Table 1. The parameter value, including the ratio of $N_{0,0,m}$ over $N_{1,0,m}$, is selected according to domain knowledge. The numbers 9 and 1 represent the minimal feasible integers that can be used to initialize $N_{0,0,m}$ and $N_{1,0,m}$. For both BFDR and SMC, the null distribution is fixed to be a zero-mean Gaussian distribution. For BFDR, the variance of the null distribution is chosen using an offline empirical approach as suggested in [40]. For both BFDR and SMC, a diffusing uniform prior of $\beta$ is selected for use.

The result of the proposed SMC sampler in terms of posterior histograms for the estimate of $\beta_{0}$ is displayed in Fig.1. It is shown that, the distributional estimate of $\beta_{0}$, given by the SMC sampler, is comparable with that by the BFDR method.

In Fig.2, we plot traces of the maximum a posterior (MAP) estimate of $\beta$ obtained from an example run of the SMC sampler. A good convergence performance of the SMC sampler in estimating $\beta$ is confirmed. The estimation results on $f_{0}$ and $f_{1}$ are presented in Table 2. The estimate of $\pi(\theta)$ given by the SMC sampler has two distribution components, while one of them has a tiny weight $0.0027$, so only parameters of the dominant component are listed here. We see that the SMC sampler gives a more accurate estimate on $\sigma_{0}$ and $\sigma_{1}$, and a worse estimate on $\mu_{1}$, as compared with the BFDR method.

The numbers of detections and errors obtained by the SMC sampler and BFDR method based on the same synthetic dataset are listed in Tables 3 and 4, respectively. For BFDR, the FDR is controlled roughly at $10\%$ and the real FDR obtained in the experiment is $36/321=11.2\%$. The FDR given by the SMC sampler is $38/331=11.5\%$, which is slightly higher than that of BFDR, while the SMC sampler detects 8 more true alternatives than BFDR. Taking both factors of FDR and detection power into consideration, we argue that the SMC sampler gives almost an equivalent testing performance as that of the BFDR method. The significant difference between the SMC sampler and BFDR lies in that the former accessed each test statistic $z_{i}$ only once, while the latter accessed each $z_{i}$ for 2200 times, during this experiment. Assume that a new data item is added to the dataset for analysis. If the BFDR method is under use, the nature of the MCMC sampling mechanism determines that the analysis has to start the whole analysis from scratch, which means that all the historical data have to be re-processed. If the proposed SMC sampler is under use, the analyst can just run one iteration of the SMC sampler, which only needs to access this new data item, without having to retracing any historical data, to update the model parameter as well as the testing result. As demonstrated by the above simulated experiment, the updated result yielded by the SMC sampler will be as accurate as that obtained by applying the BFDR method to analyze the whole updated dataset.

Finally, we conduct a simulated experiment to compare the proposed SMC sampler and the BFDR method in the computational time. Both methods are coded with MATLAB and run on a 16-core microprocessor. The comparison result is depicted in Fig. 3. We see that, the greater the amount of data to be analyzed, the more obvious advantages for the SMC sampler in terms of computation time. For the SMC sampler, the particle rejuvenation and weighting steps are straightforward to parallelise, as they require only independent operations on each particle. Hence these two steps are parallelized in our experiment. The resampling step is more difficult to parallelise, as it requires a collective operation, such as a sum across particle weights. Hence it is not parallelized. The main building block of the BFDR method is an MCMC based iterative sampling procedure, which works in a serial way, therefore there is no straightforward way to parallelise it.

VI Experimental data

The experimental data were sampled from a real experimental session where a rhesus macaque monkey was performing a visual spatial selective attention task. Experimental procedures were approved by the Institutional Animal Care and Use Committee of the University of Pittsburgh. Two task variables are involved. One is orientation. A Gabor stimulus presented in the receptive field is associated with one of two possible orientations, corresponding to 135 degrees and 45 degrees, respectively. We use ‘ori 1’ and ‘ori 2’ to denote these two orientations, respectively. The other task variable is cue. The monkey was cued to either attend towards or away from the receptive field. We use ‘cue 1’ and ‘cue 2’ to denote these two cues, respectively. These two task variables were combined into four stimulus conditions, namely $\{$ori 1, cue 1$\}$, $\{$ori 2, cue 1$\}$, $\{$ori 1, cue 2$\}$ and $\{$ori 2, cue 2$\}$. A number of repeated trials were conducted under each stimulus condition. The numbers of repeated trials are respectively 666, 759, 954, and 1077, for those four stimulus conditions. A Utah array, which consists of a 10 by 10 grid of electrodes, was implanted in area V4 of the monkey’s extrastriate visual cortex, and the spike train records were sampled from the same set of neurons across those stimulus conditions (block-randomized).

The test statistics $z_{i}$s are obtained by applying the Fisher transformation on spike count correlations [41, 42, 15, 17]. Each spike count correlation is a Pearson correlation of spike counts of a pair of neurons. For each test statistic, there are two relevant covariates: (1) inter-neuron distance, measured in micrometers; and (2) tuning-curve correlation (TCC). The dataset to be analyzed consists of four time windows, each corresponding to a specific stimulus condition.

Given the above dataset, the question is how to reveal the neuron network structure, which is determined by neuron interactions, from data. It is also expected to check if this structure varies over the change of the stimulus condition in real time.

The proposed SMC sampler is applied here to test hypotheses on whether pairs of neurons exhibit fine-time-scale ($\sim 500ms$) interactions. We adopt a standard measure termed effective sample size (ESS) to evaluate the reliability of the algorithm’s inference result [43]. This measure is defined to be $\mbox{ESS}=1/\{\sum_{m=1}^{M}(\omega_{m})^{2}\}$, which satisfies $1\leq\mbox{ESS}\leq M$. The meaning of the ESS can be stated as that, the inference given by the SMC sampler based on the $M$ weighted particles is equivalent with that obtained based on the number ESS of particles drawn directly from the target distribution [30]. So intuitively, a greater value of ESS indicates a more reliable inference result given by the SMC sampler, and vice versa. Because the ESS is a function of the particle size $M$, we employ the normalized ESS (NESS), $\mbox{NESS}=\mbox{ESS}/M$, in practice. At the transition moments of two neighboring time windows, significant declines in the NESS value appear, because the regularity assumption between the neighboring target pdfs is violated. As long as the NESS value drops below $0.1$, the SMC sampler is re-initialized and then the NESS is re-calculated. The way to initialize SMC sampler here is the same as presented in Section V. We plot the resulting NESS per time step in Fig. 4, which shows that the NESS value maintains at a relatively high level, indicating a reliable inference result of the SMC sampler.

The neuron network architecture in terms of detected neuronal interactions is inferred from data and visually displayed in Fig.5. Under four different stimulus conditions, corresponding to sub-figures (a), (b), (c), (d) of Fig.5, there are respectively 899, 723, 282 and 1722 pairwise neuronal interactions, which are detected by the SMC sampler algorithm. So it is confirmed that the neuron network structure varies long with the change in the stimulus condition. Specifically, under the first stimulus condition $\{$ori 1, cue 1$\}$, namely, when the monkey was cued to attend towards the receptive field and the stimulus orientation was 135 degrees, the number of detected neural interactions is 899. Then, when the orientation was changed to be 45 degrees, corresponding to stimulus condition $\{$ori 2, cue 1$\}$, the number of detected interactions reduces to 723. Next the monkey was cued to away from the receptive field and the stimulus orientation was changed to be 135 degrees, namely the stimulus condition $\{$ori 1, cue 2$\}$ was applied. As a result, the number of detected neural interactions reduces to 282. Lastly, the stimulus orientation was changed to be 45 degrees, corresponding to stimulus condition $\{$ori 2, cue 2$\}$, and then the number of detected neural interactions increases rapidly from 282 to 1722. As is shown, the proposed algorithm allows the analyst to inspect the evolution of the neuron interaction network structure in real time, and the correlated pattern between the stimulus condition and the network structure can be analyzed in time.

Finally, we plot the histograms of the covariate TCC associated with different declared hypotheses, in Fig. 6. We see that there is a clear correlation between the TCC and the values of the declared hypotheses, while the correlation structure varies along with the change in the stimulus condition. Specifically, under every stimulus condition, the distribution of the TCCs, associated with the declared alternative hypotheses, shrinks to the right, that is, the direction of 1; while the degrees of shrinkage associated with stimulus conditions with ‘ori 1’ are slightly greater than those with ‘ori 2’. Under all stimulus conditions except $\{$ori 2, cue 2$\}$, the distribution of the TCCs, associated with the declared null hypotheses, almost evenly covers the entire parameter space, while under stimulus condition $\{$ori 2,cue 2$\}$, this distribution shrinks to the left, that is, the direction of minus 1.

Traditional approaches have largely focused on a single scale to study the neuron network architecture. Although being illustrative, the above results indicate a potential of the proposed algorithm to expand beyond traditional approaches to capture the neuron network architecture on more scales, which include not only the time and space scales, but also the functional scale determined by the covariates.

VII Conclusions

We designed and implemented an SMC scheme for large scale multiple testing in the context of BFDR, which is based on a hierarchical parametric model that takes account of covariate information. We showed that the performance of our algorithm can be as good as the existing MCMC procedure while employing only a single pass through the data, and that, after parallelization, it greatly reduces computation time. We also illustrated the way this method could be used with neural recordings to help identify network structure over time. In other testing seŠttings the particle rejuvenation step would need to be re-designed, but the SMC scheme would still be applicable and could easily be modi€fied. ŒThe sequential scheme proposed here thus opens the door to efficient computation in many large scale multiple testing problems.

Acknowledgment

This research has been supported by National NSFC Grant 61571238 (B. Liu), NIMH Grant RO1 064537 (R. E. Kass), NIH Grants 5R90DA023426-10 (G. Vinci), K99EY025768 (A. C. Snyder) and R01EY022928 (M. A. Smith).