MC.Xmcupo.statistics: Size and Power of Several Sample RAD-Probability Mean 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 RAD-probability mean test comparisons without reference vector of proportions, using the Generalized Wald-type statistics.

Usage

	MC.Xmcupo.statistics(group.Nrs, numMC = 10, pi0, group.pi, group.theta, 
		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.
pi0	The RAD-probability mean vector.
group.pi	If "hnull": This argument is ignored. If "ha": A matrix where each row is a vector pi values for each group.
group.theta	A vector of overdispersion values for 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.

Details

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

Value

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

Examples

	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
	Nrs1 <- rep(12000, 10)
	Nrs2 <- rep(12000, 19)
	group.Nrs <- list(Nrs1, Nrs2)
	
	group.theta <- c(fit.throat$theta, fit.tonsils$theta)
	pi0 <- fit.saliva$pi
	
	### Computing size of the test statistics (Type I error)
	pval1 <- MC.Xmcupo.statistics(group.Nrs, numMC, pi0, group.theta=group.theta, type="hnull")
	pval1
	
	### Computing Power of the test statistics (Type II error)
	group.pi <- rbind(fit.throat$pi, fit.tonsils$pi)
	pval2 <- MC.Xmcupo.statistics(group.Nrs, numMC, group.pi=group.pi, group.theta=group.theta)
	pval2

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