False discovery rate regression: an application to neural synchrony detection in primary visual cortex

1.1 Multiple testing in the presence of covariates

The problem of multiple testing concerns a group of related null hypotheses h₁, …, h_n that are tested simultaneously. In its simplest form, each test yields a summary statistic z_i, and the goal is to decide which of the z_i are signals (h_i = 1) and which are null (h_i = 0). Many solutions to this problem, such as Bonferroni correction, aim to control the family-wise error rate (FWER): the probability of incorrectly rejecting at least one null hypothesis, assuming that they are all true. An alternative, which has become the dominant approach in many domains of application, is to control the false discovery rate (FDR): the proportion of false positives among those null hypotheses that are rejected (Benjamini and Hochberg, 1995). Regardless of which error rate they aim to control, however, most existing approaches obey a monotonicity property: if test statistic z_i is declared significant, and z_j is more extreme than z_i, then z_j is also declared significant. Yet in many cases, we have auxiliary covariate information about each test statistic, such as location in the brain or distance along a chromosome. If significant test statistics tend to cluster in covariate space, then monotonicity becomes undesirable, and a procedure that takes account of the covariate should perform better. In this paper, we introduce a method called false-discovery-rate regression (FDRR) that incorporates covariates directly into the multiple-testing problem.

The method we describe here builds on the two-groups model (Efron et al., 2001), a popular framework for controlling the false-discovery rate. In the two-groups model, some small fraction c of the test statistics are assumed to come from an unknown signal population, and the remainder from a known null population. Our proposal is to allow the mixing fraction c to depend upon covariates, and to estimate the form of this dependence from the data. Extensive simulation evidence shows that, by relaxing the monotonicity property in a data-dependent way, FDR regression can improve power while still controlling the global false-discovery rate. The method is implemented in the publicly available R package FDRreg (Scott, 2014).

Our motivating application is the identification of interactions among many simultaneously recorded neurons, which has become a central issue in computational neuroscience. Specifically, we use FDR regression to detect fine-time-scale neural interactions (“synchrony”) among 128 units (either single neurons or multi-unit groups) recorded simultaneously from the primary visual cortex (V1) of a rhesus macaque monkey (Kelly et al., 2010; Kelly and Kass, 2012). The experiment from which the data are drawn produced thousands of pairs of neurons, each involving a single null hypothesis of no interaction. In this case, combining all tests into a single FDR-controlling analysis would inappropriately ignore the known spatial and functional relationships among the neurons (e.g. Smith and Kohn, 2008). Our approach for false-discovery rate regression avoids this problem: it detects roughly 50% more significant neuron pairs compared with a standard analysis by exploiting the fact that spatially and functionally related neurons are more likely to exhibit synchronous firing.

1.2 The two-groups model

In the two-groups model for multiple testing, one assumes that test statistics z₁, …, z_n arise from the mixture

where c ∈ (0, 1), and where f₀ and f₁ respectively describe the null (h_i = 0) and alternative (h_i = 1) distributions of the test statistics. For each z_i, one then reports the quantity

As Efron (2008a) observed, the information contained in w_i provides a tidy methodological unification to the multiple-testing problem. Bayesians may interpret w_i as the posterior probability that z_i is a signal, while frequentists may interpret 1 − w_i as a local false-discovery rate. The global false-discovery rate of some set Z₁ of putative signals can then estimated as

Efron et al. (2001) show that this Bayesian formulation of FDR is biased upward as an estimate of frequentist FDR, and therefore conservative.

An elegant property of the local FDR approach is that it is both frequentist and fully conditional: it yields valid error rates, yet also provides a measure of significance that depends on the precise value of z_i, and not merely its inclusion in a larger set (c.f. Berger, 2003). This can be achieved, moreover, at little computational cost. To see how, observe that (2) may be re-expressed in marginalized form as

where f(z) = c · f₁(z_i) + (1 − c) · f₀(z_i) is the overall marginal density. Importantly, f(z) can be estimated from the empirical distribution of the test statistics; this is typically quite smooth, which makes estimating f(z) notably easier than a generic density-estimation problem. Therefore one may compute local FDR using cheap plug-in estimates $\widehat{f}\left(z\right)$ and $\widehat{c}$, and avoid the difficult deconvolution problem that would have to be solved in order to find f₁(z) explicitly (e.g. Efron et al., 2001; Newton, 2002; Martin and Tokdar, 2012).

1.3 FDR regression

Implicit in the two-groups model is the assumption that all tests should be combined into a single analysis with a common mixing weight c in (1). Yet for some data sets, this may be highly dubious. In our analysis of neural recordings, for example, a test statistic z_i is a measure of pairwise synchrony in the firing rates of two neurons recorded from an array of electrodes, and these z_i’s exhibit spatial dependence across the array: two nearby neurons are more likely to fire synchronously than are two neurons at a great distance. Similar considerations are likely to arise in many applications.

False-discovery-rate regression addresses this problem through a conceptually simple modification of (1), in which covariates x_i may affect the prior probability that z_i is a signal. In its most general form, the model assumes that

for an unknown regression function s(x) and known link function $G:\mathcal{R}\to (0,1)$.

This new model poses two main challenges versus the ordinary two-groups model (1). First, we must estimate a regression model for an unobserved binary outcome: whether z_i comes from f₁, and is therefore a signal. Second, because each mixing weight in (4) depends on x_i, there is no longer a common mixture distribution f(z) for all the test statistics. We therefore cannot express the Bayes probabilities in marginalized form (3), and cannot avoid estimating f₁(z) directly.

Our approach, described in detail in Section 2, is to represent f₁(z) as a location mixture of the null density, here assumed to be a Gaussian distribution:

Even in the absence of covariates, estimating the mixing density π(θ) is known to be a challenging problem, because Gaussian convolution heavily blurs out any peaks in the prior. We consider two ways of proceeding. The first is an empirical-Bayes method in which an initial plug-in estimate $\widehat{\pi }\left(\theta \right)$ is fit via predictive recursion (Newton, 2002). The regression function is then estimated by an expectation-maximization (EM) algorithm, treating $\widehat{\pi }\left(\theta \right)$ as fixed. The second is a fully Bayes method in which π(θ) and the regression function s(x) are estimated jointly using Markov-chain Monte Carlo. In simulation studies, both methods lead to better power and equally strong protection against false discoveries compared with traditional FDR-controlling approaches.

The rest of the paper proceeds as follows. The remainder of Section 1 contains a brief review of the literature on multiple testing. Section 2 describes both empirical-Bayes and fully Bayes methods for fitting the FDR regression model, and draws connections with existing approaches for controlling the false-discovery rate. It also describes how existing methods for fitting an empirical null hypothesis may be combined with the new approach (Section 2.4). Section 3 shows the results of a simulation study that validates the frequentist performance of the method. Section 4 provides background information on the neural sychrony-detection problem. Section 5 shows the results of applying FDR regression to the sychrony-detection data set. Section 6 contains discussion.

1.4 Connection with existing work

Our approach is based on the two-groups model, and therefore in the spirit of much previous work on Bayes and empirical-Bayes multiple testing, including Efron et al. (2001), Johnstone and Silverman (2004), Scott and Berger (2006), Muller et al. (2006), Efron (2008a,b), and Bogdan et al. (2008). The final reference has a comprehensive bibliography. We will make some of these connections more explicit when they arise in subsequent sections.

Other authors have considered the problem of multiple testing in the presence of correlation (e.g. Clarke and Hall, 2009; Fan et al., 2012). The focus there is on making the resulting conclusions robust to unknown correlation structure among the test statistics. Because it explicitly uses covariates to inform the outcome of each test, FDR regression is different both in aim and execution from these approaches.

On the computational side, we also draw upon a number of recent innovations. Our empirical-Bayes approach uses predictive recursion, a fast and efficient method for estimating a mixing distribution (Newton, 2002; Tokdar et al., 2009; Martin and Tokdar, 2012). Our fully Bayes approach requires drawing posterior samples from a hierarchical logistic-regression model, for which we exploit the Pólya-Gamma data-augmentation scheme introduced by Polson et al. (2013).

There is also a growing body of work on density regression, where an unknown probability distribution is allowed to change flexibly with covariates using nonparametric mixture models (e.g. Dunson et al., 2007). We do not attempt a comprehensive review of this literature, which has goals that are quite different from the present application. For example, one of the key issues that arises in multiple testing is the need to limit the flexibility of the model so that the null and alternative hypotheses are identifiable. Ensuring this property is not trivial; see Martin and Tokdar (2012). In density regression, on the other hand, only the overall density is of interest; the mixture components themselves are rarely identifiable.

Our application draws most directly on Kelly et al. (2010) and Kelly and Kass (2012). We review other relevant neuroscience literature in Section 4.

2.1 An empirical-Bayes approach

We use a version of the FDR regression model where

We have assumed a logistic link and a Gaussian error model, both of which could be modified to suit a different problem. We also assume a linear regression where s(x) = x^T β, and therefore model non-linear functions by incorporating a flexible basis set into the covariate vector x. Both μ and σ² are initially assumed to be known; in Section 2.4, we describe how to weaken this assumption by estimating an empirical null, in the spirit of Efron (2004).

The unknown parameters of the FDR regression model that must be estimated are the regression coefficients β and the mixing distribution π(θ). This section describes two

methods—one empirical-Bayes, one fully Bayes—for doing so. Both methods are implemented in the R package FDRreg.

Our empirical-Bayes approach begins with a pre-computed plug-in estimate for π(θ) in (5), ignoring the covariates. This is equivalent to assuming that c(x_i) ≡ c for all i, albeit only for the purpose of estimating π(θ). Many methods could be used for this purpose, including finite mixture models. We recommend the predictive-recursion algorithm of Newton (2002) for two reasons: speed, and the strong guarantees of accuracy proven by Tokdar et al. (2009). Predictive recursion generates a nonparametric estimate $\widehat{\pi }\left(\theta \right)$, and therefore an estimate ${\widehat{f}}_1\left(z\right)$ for the marginal density under the alternative, after a small number of passes (typically 5–10) through the data. The algorithm itself is similar to stochastic gradient descent, and is reviewed in Appendix A.

Upon fixing this empirical-Bayes estimate ${\widehat{f}}_1\left(z\right)$, and assuming that f₀(z) is known, we can fit the FDR regression model by expectation-maximization (Dempster et al., 1977). To carry this out, we introduce binary latent variables h_i such that

Marginalizing out each h_i clearly recovers the original model (4). The complete-data log-likelihood for β is

This is a smooth, concave function of β whose gradient and Hessian matrix are available in closed form. It is therefore easily maximized using standard methods, such as the Newton–Raphson algorithm. Moreover, l(β) is linear in h_i, and the conditional expected value for h_i, given β, is just the conditional probability that h_i = 1:

These facts lead to a simple EM algorithm for fitting the model (Algorithm 1; see box).

Thus the overall approach for estimating Model (5) has three steps.

• (1)

  Fix μ and σ² under the null hypothesis, or estimate an empirical null (see section 2.4), thereby defining f₀(z).

• (2)

  Use predictive recursion to estimate π(θ), and therefore f₁(z), under the two-groups model without covariate effects (see Appendix A).

• (3)

  Use f₀(z) and f₁(z) in Algorithm 1 to estimate w_i and the regression coefficients.

In principle, the estimate for π(θ) could be improved by using the information in the covariates. Despite this, our experiments show that the empirical-Bayes approach is essentially just as effective as using a full Bayes approach to estimate π(θ) and the regression function jointly. This can be explained by the fact that Gaussian deconvolution is such a poorly conditioned problem: for any finite set of observations from f₁(z) in (5), there is a large “near-equivalence” class of approximately compatible priors π(θ). Because the Bayes oracle in (6) depends on π(θ) only through f₁(z), any prior in this near-equivalence class will yield nearly the same posterior probabilities. In this sense, predictive recursion seems to provide a “good enough” estimate for π(θ), despite ignoring the covariates.

2.2 Empirical Bayes with marginalization

In Section 2.1, we ignored the covariates to estimate f₁(z). But the most direct generalization of the local-FDR approach of Efron et al. (2001) would be to ignore the covariates and estimate the overall mixture density f(z) instead. We now explain why this is a poor solution to the FDR regression problem. Let us begin with the key insight in Efron et al.’s approach to estimating local FDR, which is that the marginal f(z) is common to all test statistics, and that it can be estimated well using the empirical distribution of the z_i, without explicit deconvolution of the mixture. This motivates a simple empirical-Bayes strategy: (i) compute a nonparametric estimate $\widehat{f}\left(z\right)$ of the common marginal density, along with a likelihood- or moment-based estimate $\widehat{c}$ of the mixing fraction; and (ii) plug $\widehat{f}$ and $\widehat{c}$ into the marginalized form of the posterior probability (3) to get local FDR for each test statistic. One caveat is that $\widehat{c}$ must be chosen to ensure that (3) falls on the unit interval for all i. But this is an easy constraint to impose during estimation. See Chapter 5 of Efron (2012) for further details.

We have already remarked that the posterior probabilities in the FDR regression model (4) do not share a common mixture density f(z), and so cannot be expressed in marginalized form. Nonetheless, it is natural to wonder what happens if we simply ignore this fact, estimate a global $\widehat{f}\left(z\right)$ from the empirical distribution of the test statistics, and use

in lieu of expression (6) used in Algorithm 1. This has the seemingly desirable feature that it avoids the difficulties of explicit deconvolution. But (7) is not guaranteed to lie on the unit interval, and constraining it to do so is much more complicated than in the no-covariates case (3). Moreover, simply truncating (7) to the unit interval during the course of the estimation procedure leads to very poor answers.

In our simulation studies, we do consider the following ad-hoc modification of (7), in an attempt to mimic the original local-FDR approach as closely as possible:

Here T_u(a) is the projection of a to the unit interval, while $\widehat{f}\left(z\right)$ and $\widehat{c}$ are plug-in estimates of the marginal density and the mixing fraction using the no-covariates method described in Chapter 5 of Efron (2012). In our simulation studies, this modification (despite no longer being a valid EM algorithm) does give stable answers with qualitatively correct covariate effects. But because it zeroes out the posterior probabilities for all z_i within a neighborhood of the origin, it yields heavily biased estimates for β, and in our studies, it is less powerful than the method of Section 2.1.

2.3 Full Bayes

From a Bayesian perspective, the hierarchical model

together with priors for β and the unknown distribution π(θ), defines a joint posterior distribution over all model parameters. We use a Markov-chain Monte Carlo algorithm to sample from this posterior, drawing iteratively from three complete conditional distributions: for the mixing density π(θ); for the latent binary variables h_i that indicate whether z_i is signal or null; and for the regression coefficients β.

An important question is how to parameterize π(θ). In the no-covariates multiple-testing problem, there have been many proposals, including simple parametric families (Scott and Berger, 2006; Polson and Scott, 2012) and nonparametric priors based on mixtures of Dirichlet processes (Do et al., 2005). In principle, any of these methods could be used. In our analyses, we model π(θ) as a K-component mixture of Gaussians with unknown means, variances, and weights. We choose K via a preliminary run of the EM algorithm for deconvolution mixture models, picking the K that minimizes the Akaike information criterion (AIC). (In simulations, we found that AIC was slightly better than BIC at recovering K for the deconvolution problem, as distinct from the ordinary density-estimation problem.)

The model’s chief computational difficulty is the analytically inconvenient form of the conditional posterior distribution for β. The two major issues here are that the response h_i depends non-linearly on the parameters, and that there is no natural conjugate prior to facilitate posterior computation. These issues are present in all Bayesian analyses of the logit model, and have typically been handled using the Metropolis–Hastings algorithm. A third issue, particular to our setting, is that the binary event h_i is a latent variable, and not actually observed.

We proceed by exploiting the Pólya-Gamma data-augmentation scheme for binomial models recently proposed by Polson et al. (2013). Let x_i be the vector of covariates for test statistic i, including an intercept term. The complete conditional for β depends only h = {h_i: i = 1, …, n}, and may be written as

where p(β) is the prior. By introducing latent variables ω_i ~ PG(1, 0), each having a standard Pólya-Gamma distribution, we may re-express each term in the above product as the marginal of a more convenient joint density:

where κ_i = h_i − 1/2 and p(ω) is the density of a PG(1,0) variate. Assuming a normal prior β ~ N(c, D), it can be shown that β has a conditionally Gaussian distribution, given the diagonal matrix Ω = diag(ω₁, …, ω_n). Moreover, the conditional for each ω_i, given β, is also in the Pólya-Gamma family, and may be efficiently simulated.

Together with standard results on mixture models, the Pólya-Gamma scheme leads to a simple, efficient Gibbs sampler for the fully Bayesian FDR regression model. Further details of the Bayesian method can be found in Appendix B, including the priors we use, the conditionals needed for sampling, and the default settings implemented in FDRreg.

On both simulated and real data sets, we have observed that the empirical-Bayes and fully Bayes approaches give very similar answers for the local false-discovery rates, and thus reach similar conclusions about which cases are significant. The advantage of the Bayesian approach is that it provides a natural way to quantify uncertainty about the regression function s(x) and π(θ) jointly. This is counterbalanced by the additional computational complexity of the fully Bayesian method.

2.4 Using an empirical null

The FDR regression model (5) assumes that μ and σ² are both known, or can be derived from the distributional theory of the test statistic in question. As Efron (2004) observes, however, many data sets are poorly described by this “theoretical null,” and an “empirical null” must be estimated instead. This is a common situation in high-dimensional screening problems, where correlation among the test statistics, along with many other factors, can invalidate the theoretical null. Efron (2004) proposes two methods for estimating μ and σ² from the data: (1) maximum likelihood, and (2) central matching, whereby a quadratic function is fit to the log density of some central fraction (e.g. a third) of the data.

The estimation of an empirical null hypothesis, whether by maximum likelihood or central matching, can be incorporated into the empirical-Bayes method of Section 2.1 as a simple pre-processing step. This approach is used later in our analysis of the neural synchrony data, and is offered as an option in FDRreg. For reasons of identifiability, it is challenging to incorporate an empirical null in a fully Bayesian manner, especially without a strong prior about the null hypothesis. We refer the reader to Martin and Tokdar (2012) for a detailed discussion of this issue, as well as an alternate proposal for estimating an empirical null. When we use an empirical null in the context of the fully Bayes model described in Section 2.3, μ and σ² are always pre-computed in empirical-Bayes fashion using Efron’s method. This may be slightly less efficient than performing a full Bayesian analysis, but we believe the improved stability is worth the trade-off.

4.1 Background

The ability to record dozens or hundreds of neural spike trains simultaneously has posed many new challenges for data analysis in neuroscience (Brown et al., 2004; Buzsaki, 2004; Aertsen, 2010; Stevenson and Kording, 2011). Among these, the problem of identifying neural interactions has, since the advent of multi-unit recording, been recognized as centrally important (Perkel et al., 1967). Neural interactions may occur on sub-behavioral timescales, where two neurons may fire repeatedly within a few milliseconds of each other. It has been proposed that such fine-timescale synchrony is crucial for binding visual objects (see Gray, 1999; Shadlen and Movshon, 1999, for opposing views), enhancing the strength of communication between groups of neurons (Niebur et al., 2002), and coordinating the activity of multiple brain regions (Fries, 2009; Saalmann and Kastner, 2011). It has also been argued that the disruption of synchrony may play a role in cognitive dysfunctions and brain disorders (Uhlhaas et al., 2009).

Furthermore, there is growing recognition that synchrony and other forms of correlated spiking have an impact on population coding (Averbeck et al., 2006) and decoding (Graf et al., 2011). The proposed roles of neural synchrony in numerous computational processes and models of coding and decoding, combined with the knowledge that the amount of synchrony can depend on stimulus identity and strength (Kohn and Smith, 2005) as well as the neuronal separation (Smith and Kohn, 2008), make it particularly important that we have effective tools for measuring synchrony and determining how it varies under different experimental paradigms.

Rigorous statistical detection of synchrony in the activity of two neurons requires formulation of a statistical model that specifies the stochastic behavior of the two neural spike trains under the assumption that they are conditionally independent, given some suitable statistics or covariates (Harrison et al., 2013). When n neural spike trains are recorded there $N=\left(\begin{array}{c}\hfill n\hfill \\ \hfill 2\hfill \end{array}\right)$ null hypotheses to be tested, which raises the problem of multiplicity. In the face of this difficulty, a popular way to proceed has been to control the false-discovery rate using the Benjamini-Hochberg procedure Benjamini and Hochberg (1995), combining all N test statistic into a single analysis. Yet this omnibus approach ignores potentially useful information about the spatial and functional relationships among individual neuron pairs. We therefore use false discovery rate regression to incorporate these covariates into an investigation of synchrony in the primary visual cortex (V1). Specifically, we analyzed data from V1 neurons recorded from an anesthetized monkey in response to visual stimuli consisting of drifting sinusoidal gratings, i.e., light whose luminance was governed by a sine wave that moved along an axis having a particular orientation. Details of the experiment and recording technique may be found in Kelly et al. (2007). Drifting gratings are known to drive many V1 neurons to fire at a rate that depends on orientation. Thus, many V1 neurons will have a maximal firing rate for a grating in a particular orientation (which is the same as an orientation rotated by 180 degrees) and a minimal firing rate when the grating orientation is rotated by 90 degrees. For a given neuron, a plot of average firing rate against angle of orientation produces what is known as the neuron’s “tuning curve.” Spike trains from 128 neurons were recorded in response to gratings in 98 equally-spaced orientations, across 125 complete replications of the experiment (125 trials). Here analyze data from the first 3 seconds of each 30-second trial. The 128 neurons generated 8,128 pairs, and thus 8,128 tests of synchrony. We applied the model in Equation (4) to examine the way the probability of synchrony c(x) for a pair of neurons depends on two covariates: the distance between the neurons and the correlation of their tuning curves (i.e., the Pearson correlation between the two vectors of length 98 that contain average firing rate as a function of orientation). The idea is that when neurons are close together, or have similar tuning curves, they may be more likely to share inputs and thus more likely to produce synchronous spikes, compared to the number predicted under conditional independence. Our analysis, reported below, substantiates the observation of Smith and Kohn (2008) that the probability of fine time-scale synchrony for pairs of V1 neurons tends to decrease with the distance between the two neurons and increase with the magnitude of tuning-curve correlation.

4.2 Data pre-processing

Our analysis takes advantage of a recently-developed technique for measuring synchrony across binned spike trains (where the time bins are small, such as 5 milliseconds). For a pair of neurons labeled 1 and 2, we calculate

where $D_t^{\left(1\right)}$ and $D_t^{\left(2\right)}$ refer to relevant conditioning information for the firing activity of neurons 1 and 2 at time t, and the sum is over all time bins across all experimental trials (Kelly and Kass, 2012). The denominator of (10) is an estimate of the number of joint spikes that would be expected, by chance, if $D_t^{\left(j\right)}$ characterized the spiking activity of neuron j (with j = 1, 2) and, apart from these background effects, the neurons were independent. When $\widehat{\zeta }\approx 1$, or $\widehat{\zeta }\approx 0$, the conclusion would be that the number of observed synchronous spikes is consistent with the prediction of synchronous spiking under independence, given the background effects D^(j). Note that this conditioning information $D_t^{\left(j\right)}$ is intended to capture effects—including tuning curve information—on each neuron separately, whereas the covariates that enter the FDRR model (4) operate pairwise. Thus, for example, two independent neurons having similar tuning curves would both be driven to fire more rapidly by a grating stimulus in a particular orientation, and would therefore be likely to produce more synchronous spikes by chance than a pair of independent neurons with dissimilar tuning curves. As another example of conditioning information, under anesthesia there are pronounced periods during which most recorded neurons increase their firing rate (Brown et al., 2010). These waves of increased network activity have much lower frequency than many other physiological wave-like neural behaviors, and are called “slow waves.” One would expect slow-wave activity to account for considerable synchrony, even if, conditionally on the slow-wave activity, a pair of neurons were independent. The statistic $\widehat{\zeta }$ in formula (10) is a maximum-likelihood estimator in the continuous-time framework discussed by Kass et al. (2011), specifically their equation (22). The purpose of that framework, and of (10), is to describe the way synchrony might depend on background information. For example, Kass et al. (2011) contrasted results from two pairs of V1 neurons. Both pairs exhibited highly statistically significantly enhanced synchrony above that predicted by stimulus effects (tuning curves) alone. However, the two pairs were very different with respect to the relationship of synchrony to slow-wave network activity. In one pair, when background information characterizing the presence of slow-wave network activity was used in (10), the enhanced synchrony vanished, with $\log {\widehat{\zeta }}_1=.06\pm .15$. In the other pair, it persisted with $\log {\widehat{\zeta }}_2=.82\pm .23$, indicating the number of synchronous spikes was more than double the number predicted by slow-wave network activity together with trial-averaged firing rate.

The set of 8,128 $\widehat{\zeta }$ coefficients analyzed here, together with their standard errors, were created with a model that differed in two ways from that used by Kass et al. (2011). First, to better capture slow-wave network effects, in place of a linear model based on single count variable (for neurons i and j, Kass et al. used the total number of spikes within the past 50 ms time among all neurons other than neurons i and j) a general nonparametric function of the count was fitted, using splines. Second, a nonparametric function capturing spike history effects was used. This allows non-Poisson variability, which is important in many contexts (Kass et al., 2014), and exploratory analysis indicated that it is consequential for this data set as well.