Poisson Regression Residuals 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.
 A: 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. 


*

*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. 

*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. 

*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. 

*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. 
