Plot confidence interval (IC) in multiple glm Poisson regression

Question

I've like to plot confidence interval (IC) in multiple glm Poisson regression without sucess. If I make:

#Artificial data set
Consumption <- c(501, 502, 503, 504, 26, 27, 55, 56, 68, 69, 72, 93)
Gender <- gl(n = 2, k = 6, length = 2*6, labels = c("Male", "Female"), ordered = FALSE)
Income <- c(5010, 5020, 5030, 5040, 260, 270, 550, 560, 680, 690, 720, 930)
df3 <- data.frame(Consumption, Gender, Income)
df3

Create glm Poisson multiple regression model:

fm1 <- glm(Consumption~Gender+Income, data=df3, family=poisson)
summary(fm1)

See the significance of variables:

# ANOVA
anova(fm1,test="Chi")

Predict values and calculate IC:

df3 = cbind(df3, pred = predict(fm1, type = "response"))#Estimate values
df3 = cbind(df3, se = predict(fm1, type="link",se.fit = TRUE)) ## Confidence interval of estimated
df3 = cbind(df3, ucl=exp(df3$pred + 1.96*df3$se.fit))
df3 = cbind(df3, lcl=exp(df3$pred - 1.96*df3$se.fit))

Now if I try to plot this:

#Plot
ggplot(data=df3, mapping=aes(x=Income, y=Consumption, color=Gender)) + 
  geom_point() +  
  geom_line(mapping=aes(y=pred)) +
  geom_smooth(data=df3, aes(ymin = lcl, ymax = ucl), stat="identity") 

#

Doesn't work:

Error in if ((w[1] * sm + w[2] * cm + w[3] * dm + w[4]) < best$score) break : 
  missing value where TRUE/FALSE needed

Any member could help me?

Thanks in advance!

Answer 1

I suggest to calculate the confidence intervals of the predicted values using the following code:

pred <- predict(fm1, type="link", se.fit = TRUE)
df3 = cbind(df3, pred = pred$fit)
df3 = cbind(df3, se = pred$se.fit) 
df3 = cbind(df3, ucl=exp(df3$pred + 1.96*df3$se))
df3 = cbind(df3, lcl=exp(df3$pred - 1.96*df3$se))

or

pred <- predict(fm1, type="response", se.fit = TRUE)
df3 = cbind(df3, pred = pred$fit)
df3 = cbind(df3, se = pred$se.fit) 
df3 = cbind(df3, ucl=df3$pred + 1.96*df3$se)
df3 = cbind(df3, lcl=df3$pred - 1.96*df3$se)