Fitting a zero-inflated negative binomial regression with R

In this thread, I laid out a problem involving fitting a model that attempts to use minor league baseball statistics to predict success at the major league level (explained in full in the thread). After doing further research outside of the thread, I have come to the conclusion that a zero-inflated negative binomial model is likely the best fit given that I believe there are two processes generating the data. The first process determines whether a player will make the majors and once the player reaches the majors, a second process governs their success (as measured by WAR - also explained in the linked thread).

I ran the model and the ZINB model appears to be a reasonable fit given the following diagnostic plot of fitted values vs residuals (broken into two plots to make it easier to inspect visually).

EDIT 2: Here is a plot of the Pearson residuals.

EDIT 3: Here is a plot of the Pearson residuals vs the fitted values.

Although this is a big improvement over my previous models, there is clearly a bias for the model to underestimate a player's career WAR i.e., a majority of the residuals are greater than zero. EDIT: It turns out that the plot is misleading due to the high number of overlapping residuals. In reality, only 10% of the residuals are greater than zero. I am guessing that this may be improved by either a) a better specification of the functional form of the covariates or b) using a different yet similar model e.g., ZIP. Given that this type of regression goes well beyond what I have worked with before, I would appreciate any suggestions on further diagnostics I can use to both test this model and compare it to others and how to improve the functional form of the covariates given that the R function I have used, zeroinfl (from the pscl package), does not appear to be that flexible. Thank you!

• "a majority of the residuals >0" of itself doesn't automatically imply bias (it does if mean = median, for example); it probably does in this case because we'd actually expect a majority <0. What is the proportion of residuals >0? – Glen_b Jun 14 '13 at 1:51
• What kind of residuals are these? What code in R did you use to generate them? – Gavin Simpson Jun 14 '13 at 2:37
• @GavinSimpson I calculated these residuals myself as Y - fitted values – zgall1 Jun 14 '13 at 11:39
• @Glen_b It turns out that residual plot is misleading due to the high number of overlapping residuals. Only 10% of the residuals are actually > 0 – zgall1 Jun 14 '13 at 12:40
• So those are raw residuals. You would expect patterns in those residuals wouldn't you? The model you've fitted allows for them (variance increasing with mean). You'll also get bands with discrete data. You might like to look at the ?residuals.zeroinfl and use the resid() function to get the Pearson residuals. I would probably now predict from the model to see what sorts of predictions you get and how they relate to the observed data. Also start to think about whether you are missing key terms or functions by plotting residuals versus the predictors etc and look for pattern. – Gavin Simpson Jun 14 '13 at 15:09