Implement GLM-based firing rate estimation

[abuzarmahmood]

Closes #708

Summary

Adds utils/glm/ module for GLM-based firing rate estimation using the nemos library.

Architecture

Due to dependency conflicts between nemos and blech_clust, this module uses a two-environment architecture:

infer_glm_rates.py (orchestrator)
        │
        ├── _glm_extract_data.py  →  runs in blech_clust conda env
        │                             (uses ephys_data for spike extraction)
        │
        └── _glm_fit_models.py    →  runs in separate nemos venv
                                      (fits GLM, generates plots, writes HDF5)

Data is passed between environments via temporary numpy files in utils/glm/_temp/.

Features

• Spike history modeling: Raised cosine log basis functions for autoregressive effects
• Stimulus features: Optional stimulus timing input
• Neural coupling: Optional coupling between simultaneously recorded neurons
• Bits per spike metric: Model comparison against baseline Poisson
• Flexible grouping: Process by taste and/or brain region

Setup

1. blech_clust conda environment (existing)
2. Create nemos venv:

   python -m venv ~/nemos_env
   source ~/nemos_env/bin/activate
   pip install nemos matplotlib pandas tables

Usage

python utils/glm/infer_glm_rates.py <data_dir> [options]

The script auto-detects environments or prompts if not found.

Output

• Predicted firing rates in HDF5 under /glm_output/regions/
• Bits per spike summary CSV
• Individual neuron prediction plots
• Pickled GLM models

[abuzarmahmood]

Implementation Details

GLM Model Structure

The script fits a Poisson GLM for each neuron with the following components:

1. Spike History Filter: Captures autoregressive effects using nemos.basis.RaisedCosineLogConv. The log-spaced basis functions efficiently represent both short and long timescale dependencies with fewer parameters than raw history.

2. Stimulus Features (optional): Models the effect of stimulus onset using a convolved indicator function. This captures post-stimulus response dynamics.

3. Neural Coupling (optional): Models how activity in other neurons affects the target neuron's firing. Uses the same raised cosine basis but with reduced dimensionality.

Bits Per Spike Metric

The bits per spike metric quantifies model performance as information gain over a baseline:

bits_per_spike = (LL_model - LL_baseline) / (n_spikes * log(2))

Where:

• LL_model: Log-likelihood under the fitted GLM
• LL_baseline: Log-likelihood under a homogeneous Poisson model (mean rate)

Positive values indicate the GLM captures structure beyond the mean rate. Typical values for well-fit models range from 0.1-1.0 bits/spike.

Comparison to RNN Approach

This GLM approach complements the existing infer_rnn_rates.py:

Aspect	GLM (nemos)	RNN (blechRNN)
Interpretability	High - coefficients map to features	Low - latent dynamics
Computation	Fast - convex optimization	Slow - gradient descent
Flexibility	Structured features	Learns representations
Dependencies	nemos (pip)	blechRNN (external)

Future Extensions

Per issue #708, potential additions:

• Cross-validation for regularization strength
• Model selection across different feature sets
• Comparison plots between GLM and RNN predictions