Function for Poisson Regression

Hello,

I am attempting to perform poisson regression on an array I have that looks to follow a poisson distribution. I have been attempting to do this using glmfit() with the poisson argument, along with glmval to obtain the y values of the fitted distribution, as follows:

%Modeling spike counts as a function of time for action neuron
[histActionNeuron, bins] = histcounts(ActionNeuronSpikeTimes);
figure;
hold on;
plot(histActionNeuron)
[B, FitInfo] = glmfit(histActionNeuron,bins(1:81)+300,"poisson");
yfit = glmval(B,histActionNeuron,"probit");

histActionNeuron is a 1x81 double, which is why I used bins(1:81), and the +300 is because bins starts at -300 ms, and glmfit did not accept negative values in bins.

When I plot histActionNeuron, it looks like the following:

However, yfit turns out to be an 81x1 double of ones, indicating that there is some error in how I am performing the regression. Is there some other function I should use that would give me better results? I have been experimenting with lassoglm and poissfit, but haven't been successful. Any guidance would be appreciated.

Thanks!

Answer 1

[Saarthak Gupta]

Hi Coulter,

It seems that you are trying to fit a Poisson distribution to your data.

The “glmfit” function fits a generalized linear regression model to your data. Generalized linear regression is a statistical approach that builds on regular linear regression. This method uses a link function to relate the predicted values to the actual data, making it suitable for different kinds of data including binary, counts, and continuous.

Poisson regression is based two assumptions:

1. The response variable Y follows a Poisson distribution and,

2. The logarithm of expected value of Y (response variable) can be expressed as a linear combination of unknown parameters.

While the first assumption holds for your model, the second assumption does not.

Instead of fitting a generalized linear model, a better approach would be to fit a Poisson distribution to your ActionNeuronSpikeTimes data using “fitdist” function (Statistics and Machine Learning Toolbox). Refer to the trailing links for documentation of “fitdist” and relevant examples.

Equivalently, you can fit the Poisson probability density function to your data using the Curve Fitting Toolbox. Refer to the following code:

lambda = 4;     % parameter for Poisson distribution (to be estimated)
x = (0:20)';
y = pdf('Poisson',x,lambda);     % density values for samples in x

% create anonymous function which takes l and x as inputs and returns the
% Poisson pdf as the output. l is the parameter to be estimated.
g = fittype( @(l,x) pdf('Poisson',x,l) );
% fit the distribution
ls = fit(x,y,g,"StartPoint",1);

coeffvalues(ls)

ans = 4.0000

The estimated value for lambda is 4, as expected.

P.S. Poisson regression typically employs logarithm as its link function. However, in your use of "glmval," you have specified probit function as the link function. This function behaves differently from the logarithmic function, which could explain why your "yfit" results in an array of ones.

Refer to the following MATLAB documentation for further reference:

• Generalized Linear Regression: https://www.mathworks.com/help/stats/generalized-linear-regression.html
• glmval: https://www.mathworks.com/help/stats/glmval.html
• glmfit: https://www.mathworks.com/help/stats/glmfit.html
• fitdist: https://www.mathworks.com/help/stats/fitdist.html
• fittype: https://www.mathworks.com/help/curvefit/fittype.html
• pdf: https://www.mathworks.com/help/stats/prob.normaldistribution.pdf.html