These notes provide a function for calculating the AUC for a Poisson regression that works when the spike counts are given only up to some (unknown) normalization constant. Skip to the end for the Python code.
Poisson regression and point processes
Neuroscientists use Poisson regression (in the continuous limit: Poisson point process ) to assess the correlation between neuronal spiking and other experimental covariates. These models are well suited to the scenario in which spikes are binary all-or-nothing events, recorded based on electrical potentials.
The area under the receiver operating characteristic curve
It is common to assess the performance of a Poisson GLM using the Area Under the receiver-operator-characteristic Curve (AUC) .
Definition of the AUC: The AUC is defined as the probability that a randomly-sampled (true-negative, true-positive) pair will be correctly ranked by the predicted firing-rates. We can calculate it by integrating over all such pairs. We will work in continuous time $t\in[0,T]$ for an experiment with duration $T$.
Denote the sets of true positives and true negatives as $t_p\in\mathcal T_p$ and $t_n\in\mathcal T_n$, respectively. Consider a generic model $\lambda(t)$ which predicts the spiking probability at time $t$. Generically, the AUC can be calculated by integrating over all pairs $(t_p,t_n)$ of positive and negative examples:
\begin{equation} \begin{aligned} \text{AUC}(\lambda,\mathcal T_p,\mathcal T_n) &= \Pr\big\{ \lambda(t_n) < \lambda(t_p) \text{ for }t_n\in\mathcal T_n \text{ and }t_p\in\mathcal T_p \big\} \\ &= \frac{1}{|\mathcal T_n|} \frac{1}{|\mathcal T_p|} \int_{t_n\in\mathcal T_n} \int_{t_p\in\mathcal T_p} \delta_{\lambda(t_n)<\lambda(t_p)}, \end{aligned} \label{auc0} \end{equation}[auc0]
where $\delta_{u{>}v}$ as an indicator function, which is $1$ if $u{>}v$.
Why use the AUC?
Because it depends only on the relative ranks of predicted firing rates, AUC is insensitive to errors in predicting a cell's mean rate, firing-rate variance, or any other strictly increasing transformation of the true rates. These properties make the AUC a desirable performance metric if we only care how well a GLM predicts neuronal tuning in a very general sense. It can be viewed as a type of nonlinear (rank) correlation between the test data and model predictions.
Calcium imaging yields fractional spike counts
Calcium fluorescence imaging has emerged as a standard method for simultaneously recording from large populations of neurons. Since calcium concentrations change relatively slowly following spiking, spike times are recovered indirectly by deconvolving calcium fluorescence signals.
These time-series of deconvolved spikes are sparse, with most entries zero. However, when spikes are detected, they are often detected with a fractional (non integer) value. This introduces some subtleties when working with point-process models, which are defined in terms of all-or-nothing (integer) events.
In practice, non-integer spike counts are no issue for the Poisson regression. Although a rigorous definition of the Poisson distribution is defined only on integers, it can be evaluated for any non-negative real number. This is technically improper, but works just as well and is commonly used in quasi-Poisson regression .
However, integrating commonly used measurements of model performance with fractional spike counts is slightly more subtle.
The AUC for continuous-time spiking data
For spiking data, the "true positives" are the spike times, and the "true negatives" are times when a neuron does not spike. We will work in continuous time, in which a spike train can be treated as a sum of delta distributions at each spike time, $y(t) = \sum_{t_s\in\text{spike times}} \delta(t-t_s)$.
In this continuous case, the duration of spikes is infinitesimal. This means that the set of true positives (i.e. times containing a spike) is of measure zero, compared to the measure of the time-duration of the experiment. In other words, the probability that a given time-point $t_n$ corresponds to a true negative, $\operatorname P_{t_n\in\mathcal T_n}$ is effectively $1$. True negatives are therefore uniformly distributed, with $\operatorname P_{t_n=t | t_n\in\mathcal T_n}= 1/T$.
In contrast, the true positives are a discrete set of $K = |\mathcal T_p|$ spikes, turning the second integral in $\eqref{auc0}$ into a sum. For a continuous spiking time-series, the AUC in $\eqref{auc0}$ can then be written in terms of the spike times $t_k\in\mathcal T_p$, $k\in\{1,\dots,K\}$ as:
\begin{equation} \begin{aligned} \text{AUC}(f,\mathcal T_p) &= \frac{1}{TK} \int_{t_n\in[0,T]} \sum_{k=1}^K \delta_{\lambda(t_n)<\lambda(t_k)} \\ &= \frac{1}{TK} \sum_{k=1}^K \int_{0}^{T} \delta_{\lambda(t)<\lambda(t_k)} \,dt. \end{aligned} \label{auc0b} \end{equation}[auc0b]
The AUC for deconvolved spikes
Deconvolved spikes from fluorescence data do not recover the true spikes counts $y(t)$, only a non-negative time-series $y(t)$ that is proportional to them (and is not quantized).
This is no problem in practice. To compute the AUC, we only need the density of true positive events per unit time, not their absolute number. The density of true positives per unit time is simply $y$ divided by its sum:
\begin{equation} \begin{aligned} \Pr(\,t_p{=}t \,|\, t_p{\in}\mathcal T_p\,) = \frac {y(t)} {\int_{dt} y(t)} \end{aligned} \label{pt} \end{equation}[pt]
To calculate the AUC, we could sample $t_n\in\mathcal T_n$ uniformly over $[0,T]$, and sample $t_p\in\mathcal T_p$ using the normalized deconvolved spike counts, $p_t$. We can compute the expectation of this sampling exactly by integrating over the joint distributions of true positives and negatives.
\begin{equation} \begin{aligned} \text{AUC}(\lambda,y) &= \frac 1 {T \int_{dt} y(t)} \int_0^T y(t_p) \int_0^T \delta_{\lambda(t_p){>}\lambda(t_n)} \,dt_n\,dt_p, \end{aligned} \label{auc1} \end{equation}[auc1]
We'd like to transform this equation into a nicer form, getting rid of that $\delta$ and producing something we can type into the computer.
Define $\operatorname s(t) \in [0,1]$ as the relative rank-ordering of $\lambda(t)$, which starts at 0 for the smallest $\lambda(t)$, and counts up to $1$ for the largest $\lambda(t)$:
\begin{equation} \begin{aligned} \operatorname s(t) &:= \frac {\text{rank of $\lambda(t)$}} {T} = \frac 1 T \int_0^T \delta_{\lambda(t){>}\lambda(\tau)}\,d\tau \end{aligned} \label{s} \end{equation}[s]
Using $\operatorname s(t)$, we can write the AUC as defined in $\eqref{auc1}$ as:
\begin{equation} \begin{aligned} \text{AUC}(\lambda,y) &= \int_0^T p(t) s(t)\,dt \end{aligned} \label{auc2} \end{equation}[auc2]
In discrete time, $\eqref{auc2}$ reduces to an average of the ranks $s_t\in[0,1]$, weighted by the deconvolved spike counts $y$. Here it is in Python:
def auc(λ,y):
'''
Area Under the ROC Curve (AUC) for deconvolved spike data.
Parameters
----------
λ: 1D np.float32, length NSAMPLES
Predicted scores for each datapoint
y: 1D np.float32, length NSAMPLES
Estimated spike counts for each datapoint
Returns
-------
AUC: float
Area under the receiver-operator-characteristic curve.
'''
# Cast to numpy arrays
λ = np.array(λ).ravel()
y = np.array(y).ravel()
# Ignore NaNs
bad = ~np.isfinite(λ)|~np.isfinite(y)
λ = λ[~bad]
y = y[~bad]
# AUC is undefined if there are no positive examples
d = np.sum(y)
if d<=0: return NaN
# Density of true positives per timepoint
py = y/d
# Convert scores to fractional ranks
rs = scipy.stats.rankdata(λ)/len(λ)
AUC = np.average(rs,weights=py)
return AUC
