Given a set of Poisson-distributed data, I would expect the confidence interval for a
glm(..., family="poisson") to be asymmetrical, because of Poisson distributions being right-skewed. Instead, if I simulate the following:
# Simulate Poisson-distributed data: #sample size n <- 100 #regression coefficients beta0 <- 0.001 beta1 <- 0.2 #generate covariate values set.seed(1) x <- runif(n=n, min=0, max=15) #compute mu's mu <- exp(beta0 + beta1 * x) #generate Y-values y <- rpois(n=n, lambda=mu)
Modelling the simulated data according to:
glm1 <- glm(y~x, data=dd, family="poisson")
And predicting the expected values of y and their associated CI as follows:
# Set a sequence of hypothetic values for the explanatory variable: fake.df <- data.frame(x =seq(min(dd$x), max(dd$x), 0.01)) # This just helps keeping the following lines of code short and clean: fake.x<- fake.df$x # Predict y values and their standard errors, given fake.x: yy<- predict(glm1, data.frame(x=fake.df), type="response", se=T) # Again, this just helps keeping the following lines of code short and clean: yysep<- yy$fit+1.96*yy$se yysem<- yy$fit-1.96*yy$se # show our regression line using the bloody predict() function: lines(fake.x, predict(glm1, data.frame(x=fake.df), type="response"), col="black", lty="solid", lwd=2 ) # Confidence intervals: # lines(fake.x, yysep, lty="dashed", col="red") # lines(fake.x, yysem, lty="dashed", col="red") # Or as a shaded area: polygon(c(fake.x,rev(fake.x)),c(yysep,rev(yysem)),col="#0000ff33", border=NA)
Can anyone explain me why? Many thanks!
The following is a shortened, operational version of the answer provided by Matt Barstead, with his code modified to be consistent with mine. For details on the rationale behind it, see his answer in the Answers section.
If you wanted to include the asymmetry of the Poisson distribution into your confidence intervals you need to calculate your confidence intervals in the logged units first.
y_hat <- predict(glm1, data.frame(x=fake.df), type="link", se=T) y_hat_sep <- exp(y_hat$fit+1.96*y_hat$se) y_hat_sem <- exp(y_hat$fit-1.96*y_hat$se) polygon(c(fake.x,rev(fake.x)),c(y_hat_sep,rev(y_hat_sem)),col="#FF000050", border=NA)
Exponentiating out of the log transform after you have created a 95% CI (which should be symmetrically distributed around your linear prediction fit line), will yield an asymmetrical 95% confidence interval, though it may not be visually obvious when the error is sufficiently small. You can confirm that the interval is not symmetrical though by running
Should all be FALSE.
Here I zoom in on a section of the graph shown above to show the difference between CI based on the SE calculated by
predict(..., type="response") and the CI obtained by first computing it on the link scale and then exponentiating it (blue shade and red shade, respectively):
[Gavin Simpson addresses the issue here.]