I have made a poisson regression and then visualised the model:
library(ggplot2)
year <- 1990:2010
count <- c(29, 8, 13, 3, 20, 14, 18, 15, 10, 19, 17, 18, 24, 47, 52, 24, 25, 24, 31, 56, 48)
df <- data.frame(year, count)
my_glm <- glm(count ~ year, family = "poisson", data = df)
my_glm$model$fitted <- predict(my_glm, type = "response")
ggplot(my_glm$model) + geom_point(aes(year, count)) + geom_line(aes(year, fitted))
Now I want to add these simulated Poisson distributions to the plot:
library(plyr)
poisson_sim <- llply(my_glm$model$fitted, function(x) rpois(100, x))
The plot should look something like this (red points are photoshopped):