Adjusted regularization of cortical covariance

Abstract

It is now common to record dozens to hundreds or more neurons simultaneously, and to ask how the network activity changes across experimental conditions. A natural framework for addressing questions of functional connectivity is to apply Gaussian graphical modeling to neural data, where each edge in the graph corresponds to a non-zero partial correlation between neurons. Because the number of possible edges is large, one strategy for estimating the graph has been to apply methods that aim to identify large sparse effects using an $L_{1}$ penalty. However, the partial correlations found in neural spike count data are neither large nor sparse, so techniques that perform well in sparse settings will typically perform poorly in the context of neural spike count data. Fortunately, the correlated firing for any pair of cortical neurons depends strongly on both their distance apart and the features for which they are tuned. We introduce a method that takes advantage of these known, strong effects by allowing the penalty to depend on them: thus, for example, the connection between pairs of neurons that are close together will be penalized less than pairs that are far apart. We show through simulations that this physiologically-motivated procedure performs substantially better than off-the-shelf generic tools, and we illustrate by applying the methodology to populations of neurons recorded with multielectrode arrays implanted in macaque visual cortex areas V1 and V4.

Appendix

1.1 SAGlasso algorithm

There are several ways to build the weight matrix Q of SAGlasso. We used Gamma regression, as described in Algorithm 1, which can be implemented efficiently with standard statistical software, e.g. R packages glm (Dobson and Barnett 2008; Hastie and Pregibon 1992; McCullagh and Nelder 1989; Venables and Ripley 2002), mgcv ( Wood 2011), and gam (Hastie and Tibshirani 1990). Note that in Eq. (3), Q is typically estimated by the square rooted absolute entries of the inverse sample covariance matrix. In SAGlasso, we observed a slightly better performance without applying any transformation.

1.2 GAR algorithms

GAR Algorithms 2–4 are derived in Section 2.2, and implemented in our R package “GARggm” available in ModelDB.

In Algorithm 2, $U\sim \text {InvGaussian}(a,b)$ has p.d.f. $p(u) = \left (\frac {b}{2\pi u^{3}}\right )^{1/2}\exp \left \{-b(u-a)^{2}/(2a^{2}u)\right \}$. Moreover, given a matrix M, $M_{ij}$ is the i-th row and j-th column entry of M; $M_{-ij}$ is the j-th column of M without the i-th entry; $M_{i-j}$ is the i-th row of M without the j-th entry; and $M_{-i-j}$ is the submatrix obtained by removing the i-th row and the j-th column from M. Algorithms 3 and 4 both produce posterior samples of ${\Omega }$ whose average approximates the posterior mean of ${\Omega }$. The posterior mode of ${\Omega }$ can be obtained by solving Eq. (1) with λ∥Ω∥ replaced by $\Vert \hat {\Lambda } \odot {\Omega }\Vert _{1}$, where $\hat {\Lambda }$ is the estimated penalty matrix from either Algorithm 3 or 4, and $\odot$ denotes the entry-wise matrix multiplication. This optimization can be performed using R functions such as glasso (package glasso, Friedman et al. (2008)) with argument rho set equal to $2\hat {\Lambda } /n$; see also the R package QUIC, Hsieh et al. (2011). We solve the SPL problem in Eq. (5) by the EM algorithm of Yuan (2012) involving Glasso, and we impose the GAR penalty matrix $\hat {\Lambda }$ on S in the Maximization step to obtain the GAR-SPL estimate. For $d\sim 100$, we suggest to run the Gibbs samplers for at least $B = 2000$ iterations, including a burn-in period of 300 iterations. The Gamma regression in step 2b of Algorithm 4 can be implemented either parametrically or nonparametrically by using standard statistical software e.g. R packages glm (Dobson and Barnett 2008; Hastie and Pregibon 1992; McCullagh and Nelder 1989; Venables and Ripley 2002) and mgcv (Wood 2011); in the data analyses we used splines (Kass et al. 2014).

1.3 Computational efficiency of estimators

Table 1 contains the computation times of the graph estimators considered for $d = 50, 100$ and $n = 200, 500$, using the programming language R, CPU Quad-core 2.6 GHz Intel Core i7, and RAM 16 GB 2133 MHz DDR4. These times could be improved substantially by using a lower level language such as C++. Glasso, AGlasso, SPL, and SAGlasso are fitted with tuning parameter optimization based on ten-fold cross-validation involving 500 random splits over a fine grid of 20 values of the tuning parameter about its optimal value. GAR full Bayes (Algorithm 3; $K=\lfloor \sqrt {d} \rceil$) involved $B = 2000$ iterations, where the Gibbs sampler converged after about 300 iterations. GAR empirical Bayes (Algorithm 4; splines with 3 knots) involved 30 EM iterations, each including 500 iterations of Gibbs sampler for the E-step; the efficiency of this method may be improved by replacing the Gibbs sampler with some alternate faster approximation of Eq. (20).

Table 1 Computational time in seconds with 95% confidence intervals for $d = 50, 100$ and $n = 200, 500$