# Simulating an overdispersed negative binomial distribution

I have RNAseq data in the form of normalized counts. The counts themselves follow an over-dispersed negative binomial, and I would like to generate random data which replicate these distributions. For example, I've subsampled a set of the counts ~4200 n, and the descriptives are as follows,

            Statistic    Std. Error
Mean        46.28        32.6
Lower       -17.64
Upper       110.19
Median      0
Variance    4493154.99
STD         2119.71
Minimum     0
Maximum     136974
Skewness    63.86        0.04
Kurtosis    4123.38      0.08


I've tried the following:

data <- rnbinom(4200,46.28,0.9)


It seems to give a distribution similar to my data, but cuts off the upper bound range. So I'm not exactly sure if this is a robust simulation.

I do apologize for what is likely a rudimentary question. I'm not exactly a skilled statistician, so it's somewhat difficult to follow some of the documentation on these processes. Would anyone be kind enough to explain this to me in the way a novice would understand?

• Welcome Martin. Note that your username, identicon, & a link to your user page are automatically added to every post you make, so there is no need to sign your posts. Jul 25, 2016 at 16:06

Since you have the mean and the variance, the easiest way to do this would be to use the parameterization via the mu and size parameters of rnbinom. Quoting from the help page:

An alternative parametrization (often used in ecology) is by the mean ‘mu’ (see above), and ‘size’, the dispersion parameter, where ‘prob’ = size/(size+mu)’. The variance is ‘mu + mu^2/size’ in this parametrization.

So we can use your mean and variance estimates and solve for size:

size=46.28^2/(4493154.99-46.28))


We can simulate like this:

rnbinom(4200,mu=46.28,size=46.28^2/(4493154.99-46.28))


However...

set.seed(1)
table(rnbinom(4200,mu=46.28,size=46.28^2/(4493154.99-46.28)))

0     1     3     4     7     8     9    10    13    18    19    24    52
4174     1     1     1     1     1     1     1     1     1     1     1     1
68    69    72   115   186   197  1261  3654  3745  4212  5790 10304 14269
1     1     1     1     1     1     1     1     1     1     1     1     1
25467
1


Note especially the right tail. If you simulate this repeatedly, the maximum value and the entire right tail will fluctuate wildly. (I just sampled this five times and got maxima between 14,294 and 186,670.) This is a consequence of (a) enormous overdispersion (your SD is 45 times the mean!), and (b) a large sample size. And of course your estimate of the SD and the overdispersion is largely determined by the right tail of your observations.

Whatever you plan on doing with your data, I'd take a good hard look at it. Maybe perform some sensitivity analysis. What would your result look like if your maximum had been "only" 50,000? And so on.

It may be that whatever conclusions you draw in the course of your analysis are mainly driven by the largest ten of your observations, 0.2% of your sample of size 4200. That is not a very strong position to argue from.

• Thank you, Stephen. Very much appreciated. I believe the reason why the variance is so high is that the counts aggregate based upon the size of a genomic feature they are being assigned too. So say you have alignments which are of a size 50, and two genomic features one of length 50 and another at say 1000. By chance alone, you would expect to see 20x more alignments mapping to the feature of 1000. It's kind of an interesting question which comes from this, albeit it more bioinformatics related, specifically what's the best way of both normalizing and subsampling alignment counts. Jul 25, 2016 at 16:22
• As I alluded to before, I'm currently subsampling by 1000 bps for raw counts, so maybe it would make sense to reduce this even further? The alternative is to use a normalization like alignments per million counts which correct for feature size, though I'm not sure if that would still follow a poisson distribution. Jul 25, 2016 at 16:22
• Then again, it might just make sense to try and exclude outliers in the right tail of the distribution. Sorry I'm not very skilled with this statistics, but would you recommend a specific sensitivity test for this purpose? Jul 25, 2016 at 16:27
• If you have very different feature lengths and these directly impact on your counts, then you are comparing (genomic) apples and oranges, aren't you? I'd recommend running a negative binomial regression of each count on the corresponding length. This may be helpful. I'd recommend bootstrapping this analysis and looking at the distribution of the estimated parameters to get a feeling for how precise they are. Jul 25, 2016 at 16:30
• Yes, it certainly is apples to oranges. But seriously, thank you again. It's much appreciated. I was curious if I could pick your brain a bit more. Essentially, I have two count sets for genes and intergenic features. Unfortunately, sequencing produces noise, so to account for this we subsampled intergenic regions for specific gene features by feature size. We first normalize the counts by feature size, log transform the data to produce a normal distribution, then use the 90th percentile of the subsampled intergenic region as a cutoff for gene features. Jul 25, 2016 at 17:39