# MC.Xdc.statistics: Size and Power for the Several-Sample DM Parameter Test... In HMP: Hypothesis Testing and Power Calculations for Comparing Metagenomic Samples from HMP

## Description

This Monte-Carlo simulation procedure provides the power and size of the several sample Dirichlet-Multinomial parameter test comparison, using the likelihood-ratio-test statistics.

## Usage

 ```1 2``` ``` MC.Xdc.statistics(group.Nrs, numMC = 10, alphap, type = "ha", siglev = 0.05, est = "mom") ```

## Arguments

 `group.Nrs` A list specifying the number of reads/sequence depth for each sample in a group with one group per list entry. `numMC` Number of Monte-Carlo experiments. In practice this should be at least 1,000. `alphap` If `"hnull"`: A matrix where rows are vectors of alpha parameters for the reference group. If `"ha"`: A matrix consisting of vectors of alpha parameters for each taxa in each group. `type` If `"hnull"`: Computes the size of the test. If `"ha"`: Computes the power of the test. (default) `siglev` Significance level for size of the test / power calculation. The default is 0.05. `est` The type of parameter estimator to be used with the Likelihood-ratio-test statistics, 'mle' or 'mom'. Default value is 'mom'. (See Note 2 in details)

## Details

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

2. Note 2: 'mle' will take significantly longer time and may not be optimal for small sample sizes; 'mom' will provide a more conservative result in such a case.

3. Note 3: All components of `alphap` should be non-zero or it may result in errors and/or invalid results.

## Value

Size of the test statistics (under `"hnull"`) or power (under `"ha"`) of the test.

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29``` ``` data(saliva) data(throat) data(tonsils) ### Get a list of dirichlet-multinomial parameters for the data fit.saliva <- DM.MoM(saliva) fit.throat <- DM.MoM(throat) fit.tonsils <- DM.MoM(tonsils) ### Set up the number of Monte-Carlo experiments ### We use 1 for speed, should be at least 1,000 numMC <- 1 ### Generate the number of reads per sample ### The first number is the number of reads and the second is the number of subjects nrsGrp1 <- rep(12000, 9) nrsGrp2 <- rep(12000, 11) nrsGrp3 <- rep(12000, 12) group.Nrs <- list(nrsGrp1, nrsGrp2, nrsGrp3) ### Computing size of the test statistics (Type I error) alphap <- fit.saliva\$gamma pval1 <- MC.Xdc.statistics(group.Nrs, numMC, alphap, "hnull") pval1 ### Computing Power of the test statistics (Type II error) alphap <- rbind(fit.saliva\$gamma, fit.throat\$gamma, fit.tonsils\$gamma) pval2 <- MC.Xdc.statistics(group.Nrs, numMC, alphap) pval2 ```

### Example output

```Loading required package: dirmult

Attaching package: 'HMP'

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

weirMoM

 0.5
 1
```

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