Regarding glm.nb() and my parameter

I have been doing a negative binomial regression model using the following code My my estimate here comes out as 3.48. (the exponential of the intercept).

The data was taken randomly (with set seed) from a distribution with my 3.2 so yeah it makes sense for it to be this close. BUT When i try and find the maximum likelihood function it appears to be peaking at a my value of around 3.13, which to me doesn't make much sense. Why aren't these to my's the same. I would assume the GLM.NB is using the maximum likelihood estimator on a log-scale?

Here is the script for reproducing the likelihood function: • To me, it would be extremely suspicious if your randomly generated data always had a maximum likelihood estimate close to the original parameter: that would be an indication of a problem with the random number generation. In what way, then, do these results not "make sense" to you? – whuber Dec 30 '18 at 16:37
• The randomly generated data is from a negativ binomial distribution with my = 3.2, so ofc an MLE would give me a coefficient close to this. What i dont understand is that the coefficient of the glm regressions is 3.48 but my Maximum likelihood funktion for the same dataset is peaking at 3.12, Should the coefficient be a represitation of my MLE? Isnt glm models coefficient made from MLEs? – AlexM Dec 31 '18 at 21:44

You are forgetting that glm.nb has to estimate the size parameter of the negative binomial distribution (called theta in glm.nb) as well as the intercept coefficient.

glm.nb is maximizing the likelihood with respect to two parameters (theta and the intercept coefficient) whereas you are maximizing with respect to only one (the intercept). So naturally glm.nb will find a higher value of the likelihood than you do and a different coefficient value.

When you ran summary(fit.what) the output would have shown you that glm.nb has estimated theta to be 0.3913. If you evaluated the negative binomial likelihood with size = 0.3913 instead of your value then your MLE for the intercept would agree with that from glm.nb.

I wonder whether you may have misunderstood what the size parameter represents for the negative binomial distribution. The variance of the distribution is given by $${\rm var}(y) = \mu + \mu^2 / \theta$$ where $$\mu$$ is the mean and $$\theta$$ is the size parameter. glm.nb assumes a constant value for $$\theta$$, not one that is proportion to $$\mu$$, as you have made it in your LL code. If you try to make $$\theta$$ depend on $$\mu$$, then you will destroy the quadratic mean-variance relationship that is a hallmark of the negative binomial distribution. There is a careful discussion of negative binomial GLMs in Chapter 10 of my book with Peter Dunn (Dunn and Smyth, 2018), including some code examples using the glm.nb function.

Reference

Dunn, PK, and Smyth, GK (2018). Generalized linear models with examples in R. Springer, New York, NY. https://doi.org/10.1007/978-1-4419-0118-7

• nice. If you felt like doing the work it would be nice to generate an image/contour plot of log-likelihood as a function of (mu, theta) and show the lines corresponding to the two size values . – Ben Bolker Jan 1 at 0:32
• i.e. emdbook::curve3d(-sum(dnbinom(y,mu=mu,size=size,log=TRUE)), xlim=c(2,5),ylim=c(0.2,0.5), varnames=c("mu","size"), sys3d="contour") abline(h=3.2/7,col="red") abline(h=fit.what\$theta,col="blue") ... although doesn't look this is actually the answer - mu-hat is independent of theta-hat ... ? – Ben Bolker Jan 1 at 0:40

thanks for the answer (Gordon Smyth). Super useful. I can now reproduce the LL using a size parameter of dnbinom set to constant 0.3913, then i get the peak at aronud 3.48. However our assignment was through this expression. This yields me the likelihood function directly through this script

LL2 <- function(my)

{

R2 = gamma(y+my/7)/(factorial(y)*gamma(my/7)7^(my/7))((7/(1+7))^(y+my/7))

sum(log(R2))

}

Vectorize(LL2)

LL2(3.13425)

par(mfrow=c(1,1))

curve(logLike2,from=2,to=4,xname="my")

I was wondering if there is a way of implementing the constant theta into this equation so it will reproduce the likelihoodfunction that peaks at 3.48 instead of my previous one where the size parameter is proportional to the mu that peaks at 3.12. I'm not sure how to interpret the size parameter of p(y), since if i set the mu = 3.2 and w = 7, the function will no longer have any mu for the loglikelihoodfunction.

We were told to write the loglikelihoodfunction of the p(y) given above as a function of mu, given w = 7. Im not sure which way is better to combat this - which LL fucntion would be best to answer this? The one peaking at 3.12 or the one peaking at 3.48?

• Your p(y) has a dependence between the NB mean and size parameters. You cannot reproduce results from this likelihood by using glm.nb or by assuming theta constant or by assuming mu and theta independent. You obviously have to implement p(y) by programming it directly -- you seem to have already done that so I am unclear why you are trying to do anything else. – Gordon Smyth Jan 2 at 8:53
• That's what i guessed aswell - just wanted to make sure there was no clear relationship i was missing :D Thanks a lot – AlexM Jan 2 at 15:05