How to deal with high baseline firing rates

[vigji]

Describe the bug
Not sure whether this is a bug, or I am just being stupid with my fitting. Note: not a big expert of LNP model fitting, those are my first attempt!

I encountered this issue when trying to fit some neurons that had high (around 10-20 spk/s) firing rate.
I consistently encounter this error:

File ~/miniconda3/envs/lab-env/lib/python3.10/site-packages/nemos/glm.py:546, in GLM._initialize_parameters(self, X, y)
    [544](https://file+.vscode-resource.vscode-cdn.net/Users/vigji/code/mpa-ephys/notebooks/~/miniconda3/envs/lab-env/lib/python3.10/site-packages/nemos/glm.py:544) func_root = root(func, y.mean(axis=0, keepdims=False), method="hybr")
    [545](https://file+.vscode-resource.vscode-cdn.net/Users/vigji/code/mpa-ephys/notebooks/~/miniconda3/envs/lab-env/lib/python3.10/site-packages/nemos/glm.py:545) if not jnp.allclose(func_root.fun, 0, atol=10**-4):
--> [546](https://file+.vscode-resource.vscode-cdn.net/Users/vigji/code/mpa-ephys/notebooks/~/miniconda3/envs/lab-env/lib/python3.10/site-packages/nemos/glm.py:546)     raise ValueError(
    [547](https://file+.vscode-resource.vscode-cdn.net/Users/vigji/code/mpa-ephys/notebooks/~/miniconda3/envs/lab-env/lib/python3.10/site-packages/nemos/glm.py:547)         "Could not set the initial intercept as the inverse of the firing rate for "
    [548](https://file+.vscode-resource.vscode-cdn.net/Users/vigji/code/mpa-ephys/notebooks/~/miniconda3/envs/lab-env/lib/python3.10/site-packages/nemos/glm.py:548)         "the provided link function. "
    [549](https://file+.vscode-resource.vscode-cdn.net/Users/vigji/code/mpa-ephys/notebooks/~/miniconda3/envs/lab-env/lib/python3.10/site-packages/nemos/glm.py:549)         "Please, provide initial parameters instead!"
    [550](https://file+.vscode-resource.vscode-cdn.net/Users/vigji/code/mpa-ephys/notebooks/~/miniconda3/envs/lab-env/lib/python3.10/site-packages/nemos/glm.py:550)     )
    [551](https://file+.vscode-resource.vscode-cdn.net/Users/vigji/code/mpa-ephys/notebooks/~/miniconda3/envs/lab-env/lib/python3.10/site-packages/nemos/glm.py:551) initial_intercept = jnp.atleast_1d(func_root.x)
    [553](https://file+.vscode-resource.vscode-cdn.net/Users/vigji/code/mpa-ephys/notebooks/~/miniconda3/envs/lab-env/lib/python3.10/site-packages/nemos/glm.py:553) # Initialize parameters

ValueError: Could not set the initial intercept as the inverse of the firing rate for the provided link function. Please, provide initial parameters instead!

I do not understand what the initialisation method is doing, I would have expected intercepts to just be the log(firing_rates) from the documentation. Before keep looking at my data, I tried to reproduce the issue with the examples.

To Reproduce
This is a combination that recapitulate my issue. Code slightly adapted from the examples:

import jax.numpy as jnp
import nemos as nmo
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(123)

n_features = 18
n_neurons = 60
n_samples = 500

# random design array. Shape (n_time_points, n_features).
X = 0.5 * np.random.normal(size=(n_samples, n_features))

# log-rates & weights
b_true = np.random.uniform(size=(n_neurons,))*3  # baseline rates
w_true = np.random.uniform(size=(n_features, n_neurons))*0.1  # real weights:


# generate counts (spikes will be (n_samples, n_features)
rate = jnp.exp(jnp.dot(X, w_true) + b_true)
spikes = np.random.poisson(rate)

plt.figure()
plt.plot(spikes)
plt.show()

model = nmo.glm.PopulationGLM()
model.fit(X, spikes)

print(f"population GLM log-likelihood: {model.score(X, spikes)}")

Expected behavior
(I do not understand why the model would not fit, and I do not know what a suitable initialisation would be.)
Edit: the initialisation described in the docs - log(avg_firing) - actually works, I do not understand why it is not the default one actually implemented!

Screenshots
If applicable, add screenshots to help explain your problem.

Desktop (please complete the following information):

• OS: macOS
• nemos==0.1.6

[BalzaniEdoardo]

Hi Luigi, first of all, thanks for this issue, really appreciate all the details provided. I am looking into it now; The initialization was meant to find a numerical inverse of the link function.

For exp, we know it is log, but since we allow to pass any arbitrary non-linearity I wanted a more general approach; I had that issue already, and I thought I'd fixed it using scipy root finding to get the inverse, which was stable enough on my tests.
I'll try your synthetic data on a more recent branch, like development to see if the issue is already resolved and get back to you soon.

[BalzaniEdoardo]

@vigji I probably won't be able to merge this in by this week because we are under a deadline, but I created a branch that fixes the issue by:

1. Using known inverse link function when possible (log, and inverse soft-plus).
2. If not possible, it will try to invert the link function numerically. In the new implementation we use a loop over the neurons and call scipy.optimize.root_scalar since for each neuron can be treated independently. This makes the problem much simpler, and more stable. I tried it out on the synthetic data you generated and works well.

If you want to test it out, try this branch. If you want to check the numerical inversion on the real data too, pass as a link function lambda x: jax.numpy.exp(x) instead of the exponential directly.

Let me know,
Edo

[vigji]

Awesome, thanks! I will try out the branch.

[vigji]

Yup, this fixed the issue for me! thanks!
I'll keep the issue for ref until the branch gets merged if that is ok

[BalzaniEdoardo]

Sounds good!

…

On Mon, Oct 21, 2024, 5:08 AM Luigi Petrucco ***@***.***> wrote: Yup, this fixed the issue for me! thanks! I'll keep the issue for ref until the branch gets merged if that is ok — Reply to this email directly, view it on GitHub <#243 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AG2MOSZXDQZ64ZUF3Q6BV2TZ4TACNAVCNFSM6AAAAABPQC53H6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMRWGA3TCNJSGU> . You are receiving this because you commented.Message ID: ***@***.***>

[BalzaniEdoardo]

@vigji I merged the solution into development, it will be fixed in the next release, I'll still leave this open until then.

[davidgwyrick]

Just to jump in here, it does appear this has been fixed with the latest release, but what about low baseline firing rates? I'm running into an issue where depending on how I chunk up the data for cross-validation, some neurons will have no spikes for that window. It throws a different error, but probably for the same reason?

ValueError: Solver returned at least one NaN parameter, so solution is invalid! Try tuning optimization hyperparameters, specifically try decreasing the stepsize and/or setting acceleration=False.

Maybe the solution is to mask those neurons somehow?

[BalzaniEdoardo]

Hi David, I'll look into it, but yes, if there are no spikes for a neuron it may lead to the optimizer making the weights zero and the intercept very very large and negative, since the best possible prediction is firing rate near zero. If the negative intercept can grow unbounded it may diverge. If I were you, I would filter the neuron making sure that in each fold the mean firing rate is above a threshold (try 1hz if that makes sense for the population you are fitting). Let me know if the fit goes through with this additional step.

…

On Wed, Feb 12, 2025, 7:23 PM David Wyrick ***@***.***> wrote: Just to jump in here, it does appear this has been fixed with the latest release, but what about low baseline firing rates? I'm running into an issue where depending on how I chunk up the data for cross-validation, some neurons will have no spikes for that window. It throws a different error, but probably for the same reason? ValueError: Solver returned at least one NaN parameter, so solution is invalid! Try tuning optimization hyperparameters, specifically try decreasing the stepsize and/or setting acceleration=False. Maybe the solution is to mask those neurons somehow? — Reply to this email directly, view it on GitHub <#243 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AG2MOS7AZFSII3BAGA2JFA32PPQZLAVCNFSM6AAAAABPQC53H6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMNJVGE2DQNBQG4> . You are receiving this because you commented.Message ID: ***@***.***>

[BalzaniEdoardo]

@davidgwyrick, adding to my previous answer, if you are using pynapple and your spikes are a TsGroup, you can threshold using the rate metadata:

# assume that units is a TsGroup
high_firing_units = units[units.rate > 1]
counts = high_firing_units.count(bin_size)
# .. fit the counts

And if you must to make sure that in all folds the rate is greater than one, I would check the rate in each fold separately and grab those units only.
You can also use array indexing and mask the counts, but we suggest to use pynapple for this type of manipulation, since it is more readable and deals with gaps in your time series.

And if you can generate an example that results the same problem and give me more details on what you are trying to achieve, I can see if there is better solution.