Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

Goal: I want to visually depict the difference between experimental treatments while controlling for a continuous covariate in Poisson regression.

In my experiment, I am seeking to explain the number of species of insect on plants (richness) with a 3-level treatment factor and the size of the plant as a continuous covariate. My model is:

glm(rich ~ tmnt + size, family=poisson)

Furthermore, I have contrast coded my treatment factor to compare treatments 1&2 to treatment 3, using the following code:

1.2vs3 <- cbind(1vs2.3=c(1,1,2), 2vs3=c(-1,1,0))
contrasts(tmnt) <- 1vs2.3

(However, I don't think my contrast coding scheme will make a difference for my question here.)

I'm looking to make a bar plot of the species richness across treatments while controlling for the covariate. In multiple regression with Gaussian errors and two continuous variables, one would create a partial plot, by graphing the following residuals:

richRes <- residuals(lm(rich ~ size))
contRes <- residuals(lm(cont ~ size))
plot(contRes, richRes)

And the following two models would produce parameter estimates for the treatment effect that are identical:

lm(rich    ~ cont  +  size)
lm(richRes ~ contRes      )

However, I don't know how to extract values for plotting the effect of a factor when controlling for a covariate, and I also think that the Poisson errors make things more complicated. I thought that I could simply take the residuals of the full model lacking the tmnt and plot them against tmnt like this:

richRes <- residuals(glm(rich ~ size, family=poisson)
plot(tmnt, richRes)

which works, but the models:

glm(rich ~ tmnt + size, family=poisson) 
lm(richRes ~ tmnt, family=poisson)

give different parameter estimates, so something is not quite right.

Can anybody please provide me with a procedure that will allow me to plot richness against treatment in way that accurately reflects the difference between and variance within (error bars) treatments from the full Poisson regression, i.e., accounting for size?

Here are some made up data with which to play, if you so choose. These variables will fit into my code (contrasts excepted):

rich <- rpois(50, 10)
tmnt <- as.factor(c(rep(1, 25), rep(2, 25)))
size <- rnorm(50, 10, 3)
cont <- rnorm(50, 10, 3)
share|improve this question
+1, clear, well thought-out question. If you just want a bar chart, you can solve for the predicted mean response for each of the 3 treatment conditions at the mean of the covariate. – gung Oct 12 '12 at 22:05

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.