HMP (version 2.0.1)

MC.ZT.statistics: Size and Power of Goodness of Fit Test: Multinomial vs. Dirichlet-Multinomial

Description

This Monte-Carlo simulation procedure provides the power and size of the Multinomial vs. Dirichlet-Multinomial goodness of fit test, using the C(\(\alpha\))-optimal test statistics of Kim and Margolin (1992) (t statistics) and the C(\(\alpha\))-optimal test statistics of (Paul et al., 1989).

Usage

MC.ZT.statistics(Nrs, numMC = 10, fit, type = "ha", siglev = 0.05)

Arguments

Nrs

A vector specifying the number of reads/sequence depth for each sample.

numMC

Number of Monte-Carlo experiments. In practice this should be at least 1,000.

fit

A list (in the format of the output of dirmult function) containing the data parameters for evaluating either the size or power of the test.

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.

Value

A vector containing both the size of the test statistics (under "hnull") or power (under "ha") of the test for both the z and t statistics.

Details

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

Examples

# NOT RUN {
	data(saliva) 
	
	### Get a list of dirichlet-multinomial parameters for the data
	fit.saliva <- DM.MoM(saliva) 
	
	### 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
	nrs <- rep(15000, 25)
	
	### Computing size of the test statistics (Type I error)
	pval1 <- MC.ZT.statistics(nrs, numMC, fit.saliva, "hnull") 
	pval1
	
	### Computing Power of the test statistics (Type II error)
	pval2 <- MC.ZT.statistics(nrs, numMC, fit.saliva)
	pval2
# }