Residual diagnostics for negative binomial regression In the residual diagnostics for OLS, I understand what to look to assess any violations (e.g., normality and homoskedasticity of residuals). I was wondering what should one check for in residuals for a negative binomial regression fitted model. How could I do about getting these residual diagnostics in R? 
 A: Check out the DHARMa package in R. It uses a simulation based approach with quantile residuals to generate the type of residuals you may be interested in. And it works with glm.nb from MASS.
The essential idea is explained here and goes in three steps:


*

*Simulate plausible responses for each case. You can use the distribution of each regression coefficient (coefficient together with standard error) to generate several sets of coefficients. You can multiply each set of coefficients by the observed predictors for each case to obtain multiple simulated response values for each case.

*from the multiple response values for each case, generate the empirical cumulative density function (cdf)

*find the value of the empirical cdf that corresponds to the observed response for each case. This is your residual and it ranges from 0 to 1.


If there is no systematic misspecification in your model, all values from the empirical cdf are equally likely. And a histogram of these residuals should be uniformly distributed between 0 and 1.
The package contains additional checks.

Edit:
The steps above are not exactly correct. The biggest difference between my description and what DHARMa does is DHARMa uses the simulate() function in base R, which ignores the variability in the estimated regression coefficients. The Gelman and Hill regression text recommends taking the variability in regression coefficients into account.
A crucial step I forgot to include is: once one has generated the responses, we should place them on the response scale. For example, the predicted variables from logistic regression are logits, so one should place them on the probability scale. Next step would be to randomly generate observed scores using the predicted probabilities. Continuing with the logistic regression example, one can use rbinom() to generate 0-1 outcomes given predicted probabilities.
Additionally, when the responses are integers, such as with binary outcomes or count data models, DHARMa adds random noise to the simulated responses prior to computing the empirical cdf. This step ensures the residuals behave as one would expect in OLS if the model was not misspecified. Without it, you could have a pile up of residuals at 1 if your outcome is binary outcome.
The code for the simulate residuals function in DHARMa is relatively easy to follow for anyone looking to roll their own implementation.
