Dirichlet.multinomial: Generation of Dirichlet-Multinomial Random Samples

In HMP: Hypothesis Testing and Power Calculations for Comparing Metagenomic Samples from HMP

Description

Random generation of Dirichlet-Multinomial samples.

Usage

Dirichlet.multinomial(Nrs, shape)

Arguments

Nrs	A vector specifying the number of reads or sequence depth for each sample.
shape	A vector of Dirichlet parameters for each taxa.

Details

The Dirichlet-Multinomial distribution is given by (Mosimann, J. E. (1962); Tvedebrink, T. (2010)),

\textbf{P}≤ft ({\textbf{X}_i}=x_{i};≤ft \{ π_j \right \},θ\right )=\frac{N_{i}!}{x_{i1} !,…,x_{iK} !}\frac{∏_{j=1}^K ∏_{r=1}^{x_{ij}} ≤ft \{ π_j ≤ft ( 1-θ \right )+≤ft ( r-1 \right )θ\right \}}{∏_{r=1}^{N_i}≤ft ( 1-θ\right )+≤ft ( r-1 \right) θ}

where \textbf{x}_{i}= ≤ft [ x_{i1}, …, x_{iK} \right ] is the random vector formed by K taxa (features) counts (RAD vector), N_{i}= ∑_{j=1}^K x_{ij} is the total number of reads (sequence depth), ≤ft\{ π_j \right\} are the mean of taxa-proportions (RAD-probability mean), and θ is the overdispersion parameter.

Note: Though the test statistic supports an unequal number of reads across samples, the performance has not yet been fully tested.

Value

A data matrix of taxa counts where the rows are samples and columns are the taxa.

References

Mosimann, J. E. (1962). On the compound multinomial distribution, the multivariate β-distribution, and correlations among proportions. Biometrika 49, 65-82.
Tvedebrink, T. (2010). Overdispersion in allelic counts and theta-correction in forensic genetics. Theor Popul Biol 78, 200-210.

Examples

	data(saliva)
	
	### Generate a the number of reads per sample
	### The first number is the number of reads and the second is the number of subjects
	nrs <- rep(15000, 20) 
	
	### Get gamma from the dirichlet-multinomial parameters
	shape <- dirmult(saliva)$gamma
	
	dmData <- Dirichlet.multinomial(nrs, shape)
	dmData[1:5, 1:5]

Example output

Loading required package: dirmult

Attaching package: 'HMP'

The following object is masked from 'package:dirmult':

    weirMoM

Iteration 1: Log-likelihood value: -1426219.55743915
Iteration 2: Log-likelihood value: -1426159.69392069
Iteration 3: Log-likelihood value: -1426137.57581644
Iteration 4: Log-likelihood value: -1426134.36208055
Iteration 5: Log-likelihood value: -1426134.28420196
Iteration 6: Log-likelihood value: -1426134.28414898
         Taxa 1 Taxa 2 Taxa 3 Taxa 4 Taxa 5
Sample 1   2720   1658   1347   1875    510
Sample 2   2996   2062   1515   1389    856
Sample 3   3358   2525   1877   1272    972
Sample 4   3024   1866   1676   1168    926
Sample 5   2609   2113   1385   1455   1057

HMP documentation built on Aug. 31, 2019, 5:05 p.m.