Hi folks,

I have been doing a negative binomial regression model using the following code (note that MASS package is used for glm.nb)

set.seed(1234)

(y <- rnbinom(100, mu = 3.2, size= 3.2 / 7))

fit.what <- glm.nb(y~1)

summary(fit.what)

exp(fit.what$coefficients)

My mu-estimate here comes out as 3.48. (the exponential of the intercept).

The data was taken randomly (with set seed) from a distribution with mu at 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 mu value of around 3.13, which to me doesn't make much sense. Why aren't these to mu's the same. I would assume the GLM.NB is using the maximum likelihood estimator on a log-scale?ยจ

Here is the loglikelihoodfunction

LL <- function(my)

{

R=dnbinom(y,mu=my,size=my/7)

sum(log(R))

}

logLike <- Vectorize(LL)

LL(3.13425)

par(mfrow=c(1,1))

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