Poisson Regression Residuals

Question

I'm modeling the number of doctor visits (a count variable) on factors such as income, chronic condition, insurance, etc. I use the canned Stata command poisson y x1 x2, etc. I have the following two questions:

1) How come Stata does not allow me to execute predict resids, residuals after I've ran the above regression? I then tried to generate the residuals using the formula "y-yhat", yet I got a residual with a negative average mean (-6.96e-09...) - any thoughts here?

2) Should I be using robust SE (i.e. vce(robust)) here? I do think there can be heteroskedasticity issues since the number of doctor visits may vary highly among different values of the income variable; however, I'm not sure if robust SEs can be applied as freely in the case of nonlinear models such as Poisson.

Answer 1

There's a bundle of different questions here, but with enough statistical content to allow an answer. Note that Stata-focused questions may go better on Statalist or (if about programming, not the case here) on Stack Overflow.

1. Calculation of residuals is indeed not supported after poisson: see on-line version of the pertinent help. In the absence of an obvious reason for that, it does seem like a quirk. In addition to direct calculation, note that glm, f(poisson) fits essentially the same models and does allow calculation of residuals through predict afterwards.

2. Your concern about negative averages for residuals is unwarranted: -6.96e-09 is another way of writing -0.00000000696 which shows that the average is essentially zero (given that the variable doctor visits is presumably integer-valued). The chance of getting exactly zero for mean residual is negligible, but the non-zero part is to be thought as round-off error. If you use double as a variable type for the residuals, it is likely that you will get even closer to zero.

3. Whether you should use robust standard errors depends on the characteristics of your data. Note that the response variable varying highly among different values of a predictor is not in itself heteroscedasticity. In addition, I think you are confusing heteroscedasticity with over-dispersion. But in principle there is no reason not to try this as an easy alternative. For a positive recommendation of this as a general approach, see this blog post by William Gould.

4. Your comment about nonlinear models is a little cryptic, but Poisson models qualify as generalized linear models, so whatever the line of argument is, it should not give you cause for concern.