Sampling from under/over-dispersed count data in R I am currently working a some datasets with count data in R, in which the response is the number of activities of a given type that were performed in one day by a population. 
For each type, I build a Poisson model and test for over/underdispersion using the function dispersiontest() from package AER. Depending on the result, I switch to quasi-Poisson model when there is evidence of over- or underdispersion.
In a next step, I would like to sample and generate simulated data using the results of my models. If I have a Poisson model, I can use rpois() with lambda being the fitted value of the model. However, I have no idea how to do it in the cases of over/underdispersion. Any ideas?
 A: Simulating on the continuum from not overdispersed to over-dispersed is easy: you fix a rate $\lambda$ (no overdispersion) or you draw your rates from a gamma distribution (=overdispersion). This is one way of getting the negative binomial distribution. For example, if you want the negative binomial distribution with the parameterization that has a mean rate per time unit of $\mu$, a dispersion parameter $\kappa \in [0, \infty)$ ($\kappa=0$ = Poisson) and an observation time $t>0$, you can use functions like these:
# Simulate from NegBin(mu*t, kappa); mu>0, t>0, kappa>0
rnegbin = function(n, mu, kappa, t=1){ # Note: no protection against kappa <= 0
  rpois(n=n, mu*t*rgamma(n=n, shape=1/kappa, rate=1/kappa))
}

# Density function for NegBin(mu, kappa); mu>0, kappa>=0
dnegbin = function(x, mu, kappa, log=FALSE){
  if (kappa>0) {
    tmp = lgamma(x+1/kappa) - lgamma(x+1) - lgamma(1/kappa) + x*log(kappa*mu) - (x+1/kappa)*log1p(kappa*mu) 
  } else { # no protection against kappa <0
    tmp = x*log(mu)-mu-lgamma(x+1)
  }
  if (log==F) return(exp(tmp)) else return(tmp)
}

If you do not really mind whether your overdispersed counts are from a negative binomial, you can of course also draw your rates from a log-normal distribution (or log-rates from a normal distribution).
I don't have much experience with simulating underdispersed counts, but one I idea would be to simulate Poisson random variables and to draw a new set with increasing probability the further away from the mean you are. It's a bit tricky to keep the expected value constant though (if that matters to you).
A: I would use the predict() function with the newdata parameter.
For simulation, we would need to create a random data set of predictors. Let's say that your model has 2 explanatory variables: a factor and a continuous variable. Some pseudo code could be:
set.seed(1) # Reproducibility
simfac <- sample(c("level1","level2"),100,replace=T) # Simulate 100 factor observations with 2 levels
simcont <- runif(100,0,1) # Simulate 100 continuous uniforms between 0 and 1
dat <- data.frame(simfac, simcont)
pred <- predict(model, newdata=dat, type="response") # Predict using the new data

