Sampling from under/over-dispersed count data in R

Question

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 isn`t appropriate since it assumes that the variance is necessarily larger than the mean — Gabriel LM, Feb 13 '17 at 07: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 05 '20 at 11:24
Quasilikelihood is not a probability model, only a weighting scheme to get proper parameter estimates. There exist alternatives to the Negative Binomial which allow underdispersion like the Gamma-count model (https://arxiv.org/pdf/1312.2423.pdf) — Maximilian Aigner, Mar 06 '21 at 11:34

score 1 · Answer 1 · answered Jul 06 '21 at 07:57

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

score 0 · Answer 2 · 2017-02-10T12:04:00.237

0

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

edited Feb 10 '17 at 12:04

answered Feb 10 '17 at 11:12

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
1

@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. – 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

Sampling from under/over-dispersed count data in R

2 Answers2

Linked