The mean learned by glm.nb is not the same as the true mean?

Question

I am not able to make the predict values to be the same as the true means (mu). I'd expect that they should be the same. Does anybody see there anything wrong in my code?

library(MASS)
n = 10000
x = seq_len(n)/n
mu = exp(2*x)
y = rnegbin(n = n, mu = mu, theta = .1)
fit = glm.nb(y ~ mu)
plot(mu, predict(fit, type='response'))

Answer 1

When you fit a negative binomial regression with a log link, the mean is expressed by:

$$ \mu_i = exp( \beta_0 + \beta_1x1 +.. \beta_ix_i) $$

So you need to fit it to the log of the mean:

fit = glm.nb(y ~ log(mu))

Or you can do away with the log link:

fit = glm.nb(y ~ mu,link="identity")