What is the distribution of theta in a negative binomial model (glm.nb with R)?

Good morning every body.

My question concerns the distribution of the $\theta$ parameter in a glm with a negative binomial distribution, such as $V(X)=\mu+\theta\mu^2$.
Indeed, $\theta$ is expected to follow a Gamma distribution but how to estimate the parameters of this distribution?

In R, when I use the function glm.nb I obtain the expected value of $\theta$ and its associated SE but I am not able whith these statistics to compute the parameter of the Gamma distribution. Does anybody knows how to derive these parameters from the output of the glm.nb function?

Regards, Maxime

-
In classical statistics, $\theta$ doesn't follow any distribution, it's just a number. You are getting an estimate of that number, and a standard error associated with that estimate, from glm.nb. However, to reiterate, $\theta$ itself doesn't follow any distribution. If, on the other hand, you want to be Bayesian about it, there's no need to assume that $\theta$ is distributed Gamma, and you wouldn't use glm.nb anyway, preferring one of the MCMC packages in R. –  jbowman Jun 13 '12 at 17:07
@jbowman. I do not agree with this statement. I used a monte-carlo approach and I found that \theta follow a gamma distribution with shape and scale parameters and that these parameters could be estimated from the expected value and from the variance estimated with glm.nb. This was also corroborated by the posterior distribution of Bayesian estimates <a href="en.wikipedia.org/wiki/…; title="Wiki">Negative Binomial Distribution</a> –  Maxime Jul 9 '12 at 13:10