Fitting a GLM to discrete-time events

[AamnaLawrence]

Hello everyone! I am trying to determine whether the activity of a neuron is modulated by discrete-time events in the task (lever insertion in the box, reward delivery etc). I prepare the events variable as follows:

events = nap.TsGroup(
    {
        1: nap.Ts(LeverInsert.t),
        2: nap.Ts(Reward.t),
     },
    metadata={"event_type": ["Lever Insert", "Reward"]},time_support=nap.IntervalSet(0,np.max(LeverInsert.t)))

An example of how the Lever Insert event looks like is:

Time (s)
46.153708
107.813152
197.822303
287.831483
340.031011
430.039936
465.079573
555.088619
645.097632
735.106611
...
3625.244659
3677.16383
3733.702946
3823.711426
3913.719999
4003.728668
4093.737208
4183.745749
4273.754195
4363.76283
shape: 61

When I try to compute the features for these time inputs, I get a bunch of nans. I have been playing around with different window sizes but it did not help me. Because of the possibly faulty feature matrix, all my GLM coefficients are 0. Here is the code I am using for finding the features:

bin_size = 0.01
binned_events = events.count(bin_size)
# pass time support to make sure they span the same range
count  = spikes_real.count(bin_size, ep=binned_events.time_support)

# define a basis over 500ms
window_size = int(0.5 / bin_size)
# add the basis using the label
add_basis = sum(*(nmo.basis.RaisedCosineLogConv(5, window_size, label=label)
      for label in events.event_type))

# pass each time series individually to the basis
X = add_basis.compute_features(binned_events[:, 0], binned_events[:, 1])

Has anyone used discrete-time events as inputs to GLM? Any insights on choosing the right basis functions on NeMoS would be very helpful.

Thank you!!

[BalzaniEdoardo]

@AamnaLawrence thanks the discussion and for sharing the data!

I’ve looked into this more closely, and your approach to modeling events is sound. However, there are a few important considerations:

• The total recording time is about 1 hour and 12 minutes, during which there are 24 reward events and 61 lever presses. This event imbalance means that if both events equally influence firing rate, statistical methods will be more confident in attributing an effect to lever presses, simply due to their higher frequency. If you’re still running experiments, increasing the number of both reward and lever events would improve robustness. Alternatively, balancing the events post hoc (by undersampling) is another option before performing any model comparison.

• Standard 80/20 train-test splits (such as in 5-fold cross-validation) are not well-suited here, since each test set would contain only ~4.8 reward events (24 * 0.2), which is too few for reliable evaluation. Instead, we use a 50/50 split to ensure enough predictive power in the test set. If additional generalization testing is needed, time-based validation may also be worth considering.

• Before modeling the entire population, it's useful to select a neuron that shows clear peri-event modulation. By plotting peri-event time histograms (PETHs) for all neurons, we can visually identify one that exhibits a strong firing rate change around events. This helps in choosing appropriate basis parameters, particularly the number of basis functions and whether a causal or acausal structure is needed. In the example, unit 482 was chosen because it showed strong modulation.

• The firing rate is not stationary over time, fluctuating significantly. This is evident when smoothing the spike count (spikes_real[482].count(0.01).smooth(std=4)). A model with only an intercept and event predictors will likely struggle unless we account for non-stationarity. One effective approach is to add a self-history predictor (as we did in the model using the RaisedCosineLogConv basis function), allowing the model to capture longer timescale dependencies.

• Since classic splitting does not work well, an alternative approach is to rank variable importance using a GroupLasso-regularized model. The strategy is to incrementally increase the regularization strength and observe the order in which predictors drop out. The most important variable will remain significant the longest. This method requires balancing events before applying it, or else the model will tend to overemphasize the more frequent event.

• The first step is to fine-tune the basis function parameters by fitting an unregularized GLM and assessing its performance in terms of peri-event modulation (both qualitatively by comparing predicted and observed PETHs, and quantitatively using pseudo-R² scores). Once the model is satisfactory, we proceed to ranking variables by fitting a GroupLasso regularized model over a range of regularization strengths and tracking both the pseudo-R² score and coefficient norms.

Interpreting the results

• If the test score drops significantly compared to training, the model may be overfitting.
• If the test set is too small, the test score will have high variance, making direct train-test comparisons unreliable.
• As regularization increases, the first predictors to drop out are typically the least relevant ones, allowing us to infer the most predictive features.

This first script can be used for wrangling with the basis parameters until we are satisfied with the fit (in terms of the score we get and the prediction of the peri-event modulation).

# Fit a GLM
import pynapple as nap
import numpy as np
import matplotlib.pyplot as plt
import nemos as nmo
import jax

jax.config.update("jax_enable_x64", True)

# load data and store events in a TsGroup
data = nap.load_file("/Users/ebalzani/Downloads/NWB_phy.nwb")

lever = data["LeverInsert"]
rew = data["Reward"]

events = nap.TsGroup(
    {
        1: nap.Ts(lever.t),
        2: nap.Ts(rew.t),
    },
    metadata={"event_type": ["Lever Insert", "Reward"]}, time_support=nap.IntervalSet(0, np.max(lever.t)))

# print number of events
print(f"Num reward events: {len(rew)}")
print(f"Num lever events: {len(lever)}")

# filter for active neurons (at least > 1Hz)
spikes_real = data["units"]
spikes_real = spikes_real[spikes_real.rate>1]

# select a neuron, here to find this one I plotted all the peri-events and pick one that looked modulated by eye
# as a starting point to design a model.
# from the peri-event you guess if the rate modulation is causal (after the event you see a change in rate) or not (the rate pre and post event show some modulation)
neu_id = 482
f,axs=plt.subplots(1,2)
cc = 0
for cc, item in enumerate(events.items()):
    i, ev = item
    peth = nap.compute_perievent(
      data=spikes_real[neu_id],
      tref=ev,
      minmax=(-3, 3),
      time_unit="s")
    axs[cc].plot(np.mean(peth.count(0.05), 1)/0.05, linewidth=3, color="red")
    axs[cc].set_title(events.metadata["event_type"][i])
    axs[cc].axvline(0, *axs[cc].get_ylim())
    

# use a 50-50 split to have some predictive power in the test set
duration = lever.time_support.end[-1] - lever.time_support.start[0]
train = nap.IntervalSet(0, duration * 0.5)
test = nap.IntervalSet(train.end[-1] + 0.0001, lever.time_support.end[-1])

# define a window size 
bin_size = 0.01
window = int(2 / bin_size)

# bin events and spike times
binned_events = events.count(bin_size)
count = spikes_real.count(bin_size, ep=binned_events.time_support) # pass time support to make sure they span the same range

# define a basis:
# for the most fair the comparison use the same number of basis function and (same type, at least for the two variable that you are directly comparing)
# for this neuron, it looks like the lever press is causal (rate increases after the press), for the reward, it looks like an acausal modulation
# you can model both type of predictors ("causal" is the default and most common)
add_basis = (
    nmo.basis.RaisedCosineLinearConv(11, window, width=4, label="Lever Insert") +  # default is causal
    nmo.basis.RaisedCosineLinearConv(11, window, width=4, label="Reward", conv_kwargs={"predictor_causality": "acausal"}) +
    nmo.basis.RaisedCosineLogConv(11, window, width=4, label="Spike History")
)

# define train and test design matrix (restrict first to avoid border effects)
X_train = add_basis.compute_features(
    binned_events[:, 0].restrict(train),
    binned_events[:, 1].restrict(train),
    count.loc[neu_id].restrict(train),
)
X_test = add_basis.compute_features(
    binned_events[:, 0].restrict(test),
    binned_events[:, 1].restrict(test),
    count.loc[neu_id].restrict(test),
)

# first check the the performance of an un-regularized model
model = nmo.glm.GLM(solver_name="LBFGS", solver_kwargs={"tol": 10 ** -12}).fit(
    X_train, count.loc[neu_id].restrict(train)
)

# print the scores: if there is a massive drop in score then we are over fitting.
# if the test set is too small the test score has too much variance so the K-fold directly may not be very informative
score_train = model.score(X_train, count.loc[neu_id].restrict(train), score_type="pseudo-r2-McFadden")
score_test = model.score(X_test, count.loc[neu_id].restrict(test), score_type="pseudo-r2-McFadden")
print(f"Train score: {score_train}\nTest score: {score_test}")

# compare the peri-event from the test set based on the predicted rate and the raw spikes 
rate_test = model.predict(X_test) / bin_size
f,axs=plt.subplots(1,2)
cc = 0
for cc, item in enumerate(events.items()):
    i, ev = item
    peth = nap.compute_perievent(
      data=spikes_real[neu_id].restrict(test),
      tref=ev.restrict(test),
      minmax=(-3, 3),
      time_unit="s")
    peth_rate = nap.compute_perievent_continuous(
      data=rate_test,
      tref=ev.restrict(test),
      minmax=(-3, 3),
      time_unit="s")
    axs[cc].plot(np.mean(peth.count(0.05), 1)/0.05, linewidth=3, color="red", label="raw")
    axs[cc].plot(np.nanmean(peth_rate, 1),label="model")
    axs[cc].set_title(events.metadata["event_type"][i])
    axs[cc].axvline(0, *axs[cc].get_ylim())
    axs[cc].legend()
plt.show()

Once this is satisfactory proceed with the GroupLasso ranking of the variable. (as a note, I quickly balanced the events - not in the script below because it was not systematic - and the variables seemed to be equally contributing for this neuron).

from tqdm import tqdm  # progress bar
# create some log-spaced regularization strength for group lasso
reg_str = np.geomspace(1E-5, 0.1, 20)
# one can use the whole data for this
X = add_basis.compute_features(binned_events[:,0], binned_events[:,1], count.loc[neu_id])

# define the variable grouping by constructing a mask
mask = np.zeros((len(add_basis), X.shape[1]))
cc = 0
for k, bas in enumerate(add_basis):
    mask[k, cc:cc+bas.n_basis_funcs] = 1
    cc += bas.n_basis_funcs
print(mask)

# initalize the arrays to store the norm of the coefficients for each predictor and the score
coeff_norms = np.zeros((len(reg_str), len(add_basis)))
scores = np.zeros(len(reg_str))

# loop over reg strength
for i, reg in tqdm(enumerate(reg_str)):
    regularizer = nmo.regularizer.GroupLasso(mask=mask)
    model = nmo.glm.GLM(
        regularizer=regularizer, solver_kwargs={"tol": 10 ** -12},
        regularizer_strength=reg
    ).fit(
        X, count.loc[neu_id]
    )
    coeff_dict = add_basis.split_by_feature(model.coef_, axis=0)
    # jax funciton that applies the norm to each vector in the coeff_dict equivalent of looping over the dict items
    coeff_norms[i] = jax.tree_util.tree_leaves(jax.tree_util.tree_map(np.linalg.norm, coeff_dict))
    scores[i] = model.score(X, count.loc[neu_id], score_type="pseudo-r2-McFadden")

# plot the results (since the events are unbalance, it may be misleading to interpret the output now, but if you balance
# the number of events, you can rank the variables from most to least significant)
# in general, as one adds more and more regularization, the score drops, and the norm of the coefficient too
f, axs = plt.subplots(2, 1, sharex=True)
axs[0].plot(reg_str, scores,"-ok")
axs[0].set_xscale('log')
axs[1].plot(reg_str, coeff_norms)
axs[1].legend([b.label for b in add_basis])
plt.show()

If the basis configuration is appropriate for multiple neurons, you can use the PopulationGLM to fit multiple neuron at the same time, this would speed up the computation time.

[AamnaLawrence]

Thank you so much! I will give this a go and keep you updated

[AamnaLawrence]

Hi @BalzaniEdoardo! After running the code I noticed that my r2 scores were pretty low (~0.05ish). I wonder if training the model around a certain window near the events might be better instead of taking the entire interval from first lever insertion to nth lever insertion since there are several events that happen during that time (the rat is freely moving and may groom etc) leading to modulation in firing rate and poor fitting.

[BalzaniEdoardo]

Hello!

I see your point—I had similar concerns when I first started fitting Poisson models. Restricting the analysis to intervals around the event is a good idea. It should lead to a better score and also speed up the fit, even though it won’t change the estimated coefficients.

I have a few thoughts on pseudo-R² scores and how I interpret them that might be helpful.

Understanding Pseudo-R² in Poisson Models

First, pseudo-R² cannot be interpreted in the same way as variance explained in linear regression. In Poisson models with low counts (e.g., when binning finely or when units have a low mean firing rate), pseudo-R² values tend to be much lower. From my experience, an extremely good fit might yield a score between 0.1 and 0.3.

Unfortunately, to my knowledge, there isn’t an established reference for interpreting the absolute values of these scores. My advice is to use pseudo-R² primarily for comparing models rather than for assessing absolute goodness-of-fit. For instance, in this paper by Benjamin et al., the authors compare GLMs to more advanced nonlinear encoding models, and the reported scores are comparable to what you might observe in your case.

Validating Model Scores via Simulation

To check whether my model’s scores are reasonable given the input and count statistics, I often simulate spike trains from a GLM with fixed coefficients. Typically, I take the coefficients I learned from fitting a GLM and use them to generate synthetic data. If the pseudo-R² scores for the simulated spikes fall within the same range as my actual fits, it suggests that my observed scores are near the best I can expect.

Here’s an example of how I generate simulated spikes and compute the pseudo-R² score:

# First, fit a GLM model
model = nmo.glm.GLM(solver_name="LBFGS", solver_kwargs={"tol": 1e-12}).fit(
    X_train, count.loc[neu_id].restrict(train)
)

# Extract model coefficients
coeff_lever = add_basis.split_by_feature(model.coef_, axis=0)["Lever Insert"]
intercept = model.intercept_

# Create a new model using the estimated coefficients
model_sim = nmo.glm.GLM()
model_sim.coef_ = coeff_lever
model_sim.intercept_ = intercept

# Define the basis function for the predictor
lever_basis = nmo.basis.RaisedCosineLogConv(
    11, window, width=4, label="Lever Insert", conv_kwargs={"predictor_causality": "causal"}
)

# Extract the corresponding design matrix columns
predictor_lever = add_basis.split_by_feature(X_train, axis=1)["Lever Insert"]

# Simulate spikes (without recurrent connections, the simulate method does not allow for recurrent simulations)
simulated_spikes, simulated_rate = model_sim.simulate(jax.random.PRNGKey(123), predictor_lever)

# Compute pseudo-R² score
model_sim.score(predictor_lever, simulated_spikes, score_type="pseudo-r2-McFadden")

Running this with a single predictor gives a pseudo-R² score of 0.00080415, which is very close to zero—even though I’m literally using the true model to generate the spikes.

The reason for this is that the pseudo-R² is measuring how much better the model’s predictions are compared to a Poisson model with a constant firing rate. Since lever events are rare, most design matrix entries are zero, so the improvement over a simple mean model is small.

So your intuition is correct! 🚀

PS The fact that your score was 0.05 combined with this simulation result suggests that the self-connectivity is quite important for model predictions here.

[BalzaniEdoardo]

you can also see this by plotting the simulated rate over time, and it is basically always constant with some brief spikes when the lever event happen

[AamnaLawrence]

Hi Edoardo!
I have been trying to use different predictors for instance the cumulative number of trials (in case a neuron is tracking the time in the session) in addition to the lever insertion time and spike history. Here is how I have been doing it:

bin_size=0.01
window = int(2 / bin_size)
lever=LeverInsert.count(bin_size) #LeverInsert are the time stamps of lever insertion event
trialcum=nap.Tsd(d=np.cumsum(lever),t=lever.t) #cumulative trials



count = spikes_real.count(bin_size, ep=LeverInsert.time_support)
basis1 = nmo.basis.RaisedCosineLogConv(11, window, width=4,label='Spike History')
basis2 = nmo.basis.RaisedCosineLinearConv(11, window, width=4, label="Lever Insert")
basis3 = nmo.basis.RaisedCosineLinearConv(11, window, width=4, label='Cumulative trials')

basis=basis1+basis2+basis3

time, basis_kernel1 = basis1.evaluate_on_grid(window)
time, basis_kernel2 = basis2.evaluate_on_grid(window)
time, basis_kernel3 = basis3.evaluate_on_grid(window)

duration = LeverInsert.time_support.end[-1] - LeverInsert.time_support.start[0]
train = nap.IntervalSet(0, duration * 0.5)
test = nap.IntervalSet(train.end[-1] + 0.0001, lever.time_support.end[-1])

X_train = basis.compute_features(
   
    count.loc[neu_id].restrict(train),
    lever.restrict(train),
    trialcum.restrict(train),
   
)
X_test = basis.compute_features(

    count.loc[neu_id].restrict(test),
    lever.restrict(test),
    trialcum.restrict(test),
    
)

# first check the the performance of an un-regularized model
model = nmo.glm.GLM(solver_name="LBFGS", solver_kwargs={"tol": 10 ** -12}).fit(
    X_train, count.loc[neu_id].restrict(train)
)

score_train = model.score(X_train, count.loc[neu_id].restrict(train), score_type="pseudo-r2-McFadden")
score_test = model.score(X_test, count.loc[neu_id].restrict(test), score_type="pseudo-r2-McFadden")
print(f"Train score: {score_train}\nTest score: {score_test}")

coef=basis.split_by_feature(model.coef_, axis=0)
fig, ax = plt.subplots(1, 3, figsize=(12, 6))
temp_weights = np.einsum('b, t b -> t', coef['Spike History'], basis_kernel1)
ax[0].plot(time, temp_weights)
ax[0].set_title('Spike history kernel weight')
temp_weights = np.einsum('b, t b -> t', coef['Lever Insert'], basis_kernel2)
ax[1].plot(time, temp_weights)
ax[1].set_title('Lever insertion time kernel weight')
temp_weights = np.einsum('b, t b -> t', coef['Cumulative trials'], basis_kernel3)
ax[2].plot(time, temp_weights)
ax[2].set_title('Cumulative trial kernel weight')
fig.tight_layout(h_pad=1)

The test and train scores are around 0.005 and I get the following plots:

I want to understand a clear interpretation of the figures above and if the features make sense. I would also like to add features like cumulative rewards etc but I am not sure how to construct the input feature for that.

Thank you for your help!