0
$\begingroup$

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:

Log-likelihood formula for a binomial

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.

$\endgroup$

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$ Commented Oct 9, 2022 at 7:04
  • $\begingroup$ @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 $\endgroup$ Commented Oct 10, 2022 at 15:13
  • $\begingroup$ 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. $\endgroup$ Commented Oct 11, 2022 at 10:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.