# Size and Power of Several Sample RAD-Probability Mean Test Comparisons: Unknown Vector of Proportion 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. ```r 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. ## Returns 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 ```r 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 ```

MC.Xmcupo.statistics function