Model validation after fitting a negative binomial GLM in R Ok, I have searched and searched and just have no clue where to start.  First, what I would like to do is produce a QQ-plot (or even a readable residual plot) to look at the fit of my model.  I guess I just don't understand how the parameters that go into qnbinom() are obtained from the output of MASS::glm.nb().  I am attempting to use probplot() from package e1071, but am unsure of the inputs needed.  It would be great if someone having experience with fitting negative binomials could lend a hand.  
Secondly, I came across a residual plot here: http://www.stat.cmu.edu/~hseltman/Rclass/R8.R
I can make it to work, but I don't know how to interpret it or if I am using it correctly.  Has anyone else used this?
At the moment I am relying on the AIC and plots of fitted vs. actual values to assess the fit of my model and I would like something a little better!
Edit:
Hopefully this will clarify what I am asking.  With qnbinom(p, size, prob, mu, lower.tail = TRUE, log.p = FALSE), how do I (or is it even possible to) get the p, size, prob, mu from the output of a glm.nb fit model?  From my research, I have found that size is the dispersion parameter, but other than that I'm not sure where to go.  I know theta goes in there somehow, just not sure how to get it in the form needed. 
Edit 2: Ok, once I have a distplot(), is there a guide to interpreting it?  I am fairly positive I have a bad fit because I have a curved plot with a red line going through it (with many points on the tails far from the red line).  The prob: ML = 0.011, is this rejecting that the distribution is from the negative binomial specified?
 A: You might find distplot() from the vcd package useful either for the original data (edit: you can't use it on residuals). This plots Friendly's "negativebinomialness plots" and provides how well the negative binomial model fits
distplot(response, type = "nbinomial", ...)
To obtain the parameters: glm.nb uses the "Gamma mixture of Poisson" representation. It is actually a log-linear model that is fitted, so you should get the mean as $\exp(X\beta)$.
For example, let's say your data come from a negbin with mean 5 and theta of 1 (in the alternative representation as described above). Then you can get the mean estimate simply by  
set.seed(10)  
df <- data.frame(y=rnbinom(100,size=1,mu=5))  
m0 <- glm.nb(y~1,data=df)  
m0  
exp(coef(m0))  
m0$theta  

which are in this case 5.1 for the mean (pretty close) and 1.6 for the dispersion parameter (pretty far off). 
If you fit a model for the conditional mode, you interpret it accordingly as in every other log linear model, see this discussion on stack exchange. 
EDIT: If you want to know how to get the mean in a negbin regression model you need to sum up the linear predictor $X\beta$.
For example: I take the quine data and fit  
m1 <- glm.nb(Days~Sex,data=quine)

now males are 1 females are 0. To get the mean for males you write 
> exp(coef(m1)[1]+coef(m1)[2]*1)  
[1] 17.95455    

and for females
> exp(coef(m1)[1]+coef(m1)[2]*0)     
[1] 15.225  

Now to get the mean you must weight this with the occurence of all females and males which is 
> table(quine$Sex)  
 F  M   
80 66  

and hence the mean is
> (80/(66+80))*15.225+(66/(80+66))*17.95455  
[1] 16.45685  

This is confirmed by
> nb0 <- glm.nb(Days ~ 1, data = quine)    
> exp(coef(nb0))  
(Intercept)  
[1] 16.4589

(apart from rounding errors). 
