When do Poisson and negative binomial regressions fit the same coefficients? I’ve noticed that in R, Poisson and negative binomial (NB) regressions always seem to fit the same coefficients for categorical, but not continuous, predictors. 
For example, here's a regression with a categorical predictor:
data(warpbreaks)
library(MASS)

rs1 = glm(breaks ~ tension, data=warpbreaks, family="poisson")
rs2 = glm.nb(breaks ~ tension, data=warpbreaks)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))


Here is an example with a continuous predictor, where the Poisson and NB fit different coefficients:
data(cars)
rs1 = glm(dist ~ speed, data=cars, family="poisson")
rs2 = glm.nb(dist ~ speed, data=cars)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))


(Of course these aren't count data, and the models aren't meaningful...)
Then I recode the predictor into a factor, and the two models again fit the same coefficients:
library(Hmisc)
speedCat = cut2(cars$speed, g=5) 
#you can change g to get a different number of bins

rs1 = glm(cars$dist ~ speedCat, family="poisson")
rs2 = glm.nb(cars$dist ~ speedCat)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))


However, Joseph Hilbe’s Negative Binomial Regression gives an example (6.3, pg 118-119) where a categorical predictor, sex, is fit with slightly different coefficients by the Poisson ($b=0.883$) and NB ($b=0.881$). He says: “The incidence rate ratios between the Poisson and NB models are quite similar. This is not surprising given the proximity of $\alpha$ [corresponding to $1/\theta$ in R] to zero.”
I understand this, but in the above examples, summary(rs2) tells us that $\theta$ is estimated at 9.16 and 7.93 respectively.
So why are the coefficients exactly the same? And why only for categorical predictors?

Edit #1
Here is an example with two non-orthogonal predictors. Indeed, the coefficients are no longer the same:
data(cars)

#make random categorical predictor
set.seed(1); randomCats1 = sample( c("A","B","C"), length(cars$dist), replace=T)
set.seed(2); randomCats2 = sample( c("C","D","E"), length(cars$dist), replace=T)

rs1 = glm(dist ~ randomCats1 + randomCats2, data=cars, family="poisson")
rs2 = glm.nb(dist ~ randomCats1 + randomCats2, data=cars)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))


And, including another predictor causes the models to fit different coefficients even when the new predictor is continuous. So, it is something to do with the orthogonality of the dummy variables I created in my original example?
rs1 = glm(dist ~ randomCats1 + speed, data=cars, family="poisson")
rs2 = glm.nb(dist ~ randomCats1 + speed, data=cars)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))


 A: In order to see what is going on here, it is useful first to do the regression without the intercept, since an intercept in a categorical regression with just one predictor is meaningless:
> rs1 = glm(breaks ~ tension-1, data=warpbreaks, family="poisson")
> rs2 = glm.nb(breaks ~ tension-1, data=warpbreaks)

Since Poisson and negative binomial regressions specify the log of the mean parameter, then for categorical regression, exponentiating the coefficients will give you the actual mean parameter for each category:
>  exp(cbind("Poisson"=coef(rs1), "NB"=coef(rs2)))
          Poisson       NB
tensionL 36.38889 36.38889
tensionM 26.38889 26.38889
tensionH 21.66667 21.66667

These parameters correspond to the actual means over the different category values:
> with(warpbreaks,tapply(breaks,tension,mean))
       L        M        H 
36.38889 26.38889 21.66667 

So what is happening is that the mean parameter $\lambda$ in each case that maximizes likelihood is equal to the sample mean for each category.  
For Poisson distribution it is clear why this occurs. There is only one parameter to fit, and thus maximizing the overall likelihood of a model with a single categorical predictor is equivalent to independently finding a $\lambda$ which maximizes likelihood for the observations in each particular category.  The maximum likelihood estimator for the Poisson distribution is simply the sample mean, which is why the regression coefficients are precisely the (logs of the) sample means for each category.
For negative binomial it is not quite as simple, because there are two parameters to fit: $\lambda$ and the shape parameter $\theta$.  Moreover, the regression fits a single $\theta$ that covers the whole dataset, so in this situation a categorical regression is not simply equivalent to fitting a completely separate model for each category.  However, by examining the likelihood function, we can see that for any given theta, the likelihood function is again be maximized by setting $\lambda$ to the sample mean:
\begin{align}
L(X,\lambda,\theta) &= \prod \left(\frac{\theta}{\lambda+\theta}\right)^\theta
\frac{\Gamma(\theta + x_i)}{x_i!\Gamma(\theta)}\left(\frac{\lambda}{\lambda+\theta}\right)^{x_i}\\
\log L(X,\lambda,\theta) &= 
\sum\theta\left(\text{log}\theta-\text{log}(\lambda+\theta)\right)
+x_i\left(\text{log}\lambda-\text{log}(\lambda+\theta)\right)
+\log\left(\frac{\Gamma(\theta + x_i)}{x_i!\Gamma(\theta)}\right)\\
\frac{d}{d\lambda}\log L(X,\lambda,\theta) &= 
\sum \frac{x_i}{\lambda}-\frac{\theta+x_i}{\lambda+\theta}
=n\left(\frac{\bar{x}}{\lambda}-\frac{\bar{x}+\theta}{\lambda+\theta}\right),
\end{align}
so the maximum is attained when $\lambda=\bar{x}$.
The reason that you don't get the same coefficients for continuous data is because in a continuous regression, $\text{log}(\lambda)$ is no longer going to be a piecewise constant function of the predictor variables, but a linear one.   Maximizing the likelihood function in this case will not reduce to independently fitting a value $\lambda$ for disjoint subsets of the data, but will rather be a nontrivial problem that is solved numerically, and is likely to produce different results for different likelihood functions.
Similarly, if you have multiple categorical predictors, despite the fact that the fitted model will ultimately specify $\lambda$ as a piecewise constant function, in general there will not be enough degrees of freedom to allow $\lambda$ to be determined independently for each constant segment.  For example, suppose that you have $2$ predictors with $5$ categories each.  In this case, your model has $10$ degrees of freedom, whereas there are $5*5=25$ unique different combinations of the categories, each of which will have its own fitted value of $\lambda$.  So, assuming that the intersections of these categories are non-empty (or at least that $11$ of them are nonempty), the likelihood maximization problem again becomes nontrivial and will generally produce different outcomes for Poisson versus negative binomial or any other distribution.
