Dirichlet.multinomial: Generation of Dirichlet-Multinomial Random Samples

Description Usage Arguments Details Value Author(s) References Examples

View source: R/Dirichlet.multinomial.R

Description

Random generation of Dirichlet-Multinomial samples.

Usage

1

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.

Author(s)

Patricio S. La Rosa, Elena Deych, Berkley Shands, William D. Shannon

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
	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 Dec. 23, 2017, 5:40 p.m.