MC.Xoc.statistics: Size and Power of Several Sample-Overdispersion Test...

Description Usage Arguments Details Value Examples

View source: R/MC.Xoc.statistics.R

Description

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

Usage

1
	MC.Xoc.statistics(group.Nrs, numMC = 10, group.alphap, type = "ha", siglev = 0.05)

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.

group.alphap

If "hnull": A vector of alpha parameters for each taxa.
If "ha": A list consisting of vectors of alpha parameters for each taxa.

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.

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: All components of group.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
30
31
32
	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.tonsils$gamma
	pval1 <- MC.Xoc.statistics(group.Nrs, numMC, alphap, "hnull")
	pval1
	
	## Not run: 
		### Computing Power of the test statistics (Type II error)
		alphap <- rbind(fit.saliva$gamma, fit.throat$gamma, fit.tonsils$gamma)
		pval2 <- MC.Xoc.statistics(group.Nrs, numMC, alphap, "ha")
		pval2
	
## End(Not run)

Example output

Loading required package: dirmult

Attaching package:HMPThe following object is masked frompackage:dirmult:

    weirMoM

[1] 0.5
[1] 1

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