# Understanding the population structure correction regression

The Tien Mai(1) and Pierre Alquier(2)

(1) Oslo Centre for Biostatistics and Epidemiology, Department of
Biostatistics,
University of Oslo, Norway.
(2) RIKEN Center for Advanced Intelligence Project, Tokyo, Japan.

(email<EMAIL_ADDRESS>
###### Abstract

Although genome-wide association studies (GWAS) on complex traits have
achieved great successes, the current leading GWAS approaches simply perform
to test each genotype-phenotype association separately for each genetic
variant. Curiously, the statistical properties for using these approaches is
not known when a joint model for the whole genetic variants is considered.
Here we advance in GWAS in understanding the statistical properties of the
”population structure correction” (PSC) approach, a standard univariate
approach in GWAS. We further propose and analyse a correction to the PSC
approach, termed as ”corrected population correction” (CPC). Together with the
theoretical results, numerical simulations show that CPC is always comparable
or better than PSC, with a dramatic improvement in some special cases.

Keywords: GWAS, population structure correction, linear regression, bias,
variance.

##  1 Introduction

In high dimensional data analysis where the number of covariates $p$ is larger
than the number of samples $n$, penalized regression approaches, such as
Lasso, are one of the most popular approach [Hastie et al., 2009, Bühlmann and
Van De Geer, 2011, Giraud, 2014, Efron and Hastie, 2016]. However,
interpreting the results of the Lasso in terms of hypothesis testing or
uncertainty quantification is difficult.

Motivated by genome-wide association studies (GWAS), we focus on the following
question: given a response vector $y$ of $n$ samples and a matrix of (genetic)
covariates $X_{n\times p}$ formed by $p$ (genetic) covariates of $n$ samples,
we want to determine which covariates associate with the response. Although
the variable selection problem is a classical problem in statistics, in this
context it is still a big challenge as the number of covariates is huge
compared to the sample size, which prohibits the use of the classical methods.
Moreover, in many practical situation, the genomic data are huge and can not
even be examined on a personal laptop. For example, in human genomics, the
number of covariates (SNPs, single-nucleotide polymorphism) as well as the
number of samples are often at hundreds thousands [Bycroft et al., 2018]; or
it can be at order tens of millions covariates when using k-mers (an
alignment-free biomarker type) as in bacterial genomics [Lees et al., 2016]
with thousands of samples.

Besides the computational reason, most of the theoretical result on the Lasso
are on $\ell_{2}$ estimation of the parameter and not on variable selection.
These two objectives are known to be incompatible in general, see [Yang, 2005,
Leeb and Pötscher, 2008] (and [Zhao and Yu, 2006] illustrated this in the case
of the Lasso). Some results for variable selection were derived for
thresholded versions of the Lasso, for example in [Lounici, 2008], but are
valid under strong assumptions that are usually not satisfied in GWAS. Some
procedure leading to significance tests and confidence intervals are proposed
e.g in [Javanmard and Montanari, 2014, Van de Geer et al., 2014, Zhang and
Zhang, 2014] but they are not easy to handle and costly to compute when using
penalized regression such as Lasso.

This lead to the widely use of a univariate model for testing the association
of a trait and a covariate (say $X_{\cdot 1}$):

$\displaystyle y=\alpha_{1}X_{\cdot 1}+e$

to estimate $\alpha_{1}$ and test its significant. However, the omitted
variables have an effect, which we will model by multivariate and so the
fitted model should be:

$\displaystyle y=\beta_{1}X_{\cdot 1}+\sum_{j>1}\beta_{j}X_{\cdot
j}+\varepsilon.$

The effect of omitted variables (that is, the difference between $\alpha_{1}$
and $\beta_{1}$, depends strongly on the dependence between $X_{\cdot 1}$ and
the other covariances. In GWAS data, the covariates (often SNPs) are in some
dependent structures which is called linkage disequilibrium. This is due to
the population structure: many of SNPs have different frequencies in each
population [Price et al., 2006, Price et al., 2010].

If one uses a univariate regression and ignores the effect of the other
covariates, then they can be effectively modelled as part of the error as
$e=\sum_{j>1}\beta_{j}X_{\cdot j}+\varepsilon$. However the covariates are
correlated due to the population structure, this leads to a correlation
between the tested covariate, say $X_{\cdot 1}$, and the noise term together
with the noise of the samples. These correlations can cause in inflated type-1
error rates [Derks et al., 2017]. To handle this problem, the so-called
’population structure correction’ approach had been introduced and
successfully applied in practice, see for example [Price et al., 2006, Price
et al., 2010, Lippert et al., 2011] among others.

In principle, population structure correction is an alternative way to
implicitly model the other covariates that are not being tested at the time.
This can be done through the latent subspace of these variables. A natural way
is to use principal component analysis to extract some features that contain
most information of the other covariates $X_{\cdot-1}$ and use these features
as representatives added in the univariate regression of $X_{\cdot 1}$. In
this way, it can be seen as a dimension-reduction approach. Another way that
is also being the standard approach in GWAS is to use ’linear mixed model’
framework in which the covariates that are not directly being tested are
treated as random. However, several works had shown that inclusion of
$X_{\cdot 1}$ in calculating the principal components can lead to loss in
power [Listgarten et al., 2012, Yang et al., 2014]. This motivates and leads
to the popular usage of leave-one-chromosome-out method [Lippert et al., 2011,
Listgarten et al., 2012, Yang et al., 2014].

Although univariate regression approach with population structure correction
has become the state-of-the-art approach in GWAS, there are several numerical
works have showed that fitting a penalized multivariate regression can exceed
it, e.g [Wu et al., 2009, Visscher et al., 2017, Buzdugan et al., 2016,
Brzyski et al., 2017, Lees et al., 2020]. This can be explained as that using
population structure correction can be biased. Moreover, population structure
correction very much depends on the added latent features of untested
covariates.

In this paper, we study the statistical properties of the population structure
correction when assuming the true underlying model is a multivariate linear
regression. More specifically, we derive explicitly the bias and the variance
of the population structure correction method. Moreover, we also propose and
study a simple version of the leave-one-chromosome-out method, termed as
’corrected population correction’. We show theoretically and empirically that
’corrected population correction’ approach reduces the variance compare to the
population structure method.

The paper is organized as follow. In Section 2, the model formulation and
different methods are presented. The main results on the statistical
properties of population structure correction and corrected population
correction are given in Section 3. Some numerical simulations are conducted in
Section 4 and we conclude the paper in the final section.

##  2 Model and methods

###  2.1 Model

Given a response vector $y$ of $n$ samples and $p$ covariates $X_{\cdot j}$
with $n\ll p$, we assume that the response vector relates to the covariates by
the following linear model

$Y_{i}=\sum_{j=1}^{p}\beta_{j}X_{i,j}+\varepsilon_{i},i=1,\ldots,n$ (1)

or, summarized by $Y=X\beta+\varepsilon$ in matrix form with

$\left(\begin{array}[]{c}Y_{1}\\\ \vdots\\\
Y_{n}\end{array}\right)=\left(\begin{array}[]{c c c}X_{1,1}&\dots&X_{1,p}\\\
\vdots&\ddots&\vdots\\\
X_{n,1}&\dots&X_{n,p}\end{array}\right)\left(\begin{array}[]{c}\beta_{1}\\\
\vdots\\\
\beta_{p}\end{array}\right)+\left(\begin{array}[]{c}\varepsilon_{1}\\\
\vdots\\\ \varepsilon_{n}\end{array}\right).$

We assume that $\varepsilon_{i}\sim\mathcal{N}(0,\sigma^{2})$. With a
deterministic $X$, this leads to $\mathbb{E}(Y)=X\beta.$

The problem is to estimate the coefficient $\beta_{j}$. Up to a re-ordering of
the variables, say that the coefficient is $\beta_{1}$. For simplicity, there
is no intercept: we assume that the variables are already centered, and
normalized.

##### Notations:

For any matrix $A$, $A_{(-j)}$ denotes matrix $A$ without its $j$-th column,
and $A_{s:j}$ is the submatrix with only columns $s,s+1,\dots,j$. We use the
same convention for column vectors: $\beta_{s:j}$ means that we extract
entries from $s$ to $j$.

###  2.2 Population structure correction (PSC)

The idea is to perform a principal component analysis (PCA) on $X$ and then
use some principal components corresponding to the top leading eigenvalues. In
other words,

$X^{\top}X=\bar{W}^{\top}\left(\begin{array}[]{c c
c}\bar{\lambda}_{1}&\dots&0\\\ \vdots&\ddots&\vdots\\\
0&\dots&\bar{\lambda}_{p}\end{array}\right)\bar{W}\text{ where
}\bar{W}=\left(\begin{array}[]{c c c}\bar{W}_{1,1}&\dots&\bar{W}_{1,p}\\\
\vdots&\ddots&\vdots\\\ \bar{W}_{p,1}&\dots&\bar{W}_{p,p}\end{array}\right)$

and $\bar{\lambda}_{1}\geq\dots\geq\bar{\lambda}_{p}\geq 0$. The matrix
$X^{\top}X$ is also known as the ’kinship’ matrix in genomic research.

Here it will be more convenience to think of PCA as an SVD, that is

$X=\bar{U}\left(\begin{array}[]{c c c}\bar{\sigma}_{1}&\dots&0\\\
\vdots&\ddots&\vdots\\\
0&\dots&\bar{\sigma}_{r}\end{array}\right)\bar{V}^{\top}=\bar{U}\bar{\Sigma}\bar{V}^{\top}$

where $\bar{r}={\rm rank}(X)$ and $\bar{\sigma}_{1}>\dots>\bar{\sigma}_{r}>0$.
It is easy to see that $\bar{V}$ contains the $\bar{r}$ first columns of
$\bar{W}$, and so $\bar{V}=\bar{W}$ as soon as ${\rm rank}(\bar{X})=p$.

Then the idea is simply to estimate the model for some $k$,

$Y_{i}=\bar{\alpha}X_{1,i}+\sum_{j=1}^{k}\bar{\gamma}_{j}\bar{U}_{i,j}+\bar{e}^{(k)}_{i},$

or

$\displaystyle Y=X_{\cdot
1}\bar{\alpha}+\bar{U}_{(1:k)}\bar{\gamma}+\bar{e}^{(k)},$ (2)

hoping that $\bar{\alpha}$ is a good proxy for $\beta_{1}$ in model (1).

###  2.3 Corrected population correction (CPC)

We propose a modified procedure: first, perform a PCA on $X_{(-1)}$ ($X$
without its first column $X_{\cdot 1}$ is denoted by $X_{(-1)}$), in other
words:

$X_{(-1)}=U\left(\begin{array}[]{c c c}\sigma_{1}&\dots&0\\\
\vdots&\ddots&\vdots\\\ 0&\dots&\sigma_{r}\end{array}\right)V^{\top}=U\Sigma
V^{\top}$

where $r={\rm rank}(X_{(-1)})$ and $\sigma_{1}>\dots>\sigma_{r}>0$.

Similar to PSC, it is simply to estimate the model, for some $k$,

$Y_{i}=\alpha X_{1,i}+\sum_{j=1}^{k}\gamma_{j}U_{i,j}+e^{(k)}_{i},$

or

$\displaystyle Y=X_{\cdot 1}\,\alpha+U_{(1:k)}\gamma+e^{(k)},$ (3)

and hoping that $\alpha$ is a good proxy for $\beta_{1}$ in model (1).

We would like to note that this method is a simple version of the so-called
popular method ’leave-one-chromosome-out’ in GWAS [Yang et al., 2014].

###  2.4 Why does PSC need to be corrected?

In general, the CPC model is not “correct” in the sense that
$\mathbb{E}(Y)\neq X_{\cdot 1}\alpha+U_{1:k}\gamma$. In other words, in (3) we
don’t have $\mathbb{E}(e^{(k)})=0$, in general. However, assume that $k={\rm
rank}(X_{-1})$, then we have

$Y=X_{\cdot 1}\alpha+U\gamma+e$

and note that $X_{(-1)}=U\Sigma V^{\top}$ leads to
$X_{(-1)}V(V^{\top}V)^{-1}\Sigma^{-1}=U$. Thus, this model is equivalent to
(1) with $\alpha=\beta_{1}$, by identification:

$X_{\cdot
1}\alpha+X_{(-1)}(V(V^{\top}V)^{-1}\Sigma^{-1}\gamma)=X\beta=X_{\cdot
1}\beta_{1}+X_{(-1)}\left(\begin{array}[]{c}\beta_{2}\\\ \vdots\\\
\beta_{p}\end{array}\right)$

(and so $\varepsilon=e$ in this case). Therefore, for a well chosen $k$, the
model is actually exact. For this reason, we can reformulate the problem as:
with the true model

$Y=X_{\cdot 1}\alpha+U\gamma+\varepsilon,$ (4)

$\gamma\in\mathbb{R}^{r}$ where $r={\rm rank}(X_{(-1)})$, what is the effect
on $\alpha$ to estimate instead, for some $k$,

$Y=X_{\cdot 1}\alpha+U_{(1:k)}\gamma_{1:k}+e^{(k)}$ (5)

where we actually have $e^{(k)}=U_{(k+1):r}\gamma_{(k+1):r}+\varepsilon$. This
is simply a problem of omission of variables: what is the effect of the
omission of $U_{(k+1):r}$?

On the other hand, for $k={\rm rank}(X)$, we have

$Y=X\beta+\varepsilon=\bar{U}\bar{\Sigma}\bar{V}^{\top}\beta+\varepsilon=\bar{U}(\bar{\Sigma}\bar{V}^{\top}\beta)+\varepsilon=\bar{U}\bar{\gamma}+\varepsilon$

and so the PSC model

$\displaystyle Y=\bar{\alpha}X_{1}+\bar{U}\bar{\gamma}+\varepsilon$ (6)

is simply not identifiable (the variable $X_{\cdot 1}$ is twice in the model).
When $k<{\rm rank}(X)$, the model might be identifiable, but the fact that
$X_{\cdot 1}$ is in the first term, and “partly” in the second, will lead to a
greater bias than in CPC. This drawback has been figured out in the field of
genetic research [Yang et al., 2014]. For a formal statement see the analysis
below.

##  3 Statistical Analysis

In the following we explicitly derive the bias and the variance for each
considered methods above. These results bring insights on understanding how
the population structure correction is working practically. All technical
proofs are postponed to Section A.

###  3.1 Main theorems

We first provide some statistical properties for the CPC method.

######  Theorem 1.

Assume that model (1) or equivalently (4) holds. Then with CPC method we have

${\rm bias}(\hat{\alpha})=\frac{X_{\cdot
1}^{\top}U_{(k+1):r}\gamma_{(k+1):r}}{X_{\cdot 1}^{\top}X_{\cdot
1}-\left\|X_{\cdot 1}^{\top}U_{1:k}\right\|^{2}},$ $\displaystyle{\rm
Var}(\hat{\alpha})=\frac{\sigma^{2}}{X_{\cdot 1}^{\top}X_{\cdot
1}-\left\|X_{\cdot 1}^{\top}U_{1:k}\right\|^{2}}.$ (7)

It can be seen that $X_{\cdot 1}^{\top}U_{(k+1):r}\gamma_{(k+1):r}$ measures
the correlation between $X_{\cdot 1}$ and the other part of $X$ (as in the
true model (4)) which was not included in the wrong model (5). Obviously, if
the wrong model is actually not “too wrong” in the sense that
$\|U_{(k+1):r}\gamma_{(k+1):r}\|\simeq 0$ then the bias would be small. But
when this is not the case, the term is problematic only if $X_{\cdot 1}$ is
correlated with this quantity.

The denominator $X_{\cdot 1}^{\top}X_{\cdot 1}-\left\|X_{\cdot
1}^{\top}U_{1:k}\right\|^{2}$ is just an identifiability term: if $X_{\cdot
1}$ is too correlated with the other variables used in model (5), then the
variance of $\hat{\alpha}$ will increase (as usual) but also the bias due to
misspecification.

######  Remark 1.

As $\gamma$ is unknown in practice, but the $c_{j}:=X_{\cdot 1}^{\top}U_{(j)}$
are observed. So we can give a result under an assumption that depends only on
$\gamma$. For example if we assume that $\|\gamma\|_{1}\leq B$ (as in the
Lasso) then

$|{\rm
bias}(\hat{\alpha})|\leq\frac{\|\gamma\|_{1}\sum_{j=k+1}^{r}|c_{j}|}{X_{\cdot
1}^{\top}X_{\cdot
1}-\sum_{j=1}^{k}c_{j}^{2}}\leq\frac{B\sum_{j=k+1}^{r}|c_{j}|}{X_{\cdot
1}^{\top}X_{\cdot 1}-\sum_{j=1}^{k}c_{j}^{2}}.$

Statistical properties of the PSC method are given in the following theorem.

######  Theorem 2.

Assume that model (1) holds. For the model (6), with PSC method, we have

${\rm bias}(\hat{\bar{\alpha}})=\frac{X_{\cdot
1}^{\top}(\bar{U}_{(k+1):\bar{r}}\bar{\gamma}_{(k+1):\bar{r}}-X_{\cdot
1}\beta_{1})-\beta_{1}\left\|X_{\cdot
1}^{\top}\bar{U}_{1:k}\right\|^{2}}{X_{\cdot 1}^{\top}X_{\cdot
1}-\left\|X_{\cdot 1}^{\top}\bar{U}_{1:k}\right\|^{2}};$ $\displaystyle{\rm
Var}(\hat{\bar{\alpha}})=\frac{\sigma^{2}}{X_{\cdot 1}^{\top}X_{\cdot
1}-\left\|X_{\cdot 1}^{\top}\bar{U}_{1:k}\right\|^{2}}.$ (8)

In the following theorem, we derive the relationship between the variances of
these two methods.

######  Theorem 3.

The corrected population correction method reduces the variance of the
original population structure correction, i.e ${\rm Var}(\hat{\alpha})\leq{\rm
Var}(\hat{\bar{\alpha}}).$

######  Remark 2.

From Theorem 3, it states that the corrected population correction (CPC)
always returns estimate with smaller variance comparing to the structured
population correction (SPC). From simulations, we conjecture that the biasness
of CPC method is also smaller than those from SPC method, however this is not
easy to show from our analysis.

###  3.2 Reliable implication check

Assuming that $X$ is normalized, that is $X_{\cdot 1}^{\top}X_{\cdot 1}=1$. As
we have that $|X_{\cdot
1}^{\top}U_{(k+1):r}\gamma_{(k+1):r}|\leq\sum_{s=k+1}^{r}|X_{\cdot
1}^{\top}U_{s}\gamma_{s}|,$ and using Cauchy-Schwarz inequality yields
$|X_{\cdot 1}^{\top}U_{s}|^{2}\leq\|X_{\cdot 1}\|^{2}\|U_{s}\|^{2}\leq 1$ and
we obtain $|X_{\cdot
1}^{\top}U_{(k+1):r}\gamma_{(k+1):r}|\leq\sum_{s=k+1}^{r}|\gamma_{s}|.$ Thus,
we have

$\left|{\rm
bias}(\hat{\alpha})\right|\leq\frac{\left\|\gamma_{(k+1):r}\right\|_{1}}{\left|1-\left\|X_{\cdot
1}^{\top}U_{1:k}\right\|^{2}\right|}:=\frac{N}{D},\quad{\rm
Var}(\hat{\alpha})\leq\frac{\sigma^{2}}{\left|1-\left\|X_{\cdot
1}^{\top}U_{1:k}\right\|^{2}\right|}:=\frac{\sigma^{2}}{D}.$

Now, CLT yields
$\hat{\alpha}\sim\mathcal{N}(\mathbb{E}\left(\hat{\alpha}),{\rm
Var}(\hat{\alpha})\right)$. We want to test the null hypothesis that

$H_{0}:\alpha=0.$

Under this null hypothesis, we have
$|\mathbb{E}(\hat{\alpha})|\leq\frac{N}{D}$ and $(1-a)100\%$ confidence
interval is
$\mathbb{E}(\hat{\alpha})\in\left[\hat{\alpha}-z_{\frac{a}{2}}\sqrt{{\rm
Var}(\hat{\alpha})},\hat{\alpha}+z_{\frac{a}{2}}\sqrt{{\rm
Var}(\hat{\alpha})}\right]$ or
$\hat{\alpha}\in\left[-\mathbb{E}(\hat{\alpha})-z_{\frac{a}{2}}\sqrt{{\rm
Var}(\hat{\alpha})},\mathbb{E}(\hat{\alpha})+z_{\frac{a}{2}}\sqrt{{\rm
Var}(\hat{\alpha})}\right]$. Thus we obtain

$\displaystyle\hat{\alpha}\in\left[-\frac{N}{D}-z_{\frac{a}{2}}\sqrt{\frac{\sigma^{2}}{D}},\frac{N}{D}+z_{\frac{a}{2}}\sqrt{\frac{\sigma^{2}}{D}}\right],$

where $z_{\frac{a}{2}}=\Phi^{-1}(1-\frac{a}{2})$ and $\Phi^{-1}(\cdot)$ is the
normal cumulative distribution function.

The above analysis lead to the following tests for the null hypothesis:

$\displaystyle{\rm
If}\quad\hat{\alpha}\notin\left[-\frac{N}{D}-z_{\frac{a}{2}}\sqrt{\frac{\sigma^{2}}{D}},\frac{N}{D}+z_{\frac{a}{2}}\sqrt{\frac{\sigma^{2}}{D}}\right]\quad{\rm
then\,\,reject}\,\,H_{0}.$

If the noise variance $\sigma$ is not known, one can use a consistent estimate
$\hat{\sigma}$, for example as in [Dicker, 2014, Sun and Zhang, 2012, Janson
et al., 2017] and the confidence interval would become

$\displaystyle\hat{\alpha}\in\left[-\frac{N}{D}-t_{(n-1,\frac{a}{2})}\sqrt{\frac{\hat{\sigma}^{2}}{D}},\frac{N}{D}+t_{(n-1,\frac{a}{2})}\sqrt{\frac{\hat{\sigma}^{2}}{D}}\right],$

where the normal cumulative distribution function is replaced by the student
distribution.

##  4 Numerical simulations

### Setup

In this section, we investigate basic properties of the PSC and CPC methods
studied above. We fix $p=100,n=1000$ for low dimension setting and
$p=1000,n=600$ for high dimension setting. The noise variance is fixed at
$\sigma^{2}=1$.

We generate the parameter $\beta\in\mathbb{R}^{p}$ such that its first
component $\beta_{1}$ is fixed to 1, and other non-zero components was sampled
uniformly at random from $\\{\pm 1\\}$. The sparsity of $\beta$ will be
changed in each setting corresponding to $\|\beta\|_{0}=20,100$. The response
$Y$ is simulated as in linear model (1).

For each setting, we simulated 100 independent datasets and report the average
results together with their standard deviations. The number of principal
components $k$ added in models (3) and (2) are varied from 1 to 30.

### Example: worst case scenario for PSC

Here we show cases that PSC does not work well while CPC performs superior
results.

We consider the structured $X$ such that its first two columns $X_{\cdot 1}$
and $X_{\cdot 2}$ are corresponding to its first two leading principal
components. A brief summary of the data can be found in the Figure 1.

Figure 1: Summary the structured X: When $X_{\cdot-1}$ is removed, a principal
component is also removed.

In this case, it is clear to see that PSC can actually be very biased whereas
CPC is very stable and accurate, see Figure 2. This example demonstrates that
including $X_{1}$ (the covariate being tested) in the calculation of the
principal components can be very harmful.

Figure 2: Structured $X$. Estimates of $\beta_{1}$ with different number of
Principal components (PCs)

### Behaviour in other cases

We further consider the following settings for the design matrix:

  * •

Independent $X$: In this setting, $X_{ij}\sim\mathcal{N}(0,1)$.

  * •

Dependent $X$: We consider $X_{i\cdot}\sim\mathcal{N}(0,\Sigma)$, where
$\Sigma_{ij}=0.5^{|i-j|}$.

  * •

Binary $X$: the $X_{ij}$ is simulated from the set $\\{\pm 1\\}$ with equal
probability.

Results from simulations, Figures 6,6 and 6 confirm our theoretical results
above. In general, the CPC method performs similarly to PSC method. However,
CPC return the results with less variation than the PSC method. Moreover, the
PSC method is very much depending on $k$, the number of principal components
added in the model.

### Real data assessment in a wheat GWAS data

We apply two methods to a real wheat GWAS data which is available in the R
package ’BGLR’ [Pérez and de Los Campos, 2014]. The data consists of 599 wheat
lines: lines (responses) were evaluated for grain yield and each line has been
genotyped with 1279 markers.

We run CPC and PSC across 1279 covariates and report the absolute errors

$|\hat{\beta}_{j}^{CPS}-\hat{\beta}_{j}^{PSC}|$

and the relative errors

$\frac{|\hat{\beta}_{j}^{CPC}-\hat{\beta}_{j}^{PSC}|}{|\hat{\beta}_{j}^{PSC}|}.$

These results are given in Figure 3.

Regarding the histogram in Figure 3, the conclusion is clear: for most
coefficients, PSC and CPC lead to similar estimation, but for some of them,
the deviation is extremely high. There are in total 55 covariates such that
their relative errors are greater than 0.5 (and there are 33 covariates such
that their relative errors are greater than 1). Therefore, including the
covariate being tested in the calculation of the principal components could
create a huge difference.

Figure 3: Histogram of the absolute errors and relative errors of 1279
covariates in wheat data with 10 principal components.

##  5 Closing Discussion

In this paper, we have dicussed the statistical properties of the widely used
method, structure population correction method, in genome-wide association
studies. We have also proposed and studied a simple version of the ’leave-one-
chromosome-out’ in GWAS, termed as Corrected population correction method. Our
theoretical analysis and simulations show that the structure population
correction method (although efficient computationally) should be used with
more careful as it comes with higher variance due to model-mispecification.
The corrected population correction method, which requires higher
computational cost, returns better results as it avoids model-mispecification.

## Acknowledgments

T.T.M would like to thank Jukka Corander and John A Lees for useful discussion
on GWAS. The research of T.T.M was supported by the European Research Council
(SCARABEE project) no. 742158.

## Conflict of interest

The authors declare no potential conflict of interests.

## Availability of data and materials

The R codes and data used in the numerical experiments are available at:
https://github.com/tienmt/understand_SPC .

Figure 4: Independent $X$. Estimates of $\beta_{1}$ with different number of
Principal components (PCs)

Figure 5: Dependent $X$. Estimates of $\beta_{1}$ with different number of
Principal components (PCs)

Figure 6: Binary $X$. Estimates of $\beta_{1}$ with different number of
Principal components (PCs)

##  Appendix A Proofs

###  A.1 Proof of Theorem 1 for CPC

#### Analysis of the bias

In (5), let’s define the matrix $Z=\left(X_{\cdot 1}|U_{1:k}\right)$, put
$\zeta=(\alpha|\gamma_{i:k}^{\top})^{\top}$, then it becomes:

$Y=Z\zeta+e^{(k)}$

and the least square estimator is given by

$\displaystyle\hat{\zeta}=(Z^{\top}Z)^{-1}Z^{\top}Y=(Z^{\top}Z)^{-1}Z^{\top}[Z\zeta+e^{(k)}]=\zeta+(Z^{\top}Z)^{-1}Z^{\top}e^{(k)}.$
(9)

So the bias of this estimator is simply

${\rm
bias}(\hat{\zeta})=\mathbb{E}\left(\hat{\zeta}\right)-\zeta=(Z^{\top}Z)^{-1}Z^{\top}\mathbb{E}(e^{(k)}),$

that is

$\displaystyle{\rm
bias}(\hat{\zeta})=(Z^{\top}Z)^{-1}Z^{\top}U_{(k+1):r}\gamma_{(k+1):r}.$ (10)

Let us now make this more explicit. First,

$Z^{\top}U_{(k+1):r}\gamma_{(k+1):r}=\left(\begin{array}[]{c}X_{\cdot
1}^{\top}U_{(k+1):r}\gamma_{(k+1):r}\\\ 0\\\ \vdots\\\ 0\end{array}\right),$

and

$Z^{T}Z=\left(\begin{array}[]{c c c c}X_{\cdot 1}^{\top}X_{\cdot 1}&X_{\cdot
1}^{\top}U_{(1)}&\dots&X_{\cdot 1}^{\top}U_{(k)}\\\ X_{\cdot
1}^{\top}U_{(1)}&1&&\\\ \vdots&&\ddots&\vdots\\\ X_{\cdot
1}^{\top}U_{(k)}&0&\dots&1\end{array}\right).$

Now we know that the bias of $\hat{\zeta}$ in (10) satisfies

$(Z^{T}Z){\rm bias}(\hat{\zeta})=Z^{\top}U_{(k+1):r}\gamma_{(k+1):r},$

and it can be written explicitly as

$\displaystyle\left(\begin{array}[]{c c c c}X_{\cdot 1}^{\top}X_{\cdot
1}&X_{\cdot 1}^{\top}U_{(1)}&\dots&X_{\cdot 1}^{\top}U_{(k)}\\\ X_{\cdot
1}^{\top}U_{(1)}&1&\dots&0\\\ \vdots&&\ddots&\vdots\\\ X_{\cdot
1}^{\top}U_{(k)}&0&\dots&1\end{array}\right)\left(\begin{array}[]{c}{\rm
bias}(\hat{\alpha})\\\ {\rm bias}(\hat{\gamma}_{1})\\\ \vdots\\\ {\rm
bias}(\hat{\gamma}_{k})\end{array}\right)=\left(\begin{array}[]{c}X_{\cdot
1}^{\top}U_{(k+1):r}\gamma_{(k+1):r}\\\ 0\\\ \vdots\\\ 0\end{array}\right).$

The generic equation in the second part of the system is

$X_{\cdot 1}^{\top}U_{(j)}{\rm bias}(\hat{\alpha})+{\rm
bias}(\hat{\gamma}_{j})=0$

yielding

${\rm bias}(\hat{\gamma}_{j})=-X_{\cdot 1}^{\top}U_{(j)}{\rm
bias}(\hat{\alpha}).$

Plugging this into the first equation:

$\displaystyle X_{\cdot 1}^{\top}X_{\cdot 1}{\rm
bias}(\hat{\alpha})+\sum_{j=1}^{k}X_{\cdot 1}^{\top}U_{(j)}{\rm
bias}(\hat{\gamma}_{j})=X_{\cdot 1}^{\top}U_{(k+1):r}\gamma_{(k+1):r}$

gives

$\left[X_{\cdot 1}^{\top}X_{\cdot 1}-\sum_{j=1}^{k}\left(X_{\cdot
1}^{\top}U_{(j)}\right)^{2}\right]{\rm bias}(\hat{\alpha})=X_{\cdot
1}^{\top}U_{(k+1):r}\gamma_{(k+1):r}.$

Thus, we obtain

$\displaystyle{\rm bias}(\hat{\alpha})=\frac{X_{\cdot
1}^{\top}U_{(k+1):r}\gamma_{(k+1):r}}{X_{\cdot 1}^{\top}X_{\cdot
1}-\sum_{j=1}^{k}\left(X_{\cdot 1}^{\top}U_{(j)}\right)^{2}}=\frac{X_{\cdot
1}^{\top}U_{(k+1):r}\gamma_{(k+1):r}}{X_{\cdot 1}^{\top}X_{\cdot
1}-\left\|X_{\cdot 1}^{\top}U_{1:k}\right\|^{2}}.$

#### Variance analysis of CPC

From (9), assuming that ${\rm Var}(\epsilon)=\sigma^{2}$, we have

$\displaystyle{\rm Var}(\hat{\zeta})=(Z^{\top}Z)^{-1}Z^{\top}{\rm
Var}(u^{k})Z(Z^{\top}Z)^{-1}=\sigma^{2}(Z^{\top}Z)^{-1},$

or

$\displaystyle Z^{\top}Z{\rm Var}(\hat{\zeta})=\sigma^{2}\mathbf{I}.$

As we are only interested in estimating the variance of $\hat{\alpha}$, from
the above formula we obtain

$\displaystyle\begin{cases}(X_{\cdot 1}^{\top}X_{\cdot 1}){\rm
Var}(\hat{\alpha})+\sum_{j=1}^{k}Cov(\hat{\alpha},\hat{\gamma_{j}})X_{\cdot
1}^{\top}U_{(j)}=\sigma^{2},\\\ (X_{\cdot 1}^{\top}U_{(j)}){\rm
Var}(\hat{\alpha})+Cov(\hat{\alpha},\hat{\gamma_{j}})=0.\end{cases}$

Substituting the second equation into the first one to get

$\displaystyle{\rm Var}(\hat{\alpha})=\frac{\sigma^{2}}{X_{\cdot
1}^{\top}X_{\cdot 1}-\sum_{j=1}^{k}(X_{\cdot 1}^{\top}U_{(j)})^{2}}.$

###  A.2 Proof of Theorem 2 for PSC

#### Analysis of the bias

In (2), let’s define the matrix $\bar{Z}=\left(X_{\cdot
1}|\bar{U}_{1:k}\right)$, put
$\bar{\zeta}=(\bar{\alpha}|\bar{\gamma}_{i:k}^{\top})^{\top}$, and use the
least square estimator:

$\displaystyle\hat{\bar{\zeta}}=(\bar{Z}^{\top}\bar{Z})^{-1}\bar{Z}^{\top}Y$
(11)

and all we have is that

$Y=\bar{U}\bar{\gamma}+\varepsilon.$

Thus, we have

$\mathbb{E}(\hat{\bar{\zeta}})=\mathbb{E}\left[(\bar{Z}^{\top}\bar{Z})^{-1}\bar{Z}^{\top}Y\right]=(\bar{Z}^{\top}\bar{Z})^{-1}\bar{Z}^{\top}\bar{U}\bar{\gamma}$

or,

$(\bar{Z}^{\top}\bar{Z})\mathbb{E}(\bar{\zeta})=\bar{Z}^{\top}\bar{U}\bar{\gamma}.$
(12)

First,

$\bar{Z}^{T}\bar{Z}=\left(\begin{array}[]{c c c c}X_{\cdot 1}^{\top}X_{\cdot
1}&X_{\cdot 1}^{\top}\bar{U}_{(1)}&\dots&X_{\cdot 1}^{\top}\bar{U}_{(k)}\\\
X_{\cdot 1}^{\top}\bar{U}_{(1)}&1&&\\\ \vdots&&\ddots&\vdots\\\ X_{\cdot
1}^{\top}\bar{U}_{(k)}&0&\dots&1\end{array}\right)$

and then

$\bar{Z}^{T}\bar{U}=\left(\begin{array}[]{c c c c c c}&&X_{\cdot
1}^{\top}\bar{U}&&&\\\ \hline\cr 1&\dots&0&0&\dots&0\\\
\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
0&\dots&1&0&\dots&0\end{array}\right)\Rightarrow\bar{Z}^{T}\bar{U}\bar{\gamma}=\left(\begin{array}[]{c}X_{\cdot
1}^{\top}\bar{U}\bar{\gamma}\\\ \bar{\gamma}_{1}\\\ \vdots\\\
\bar{\gamma}_{k}\end{array}\right).$

Equation (12) above becomes:

$\left\\{\begin{array}[]{l l}X_{\cdot 1}^{\top}X_{\cdot
1}\mathbb{E}(\hat{\bar{\alpha}})+\sum_{i=1}^{k}X_{\cdot
1}^{\top}\bar{U}_{(i)}\mathbb{E}(\hat{\bar{\gamma}}_{i})&=X_{\cdot
1}^{\top}\bar{U}\bar{\gamma}\\\ X_{\cdot
1}^{\top}\bar{U}_{(1)}\mathbb{E}(\hat{\bar{\alpha}})+\mathbb{E}(\hat{\bar{\gamma}}_{1})&=\bar{\gamma}_{1}\\\
\quad\quad\vdots&\vdots\\\ X_{\cdot
1}^{\top}\bar{U}_{(k)}\mathbb{E}(\hat{\bar{\alpha}})+\mathbb{E}(\hat{\bar{\gamma}}_{k})&=\bar{\gamma}_{k}.\end{array}\right.$

The generic equations, for $1\leq j\leq k$, can be rewritten as

$X_{\cdot
1}^{\top}\bar{U}_{(j)}\mathbb{E}(\hat{\bar{\alpha}})+\mathbb{E}(\hat{\bar{\gamma}}_{j})=\bar{\gamma}_{j},$

giving

$\mathbb{E}(\hat{\bar{\gamma}}_{j})=\bar{\gamma}_{j}-X_{\cdot
1}^{\top}\bar{U}_{(j)}\mathbb{E}(\hat{\bar{\alpha}}).$

Plugging this into the first equation gives

$\mathbb{E}(\hat{\bar{\alpha}})\left[X_{\cdot 1}^{\top}X_{\cdot
1}-\sum_{i=1}^{k}\left(X_{\cdot
1}^{\top}\bar{U}_{(i)}\right)^{2}\right]+\sum_{i=1}^{k}X_{\cdot
1}^{\top}\bar{U}_{(i)}\bar{\gamma}_{i}=X_{\cdot 1}^{\top}\bar{U}\bar{\gamma},$

that is

$\mathbb{E}(\hat{\bar{\alpha}})\left[X_{\cdot 1}^{\top}X_{\cdot
1}-\sum_{i=1}^{k}\left(X_{\cdot
1}^{\top}\bar{U}_{(i)}\right)^{2}\right]+X_{\cdot
1}^{\top}\bar{U}_{1:k}\bar{\gamma}_{1:k}=X_{\cdot
1}^{\top}\bar{U}\bar{\gamma}$

and thus

$\mathbb{E}(\hat{\bar{\alpha}})\left[X_{\cdot 1}^{\top}X_{\cdot
1}-\sum_{i=1}^{k}\left(X_{\cdot
1}^{\top}\bar{U}_{(i)}\right)^{2}\right]=X_{\cdot
1}^{\top}\bar{U}_{(k+1):\bar{r}}\bar{\gamma}_{(k+1):\bar{r}}.$

Finally,

$\mathbb{E}(\hat{\bar{\alpha}})=\frac{X_{\cdot
1}^{\top}\bar{U}_{(k+1):\bar{r}}\bar{\gamma}_{(k+1):\bar{r}}}{X_{\cdot
1}^{\top}X_{\cdot 1}-\left\|X_{\cdot 1}^{\top}\bar{U}_{1:k}\right\|^{2}}.$

Now for the bias,

$\displaystyle{\rm bias}(\hat{\bar{\alpha}})$
$\displaystyle=\mathbb{E}(\hat{\bar{\alpha}})-\beta_{1}$
$\displaystyle=\frac{X_{\cdot
1}^{\top}(\bar{U}_{(k+1):\bar{r}}\bar{\gamma}_{(k+1):\bar{r}}-X_{\cdot
1}\beta_{1})+\beta_{1}\left\|X_{\cdot
1}^{\top}\bar{U}_{1:k}\right\|^{2}}{X_{\cdot 1}^{\top}X_{\cdot
1}-\left\|X_{\cdot 1}^{\top}\bar{U}_{1:k}\right\|^{2}}.$

#### Variance of PSC

From (11), we have

$\displaystyle{\rm
Var}(\hat{\bar{\zeta}})=(\bar{Z}^{\top}\bar{Z})^{-1}\bar{Z}^{\top}{\rm
Var}(Y)\bar{Z}(\bar{Z}^{\top}\bar{Z})^{-1}=(\bar{Z}^{\top}\bar{Z})^{-1}\sigma^{2}$

or

$\displaystyle(\bar{Z}^{\top}\bar{Z}){\rm
Var}(\hat{\bar{\zeta}})=\sigma^{2}\mathbf{I}.$

As we are only interested in estimating the variance of $\hat{\bar{\alpha}}$,
from the above formula we obtain

$\displaystyle\begin{cases}(X_{\cdot 1}^{\top}X_{\cdot 1}){\rm
Var}(\hat{\bar{\alpha}})+\sum_{j=1}^{k}Cov(\hat{\bar{\alpha}},\hat{\bar{\gamma}}_{j})X_{\cdot
1}^{\top}U_{(j)}=\sigma^{2},\\\ (X_{\cdot 1}^{\top}U_{(j)}){\rm
Var}(\hat{\bar{\alpha}})+Cov(\hat{\bar{\alpha}},\hat{\bar{\gamma}}_{j})=0.\end{cases}$

Substituting the second equation into the first to get

$\displaystyle{\rm Var}(\hat{\bar{\alpha}})=\frac{\sigma^{2}}{X_{\cdot
1}^{\top}X_{\cdot 1}-\sum_{j=1}^{k}(X_{\cdot 1}^{\top}\bar{U}_{(j)})^{2}}.$

###  A.3 Proof for Theorem 3

###### Proof for Theorem 3..

We have that

$\displaystyle\left\|X_{(-1)}^{\top}u\right\|^{2}=\left\|(\mathbf{0},X_{2},\ldots,X_{p})^{\top}u\right\|^{2}\leq\left\|X^{\top}u\right\|^{2}.$

And from the definition that
$\bar{U_{j}}=\arg\max_{\|u\|^{2}=1}\left\|X^{\top}u\right\|^{2}$, we get

$\displaystyle\left\|X^{\top}U_{1:k}\right\|^{2}\leq\left\|X^{\top}\bar{U}_{1:k}\right\|^{2}$

or

$\displaystyle\left\|X_{\cdot
1}^{\top}U_{1:k}\right\|^{2}+\sum_{j=2}^{p}\left\|X_{(j)}^{\top}U_{1:k}\right\|^{2}$
$\displaystyle\leq\left\|X_{\cdot
1}^{\top}\bar{U}_{1:k}\right\|^{2}+\sum_{j=2}^{p}\left\|X_{j}^{\top}\bar{U}_{1:k}\right\|^{2},$
$\displaystyle\left\|X_{\cdot
1}^{\top}U_{1:k}\right\|^{2}+\max_{\|x\|_{2}^{2}=1}\sum_{j=2}^{p}\left\|X_{(j)}^{\top}x_{1:k}\right\|^{2}$
$\displaystyle\leq\left\|X_{\cdot
1}^{\top}\bar{U}_{1:k}\right\|^{2}+\sum_{j=2}^{p}\left\|X_{j}^{\top}\bar{U}_{1:k}\right\|^{2}.$

By taking $x=\bar{U}_{1:k}$, we obtain

$\displaystyle\left\|X_{\cdot
1}^{\top}U_{1:k}\right\|^{2}+\sum_{j=2}^{p}\left\|X_{(j)}^{\top}\bar{U}_{1:k}\right\|^{2}\leq\left\|X_{\cdot
1}^{\top}\bar{U}_{1:k}\right\|^{2}+\sum_{j=2}^{p}\left\|X_{j}^{\top}\bar{U}_{1:k}\right\|^{2}.$

And thus

$\displaystyle\left\|X_{\cdot 1}^{\top}U_{1:k}\right\|^{2}\leq\left\|X_{\cdot
1}^{\top}\bar{U}_{1:k}\right\|^{2}.$

This yields the conclusion of the theorem. ∎

##  Appendix B Simulation results in other cases

Figure 7: Binary $X$. Estimates of $\beta_{1}$ with different number of
Principal components (PCs)

Figure 8: Structured $X$. Estimates of $\beta_{1}$ with different number of
Principal components (PCs)

Figure 9: Independent $X$. Estimates of $\beta_{1}$ with different number of
Principal components (PCs)

Figure 10: Dependent $X$. Estimates of $\beta_{1}$ with different number of
Principal components (PCs)

## References

  * [Brzyski et al., 2017] Brzyski, D., Peterson, C. B., Sobczyk, P., Candès, E. J., Bogdan, M., and Sabatti, C. (2017).  Controlling the rate of gwas false discoveries.  Genetics, 205(1):61–75.
  * [Bühlmann and Van De Geer, 2011] Bühlmann, P. and Van De Geer, S. (2011).  Statistics for high-dimensional data: methods, theory and applications.  Springer Science & Business Media.
  * [Buzdugan et al., 2016] Buzdugan, L., Kalisch, M., Navarro, A., Schunk, D., Fehr, E., and Bühlmann, P. (2016).  Assessing statistical significance in multivariable genome wide association analysis.  Bioinformatics, 32(13):1990–2000.
  * [Bycroft et al., 2018] Bycroft, C., Freeman, C., Petkova, D., Band, G., Elliott, L. T., Sharp, K., Motyer, A., Vukcevic, D., Delaneau, O., O’Connell, J., et al. (2018).  The UK biobank resource with deep phenotyping and genomic data.  Nature, 562(7726):203.
  * [Derks et al., 2017] Derks, E., Zwinderman, A., and Gamazon, E. (2017).  The relation between inflation in type-i and type-ii error rate and population divergence in genome-wide association analysis of multi-ethnic populations.  Behavior genetics, 47(3):360–368.
  * [Dicker, 2014] Dicker, L. H. (2014).  Variance estimation in high-dimensional linear models.  Biometrika, 101(2):269–284.
  * [Efron and Hastie, 2016] Efron, B. and Hastie, T. (2016).  Computer age statistical inference, volume 5.  Cambridge University Press.
  * [Giraud, 2014] Giraud, C. (2014).  Introduction to high-dimensional statistics.  Chapman and Hall/CRC.
  * [Hastie et al., 2009] Hastie, T., Tibshirani, R., and Friedman, J. (2009).  The elements of statistical learning: data mining, inference, and prediction.  Springer Science & Business Media.
  * [Janson et al., 2017] Janson, L., Barber, R. F., and Candes, E. (2017).  Eigenprism: inference for high dimensional signal-to-noise ratios.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(4):1037–1065.
  * [Javanmard and Montanari, 2014] Javanmard, A. and Montanari, A. (2014).  Confidence intervals and hypothesis testing for high-dimensional regression.  The Journal of Machine Learning Research, 15(1):2869–2909.
  * [Leeb and Pötscher, 2008] Leeb, H. and Pötscher, B. M. (2008).  Sparse estimators and the oracle property, or the return of hodges’ estimator.  Journal of Econometrics, 142(1):201–211.
  * [Lees et al., 2020] Lees, J. A., Mai, T. T., Galardini, M., Wheeler, N. E., Horsfield, S. T., Parkhill, J., and Corander, J. (2020).  Improved prediction of bacterial genotype-phenotype associations using interpretable pangenome-spanning regressions.  MBio, 11(4):e01344–20.
  * [Lees et al., 2016] Lees, J. A., Vehkala, M., Välimäki, N., Harris, S. R., Chewapreecha, C., Croucher, N. J., Marttinen, P., Davies, M. R., Steer, A. C., Tong, S. Y., et al. (2016).  Sequence element enrichment analysis to determine the genetic basis of bacterial phenotypes.  Nature communications, 7:12797.
  * [Lippert et al., 2011] Lippert, C., Listgarten, J., Liu, Y., Kadie, C. M., Davidson, R. I., and Heckerman, D. (2011).  Fast linear mixed models for genome-wide association studies.  Nature methods, 8(10):833.
  * [Listgarten et al., 2012] Listgarten, J., Lippert, C., Kadie, C. M., Davidson, R. I., Eskin, E., and Heckerman, D. (2012).  Improved linear mixed models for genome-wide association studies.  Nature methods, 9(6):525.
  * [Lounici, 2008] Lounici, K. (2008).  Sup-norm convergence rate and sign concentration property of lasso and dantzig estimators.  Electronic Journal of statistics, 2:90–102.
  * [Pérez and de Los Campos, 2014] Pérez, P. and de Los Campos, G. (2014).  Genome-wide regression and prediction with the bglr statistical package.  Genetics, 198(2):483–495.
  * [Price et al., 2006] Price, A. L., Patterson, N. J., Plenge, R. M., Weinblatt, M. E., Shadick, N. A., and Reich, D. (2006).  Principal components analysis corrects for stratification in genome-wide association studies.  Nature genetics, 38(8):904.
  * [Price et al., 2010] Price, A. L., Zaitlen, N. A., Reich, D., and Patterson, N. (2010).  New approaches to population stratification in genome-wide association studies.  Nature Reviews Genetics, 11(7):459.
  * [Sun and Zhang, 2012] Sun, T. and Zhang, C.-H. (2012).  Scaled sparse linear regression.  Biometrika, 99(4):879–898.
  * [Van de Geer et al., 2014] Van de Geer, S., Bühlmann, P., Ritov, Y., and Dezeure, R. (2014).  On asymptotically optimal confidence regions and tests for high-dimensional models.  The Annals of Statistics, 42(3):1166–1202.
  * [Visscher et al., 2017] Visscher, P. M., Wray, N. R., Zhang, Q., Sklar, P., McCarthy, M. I., Brown, M. A., and Yang, J. (2017).  10 years of gwas discovery: biology, function, and translation.  The American Journal of Human Genetics, 101(1):5–22.
  * [Wu et al., 2009] Wu, T. T., Chen, Y. F., Hastie, T., Sobel, E., and Lange, K. (2009).  Genome-wide association analysis by lasso penalized logistic regression.  Bioinformatics, 25(6):714–721.
  * [Yang et al., 2014] Yang, J., Zaitlen, N. A., Goddard, M. E., Visscher, P. M., and Price, A. L. (2014).  Advantages and pitfalls in the application of mixed-model association methods.  Nature genetics, 46(2):100.
  * [Yang, 2005] Yang, Y. (2005).  Can the strengths of aic and bic be shared? a conflict between model indentification and regression estimation.  Biometrika, 92(4):937–950.
  * [Zhang and Zhang, 2014] Zhang, C.-H. and Zhang, S. S. (2014).  Confidence intervals for low dimensional parameters in high dimensional linear models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76(1):217–242.
  * [Zhao and Yu, 2006] Zhao, P. and Yu, B. (2006).  On model selection consistency of lasso.  The Journal of Machine Learning Research, 7:2541–2563.


# Building population models for large-scale neural recordings: opportunities
and pitfalls

Cole Hurwitz*, Nina Kudryashova*, Arno Onken and Matthias H. Hennig†
University of Edinburgh, Institute for Adaptive and Neural Computation
Edinburgh, EH8 9AB, United Kingdom
$\dagger$ Correspondence<EMAIL_ADDRESS>
* equal contribution

(April, 2021)

###### Abstract

Modern recording technologies now enable simultaneous recording from large
numbers of neurons. This has driven the development of new statistical models
for analyzing and interpreting neural population activity. Here we provide a
broad overview of recent developments in this area. We compare and contrast
different approaches, highlight strengths and limitations, and discuss
biological and mechanistic insights that these methods provide.

## Introduction

Large scale recordings from neurons with dense, high channel count probes now
yield sufficient data to understand the precise covariation in neural
populations [35, 65, 79]. This departure from analyzing the tuning properties
of individual neuron, often pooled across experiments, requires model-based
approaches that capture the complexity of high-dimensional population activity
while remaining computationally tractable.

Statistical models have emerged as an essential tool for overcoming this
challenge. Informally, these models can be grouped into two categories: fully
observed models and latent variable models. Fully observed models try to
explicitly capture the interactions between neurons by directly modeling the
joint activity of the population. Since modeling the full space of activity
patterns is intractable, a key ingredient for tractable fully observed models
is an efficient description of the interactions between neurons that retains
essential characteristics. Latent variable models, on the other hand, assume
that population activity can be summarized by a small set of variables called
latent factors. Rather than directly modeling interactions, latent variable
models capture neuronal interactions through activation of these factors.

This conceptual difference means that the specific hypotheses each of these
model classes can address are also different and that conclusions drawn from
such analysis have to be considered in the light of the main assumptions the
model is based upon. In this review, we first briefly discuss how single
neuron activity is extracted from high-density recordings. Next, we provide an
overview of current developments and show which questions can and cannot be
addressed with each model class. We close with a discussion of limitations and
future directions.

## Extracting modeling data

Figure 1: The process of extracting the activity of single neurons with a
spike sorter is error-prone with different sorters having different failure
modes. It is, therefore, recommended to make use of different sorters before
performing model-based analysis. A first method is to perform the same
analysis on the output of different sorters to see if functional insights are
independent of the chosen spike sorter (bottom). Alternatively, one could
obtain a consensus sorting (top) where false positive units are largely
removed by taking the agreement between different sorters. Any analyses could
then be performed on this consensus sorting. These pipelines, together with
various pre- and postprocessing steps as well as automated curation, can be
easily implemented using the SpikeInterface framework [9] which supports all
commonly used spike sorters.

Large scale recordings can be obtained with high density extracellular probes
or using calcium imaging. From extracellular recordings, spike trains of
individual neurons are extracted through a process called spike sorting. There
are several algorithmic approaches for spike sorting of data from dense probes
which are reviewed elsewhere [68, 30, 12, 41]. Accurately extracting spike
trains is essential since mistakes, including incorrectly merged or split
units, unresolved overlapping spikes, or incorrect clustering due to probe
drift, have been shown to bias subsequent analysis [4, 91, 92]. For example,
in regions with high-firing rates and synchronicity, spike sorting mistakes
can lead to artificial correlations between independently firing neurons [4,
17].

It is important to acknowledge that none of the existing spike sorting
algorithms are error-free, and that the scale of these recordings now makes
manual curation difficult. To aid selection of an appropriate spike sorter, a
large collection of ground truth comparisons is available on the SpikeForest
website [45]. Moreover, the comparison of spike sorter outputs has shown
surprisingly low agreement in part due to a large number of false positive
units [9]. A viable approach to mitigate against spurious findings is,
therefore, to perform the same analysis on multiple sorter outputs or on an
ensemble agreement among multiple sorters; this can be done using a software
framework such as SpikeInterface (Figure 1) [9].

While calcium imaging allows unequivocal identification of single neurons, the
much lower sampling rate and the binding kinetics of calcium indicators mean
the recorded signal may be an incomplete record of the spiking activity of
these neurons [32]. A first approach is to recover neural spike trains using
deconvolution methods. Similar to spike sorting however, different algorithms
have been shown to produce conflicting results when analyzing correlations
between neurons [19]. It is therefore prudent to avoid analysis based only on
a single convolution algorithm. Second, recent work aims to to include a model
of the dynamics giving rise to the calcium signal directly into statistical
models [63]. This approach is promising in particular for latent variable
models if the different sources of variation (neural and calcium dynamics) can
be effectively disentangled.

## Fully observed models

Figure 2: Fully observed models construct a joint distribution of population
activity with various levels of detail. The horizontal arrow orders the
methods according to their expressivity (i.e. how much dependence structure
the method can account for).

When recording from just 20 neurons, there are over one million possible
instantaneous ON/OFF patterns of spiking ($2^{20}$) for a small time bin; this
number grows to one billion for 30 neurons ($2^{30}$). While the number of
spike patterns that are actually observed is much smaller and constrained by
the connectivity between neurons, estimating all pattern frequencies reliably
from typical recordings is impossible even for small populations. Fully
observed statistical models are a tool to achieve precisely this by
approximating this full distribution using as few parameters as possible
together with suitable constraints [3, 20, 100]. As always, different
approaches exist with different trade-offs between complexity, tractability
and interpretability (see Fig. 2 and Tab. 1). A main purpose of these models
is to ask how much information recorded activity contains about stimuli or
behavior, which can be addressed using information-theoretical or decoding
approaches.

#### Models

Maximum Entropy models (MaxEnt) provide a principled way to construct
probability distributions for spiking activity by finding the least structured
distribution that fulfills specific constraints obtained from the data [78].
Typical constraints for MaxEnt models are the firing rates of and pairwise
correlations between the neurons [80, 83, 88] and the population synchrony
[88]. While MaxEnt models have been shown to model spiking activity well, the
computational cost of fitting them is prohibitive for large neuronal
populations. This has led to the development of analytically tractable models
constrained by the dependencies between individual neurons and the aggregate
population activity [88, 25, 53, 51]. Related and more scalable is the
Dichotomized Gaussian model [2, 43] which can be used to model spike count
variables with arbitrary marginal distributions [43].

MaxEnt models cannot easily take external inputs into account. If these are
well specified, generalized linear models (GLMs) can successfully describe the
activity of large populations [90, 60, 61]. In these models, external inputs,
inputs from other neurons, and each neuron’s spiking history are weighted,
summed up and transformed by a rate function, which drives a stochastic
process to model spiking activity. For typical rate functions, fitting GLMs is
a convex optimization problem which can be solved efficiently, making GLMs
scalable to large-scale neural recordings [90, 61]. Along with modeling the
observed neural population, the activity of unobserved neurons can be inferred
[5], and knowing the topology of the external covariates can further aid in
inferring the activity of the unseen neurons [86]. Despite their flexibility,
GLMs cannot model networks with neurons that perform non-linear integration of
multiple external inputs [6].

A more direct representation of statistical dependencies is provided by copula
models, which decompose neural population activity into the distributions of
individual neuron activities and the dependence between these neurons,
represented by a copula [34]. Copula models can represent inter-neuronal
dependencies that go beyond linear (Gaussian) correlations [7], as well as
non-linear stimulus-response relationships [37]. These highly flexible models
can for example be used to jointly model different recording modalities such
as local field potentials and individual neurons’ spiking activities [1, 54].
They are particularly promising for analyzing recordings of neural data and
exploring its relation to behavior and stimuli which typically follow
different statistics and tend to have complex dependencies [37, 33].

#### Interpretation and biological insights

Model |  References |  Number of parameters |  Closed-form pattern probabilities? |  Fit for large $N$? |
---|---|---|---|---|---
Dichotomized Gaussian |  Amari et al., [2] Macke et al., [43] | $\sim N^{2}$ | No | Yes |
Pairwise MaxEnt |  Schneidman et al., [80] Shlens et al., [83] | $\sim N^{2}$ | No1 | Difficult |
Tractable MaxEnt |  Tkačik et al., [88] Gardella et al., [25] O’Donnell et al., [51] | $\sim N^{2}$ | Yes1 | Yes |
GLM |  Pillow et al., [61] | $\sim DN^{2}$ | No | Difficult |
Vine Copula |  Aas et al., [1] Onken and Panzeri, [54] | $\sim DN^{2}$ | Yes | Yes |
Table 1: Fully observed models: A table (adapted from O’Donnell et al., [51])
characterizing model properties and limitations. Here, $N$ is the number of
neurons and $D$ is the number of coefficients per interaction term, such as
filter sizes for GLM or number of parameters for parametric copula families.
Sampling is possible from all of these models. 1. For a more detailed
comparison of MaxEnt models, see Table 1 in [78].

Fully observed models can be used in two main areas. First, they provide
access to the full joint activity distribution of the population in place of
direct, typically biased estimates from the recordings for information
theoretic analyses and for decoding [33, 37, 71, 59]. GLMs and copula models
that were fit to population activity have been shown to have high decoding
performance of external variables [90, 61, 38, 39, 33, 37]. For instance, GLMs
with coupling filters were shown to capture 40% more visual information from
the retina than optimal linear decoding [61], indicating that GLMs can model
additional details in the activity that are relevant for representing the
stimulus. In large-scale datasets with complex statistics, GLMs may however
not be appropriate. In these cases, copula models can be useful for measuring
the amount of information in the population response about external covariates
[54, 33, 37]. Additionally, for complex stimuli with spatial or temporal
correlations, one must distinguish ‘stimulus correlations’ that arise from
stimulus statistics from ‘noise correlations’ due to circuit interactions.
Including this yields better performance in MaxEnt models [26], and prevents
spurious effects in GLMs [46]. More generally, comparisons of these models
with linear decoders can quantify the information that emerges from the
network interactions, i.e. how the whole network is different from the sum of
its parts [60, 61, 80, 52].

Second, the parameters in many fully observed models can be interpreted as the
strength of the interactions between neurons. For example, the coupling
parameters in dichotomized Gaussian and in pairwise MaxEnt models describe the
strength of pairwise interactions between single neurons [43, 80]. As these
connections are inferred from the population activity alone, they are usually
referred to as functional interactions or connections, to emphasize the
difference to synaptic connections. This analysis has been used to uncover the
spatial extent of interactions in the retina [83, 52] and to characterize
interactions between neurons with different selectivity in the entorhinal
cortex [14].

Whether or not inferred functional connectivity is biologically interpretable
critically depends on the identifiability of the relevant parameters, and here
caution should be exercised. For instance, while the coupling filters in GLMs
have unique solutions, filter parameter changes do not necessarily change the
GLM responses much [8] suggesting that, in practice, coupling filters may not
be properly identifiable. Similarly, in MaxEnt models inferred functional
connectivity is not well constrained by the recorded spike trains. As a
result, many coupling parameters, even those with large absolute values, can
be altered without significantly changing the network activity predicted by
the model [55, 62].

Limited data can be a main contributor to poor parameter identifiability. To
address this, the uncertainty of the model parameters can be quantified with
approximate Bayesian inference methods [99]. A recently developed LR-GLM
utilized low-rank data approximations to scale approximate Bayesian inference
methods to high-dimensional real datasets making these methods more applicable
to large-scale neural recordings [89]. However, closer inspection of MaxEnt
models suggests weakly constrained couplings could also be an intrinsic
property of neural circuits. Simulations have shown that re-wiring even a
large number of synapses may only have minor functional consequences [48],
hence the corresponding inferred connections are also weakly specified.
Analysis of long term recordings has shown that these less important
connections are also subject to more intensive re-modeling over time [55, 62].
Therefore, taking into account the reliability of the inferred couplings in
functional connectivity models allows for the identification of the neurons,
connections and circuit motifs that most determine the population activity
statistics [31, 55, 62].

Figure 3: Latent variable models assume that neural activity is generated by a
set of low-dimensional latent states. (a) When modeling the temporal evolution
of the latent states, one can explicitly model the dynamics with a state-space
model or can instead model the statistics of the latent trajectories with a
Gaussian Process (GP). (b) The latent trajectories can be related back to
neural activity (e.g. the firing rates) with a variety of mapping functions
including linear functions, neural networks, and GPs. (c) An additional step
to further constrain the dynamics is to reconstruct the spiking activity using
a point-process distribution such as the Poisson distribution. (d) Recently,
it has been shown that explicitly modeling the relationship between the latent
factors and the observed behavioral variables can further constrain the
factors to be low-dimensional and ’behaviorally relevant’ [76]

## Latent variable models

Despite the apparent high-dimensionality of population activity, recent
studies suggest that the activity patterns underlying neural function are
confined to low-dimensional manifolds [47, 13, 74, 29, 22, 18]. In other
words, neural computation results from the activation of specific population-
wide activity patterns rather than the independent modulation of single
neurons [21]. In this regime, latent variable models are a natural choice as
they aim to capture arbitrarily complicated response structure in the neural
activity with only a few latent variables.

### Models

Latent variable models aim to capture the population activity structure
through the temporal evolution of their learned latent variables. There are
two main approaches for modeling the latent trajectories (shown in Figure 3
and summarized in Table 2): State-space representations, where the dynamics of
the latent trajectories are modeled explicitly with an evolution function that
relates past latent states to future latent states, and moment representations
where, rather than having an explicit model of the dynamics, the statistics of
the latent trajectories over time are modeled. With both approaches, there is
an observation model which relates the neural trajectories back to the
observed activity.

Traditional dimensionality reduction techniques such as principal component
analysis (PCA) and factor analysis (FA) can be re-interpreted as static state-
space models without temporal dependence between states [69]. Extending these,
linear dynamical state-space models such as linear dynamical systems (LDS)
[85, 38, 10, 58] and jPCA [13, 82, 50] can capture simple temporal
dependencies between latent states. While efficient, these methods cannot
model non-linear dynamics which are thought to underlie rhythmic motor
patterns [72, 28], decision making [64], and pathologies such as epilepsy
[27]. To capture these dynamics, non-linear dynamical state-space models have
been introduced including approaches that model dynamics with recurrent neural
networks (RNNs) [56, 81], piecewise-linear RNNs [15], switching linear
dynamics [57], and recurrent switching linear dynamics [42, 98]. A
particularly efficient RNN-based model is the recently introduced latent
factor analysis via dynamical systems (LFADS) which utilizes neural networks
to reduce the computational cost of inferring dynamics, making it more
practical for large-scale neural recordings [56]. Despite these recent
developments, non-linear state-space models are typically computationally
expensive and parameter inference is difficult [73].

An alternative to explicitly learning a dynamical model for the neural
population is to model a statistical representation of the trajectories,
typically as a Gaussian process (GP). GPs define a distribution over functions
and allow to constrain certain features of the function such as its
smoothness. Unlike state-space models, GP-based models provide both
uncertainty quantification and principled model selection. The most common GP-
based latent variable model is the Gaussian-process factor analysis (GPFA)
[11]. In this model, each trajectory of the latent state is sampled from an
independent GP. This independence constraint, however, limits the ability of
GPFA to model joint temporal dependencies making them less useful for
uncovering temporal structure in neural data; for instance, GPFA cannot
disentangle motor cortex trajectories during hand reach [73]. Gaussian Process
Factor Analysis with Dynamical Structure (GPFADS) seeks to remedy this by
constraining the learned trajectories to have lower probabilities of occurring
in reverse [73]. As a result, GPFADS can disentangle the trajectories
generated by non-linear non-reversible dynamical systems such as primary motor
cortex arm reach data. Overall, while GP-based latent variable models provide
useful analytical properties and are generally quite efficient, they still
need to demonstrate that they can capture the complex temporal dynamics that
underlie neural activity as well as state-space models can.

An important aspect of all latent variable models is the observation model
which relates the trajectories back to the neural activity. Common approaches
include modeling neural activity with a conditional Gaussian model or the more
plausible conditional Poisson model [44]. The mapping between the latent
trajectories and neural activity distribution can take on a number of
functional forms including linear maps [85, 11, 56, 96], GPs [40, 94, 81], or
neural networks [24]. Most models, however, rely on a simple linear mapping
function as it forces the latent variables to capture more information about
the joint neural activity. Recently, observation models that take into account
external behavioral variables have been introduced for capturing ‘behaviorally
relevant’ latent trajectories. Sani et al., [76] introduced a novel linear
Gaussian state-space model (PSID) that jointly models the neural activity and
the recorded behavioral variables. Intriguingly, PSID found that behaviorally
relevant latent trajectories are significantly lower-dimensional than
previously thought.

#### Interpretation and biological insights

Unlike fully observed models, latent variable models provide a very succinct
description of joint population activity in the form of latent trajectories. A
simple yet effective approach to interpreting these latent trajectories is
visualization; reducing the number of latent variables down to just two or
three allows for insightful visualizations. For example, Santhanam et al.,
[77] used factor analysis to reduce preparatory neural activity in the
premotor cortex to just three dimensions that, when visualized, provided
evidence that the population activity contained information that could
discriminate between goal-related target conditions in a delayed center-out
reach task.

It may however not be appropriate to reduce the number of latent variables to
two or three for visualization without sacrificing model fit111Model fit can
be measured in different ways including the amount of variance explained or
the cross-validated likelihood of observed activity given the model. GP-based
latent variable models can also provide uncertainty estimates for the latent
trajectories.. A more principled and flexible way to interpret the learned
latent trajectories is to relate them back to a measured external variable
[93]. For example, one can measure how informative a latent variable is about
an external variable by decoding the external variable at each time step using
a simple method such as a least squares regression. While this approach has
been used to interpret population response structure in the olfactory and
visual system [75, 16, 84, 95, 87], its value is best demonstrated in studies
that correlate latent trajectories in the premotor and motor cortex to arm
reaching tasks. Churchland et al., [13] found that low-dimensional rotational
trajectories in the primary motor cortex were correlated with arm reaches.
Gallego et al., [21] extended this finding, demonstrating that these low-
dimensional trajectories are a stable neural correlate for consistent
execution of arm reaches over the course of many years. Moreover, Ramanathan
et al., [66] showed that following a stroke, diminished reaching function in
the contralesional arm was correlated with a loss of motor cortical neural
trajectories. The trajectories only reemerged after motor recovery and proved
to be a useful neuromodulatory target for therapeutic electrical stimulation.

Despite the success of latent variable models, their usefulness can be
dependent on the properties of the neuronal population and the specific
experimental setting. Firstly, latent variable models are most useful when the
neural population activity is, in fact, low-dimensional. While a number of
studies have independently found this to be the case for their datasets,
recent work by Stringer et al., [87] suggests that the manifold dimensionality
of stimulus-evoked activity in the visual cortex is actually as high as it can
be without becoming non-differentiable. Moreover, Gao et al., [23] suggests
that population activity is as high-dimensional as possible given the
simplicity of the given stimuli or tasks. While latent variable models can
still be used to model this data (by increasing the number of latent
variables), they may no longer be as interpretable. Another assumption in some
models is that the latent trajectories are autonomous222Pandarinath et al.,
[56] showed that LFADS could capture task-related inputs to a neural
population by incorporating them into the generative process. It is still an
open question how generalizable this approach is to other types of inputs.,
i.e. they only depend on an initial condition (usually trial-specific) and an
evolution function. This limits the ability of these models to capture
unmeasured inputs to the neural population from other brain regions. Finally,
and most importantly, latent variable models, by design, will mix all sources
of neural variability in the latent space [93]. This can make interpreting the
latent variables challenging in particular for complex behaviors or diverse
stimulus ensembles. These constraints may be the reason why many of the
conceptual breakthroughs provided by latent variable models are for motor
cortical activity where the dimensionality appears low, the trajectories are
largely autonomous given the trial conditions, and the variation is largely
explained by the measured behavior. A possible way forward to improve the
interpretability of these models in more complex experimental regimes is to
explicitly model external variables either in the observation model [36, 76]
or in the latent space [97].

Model | Trajectories | Mapping Function | Observation | Single-trial1
---|---|---|---|---
PCA/FA [69] | Static | Linear | Gaussian | No
dPCA [36] | Static | Linear | Gaussian2 | No
jPCA [13] | Linear3 | Linear | Gaussian | No
LDS [85] | Linear | Linear | Gaussian | Yes
PLDS [44] | Linear | Linear | Poisson | Yes
PSID [76] | Linear | Linear | Gaussian4 | Yes
PfLDS [24] | Linear | Neural Network | Poisson | Yes
SLDS [57] | Switching Linear | Linear | Gaussian | Yes
RSLDS [42, 98] | Recurrent Switching Linear5 | Linear | Gaussian | Yes
PLRNN-SSM [15] | Piecewise-linear RNN | Linear | Gaussian | Yes
LFADS [56] | RNN | Linear | Poisson | Yes
GP-RNN [81] | RNN | GP | Poisson/Gaussian | Yes
GPFA [11] | GP | Linear | Gaussian | Yes
GPFADS [73] | GP6 | Linear | Gaussian | Yes
vLGP [96] | GP | Linear | Poisson7 | Yes
P-GPLVM [94] | GP | GP | Poisson | Yes
Table 2: Latent Variable Models: A table (adapted from [81]) characterizing
the latent variable models referenced in this review. 1. The ability to
extract dynamics from single-trials is important for capturing important
variability that can be obscured with trial-averaging [49]. 2. dPCA models
trial-averaged and task averaged neural activity to demix different sources of
variability in the data. 3. jPCA is equivalent to a linear dynamical model
with specific constraints for modeling rotational dynamics [50]. 4. PSID
models both the observed neural activity and behaviour with separate Gaussian
observation models. 5. RSLDS can model multiple neural populations and their
time-varying interactions using switching linear dynamics that is governed by
a discrete latent state. 6. GPFADS introduces a new GP kernel which allows for
modeling the dynamical behaviour of temporal-irreversibility. 7. The activity
of each neuron in vLGP is linearly dependent on the latent variable and the
self-history of the neuron.

## Outlook

While both fully observed and latent variable models can be used to overcome
the curse of dimensionality in large-scale recordings, the two classes yield
distinct yet complementary insights into circuit function. In addition, they
each have important limitations that are not easily resolved. Fully observed
models can suffer from poor parameter identifiability because data sets have a
limited size and potentially because many parameters are actually not well
constrained in real neural circuits. Possible future direction are to either
simplify the description of the interactions to reduce the number of
parameters [67] or to augment a model with biologically informed latent
dynamics [70].

Latent variable models assume that circuit dynamics are inherently low-
dimensional and place specific constraints on how these dynamics unfold. It is
perhaps for this reason that we have, so far, few insights into neural
dynamics during complex naturalistic behavior; experimental design still
favors simple, repeated stimuli and behaviors where sources of neural
variability are easy to de-mix. In more complex experimental regimes, it may
no longer be sufficient to generically model the response structure of a
neural population and perform ex post facto analyses to relate it back to the
variables of interest. Instead, there is increasing evidence that approaches
that explicitly model the measured external variables have improved
interpretability and can reveal novel insights about neural function [76, 97].

## Acknowledgments

AO is supported by the Engineering and Physical Sciences Research Council (EP/
S005692/1). CLH is supported by the Thouron Award and by the School of
Informatics, University of Edinburgh.

## References

  * Aas et al., [2009] Aas, K., Czado, C., Frigessi, A., and Bakken, H. (2009).  Pair-copula constructions of multiple dependence.  Insurance: Mathematics and economics, 44(2):182–198.
  * Amari et al., [2003] Amari, S. I., Nakahara, H., Wu, S., and Sakai, Y. (2003).  Synchronous firing and higher-order interactions in neuron pool.  Neural Comput, 15.
  * Averbeck et al., [2006] Averbeck, B. B., Latham, P. E., and Pouget, A. (2006).  Neural correlations, population coding and computation.  Nat Rev Neurosci, 7(5):358–366.
  * Bar-Gad et al., [2001] Bar-Gad, I., Ritov, Y., Vaadia, E., and Bergman, H. (2001).  Failure in identification of overlapping spikes from multiple neuron activity causes artificial correlations.  J Neurosci Methods, 107(1-2):1–13.
  * Battistin et al., [2015] Battistin, C., Hertz, J., Tyrcha, J., and Roudi, Y. (2015).  Belief propagation and replicas for inference and learning in a kinetic ising model with hidden spins.  J Stat Mech Theory Exp, 2015(5):P05021.
  * Benjamin et al., [2018] Benjamin, A. S., Fernandes, H. L., Tomlinson, T., Ramkumar, P., VerSteeg, C., Chowdhury, R. H., Miller, L. E., and Kording, K. P. (2018).  Modern machine learning as a benchmark for fitting neural responses.  Front Comput Neurosc, 12:56.
  * Berkes et al., [2008] Berkes, P., Wood, F., and Pillow, J. (2008).  Characterizing neural dependencies with copula models.  In Advances in Neural Information Processing Systems, volume 21, pages 129–136.
  * Brinkman et al., [2018] Brinkman, B. A., Rieke, F., Shea-Brown, E., and Buice, M. A. (2018).  Predicting how and when hidden neurons skew measured synaptic interactions.  PLoS Comput Biol, 14(10):e1006490.
  * Buccino et al., [2020] Buccino, A. P., Hurwitz, C. L., Garcia, S., Magland, J., Siegle, J. H., Hurwitz, R., and Hennig, M. H. (2020).  SpikeInterface, a unified framework for spike sorting.  eLife, 9:e61834.  * This paper describes an open source framework to create fully reproducible analysis pipelines for extracellular recordings. It allows runnning all major current spike sorting algorithms on the same data sets, includes common methods for pre- and postprocessing and automated curation, and a systematic comparison between sorter outputs. Use of this toolkit can both help to avoid failure modes of common spike sorters, and aid more reproducible analysis.
  * Buesing et al., [2013] Buesing, L., Macke, J. H., and Sahani, M. (2013).  Spectral learning of linear dynamics from generalised-linear observations with application to neural population data.  In Advances in Neural Information Processing Systems 25: 26th Conference on Neural Information Processing Systems (NIPS 2012), pages 1691–1699.
  * Byron et al., [2009] Byron, M. Y., Cunningham, J. P., Santhanam, G., Ryu, S. I., Shenoy, K. V., and Sahani, M. (2009).  Gaussian-process factor analysis for low-dimensional single-trial analysis of neural population activity.  In Advances in Neural Information Processing Systems, pages 1881–1888.
  * Carlson and Carin, [2019] Carlson, D. and Carin, L. (2019).  Continuing progress of spike sorting in the era of big data.  Curr Opin Neurobiol, 55:90–96.
  * Churchland et al., [2012] Churchland, M. M., Cunningham, J. P., Kaufman, M. T., Foster, J. D., Nuyujukian, P., Ryu, S. I., and Shenoy, K. V. (2012).  Neural population dynamics during reaching.  Nature, 487(7405):51–56.
  * Dunn et al., [2015] Dunn, B., Mørreaunet, M., and Roudi, Y. (2015).  Correlations and Functional Connections in a Population of Grid Cells.  PLoS Comput Biol, 11(2):e1004052.
  * Durstewitz, [2017] Durstewitz, D. (2017).  A state space approach for piecewise-linear recurrent neural networks for identifying computational dynamics from neural measurements.  PLoS computational biology, 13(6):e1005542.
  * Ecker et al., [2014] Ecker, A. S., Berens, P., Cotton, R. J., Subramaniyan, M., Denfield, G. H., Cadwell, C. R., Smirnakis, S. M., Bethge, M., and Tolias, A. S. (2014).  State dependence of noise correlations in macaque primary visual cortex.  Neuron, 82(1):235–248.
  * Ecker et al., [2010] Ecker, A. S., Berens, P., Keliris, G. A., Bethge, M., Logothetis, N. K., and Tolias, A. S. (2010).  Decorrelated neuronal firing in cortical microcircuits.  science, 327(5965):584–587.
  * Elsayed and Cunningham, [2017] Elsayed, G. F. and Cunningham, J. P. (2017).  Structure in neural population recordings: an expected byproduct of simpler phenomena?  Nat Neurosci, 20(9):1310.
  * Evans et al., [2019] Evans, M. H., Petersen, R. S., and Humphries, M. D. (2019).  On the use of calcium deconvolution algorithms in practical contexts.  bioRxiv, page 871137.
  * Franke et al., [2016] Franke, F., Fiscella, M., Sevelev, M., Roska, B., Hierlemann, A., and da Silveira, R. A. (2016).  Structures of neural correlation and how they favor coding.  Neuron, 89.
  * Gallego et al., [2020] Gallego, J. A., Perich, M. G., Chowdhury, R. H., Solla, S. A., and Miller, L. E. (2020).  Long-term stability of cortical population dynamics underlying consistent behavior.  Nat Neurosci, 23(2):260–270.  ** The paper demonstrates that latent dynamics are a stable neural correlate for the execution of learned behaviors. They show that even across years of recordings, latent dynamics in the primary motor cortex during reach remain surprisingly consistent despite significant changes in single neuron activities.
  * Gallego et al., [2017] Gallego, J. A., Perich, M. G., Miller, L. E., and Solla, S. A. (2017).  Neural manifolds for the control of movement.  Neuron, 94(5):978–984.
  * Gao et al., [2017] Gao, P., Trautmann, E., Yu, B., Santhanam, G., Ryu, S., Shenoy, K., and Ganguli, S. (2017).  A theory of multineuronal dimensionality, dynamics and measurement.  BioRxiv.
  * Gao et al., [2016] Gao, Y., Archer, E., Paninski, L., and Cunningham, J. P. (2016).  Linear dynamical neural population models through nonlinear embeddings.  arXiv preprint arXiv:1605.08454.
  * Gardella et al., [2016] Gardella, C., Marre, O., and Mora, T. (2016).  A tractable method for describing complex couplings between neurons and population rate.  eNeuro, 3.
  * Granot-Atedgi et al., [2013] Granot-Atedgi, E., Tkačik, G., Segev, R., and Schneidman, E. (2013).  Stimulus-dependent maximum entropy models of neural population codes.  PLoS Comput Biol, 9(3):e1002922.
  * Haghighi and Markazi, [2017] Haghighi, H. S. and Markazi, A. (2017).  A new description of epileptic seizures based on dynamic analysis of a thalamocortical model.  Sci Rep, 7(1):1–10.
  * [28] Hall, T. M., de Carvalho, F., and Jackson, A. (2014a).  A common structure underlies low-frequency cortical dynamics in movement, sleep, and sedation.  Neuron, 83(5):1185–1199.
  * [29] Hall, T. M., Nazarpour, K., and Jackson, A. (2014b).  Real-time estimation and biofeedback of single-neuron firing rates using local field potentials.  Nat Commun, 5(1):1–12.
  * Hennig et al., [2019] Hennig, M. H., Hurwitz, C., and Sorbaro, M. (2019).  Scaling Spike Detection and Sorting for Next-Generation Electrophysiology.  In Chiappalone, M., Pasquale, V., and Frega, M., editors, In Vitro Neuronal Networks: From Culturing Methods to Neuro-Technological Applications, Advances in Neurobiology, pages 171–184. Springer International Publishing, Cham.
  * Herzog et al., [2018] Herzog, R., Escobar, M.-J., Cofre, R., Palacios, A. G., and Cessac, B. (2018).  Dimensionality Reduction on Spatio-Temporal Maximum Entropy Models of Spiking Networks.  bioRxiv, page 278606.
  * Huang et al., [2021] Huang, L., Ledochowitsch, P., Knoblich, U., Lecoq, J., Murphy, G. J., Reid, R. C., de Vries, S. E., Koch, C., Zeng, H., Buice, M. A., et al. (2021).  Relationship between simultaneously recorded spiking activity and fluorescence signal in gcamp6 transgenic mice.  Elife, 10:e51675.
  * Ince et al., [2017] Ince, R. A., Giordano, B. L., Kayser, C., Rousselet, G. A., Gross, J., and Schyns, P. G. (2017).  A statistical framework for neuroimaging data analysis based on mutual information estimated via a gaussian copula.  Hum Brain Mapp, 38(3):1541–1573.
  * Jenison and Reale, [2004] Jenison, R. L. and Reale, R. A. (2004).  The shape of neural dependence.  Neural Comput, 16.
  * Jun et al., [2017] Jun, J. J., Steinmetz, N. A., Siegle, J. H., Denman, D. J., Bauza, M., Barbarits, B., Lee, A. K., Anastassiou, C. A., Andrei, A., Aydın, Ç., Barbic, M., Blanche, T. J., Bonin, V., Couto, J., Dutta, B., Gratiy, S. L., Gutnisky, D. A., Häusser, M., Karsh, B., Ledochowitsch, P., Lopez, C. M., Mitelut, C., Musa, S., Okun, M., Pachitariu, M., Putzeys, J., Rich, P. D., Rossant, C., Sun, W.-l., Svoboda, K., Carandini, M., Harris, K. D., Koch, C., O’Keefe, J., and Harris, T. D. (2017).  Fully integrated silicon probes for high-density recording of neural activity.  Nature, 551(7679):232–236.
  * Kobak et al., [2016] Kobak, D., Brendel, W., Constantinidis, C., Feierstein, C. E., Kepecs, A., Mainen, Z. F., Qi, X.-L., Romo, R., Uchida, N., and Machens, C. K. (2016).  Demixed principal component analysis of neural population data.  Elife, 5:e10989.
  * Kudryashova et al., [2020] Kudryashova, N., Amvrosiadis, T., Dupuy, N., Rochefort, N., and Onken, A. (2020).  Parametric Copula-GP model for analyzing multidimensional neuronal and behavioral relationships.  arXiv preprint arXiv:2008.01007.  * The paper proposes a scalable framework for modeling neuronal population activity and external covariates (e.g. behavior). The advantage of the method is the mutual information estimation: it provides more accurate estimates than non-copula methods, especially for large datasets. The GP parametrization allows for conditioning on continuous variables, e.g. on time or position, making the method suitable for studying dynamic trajectories.
  * Kulkarni and Paninski, [2007] Kulkarni, J. E. and Paninski, L. (2007).  Common-input models for multiple neural spike-train data.  Network: Computation in Neural Systems, 18(4):375–407.
  * Lawhern et al., [2010] Lawhern, V., Wu, W., Hatsopoulos, N., and Paninski, L. (2010).  Population decoding of motor cortical activity using a generalized linear model with hidden states.  Journal of neuroscience methods, 189(2):267–280.
  * Lawrence and Hyvärinen, [2005] Lawrence, N. and Hyvärinen, A. (2005).  Probabilistic non-linear principal component analysis with gaussian process latent variable models.  J Mach Learn Res, 6(11).
  * Lefebvre et al., [2016] Lefebvre, B., Yger, P., and Marre, O. (2016).  Recent progress in multi-electrode spike sorting methods.  Journal of Physiology-Paris, 110(4, Part A):327–335.
  * Linderman et al., [2017] Linderman, S., Johnson, M., Miller, A., Adams, R., Blei, D., and Paninski, L. (2017).  Bayesian learning and inference in recurrent switching linear dynamical systems.  In Artificial Intelligence and Statistics, pages 914–922. PMLR.
  * Macke et al., [2009] Macke, J. H., Berens, P., Ecker, A. S., Tolias, A. S., and Bethge, M. (2009).  Generating spike trains with specified correlation coefficients.  Neural Comput, 21.
  * Macke et al., [2012] Macke, J. H., Buesing, L., Cunningham, J. P., Yu, B. M., Shenoy, K. V., and Sahani, M. (2012).  Empirical models of spiking in neural populations.  In Advances in Neural Information Processing Systems 24: 25th conference on Neural Information Processing Systems (NIPS 2011), pages 1350–1358.
  * Magland et al., [2020] Magland, J., Jun, J. J., Lovero, E., Morley, A. J., Hurwitz, C. L., Buccino, A. P., Garcia, S., and Barnett, A. H. (2020).  SpikeForest, reproducible web-facing ground-truth validation of automated neural spike sorters.  eLife, 9:e55167.
  * Mahuas et al., [2020] Mahuas, G., Isacchini, G., Marre, O., Ferrari, U., and Mora, T. (2020).  A new inference approach for training shallow and deep generalized linear models of noisy interacting neurons.  arXiv preprint arXiv:2006.06497.  ** The paper proposes a method for training GLMs with complex correlated stimuli. It disentangles stimulus correlations from couplings, which improves generalization and prevents self-excitation in the network dynamics.
  * Mante et al., [2013] Mante, V., Sussillo, D., Shenoy, K. V., and Newsome, W. T. (2013).  Context-dependent computation by recurrent dynamics in prefrontal cortex.  Nature, 503(7474):78–84.
  * Mongillo et al., [2018] Mongillo, G., Rumpel, S., and Loewenstein, Y. (2018).  Inhibitory connectivity defines the realm of excitatory plasticity.  Nat Neurosci, 21(10):1463–1470.
  * Nawrot et al., [1999] Nawrot, M., Aertsen, A., and Rotter, S. (1999).  Single-trial estimation of neuronal firing rates: from single-neuron spike trains to population activity.  Journal of neuroscience methods, 94(1):81–92.
  * Nemati et al., [2014] Nemati, S., Linderman, S. W., Chen, Z., et al. (2014).  A probabilistic modeling approach for uncovering neural population rotational dynamics.  In Computational and Systems Neuroscience (Cosyne).
  * O’Donnell et al., [2017] O’Donnell, C., alves, J. T. G., Whiteley, N., Portera-Cailliau, C., and Sejnowski, T. J. (2017).  The Population Tracking Model: A Simple, Scalable Statistical Model for Neural Population Data.  Neural Comput, 29(1):50–93.
  * Ohiorhenuan et al., [2010] Ohiorhenuan, I. E., Mechler, F., Purpura, K. P., Schmid, A. M., Hu, Q., and Victor, J. D. (2010).  Sparse coding and high-order correlations in fine-scale cortical networks.  Nature, pages 1–6.
  * Okun et al., [2015] Okun, M., Steinmetz, N. A., Cossell, L., Iacaruso, M. F., Ko, H., Barthó, P., Moore, T., Hofer, S. B., Mrsic-Flogel, T. D., Carandini, M., et al. (2015).  Diverse coupling of neurons to populations in sensory cortex.  Nature, 521(7553):511–515.
  * Onken and Panzeri, [2016] Onken, A. and Panzeri, S. (2016).  Mixed vine copulas as joint models of spike counts and local field potentials.  In Advances in Neural Information Processing Systems.
  * Panas et al., [2015] Panas, D., Amin, H., Maccione, A., Muthmann, O., van Rossum, M., Berdondini, L., and Hennig, M. H. (2015).  Sloppiness in spontaneously active neuronal networks.  J Neurosci, 35(22):8480–8492.
  * Pandarinath et al., [2018] Pandarinath, C., O’Shea, D. J., Collins, J., Jozefowicz, R., Stavisky, S. D., Kao, J. C., Trautmann, E. M., Kaufman, M. T., Ryu, S. I., Hochberg, L. R., Henderson, J. M., Shenoy, K. V., Abbott, L. F., and Sussillo, D. (2018).  Inferring single-trial neural population dynamics using sequential auto-encoders.  Nat Methods, 15(10):805–815.  ** This paper introduces a novel RNN-based, state-space model termed latent factor analysis via dynamical systems (LFADS). LFADS utilizes a variational autoencoder to amortize the cost of inference over the dynamics making it more practical for large-scale neural recordings.
  * Petreska et al., [2011] Petreska, B., Yu, B. M., Cunningham, J. P., Santhanam, G., Ryu, S. I., Shenoy, K. V., and Sahani, M. (2011).  Dynamical segmentation of single trials from population neural data.  In Proceedings of the 24th International Conference on Neural Information Processing Systems, NIPS’11, pages 756–764, Red Hook, NY, USA. Curran Associates Inc.
  * Pfau et al., [2013] Pfau, D., Pnevmatikakis, E. A., and Paninski, L. (2013).  Robust learning of low-dimensional dynamics from large neural ensembles.  In Advances in Neural Information Processing Systems.
  * Pillow et al., [2011] Pillow, J. W., Ahmadian, Y., and Paninski, L. (2011).  Model-based decoding, information estimation, and change-point detection techniques for multineuron spike trains.  Neural Comput, 23(1):1–45.
  * Pillow et al., [2005] Pillow, J. W., Paninski, L., Uzzell, V. J., Simoncelli, E. P., and Chichilnisky, E. J. (2005).  Prediction and decoding of retinal ganglion cell responses with a probabilistic spiking model.  J Neurosci, 25(47):11003–13.
  * Pillow et al., [2008] Pillow, J. W., Shlens, J., Paninski, L., Sher, A., Litke, A. M., Chichilnisky, E. J., and Simoncelli, E. P. (2008).  Spatio-temporal correlations and visual signalling in a complete neuronal population.  Nature, 454(7207):995–999.
  * Ponce-Alvarez et al., [2020] Ponce-Alvarez, A., Mochol, G., Hermoso-Mendizabal, A., de la Rocha, J., and Deco, G. (2020).  Cortical state transitions and stimulus response evolve along stiff and sloppy parameter dimensions, respectively.  eLife, 9:e53268.  * Using maximum entropy models, this paper confirms that population activity in cortical networks is constrained only by a small fraction of highly sensitive coupling parameters. In in vivo recordings from auditory cortex, the authors show that transitions between synchronised and de-synchronised states reflect changes in sensitive parameters. In contrast, auditory stimuli affect primarily the insensitive couplings, hence stimulation has only a minor effect on the overall statistics of cortical activity.
  * Prince et al., [2021] Prince, L. Y., Bakhtiari, S., Gillon, C. J., and Richards, B. A. (2021).  Parallel inference of hierarchical latent dynamics in two-photon calcium imaging of neuronal populations.  bioRxiv.
  * Rabinovich et al., [2008] Rabinovich, M. I., Huerta, R., Varona, P., and Afraimovich, V. S. (2008).  Transient cognitive dynamics, metastability, and decision making.  PLoS Comput Biol, 4(5):e1000072.
  * Raducanu et al., [2017] Raducanu, B. C., Yazicioglu, R. F., Lopez, C. M., Ballini, M., Putzeys, J., Wang, S., Andrei, A., Rochus, V., Welkenhuysen, M., Helleputte, N. v., et al. (2017).  Time multiplexed active neural probe with 1356 parallel recording sites.  Sens, 17(10):2388.
  * Ramanathan et al., [2018] Ramanathan, D. S., Guo, L., Gulati, T., Davidson, G., Hishinuma, A. K., Won, S.-J., Knight, R. T., Chang, E. F., Swanson, R. A., and Ganguly, K. (2018).  Low-frequency cortical activity is a neuromodulatory target that tracks recovery after stroke.  Nat Med, 24(8):1257–1267.
  * René et al., [2020] René, A., Longtin, A., and Macke, J. H. (2020).  Inference of a mesoscopic population model from population spike trains.  Neural Comput, 32(8):1448–1498.
  * Rey et al., [2015] Rey, H. G., Pedreira, C., and Quian Quiroga, R. (2015).  Past, present and future of spike sorting techniques.  Brain Res Bull, 119:106–117.
  * Roweis and Ghahramani, [1999] Roweis, S. and Ghahramani, Z. (1999).  A unifying review of linear gaussian models.  Neural Comput, 11(2):305–345.
  * Rule et al., [2019] Rule, M. E., Schnoerr, D., Hennig, M. H., and Sanguinetti, G. (2019).  Neural field models for latent state inference: Application to large-scale neuronal recordings.  PLoS Comput Biol, 15(11):e1007442.
  * Runyan et al., [2017] Runyan, C. A., Piasini, E., Panzeri, S., and Harvey, C. D. (2017).  Distinct timescales of population coding across cortex.  Nature, 548.
  * Russo et al., [2020] Russo, A. A., Khajeh, R., Bittner, S. R., Perkins, S. M., Cunningham, J. P., Abbott, L. F., and Churchland, M. M. (2020).  Neural trajectories in the supplementary motor area and motor cortex exhibit distinct geometries, compatible with different classes of computation.  Neuron, 107(4):745–758.  * The paper explores a new non-linear pattern of activity in supplementary motor area during multi-cycle rhythmic movement: helical population-trajectory geometry. The separation between single-cycles in these trajectories allows for counting the number of revolutions.
  * Rutten et al., [2020] Rutten, V., Bernacchia, A., Sahani, M., and Hennequin, G. (2020).  Non-reversible gaussian processes for identifying latent dynamical structure in neural data.  In Advances in Neural Information Processing Systems, volume 33.  ** The first paper to consider time-irreversibility in GPFA models by constructing a new kernel function. This method, termed GPFADS, assigns lower probabilities to the time-reversed variants of the true trajectories. Unlike traditional GPFA, it can disentangle trajectories in the arm reach dataset [13].
  * Sadtler et al., [2014] Sadtler, P. T., Quick, K. M., Golub, M. D., Chase, S. M., Ryu, S. I., Tyler-Kabara, E. C., Byron, M. Y., and Batista, A. P. (2014).  Neural constraints on learning.  Nature, 512(7515):423–426.
  * Saha et al., [2013] Saha, D., Leong, K., Li, C., Peterson, S., Siegel, G., and Raman, B. (2013).  A spatiotemporal coding mechanism for background-invariant odor recognition.  Nat. Neurosci., 16(12):1830–1839.
  * Sani et al., [2020] Sani, O. G., Abbaspourazad, H., Wong, Y. T., Pesaran, B., and Shanechi, M. M. (2020).  Modeling behaviorally relevant neural dynamics enabled by preferential subspace identification.  Nat Neurosci, pages 1–10.  ** This paper introduces a novel linear state-space model (PSID) that jointly models the neural activity and the recorded behavioural variables as generated by the behaviourally relevant and irrelevant latent dynamics. The authors show that PSID discovers significantly lower-dimensional latent dynamics than other models and is the only one that can capture reverse directional rotational dynamics in the primary motor cortex when the arm is returning versus when the arm is reaching.
  * Santhanam et al., [2009] Santhanam, G., Yu, B. M., Gilja, V., Ryu, S. I., Afshar, A., Sahani, M., and Shenoy, K. V. (2009).  Factor-analysis methods for higher-performance neural prostheses.  J Neurophysiol, 102(2):1315–1330.
  * Savin and Tkačik, [2017] Savin, C. and Tkačik, G. (2017).  Maximum entropy models as a tool for building precise neural controls.  Curr Opin Neurobiol, 46:120–126.
  * Saxena and Cunningham, [2019] Saxena, S. and Cunningham, J. P. (2019).  Towards the neural population doctrine.  Curr Opin Neurobiol, 55:103–111.
  * Schneidman et al., [2006] Schneidman, E., Berry, M. J., Segev, R., and Bialek, W. (2006).  Weak pairwise correlations imply strongly correlated network states in a neural population.  Nature, 440.
  * She and Wu, [2020] She, Q. and Wu, A. (2020).  Neural dynamics discovery via gaussian process recurrent neural networks.  In Uncertainty in Artificial Intelligence, pages 454–464. PMLR.  * This paper extends RNN-based dynamical modelling approaches with a Gaussian Process observation model. They show that their more expressive observation model is able to extract insightful low-dimensional dynamics and better fit neural datasets.
  * Shenoy et al., [2013] Shenoy, K. V., Sahani, M., and Churchland, M. M. (2013).  Cortical control of arm movements: a dynamical systems perspective.  Annu Rev Neurosci, 36.
  * Shlens et al., [2006] Shlens, J., Field, G. D., Gauthier, J. L., Grivich, M. I., Petrusca, D., Sher, A., Litke, A. M., and Chichilnisky, E. J. (2006).  The structure of multi-neuron firing patterns in primate retina.  J Neurosci, 26(32):8254–66.
  * Si et al., [2019] Si, G., Kanwal, J. K., Hu, Y., Tabone, C. J., Baron, J., Berck, M., Vignoud, G., and Samuel, A. D. (2019).  Structured odorant response patterns across a complete olfactory receptor neuron population.  Neuron, 101(5):950–962.
  * Smith and Brown, [2003] Smith, A. C. and Brown, E. N. (2003).  Estimating a state-space model from point process observations.  Neural Comput, 15(5):965–991.
  * Spreemann et al., [2018] Spreemann, G., Dunn, B., Botnan, M. B., and Baas, N. A. (2018).  Using persistent homology to reveal hidden covariates in systems governed by the kinetic ising model.  Phys Rev E, 97(3):032313.
  * Stringer et al., [2019] Stringer, C., Pachitariu, M., Steinmetz, N., Carandini, M., and Harris, K. D. (2019).  High-dimensional geometry of population responses in visual cortex.  Nature, 571(7765):361–365.
  * Tkačik et al., [2013] Tkačik, G., Marre, O., Mora, T., Amodei, D., Berry, M. J., and Bialek, W. (2013).  The simplest maximum entropy model for collective behavior in a neural network.  J Stat Mech Theory Exp, 2013.
  * Trippe et al., [2019] Trippe, B., Huggins, J., Agrawal, R., and Broderick, T. (2019).  Lr-glm: High-dimensional Bayesian inference using low-rank data approximations.  In International Conference on Machine Learning, pages 6315–6324. PMLR.
  * Truccolo et al., [2005] Truccolo, W., Eden, U. T., Fellows, M. R., Donoghue, J. P., and Brown, E. N. (2005).  A point process framework for relating neural spiking activity to spiking history, neural ensemble, and extrinsic covariate effects.  J Neurophysiol, 93(2):1074–1089.
  * Ventura, [2009] Ventura, V. (2009).  Traditional waveform based spike sorting yields biased rate code estimates.  Proc Natl Acad Sci, 106(17):6921–6926.
  * Ventura and Gerkin, [2012] Ventura, V. and Gerkin, R. C. (2012).  Accurately estimating neuronal correlation requires a new spike-sorting paradigm.  Proc Natl Acad Sci, 109(19):7230–7235.
  * Whiteway and Butts, [2019] Whiteway, M. R. and Butts, D. A. (2019).  The quest for interpretable models of neural population activity.  Curr Opin Neurobiol, 58:86–93.
  * Wu et al., [2017] Wu, A., Roy, N. A., Keeley, S., and Pillow, J. W. (2017).  Gaussian process based nonlinear latent structure discovery in multivariate spike train data.  Adv Neur In, 30:3496.
  * Zhao and Park, [2016] Zhao, Y. and Park, I. M. (2016).  Interpretable nonlinear dynamic modeling of neural trajectories.  arXiv preprint arXiv:1608.06546.
  * Zhao and Park, [2017] Zhao, Y. and Park, I. M. (2017).  Variational latent gaussian process for recovering single-trial dynamics from population spike trains.  Neural Comput, 29(5):1293–1316.
  * Zhou and Wei, [2020] Zhou, D. and Wei, X.-X. (2020).  Learning identifiable and interpretable latent models of high-dimensional neural activity using pi-vae.  arXiv preprint arXiv:2011.04798.
  * Zoltowski et al., [2020] Zoltowski, D., Pillow, J., and Linderman, S. (2020).  A general recurrent state space framework for modeling neural dynamics during decision-making.  In International Conference on Machine Learning, pages 11680–11691. PMLR.
  * Zoltowski and Pillow, [2018] Zoltowski, D. M. and Pillow, J. W. (2018).  Scaling the poisson glm to massive neural datasets through polynomial approximations.  In Advances in Neural Information Processing Systems, volume 31, page 3517. NIH Public Access.
  * Zylberberg et al., [2016] Zylberberg, J., Cafaro, J., Turner, M. H., Shea-Brown, E., and Rieke, F. (2016).  Direction-selective circuits shape noise to ensure a precise population code.  Neuron, 89.