# Why does my model predict symmetrical CI for Poisson-distributed data?

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)  Produces: Can anyone explain me why? Many thanks! EDIT 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 abs(y_hat$$fit-y_hat_sem)==abs(y_hat$$fit-y_hat_sep)  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.] ## 1 Answer I think the answer to your question is that se=T, as you seem to know returns the standard error. That is a fixed value, whether you keep the value in logged units or convert it back to the original y-scale. You ensure that it symmetrically surrounds your predicted line when you add and subtract the product of the standard error and 1.96 to create your intervals. 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. yy<- predict(glm1, data.frame(x=fake.df), type="link", se=T) yysep<- exp(yy$$fit+1.96*yy$$se) yysem<- exp(yy$$fit-1.96*yy$$se) y_hat<-predict(glm1, data.frame(x=fake.df), type="response", se=T)$fit


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 abs(y_hat-ysem)==abs(y_hat-yysep). Should all be FALSE.

• Thanks! I added an abridged and slightly edited version of your answer below my question. – Marco Plebani Dec 11 '18 at 9:14
• – Marco Plebani Apr 28 '19 at 8:42