# Statsmodels GLM binomial distribution: manual calculation of the null log-likelihood

For my better understanding, I'm trying to calculate manually the log-likelihood that statsmodels GLM is reporting in the summary results. However, my manual calculation does not match the result that I see in the summary.

Let's take for example the null log-likelihood from the GLM model applied to the the star98 dataset used in the statsmodel example pages in which the binomial distribution was applied.

import statsmodels.api as sm
glm_binom = sm.GLM(data.endog, data.exog, family=sm.families.Binomial())
res = glm_binom.fit()

# Null log-likelihood given by the statsmodels GLM results
res.llnull


gives a value of -18131.91.

To calculate that manually I used the formula that I found at slide 23 of this presentation and that I copy here:

I then replaced Xi*Beta with res.null because in the case of the null-model this is the constant prediction.

The formula translated into python is then the following:

np.sum((yi)*res.null) - np.sum(np.log(1+np.exp(res.null)))


Then the question is: what is exactly yi? I would assume that's the ratio between NABOVE and total number of cases i.e. NABOVE + NBELOW.

yi = data.endog['NABOVE']/(data.endog['NABOVE']+data.endog['NBELOW'])


In other words the success/fail rate that the model is trying to predict. But if I apply that, then the result is very far from what I read out of the res.llnull printout.

What's wrong in what I do?

Thanks.

The null likelihood is computed using the marginal estimate of the mean of the outcome. In other words, it is the negative log likelihood of a model with a single parameter (namely, the intercept).

Here is how to compute the null likelihood using python and a comparison with statsmodels

import numpy as np
import statsmodels.api as sm
from scipy.special import expit

# Simulate some data
x = np.random.normal(size = 100)
eta = -2 + 0.8*x
p = expit(eta)
y = np.random.binomial(p=p, n=1)

# Fit the model
fit = sm.Logit(y, X).fit()
fit.summary()

marginal_p = y.mean()

# Here I am computing the likelihood using a model with a single parameter
# I know that p will just be the average of the ys in this model
nLL = (np.log(marginal_p)*y + np.log(1-marginal_p)*(1-y))

# Does not yield an exception, so this must be true.
np.testing.assert_allclose(nLL.sum(), fit.llnull)


I don't see you computing the average if the y in your solution, which si probably the source of the error.

• While I understand and follow your example, I still don't get the numbers match with my dataset. I added the code I'm using in my original question. I assume the problem lies in what I assign to yi but I can't figure it out. Commented Oct 9, 2022 at 7:04
• @GianlucaColangelo Well to start, I don't see any methods in the fitted mode corresponding to null. All I see is res.llnull, and that isn't the marginal probability its the log likelihood of the null value. So let's start with this: What is res.null supposed to give you Commented Oct 10, 2022 at 15:13
• res.null it's an array representing the fitted values of the null model. So all elements of the array are equal and basically representing y.mean(). So what you call marginal_p in your example. Actually, despite being close the two quantities are not exactly the same so I could replace the above formula with this one: np.sum((y)*y.mean()) - len(y)*(np.log(1+np.exp(y.mean()))) but the result is still pretty much the same. Commented Oct 11, 2022 at 10:10