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

Question

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
data = sm.datasets.star98.load()
data.exog = sm.add_constant(data.exog, prepend=False)
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.

Answer 1

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 statsmodels.tools import add_constant
from scipy.special import expit

# Simulate some data
x = np.random.normal(size = 100)
X = add_constant(x)
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.