Sometimes the spike-triggered average is just as good as Poisson GLMs for spike-train analysis

(hopefully no major mistakes in this one; get PDF here)

The Poisson GLM for spiking data

Generalized Linear Models (GLMs) are similar to linear regression, but account for nonlinearities and non-uniform noise in the observations. In neuroscience, it is common to predict a sequence of spikes $Y=\{y_1,..,y_T\}$, $y_i\in \{0,1\}$, from a series of observations $X=\{{\mathbf{x}}_1,..,{\mathbf{x}}_T\}$, using a Poisson GLM:

$$\begin{array}{rl}{y}_i& \sim \mathrm{Poisson}({\lambda }_i\cdot \mathrm{\Delta }t)\\ {\lambda }_i& =\exp ({\mathbf{a}}^{\mathrm{\top }}{\mathbf{x}}_i+b)\end{array}$$

These models are fit by minimizing the negative log-likelihood of the observations, given the vector of regression weights $\mathbf{a}$ and mean offset parameter $b$:

$$\underset{a,b}{\mathrm{argmin}}⟨\lambda -y\ln \lambda ⟩=\underset{\mathbf{a},b}{\mathrm{argmin}}⟨e^{{\mathbf{a}}^{\mathrm{\top }}{\mathbf{x}}_i+b}-y[{\mathbf{a}}^{\mathrm{\top }}{\mathbf{x}}_i+b]⟩$$

Such models are usually optimized via gradient descent, or iterated re-weighted least squares. The gradient

$$\begin{array}{rl}{\partial }_{\mathbf{b}}(...)& =⟨\lambda -y⟩\\ {\partial }_{\mathbf{a}}(...)& =⟨(\lambda -y){\mathbf{x}}^{\mathrm{\top }}⟩\end{array}$$

This is convex, and the optimum occurs at an inflection point where the gradient is zero, i.e.:

$$\begin{array}{rl}⟨\lambda ⟩& =⟨y⟩\\ ⟨\lambda {\mathbf{x}}^{\mathrm{\top }}⟩& =⟨y{\mathbf{x}}^{\mathrm{\top }}⟩\end{array}$$

For simplicity, assume that $\mathbf{x}$ are a collection of uncorrelated zero-mean Gaussian variables. (If $\mathbf{x}$ has a multivariate Gaussian distribution, we an always transform into coordinates so that this is true by subtracting the mean and whitening the data.)

Linear models are often fine

What happens if we expand our likelihood at first order, and solve for a zero of the gradient? We start by expanding around the point $b_0=\ln ⟨y⟩$:

$$\begin{array}{rl}\lambda & =e^{{\mathbf{a}}^{\mathrm{\top }}\mathbf{x}+b}=e^{{\mathbf{a}}^{\mathrm{\top }}\mathbf{x}+b-b_0}{e}^{b_0}\\ & \approx [1+{\mathbf{a}}^{\mathrm{\top }}\mathbf{x}+b-b_0]e^{b_0}\\ & =[1+{\mathbf{a}}^{\mathrm{\top }}\mathbf{x}+b-b_0]⟨y⟩\end{array}$$

Substituting this approximation in to the identity for the optimal coefficients:

$$\begin{array}{rl}& ⟨[1+{\mathbf{a}}^{\mathrm{\top }}\mathbf{x}+b-b_0]⟩=1\\ & ⟨[1+{\mathbf{a}}^{\mathrm{\top }}\mathbf{x}+b-b_0]{\mathbf{x}}^{\mathrm{\top }}⟩=⟨y{\mathbf{x}}^{\mathrm{\top }}⟩/⟨y⟩\end{array}$$

If we have chosen $\mathrm{\Delta }t\ll 1$, so that $y$ is indeed binary, then the right hand side of the second equation is simply the spike-triggered average. Denote this as ${\mathbf{x}}_1^{\mathrm{\top }}=⟨y{\mathbf{x}}^{\mathrm{\top }}⟩/⟨y⟩$. Solving for $(\mathbf{a},b)$ in the above yields:

$$\begin{array}{rl}b& =b_0\\ \mathbf{a}& ={\mathbf{x}}_1\end{array}$$

Which is to say: to a first approximation, the coefficients $\mathbf{a}$ of the Poisson GLM are simply the spike triggered average (of the mean-centered and whitened $\mathbf{x}$).

Intuitively, this means if you're just trying to demonstrate some basic correlation between a variable and neuronal spiking, linear methods like the spike-triggered-average are fine. Moving to a GLM might give you more accuracy, since it has a better model of the observation process, but it often doesn't tell you much more, qualitatively.

For Gaussian covariates, GLMs capture how noise and nonlinearity interact to influence firing rate

Let's see if we can get further intuition if $\mathbf{x}$ is Gaussian, beyond first order. Work in normalized firing rate by defining $\overline{\lambda }=\lambda /⟨y⟩$. The optimum can then be written as:

$$\begin{array}{rl}⟨\overline{\lambda }⟩& =1\\ ⟨\overline{\lambda }\mathbf{x}⟩& ={\mathbf{x}}_1\end{array}$$

If $\mathbf{x}$ is Gaussian, the expectation $⟨\overline{\lambda }{\mathbf{x}}^{\mathrm{\top }}⟩$ has a closed form. Defining $\overline{b}=b-b_0=b-\ln ⟨y⟩$:

$$⟨\overline{\lambda }\mathbf{x}⟩={\int }_{d\mathbf{x}}\mathbf{x}{e}^{{\mathbf{a}}^{\mathrm{\top }}\mathbf{x}+\overline{b}}(2\pi {)}^{-D/2}|{\mathrm{\Sigma }}_{\mathbf{x}}{|}^{-1/2}\exp (-\frac{1}{2}(\mathbf{x}-{\mu }_{\mathbf{x}}{)}^{\mathrm{\top }}{\mathrm{\Sigma }}_{\mathbf{x}}^{-1}(\mathbf{x}-{\mu }_{\mathbf{x}})),$$

where $D$ is the dimensionality of $\mathbf{x}$. Working in coordinates where ${\mathrm{\Sigma }}_{\mathbf{x}}=I$ and ${\mu }_{\mathbf{x}}=0$, this simplifies to:

$$⟨\overline{\lambda }\mathbf{x}⟩=(2\pi {)}^{-D/2}{\int }_{d\mathbf{x}}\mathbf{x}{e}^{{\mathbf{a}}^{\mathrm{\top }}\mathbf{x}+\overline{b}-\frac{1}{2}{\mathbf{x}}^{\mathrm{\top }}\mathbf{x}}.$$

This integral can be solved by completing the square:

$$\begin{array}{rl}{\mathbf{a}}^{\mathrm{\top }}\mathbf{x}+\overline{b}-\frac{1}{2}{\mathbf{x}}^{\mathrm{\top }}\mathbf{x}& =-\frac{1}{2}[-2{\mathbf{a}}^{\mathrm{\top }}\mathbf{x}-2\overline{b}+{\mathbf{x}}^{\mathrm{\top }}\mathbf{x}]\\ & =-\frac{1}{2}[{\mathbf{x}}^{\mathrm{\top }}\mathbf{x}-2{\mathbf{a}}^{\mathrm{\top }}\mathbf{x}+{\mathbf{a}}^{\mathrm{\top }}\mathbf{a}-{\mathbf{a}}^{\mathrm{\top }}\mathbf{a}-2\overline{b}]\\ & =-\frac{1}{2}[(\mathbf{x}-\mathbf{a}{)}^{\mathrm{\top }}(\mathbf{x}-\mathbf{a})-{\mathbf{a}}^{\mathrm{\top }}\mathbf{a}-2\overline{b}]\\ & =-\frac{1}{2}[(\mathbf{x}-\mathbf{a}{)}^{\mathrm{\top }}(\mathbf{x}-\mathbf{a})]+[\frac{1}{2}{\mathbf{a}}^{\mathrm{\top }}\mathbf{a}+\overline{b}]\end{array}$$$$⟨\overline{\lambda }{\mathbf{x}}^{\mathrm{\top }}⟩=(2\pi {)}^{-D/2}{\int }_{d\mathbf{x}}\mathbf{x}{e}^{-\frac{1}{2}[(\mathbf{x}-\mathbf{a}{)}^{\mathrm{\top }}(\mathbf{x}-\mathbf{a})]+[\frac{1}{2}{\mathbf{a}}^{\mathrm{\top }}\mathbf{a}+\overline{b}]}.$$

This is a Gaussian integral. It evaluates to:

$$⟨\overline{\lambda }\mathbf{x}⟩=\mathbf{a}{e}^{\frac{1}{2}{\mathbf{a}}^{\mathrm{\top }}\mathbf{a}+\overline{b}}.$$

At the optimum, we therefore have the condition that the optimal decoding weights $\mathbf{a}$ should be proportional to the spike-triggered average ${\mathbf{x}}_1$:

$$\begin{array}{rl}{\mathbf{x}}_1& =\mathbf{a}{e}^{\frac{1}{2}{\mathbf{a}}^{\mathrm{\top }}\mathbf{a}+\overline{b}}\\ & =\mathbf{a}{e}^{\frac{1}{2}\Vert \mathbf{a}{\Vert }^2+\overline{b}}\\ \Rightarrow \mathbf{a}& \propto {\mathbf{x}}_1\end{array}$$

Solving similarly for $⟨\overline{\lambda }⟩$ gives $⟨\overline{\lambda }⟩=\exp (\frac{1}{2}{\mathbf{a}}^{\mathrm{\top }}\mathbf{a}+\overline{b})$, which seems to imply that the Poisson GLM has a closed-form solution in the limit where is $y$ sparse and $\mathbf{x}$ is Gaussian:

$$\begin{array}{rl}\mathbf{a}& ={\mathbf{x}}_1\\ b& =\ln ⟨y⟩-\frac{1}{2}\Vert {\mathbf{x}}_1{\Vert }^2\end{array}$$

This further supports the idea that the spike-triggered average is an fast, adequate approach to statistics on spike trains. In some cases, it may be unnecessary to train a GLM.

Of course, Poisson GLMs can be fit almost as quickly as linear regressions or spike-triggered averages, using iteratively re-weighted least squares. They can also gracefully and automatically handle cases where the covariates $\mathbf{x}$ are not jointly Gaussian. So perhaps the question should be why not use a GLM.