  Accurate statistical models of neural spike responses can characterize the
information carried by neural populations. But the limited samples of spike
counts during recording usually result in model overfitting. Besides, current
models assume spike counts to be Poisson-distributed, which ignores the fact
that many neurons demonstrate over-dispersed spiking behaviour. Although the
Negative Binomial Generalized Linear Model (NB-GLM) provides a powerful tool
for modeling over-dispersed spike counts, the maximum likelihood-based standard
NB-GLM leads to highly variable and inaccurate parameter estimates. Thus, we
propose a hierarchical parametric empirical Bayes method to estimate the neural
spike responses among neuronal population. Our method integrates both
Generalized Linear Models (GLMs) and empirical Bayes theory, which aims to (1)
improve the accuracy and reliability of parameter estimation, compared to the
maximum likelihood-based method for NB-GLM and Poisson-GLM; (2) effectively
capture the over-dispersion nature of spike counts from both simulated data and
experimental data; and (3) provide insight into both neural interactions and
spiking behaviours of the neuronal populations. We apply our approach to study
both simulated data and experimental neural data. The estimation of simulation
data indicates that the new framework can accurately predict mean spike counts
simulated from different models and recover the connectivity weights among
neural populations. The estimation based on retinal neurons demonstrate the
proposed method outperforms both NB-GLM and Poisson-GLM in terms of the
predictive log-likelihood of held-out data. Codes are available in
https://doi.org/10.5281/zenodo.4704423
