3
$\begingroup$

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:

enter image description here

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):

enter image description here

[Gavin Simpson addresses the issue here.]

$\endgroup$
3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.