# 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?

• Have you thought of swapping to a negative binomial model and then sampling using that distribution? – mdewey Feb 10 '17 at 13:45
• @mdewey this works fine indeed in the cases of overdispersion. However, sometimes I have underdispersion as well, and then the negative binomial isnt appropriate since it assumes that the variance is necessarily larger than the mean – Gabriel LM Feb 13 '17 at 7:30
• It sounds like this problem could be approached fully non-parametrically. This way you wouldn't have to make any assumptions about the superpopulation distribution, and you'd get error bars aroud the simulated data. – pglpm Jul 5 at 11:24

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
`
• Thanks for the reply! The prediction does work for the model, no matter what strategy I use it. However, the problem is that the result is a continuous variable, which will be equal to E[Y|covariates]. However, since I have count data, I would to draw only discrete values so it would make sense. An analogy: if I draw from a Poisson distribution with a non-integer Lambda, I will still get integer values – Gabriel LM Feb 10 '17 at 11:46
• @GabrielLM: The issue might be that in the predict function, you need to use type="response" (I've updated my answer). This applies the link function to the predictor. – user64106 Feb 10 '17 at 12:05
• the type="response" will apply the link function indeed, but the result will be the expected response given the value for the covariates in the new data. However, it still does not solve the problem of drawing from the discrete distribution, as the outcome will be a continuous. Let me give an example: My mean is 1.1 and variance is 0.1, thus most responses are 1. Using the predict function, all resulting values will be close to 1.1. What I would like is to be able to draw plausible data from it, i.e., most individuals would have a response of 1, but a few would have 0 or 2, etc. – Gabriel LM Feb 10 '17 at 12:51