Interpreting plot of residuals vs. fitted values from Poisson regression

I am trying to fit data with a GLM (poisson regression) in R. When I plotted the residuals vs the fitted values, the plot created multiple (almost linear with a slight concave curve) "lines". What does this mean?

library(faraway)
modl <- glm(doctorco ~ sex + age + agesq + income + levyplus + freepoor +
freerepa + illness + actdays + hscore + chcond1 + chcond2,
family=poisson, data=dvisits)
plot(modl)


-
 I don't know if you can upload the plot (sometimes newcomers can't), but if not, could you at least add some data & R code to your question so people can evaluate it? – gung Mar 22 '12 at 15:46 Jocelyn, I've updated your post with information you put in a comment. I also tagged this as homework since you talked about an assignment. – chl♦ Mar 22 '12 at 16:20 Thanks so much! – jocelyn Mar 22 '12 at 16:22 try plot(jitter(mod1)) to see if the graph is a bit more readable. Why don't you define residuals for us and give us your best guess as interpreting the graph yourself. – Michael Bishop Mar 22 '12 at 16:45 Also, look up the poisson distribution and then graph your outcome variable. plot(doctorco) - it should be strictly positive. – Michael Bishop Mar 22 '12 at 16:49

This is the appearance you expect of such a plot when the dependent variable is discrete.

Each curvilinear trace of points on the plot corresponds to a fixed value $k$ of the dependent variable $y$. Every case where $y=k$ has a prediction $\hat{y}$; its residual--by definition--equals $k-\hat{y}$. The plot of $k-\hat{y}$ versus $\hat{y}$ is obviously a line with slope $-1$. In Poisson regression, the x-axis is shown on a log scale: it is $\log(\hat{y})$. The curves now bend down exponentially. As $k$ varies, these curves rise by integral amounts. Exponentiating them gives a set of quasi-parallel curves. (To prove this, the plot will be explicitly constructed below, separately coloring the points by the values of $y$.)

We can reproduce the plot in question quite closely by means of a similar but arbitrary model (using small random coefficients):

# Create random data for a random model.
set.seed(17)
n <- 2^12                       # Number of cases
k <- 12                         # Number of variables
beta = rnorm(k, sd=0.2)         # Model coefficients
x <- matrix(rnorm(n*k), ncol=k) # Independent values
y <- rpois(n, lambda=exp(-0.5 + x %*% beta + 0.1*rnorm(n)))

# Wrap the data into a data frame, create a formula, and run the model.
df <- data.frame(cbind(y,x))
vars <- apply(matrix(1:k, nrow=1), 1, function(i) paste("V", i+1, sep=""))
s.formula <- paste("y ~", paste(s.formula, collapse="+"))
modl <- glm(as.formula(s.formula), family=poisson, data=df)

# Construct a residual vs. prediction plot.
b <- coefficients(modl)
y.hat <- x %*% b[-1] + b[1]     # *Logs* of the predicted values
y.res <- y - exp(y.hat)         # Residuals
colors <- 1:(max(y)+1)          # One color for each possible value of y
plot(y.hat, y.res, col=colors[y+1], main="Residuals v. Fitted")


-
(+1) The color goes a long way in showing what is happening. – cardinal Jun 11 '12 at 23:11

Sometimes stripes like these in residual plots represent points with (almost) identical observed values that get different predictions. Look at your target values: how many unique values are they? If my suggestion is correct there should be 9 unique values in your training data set.

-
 +1. (There's actually a tenth one represented as a single point near the upper right corner.) The values, of course, are $0, 1, \ldots, 9$. – whuber♦ Jun 11 '12 at 20:54

This pattern is characteristic of an incorrect match of the family and/or link. If you have overdispersed data then perhaps you should consider the negative binomial (count) or gamma (continuous) distributions. Also you should be plotting your residuals against the transformed linear predictor, not the predictors when using generalized linear models. To transform the Poisson predictor you need to take 2 times the square root of the linear predictor and plot your residuals against that. The residuals further more should not be exclusively pearson residuals, try deviance residuals and studentized resids.

-
Why 2 times the square root, when the canonical link of the poisson family in a glm is log? Shouldn't it be exp() of the linear predictor? But I don't see what the problem is with plotting residuals against the linear predictor itself, which I think is what is being done here - perhaps you could expand on that. – Peter Ellis Jun 11 '12 at 19:40
Would you mind explaining just what aspect of the "pattern" is drawing your attention to a possible model mis-specification, Ryan? It seems to be a subtle thing, but is potentially an important insight. – whuber Jun 12 '12 at 11:49