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)
  • 2
    $\begingroup$ +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. $\endgroup$ Oct 12, 2012 at 22:05


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.