HMP (version 2.0.1)

MC.Xmcupo.statistics: Size and Power of Several Sample RAD-Probability Mean Test Comparisons: Unknown Vector of Proportion

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.

Value

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

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) 
	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
# }