| Title: | Power Analysis and Sample Size Calculation |
|---|---|
| Description: | Power analysis and sample size calculation for Welch, Hsu (Hedderich and Sachs (2018), ISBN:978-3-662-56657-2) and Xiao (Xiao (2018), <doi:10.17654/TS054010021>) as well as permutation (Janssen (1997), <doi:10.1016/S0167-7152(97)00043-6>) and bootstrap (Davison and Hinkley (1997), ISBN:978-0-511-80284-3) t-tests including Monte-Carlo simulations of empirical power and type-I-error. Monte-Carlo simulations of empirical power and type-I-error for conditional two-sample tests. Power and sample size calculation for Wilcoxon rank sum and signed rank tests via Monte-Carlo simulations. Power and sample size required for the evaluation of a diagnostic test(-system) (Flahault et al. (2005), <doi:10.1016/j.jclinepi.2004.12.009>; Dobbin and Simon (2007), <doi:10.1093/biostatistics/kxj036>) as well as for a single proportion (Fleiss et al. (2003), ISBN:978-0-471-52629-2; Piegorsch (2004), <doi:10.1016/j.csda.2003.10.002>; Thulin (2014), <doi:10.1214/14-ejs909>), comparing two negative binomial rates (Zhu and Lakkis (2014), <doi:10.1002/sim.5947>), ANCOVA (Shieh (2020), <doi:10.1007/s11336-019-09692-3>), reference ranges (Jennen-Steinmetz and Wellek (2005), <doi:10.1002/sim.2177>), multiple primary endpoints (Sozu et al. (2015), ISBN:978-3-319-22005-5), and AUC (Hanley and McNeil (1982), <doi:10.1148/radiology.143.1.7063747>). |
| Authors: | Matthias Kohl [aut, cre] (0000-0001-9514-8910) |
| Maintainer: | Matthias Kohl <[email protected]> |
| License: | LGPL-3 |
| Version: | 1.2 |
| Built: | 2026-06-04 08:11:13 UTC |
| Source: | https://github.com/stamats/mkpower |
Power analysis and sample size calculation for Welch and Hsu (Hedderich and Sachs (2018), ISBN:978-3-662-56657-2) t-tests including Monte-Carlo simulations of empirical power and type-I-error. Power and sample size calculation for Wilcoxon rank sum and signed rank tests via Monte-Carlo simulations. Power and sample size required for the evaluation of a diagnostic test(-system) (Flahault et al. (2005), <doi:10.1016/j.jclinepi.2004.12.009>; Dobbin and Simon (2007), <doi:10.1093/biostatistics/kxj036>) as well as for a single proportion (Fleiss et al. (2003), ISBN:978-0-471-52629-2; Piegorsch (2004), <doi:10.1016/j.csda.2003.10.002>; Thulin (2014), <doi:10.1214/14-ejs909>), comparing two negative binomial rates (Zhu and Lakkis (2014), <doi:10.1002/sim.5947>), ANCOVA (Shieh (2020), <doi:10.1007/s11336-019-09692-3>), reference ranges (Jennen-Steinmetz and Wellek (2005), <doi:10.1002/sim.2177>), multiple primary endpoints (Sozu et al. (2015), ISBN:978-3-319-22005-5), and AUC (Hanley and McNeil (1982), <doi:10.1148/radiology.143.1.7063747>).
library(MKpower)
Matthias Kohl https://www.stamats.de
Maintainer: Matthias Kohl [email protected]
Produce histograms for simulations of power and type-I-error of tests.
## S3 method for class 'sim.power.ttest' hist(x, color.hline = "orange", ...) ## S3 method for class 'sim.power.wtest' hist(x, color.hline = "orange", ...)## S3 method for class 'sim.power.ttest' hist(x, color.hline = "orange", ...) ## S3 method for class 'sim.power.wtest' hist(x, color.hline = "orange", ...)
x |
object of class |
color.hline |
color of horizontal line indicating uniform distribution of p values. |
... |
further arguments that may be passed through). |
The plot generates a ggplot2 object that is shown.
Missing values are handled by the ggplot2 functions.
Object of class gg and ggplot.
Matthias Kohl [email protected]
res1 <- sim.power.t.test(nx = 5, rx = rnorm, rx.H0 = rnorm, ny = 10, ry = function(x) rnorm(x, mean = 3, sd = 3), ry.H0 = function(x) rnorm(x, sd = 3), methods = c("student", "welch", "hsu")) hist(res1) res2 <- sim.power.wilcox.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm) hist(res2) res3 <- sim.power.cond.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm) hist(res3)res1 <- sim.power.t.test(nx = 5, rx = rnorm, rx.H0 = rnorm, ny = 10, ry = function(x) rnorm(x, mean = 3, sd = 3), ry.H0 = function(x) rnorm(x, sd = 3), methods = c("student", "welch", "hsu")) hist(res1) res2 <- sim.power.wilcox.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm) hist(res2) res3 <- sim.power.cond.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm) hist(res3)
Compute sample size for ANCOVA.
power.ancova(n = NULL, mu = NULL, var = 1, nr.covs = 1L, group.ratio = NULL, contr.mat = NULL, sig.level = 0.05, power = NULL, n.max = 1000L, rel.tol = .Machine$double.eps^0.25)power.ancova(n = NULL, mu = NULL, var = 1, nr.covs = 1L, group.ratio = NULL, contr.mat = NULL, sig.level = 0.05, power = NULL, n.max = 1000L, rel.tol = .Machine$double.eps^0.25)
n |
vector of sample sizes per groups. |
mu |
vector of mean values of the groups. |
var |
error variance. |
nr.covs |
number of covariates (larger or equal than 1). |
group.ratio |
vector of group sizes relative to group 1; i.e., first
entry should always be one. If |
contr.mat |
matrix of contrasts (number of columns must be idential to
number of groups). If |
sig.level |
significance level (type I error probability) |
power |
power of test (1 minus type II error probability) |
n.max |
maximum sample size considered in the computations. |
rel.tol |
relative tolerance passed to function |
Exactly one of the parameters n and power must be passed as
NULL, and that parameter is determined from the other.
The function includes an implementation of the exact approach of Shieh (2020).
It is based on the code provided in the supplement of Shieh (2020), but
uses integrate instead of the trapezoid rule and uniroot for
finding the required sample size.
Object of class "power.htest", a list of the arguments
(including the computed one) augmented with a note element.
Matthias Kohl [email protected]
G. Shieh (2020). Power Analysis and Sample Size Planning in ANCOVA Designs. Psychometrika 85:101-120. doi:10.1007/s11336-019-09692-3.
S.E. Maxwell and H.D. Delaney (2004). Designing experiments and analyzing data: A model comparison perspective (2nd ed.). Mahwah, NJ: Lawrence Erlbaum Associates.
power.anova.test, power.t.test
## Default matrix of contrasts ## 3 groups cbind(rep(1,2), -diag(2)) ## 4 groups cbind(rep(1,3), -diag(3)) ## Table 1 in Shieh (2020) power.ancova(mu=c(400, 450, 500), var = 9900, power = 0.8) power.ancova(n = rep(63/3, 3), mu=c(400, 450, 500), var = 9900) power.ancova(mu=c(400, 450, 500), var = 9900, power = 0.8, nr.covs = 10) power.ancova(n = rep(72/3, 3), mu=c(400, 450, 500), var = 9900, nr.covs = 10) ## Table 2 in Shieh (2020) power.ancova(mu=c(400, 450, 500), var = 7500, power = 0.8) power.ancova(n = rep(48/3, 3), mu=c(400, 450, 500), var = 7500) power.ancova(mu=c(400, 450, 500), var = 7500, power = 0.8, nr.covs = 10) power.ancova(n = rep(60/3, 3), mu=c(400, 450, 500), var = 7500, nr.covs = 10) ## Table 3 in Shieh (2020) power.ancova(mu=c(400, 450, 500), var = 1900, power = 0.8) power.ancova(n = rep(18/3, 3), mu=c(400, 450, 500), var = 1900) power.ancova(mu=c(400, 450, 500), var = 1900, power = 0.8, nr.covs = 10) power.ancova(n = rep(27/3, 3), mu=c(400, 450, 500), var = 1900, nr.covs = 10) ## ANOVA approach for Table 1-3 power.anova.test(groups = 3, between.var = var(c(400, 450, 500)), within.var = 10000, power = 0.8) power.anova.test(n = 63/3, groups = 3, between.var = var(c(400, 450, 500)), within.var = 10000) ## Table 4 in Shieh (2020) power.ancova(mu=c(410, 450, 490), var = 9900, power = 0.8) power.ancova(n = rep(96/3, 3), mu=c(410, 450, 490), var = 9900) power.ancova(mu=c(410, 450, 490), var = 9900, power = 0.8, nr.covs = 10) power.ancova(n = rep(105/3, 3), mu=c(410, 450, 490), var = 9900, nr.covs = 10) ## Table 5 in Shieh (2020) power.ancova(mu=c(410, 450, 490), var = 7500, power = 0.8) power.ancova(n = rep(72/3, 3), mu=c(410, 450, 490), var = 7500) power.ancova(mu=c(410, 450, 490), var = 7500, power = 0.8, nr.covs = 10) power.ancova(n = rep(84/3, 3), mu=c(410, 450, 490), var = 7500, nr.covs = 10) ## Table 6 in Shieh (2020) power.ancova(mu=c(410, 450, 490), var = 1900, power = 0.8) power.ancova(n = rep(24/3, 3), mu=c(410, 450, 490), var = 1900) power.ancova(mu=c(410, 450, 490), var = 1900, power = 0.8, nr.covs = 10) power.ancova(n = rep(33/3, 3), mu=c(410, 450, 490), var = 1900, nr.covs = 10) ## ANOVA approach for Table 4-6 power.anova.test(groups = 3, between.var = var(c(410, 450, 490)), within.var = 10000, power = 0.8) power.anova.test(n = 96/3, groups = 3, between.var = var(c(410, 450, 490)), within.var = 10000) ############################################################################### ## Example from Maxwell and Delaney (2004) according to Shieh (2020) ############################################################################### ## ANCOVA (balanced design) power.ancova(n = rep(30/3, 3), mu=c(7.5366, 11.9849, 13.9785), var = 29.0898) power.ancova(mu=c(7.5366, 11.9849, 13.9785), var = 29.0898, power = 0.8) power.ancova(mu=c(7.5366, 11.9849, 13.9785), var = 29.0898, power = 0.9) ## ANOVA power.anova.test(n = 30/3, groups = 3, between.var = var(c(7.5366, 11.9849, 13.9785)), within.var = 29.0898) power.anova.test(groups = 3, between.var = var(c(7.5366, 11.9849, 13.9785)), within.var = 29.0898, power = 0.8) power.anova.test(groups = 3, between.var = var(c(7.5366, 11.9849, 13.9785)), within.var = 29.0898, power = 0.9) ## ANCOVA - imbalanced design power.ancova(mu=c(7.5366, 11.9849, 13.9785), var = 29.0898, power = 0.8, group.ratio = c(1, 1.25, 1.5)) power.ancova(n = c(13, 16, 19), mu=c(7.5366, 11.9849, 13.9785), var = 29.0898, group.ratio = c(1, 1.25, 1.5)) power.ancova(mu=c(7.5366, 11.9849, 13.9785), var = 29.0898, power = 0.8, group.ratio = c(1, 0.8, 2/3)) power.ancova(n = c(17, 14, 12), mu=c(7.5366, 11.9849, 13.9785), var = 29.0898, group.ratio = c(1, 0.8, 2/3))## Default matrix of contrasts ## 3 groups cbind(rep(1,2), -diag(2)) ## 4 groups cbind(rep(1,3), -diag(3)) ## Table 1 in Shieh (2020) power.ancova(mu=c(400, 450, 500), var = 9900, power = 0.8) power.ancova(n = rep(63/3, 3), mu=c(400, 450, 500), var = 9900) power.ancova(mu=c(400, 450, 500), var = 9900, power = 0.8, nr.covs = 10) power.ancova(n = rep(72/3, 3), mu=c(400, 450, 500), var = 9900, nr.covs = 10) ## Table 2 in Shieh (2020) power.ancova(mu=c(400, 450, 500), var = 7500, power = 0.8) power.ancova(n = rep(48/3, 3), mu=c(400, 450, 500), var = 7500) power.ancova(mu=c(400, 450, 500), var = 7500, power = 0.8, nr.covs = 10) power.ancova(n = rep(60/3, 3), mu=c(400, 450, 500), var = 7500, nr.covs = 10) ## Table 3 in Shieh (2020) power.ancova(mu=c(400, 450, 500), var = 1900, power = 0.8) power.ancova(n = rep(18/3, 3), mu=c(400, 450, 500), var = 1900) power.ancova(mu=c(400, 450, 500), var = 1900, power = 0.8, nr.covs = 10) power.ancova(n = rep(27/3, 3), mu=c(400, 450, 500), var = 1900, nr.covs = 10) ## ANOVA approach for Table 1-3 power.anova.test(groups = 3, between.var = var(c(400, 450, 500)), within.var = 10000, power = 0.8) power.anova.test(n = 63/3, groups = 3, between.var = var(c(400, 450, 500)), within.var = 10000) ## Table 4 in Shieh (2020) power.ancova(mu=c(410, 450, 490), var = 9900, power = 0.8) power.ancova(n = rep(96/3, 3), mu=c(410, 450, 490), var = 9900) power.ancova(mu=c(410, 450, 490), var = 9900, power = 0.8, nr.covs = 10) power.ancova(n = rep(105/3, 3), mu=c(410, 450, 490), var = 9900, nr.covs = 10) ## Table 5 in Shieh (2020) power.ancova(mu=c(410, 450, 490), var = 7500, power = 0.8) power.ancova(n = rep(72/3, 3), mu=c(410, 450, 490), var = 7500) power.ancova(mu=c(410, 450, 490), var = 7500, power = 0.8, nr.covs = 10) power.ancova(n = rep(84/3, 3), mu=c(410, 450, 490), var = 7500, nr.covs = 10) ## Table 6 in Shieh (2020) power.ancova(mu=c(410, 450, 490), var = 1900, power = 0.8) power.ancova(n = rep(24/3, 3), mu=c(410, 450, 490), var = 1900) power.ancova(mu=c(410, 450, 490), var = 1900, power = 0.8, nr.covs = 10) power.ancova(n = rep(33/3, 3), mu=c(410, 450, 490), var = 1900, nr.covs = 10) ## ANOVA approach for Table 4-6 power.anova.test(groups = 3, between.var = var(c(410, 450, 490)), within.var = 10000, power = 0.8) power.anova.test(n = 96/3, groups = 3, between.var = var(c(410, 450, 490)), within.var = 10000) ############################################################################### ## Example from Maxwell and Delaney (2004) according to Shieh (2020) ############################################################################### ## ANCOVA (balanced design) power.ancova(n = rep(30/3, 3), mu=c(7.5366, 11.9849, 13.9785), var = 29.0898) power.ancova(mu=c(7.5366, 11.9849, 13.9785), var = 29.0898, power = 0.8) power.ancova(mu=c(7.5366, 11.9849, 13.9785), var = 29.0898, power = 0.9) ## ANOVA power.anova.test(n = 30/3, groups = 3, between.var = var(c(7.5366, 11.9849, 13.9785)), within.var = 29.0898) power.anova.test(groups = 3, between.var = var(c(7.5366, 11.9849, 13.9785)), within.var = 29.0898, power = 0.8) power.anova.test(groups = 3, between.var = var(c(7.5366, 11.9849, 13.9785)), within.var = 29.0898, power = 0.9) ## ANCOVA - imbalanced design power.ancova(mu=c(7.5366, 11.9849, 13.9785), var = 29.0898, power = 0.8, group.ratio = c(1, 1.25, 1.5)) power.ancova(n = c(13, 16, 19), mu=c(7.5366, 11.9849, 13.9785), var = 29.0898, group.ratio = c(1, 1.25, 1.5)) power.ancova(mu=c(7.5366, 11.9849, 13.9785), var = 29.0898, power = 0.8, group.ratio = c(1, 0.8, 2/3)) power.ancova(n = c(17, 14, 12), mu=c(7.5366, 11.9849, 13.9785), var = 29.0898, group.ratio = c(1, 0.8, 2/3))
Compute sample size, power, delta, or significance level of a diagnostic test for an expected sensititivy or specificity.
power.diagnostic.test(sens = NULL, spec = NULL, n = NULL, delta = NULL, sig.level = 0.05, power = NULL, prev = NULL, method = c("exact", "asymptotic"), NMAX = 1e4) ssize.sens.ci(sens = NULL, n = NULL, delta = NULL, sig.level = 0.05, power = NULL, prev = NULL, method = c("exact", "asymptotic"), NMAX = 1e4) ssize.spec.ci(spec = NULL, n = NULL, delta = NULL, sig.level = 0.05, power = NULL, prev = NULL, method = c("exact", "asymptotic"), NMAX = 1e4)power.diagnostic.test(sens = NULL, spec = NULL, n = NULL, delta = NULL, sig.level = 0.05, power = NULL, prev = NULL, method = c("exact", "asymptotic"), NMAX = 1e4) ssize.sens.ci(sens = NULL, n = NULL, delta = NULL, sig.level = 0.05, power = NULL, prev = NULL, method = c("exact", "asymptotic"), NMAX = 1e4) ssize.spec.ci(spec = NULL, n = NULL, delta = NULL, sig.level = 0.05, power = NULL, prev = NULL, method = c("exact", "asymptotic"), NMAX = 1e4)
sens |
Expected sensitivity; either |
spec |
Expected specificity; either |
n |
Number of cases if |
delta |
|
sig.level |
Significance level (Type I error probability) |
power |
Power of test (1 minus Type II error probability) |
prev |
Expected prevalence, if |
method |
exact or asymptotic formula; default |
NMAX |
Maximum sample size considered in case |
Either sens or spec has to be specified which leads to
computations for either cases or controls.
Exactly one of the parameters n, delta, sig.level,
and power must be passed as NULL, and that parameter is determined
from the others. Notice that sig.level has a non-NULL default
so NULL must be explicitly passed if you want to compute it.
The computations are based on the formulas given in the Appendix of Flahault et al. (2005). Please be careful, in Equation (A1) the numerator should be squared, in equation (A2) and (A3) the second exponent should be n-i and not i.
As noted in Chu and Cole (2007) power is not a monotonically increasing
function in n but rather saw toothed (see also Chernick and Liu (2002)).
Hence, in our calculations we use the more conservative approach II);
i.e., the minimum sample size n such that the actual power is
larger or equal power andsuch that for any sample size larger
than n it also holds that the actual power is larger or equal
power.
Object of class "power.htest", a list of the arguments
(including the computed one) augmented with method and
note elements.
uniroot is used to solve the power equation if delta is unknown,
so you may see errors from it, notably about inability to bracket the
root when invalid arguments are given.
Matthias Kohl [email protected]
A. Flahault, M. Cadilhac, and G. Thomas (2005). Sample size calculation should be performed for design accuracy in diagnostic test studies. Journal of Clinical Epidemiology, 58(8):859-862.
H. Chu and S.R. Cole (2007). Sample size calculation using exact methods in diagnostic test studies. Journal of Clinical Epidemiology, 60(11):1201-1202.
M.R. Chernick amd C.Y. Liu (2002). The saw-toothed behavior of power versus sample size and software solutions: single binomial proportion using exact methods. Am Stat, 56:149-155.
## see n2 on page 1202 of Chu and Cole (2007) power.diagnostic.test(sens = 0.99, delta = 0.14, power = 0.95) # 40 power.diagnostic.test(sens = 0.99, delta = 0.13, power = 0.95) # 43 power.diagnostic.test(sens = 0.99, delta = 0.12, power = 0.95) # 47 power.diagnostic.test(sens = 0.98, delta = 0.13, power = 0.95) # 50 power.diagnostic.test(sens = 0.98, delta = 0.11, power = 0.95) # 58 ## see page 1201 of Chu and Cole (2007) power.diagnostic.test(sens = 0.95, delta = 0.1, n = 93) ## 0.957 power.diagnostic.test(sens = 0.95, delta = 0.1, n = 93, power = 0.95, sig.level = NULL) ## 0.0496 power.diagnostic.test(sens = 0.95, delta = 0.1, n = 102) ## 0.968 power.diagnostic.test(sens = 0.95, delta = 0.1, n = 102, power = 0.95, sig.level = NULL) ## 0.0471 ## yields 102 not 93! power.diagnostic.test(sens = 0.95, delta = 0.1, power = 0.95) ## function only for sensitivity ssize.sens.ci(sens = 0.99, delta = 0.14, power = 0.95) # 40 ## function only for specificity ssize.spec.ci(spec = 0.99, delta = 0.13, power = 0.95) # 43## see n2 on page 1202 of Chu and Cole (2007) power.diagnostic.test(sens = 0.99, delta = 0.14, power = 0.95) # 40 power.diagnostic.test(sens = 0.99, delta = 0.13, power = 0.95) # 43 power.diagnostic.test(sens = 0.99, delta = 0.12, power = 0.95) # 47 power.diagnostic.test(sens = 0.98, delta = 0.13, power = 0.95) # 50 power.diagnostic.test(sens = 0.98, delta = 0.11, power = 0.95) # 58 ## see page 1201 of Chu and Cole (2007) power.diagnostic.test(sens = 0.95, delta = 0.1, n = 93) ## 0.957 power.diagnostic.test(sens = 0.95, delta = 0.1, n = 93, power = 0.95, sig.level = NULL) ## 0.0496 power.diagnostic.test(sens = 0.95, delta = 0.1, n = 102) ## 0.968 power.diagnostic.test(sens = 0.95, delta = 0.1, n = 102, power = 0.95, sig.level = NULL) ## 0.0471 ## yields 102 not 93! power.diagnostic.test(sens = 0.95, delta = 0.1, power = 0.95) ## function only for sensitivity ssize.sens.ci(sens = 0.99, delta = 0.14, power = 0.95) # 40 ## function only for specificity ssize.spec.ci(spec = 0.99, delta = 0.13, power = 0.95) # 43
Compute the power of the two-sample Hsu t test, or determine parameters to obtain a target power; see Section 7.4.4 in Hedderich and Sachs (2016).
power.hsu.t.test(n = NULL, delta = NULL, sd1 = 1, sd2 = 1, sig.level = 0.05, power = NULL, alternative = c("two.sided", "one.sided"), strict = FALSE, tol = .Machine$double.eps^0.25) power.hsu.t.test2(n = NULL, delta = NULL, sd1 = 1, sd2 = 1, sig.level = 0.05, power = NULL, alternative = c("two.sided", "one.sided"), strict = FALSE, tol = .Machine$double.eps^0.25, abstol = .Machine$double.eps^0.5) ssize.hsu.t.test(delta = 1, sd1 = 1, sd2 = 1, sig.level = 0.05, power = 0.8, alternative = c("two.sided", "one.sided"), strict = FALSE)power.hsu.t.test(n = NULL, delta = NULL, sd1 = 1, sd2 = 1, sig.level = 0.05, power = NULL, alternative = c("two.sided", "one.sided"), strict = FALSE, tol = .Machine$double.eps^0.25) power.hsu.t.test2(n = NULL, delta = NULL, sd1 = 1, sd2 = 1, sig.level = 0.05, power = NULL, alternative = c("two.sided", "one.sided"), strict = FALSE, tol = .Machine$double.eps^0.25, abstol = .Machine$double.eps^0.5) ssize.hsu.t.test(delta = 1, sd1 = 1, sd2 = 1, sig.level = 0.05, power = 0.8, alternative = c("two.sided", "one.sided"), strict = FALSE)
n |
number of observations (per group) |
delta |
(expected) true difference in means |
sd1 |
(expected) standard deviation of group 1 |
sd2 |
(expected) standard deviation of group 2 |
sig.level |
significance level (Type I error probability) |
power |
power of test (1 minus Type II error probability) |
alternative |
one- or two-sided test. Can be abbreviated. |
strict |
use strict interpretation in two-sided case |
tol |
numerical tolerance used in root finding, the default providing (at least) four significant digits. |
abstol |
controls the absolute convergence tolerance of the
|
In case of power.hsu.t.test and power.hsu.t.test2,
exactly one of the parameters n, delta, power,
sd1, sd2 and sig.level must be passed as NULL,
and that parameter is determined from the others. Notice that the last three
have non-NULL defaults, so NULL must be explicitly passed if you want to
compute them.
In the case of the power.hsu.t.test function, both groups are assumed
to be of equal size (i.e. balanced design), whereas the power.hsu.t.test2
function allows for unequal sample sizes per group (i.e. unbalanced design).
For unequal variances, the balanced design does not give the minimum total
sample size (sum of the group sample sizes).
The Function ssize.hsu.t.test returns the optimal group sample sizes
by minimizing the total sample size (sum of the group sample sizes) that achieves
the given power. The solution of function power.hsu.t.test2 might not
represent the global optimum in terms of total sample size.
If strict = TRUE is used, the power will include the probability of
rejection in the opposite direction of the true effect, in the two-sided
case. Without this the power will be half the significance level if the
true difference is zero.
Object of class "power.htest", a list of the arguments
(including the computed one) augmented with method and
note elements.
The function and its documentation was adapted from power.t.test
implemented by Peter Dalgaard and based on previous work by Claus Ekstroem.
uniroot is used to solve the power equation for unknowns, so
you may see errors from it, notably about inability to bracket the
root when invalid arguments are given.
Matthias Kohl [email protected]
J. Hedderich, L. Sachs. Angewandte Statistik: Methodensammlung mit R. Springer 2016.
power.welch.t.test, power.t.test,
t.test, uniroot
## more conservative than classical or Welch t-test power.t.test(n = 20, delta = 1) power.welch.t.test(n = 20, delta = 1) power.hsu.t.test(n = 20, delta = 1) power.hsu.t.test(power = .90, delta = 1) power.hsu.t.test(power = .90, delta = 1, alternative = "one.sided") ## sd1 = 0.5, sd2 = 1 power.welch.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) power.hsu.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) power.welch.t.test2(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) # 22+87 = 109 power.hsu.t.test2(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) # 44+59 = 103 power.hsu.t.test2(n = c(44, 59), delta = 0.5, sd1 = 0.5, sd2 = 1) ssize.welch.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) # 35+62 = 97 ssize.hsu.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) # 39+62 = 101 power.hsu.t.test2(n = c(39, 62), delta = 0.5, sd1 = 0.5, sd2 = 1) power.welch.t.test2(n = c(39, 62), delta = 0.5, sd1 = 0.5, sd2 = 1)## more conservative than classical or Welch t-test power.t.test(n = 20, delta = 1) power.welch.t.test(n = 20, delta = 1) power.hsu.t.test(n = 20, delta = 1) power.hsu.t.test(power = .90, delta = 1) power.hsu.t.test(power = .90, delta = 1, alternative = "one.sided") ## sd1 = 0.5, sd2 = 1 power.welch.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) power.hsu.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) power.welch.t.test2(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) # 22+87 = 109 power.hsu.t.test2(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) # 44+59 = 103 power.hsu.t.test2(n = c(44, 59), delta = 0.5, sd1 = 0.5, sd2 = 1) ssize.welch.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) # 35+62 = 97 ssize.hsu.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) # 39+62 = 101 power.hsu.t.test2(n = c(39, 62), delta = 0.5, sd1 = 0.5, sd2 = 1) power.welch.t.test2(n = c(39, 62), delta = 0.5, sd1 = 0.5, sd2 = 1)
The function calculates either sample size or power for continuous multiple primary endpoints for at least one endpoint with known covariance.
power.mpe.atleast.one(K, n = NULL, delta = NULL, Sigma, SD, rho, sig.level = 0.05/K, power = NULL, n.max = 1e5, tol = .Machine$double.eps^0.25)power.mpe.atleast.one(K, n = NULL, delta = NULL, Sigma, SD, rho, sig.level = 0.05/K, power = NULL, n.max = 1e5, tol = .Machine$double.eps^0.25)
K |
number of endpoints |
n |
optional: sample size |
delta |
expected effect size |
Sigma |
A covariance of known matrix |
SD |
known standard deviations (length |
rho |
known correlations (length |
sig.level |
Significance level (Type I error probability) |
power |
optional: Power of test (1 minus Type II error probability) |
n.max |
upper end of the interval to be search for |
tol |
The desired accuracy |
The function can be used to either compute sample size or power for continuous multiple primary endpoints with known covariance where a significant difference for at least one endpoint is expected. The implementation is based on the formulas given in the references below.
The null hypothesis reads for
all where Tk is treatment k,
Ck is control k and K is the number of co-primary endpoints.
One has to specify either n or power, the other parameter is
determined. Moreover, either covariance matrix Sigma or standard
deviations SD and correlations rho must be given.
Object of class power.mpe.test, a list of arguments (including the
computed one) augmented with method and note elements.
The function first appeared in package mpe, which is now archived on CRAN.
Srinath Kolampally, Matthias Kohl [email protected]
Sugimoto, T. and Sozu, T. and Hamasaki, T. (2012). A convenient formula for sample size calculations in clinical trials with multiple co-primary continuous endpoints. Pharmaceut. Statist., 11: 118-128. doi:10.1002/pst.505
Sozu, T. and Sugimoto, T. and Hamasaki, T. and Evans, S.R. (2015). Sample Size Determination in Clinical Trials with Multiple Endpoints. Springer Briefs in Statistics, ISBN 978-3-319-22005-5.
## compute power power.mpe.atleast.one(K = 2, delta = c(0.2,0.2), Sigma = diag(c(1,1)), power = 0.8) ## compute sample size power.mpe.atleast.one(K = 2, delta = c(0.2,0.2), Sigma = diag(c(2,2)), power = 0.9) ## known covariance matrix Sigma <- matrix(c(1.440, 0.840, 1.296, 0.840, 0.840, 1.960, 0.168, 1.568, 1.296, 0.168, 1.440, 0.420, 0.840, 1.568, 0.420, 1.960), ncol = 4) ## compute power power.mpe.atleast.one(K = 4, n = 60, delta = c(0.5, 0.75, 0.5, 0.75), Sigma = Sigma) ## equivalent: known SDs and correlation rho power.mpe.atleast.one(K = 4, n = 60, delta = c(0.5, 0.75, 0.5, 0.75), SD = c(1.2, 1.4, 1.2, 1.4), rho = c(0.5, 0.9, 0.5, 0.1, 0.8, 0.25))## compute power power.mpe.atleast.one(K = 2, delta = c(0.2,0.2), Sigma = diag(c(1,1)), power = 0.8) ## compute sample size power.mpe.atleast.one(K = 2, delta = c(0.2,0.2), Sigma = diag(c(2,2)), power = 0.9) ## known covariance matrix Sigma <- matrix(c(1.440, 0.840, 1.296, 0.840, 0.840, 1.960, 0.168, 1.568, 1.296, 0.168, 1.440, 0.420, 0.840, 1.568, 0.420, 1.960), ncol = 4) ## compute power power.mpe.atleast.one(K = 4, n = 60, delta = c(0.5, 0.75, 0.5, 0.75), Sigma = Sigma) ## equivalent: known SDs and correlation rho power.mpe.atleast.one(K = 4, n = 60, delta = c(0.5, 0.75, 0.5, 0.75), SD = c(1.2, 1.4, 1.2, 1.4), rho = c(0.5, 0.9, 0.5, 0.1, 0.8, 0.25))
The function calculates either sample size or power for continuous multiple co-primary endpoints with known covariance.
power.mpe.known.var(K, n = NULL, delta = NULL, Sigma, SD, rho, sig.level = 0.05, power = NULL, n.max = 1e5, tol = .Machine$double.eps^0.25)power.mpe.known.var(K, n = NULL, delta = NULL, Sigma, SD, rho, sig.level = 0.05, power = NULL, n.max = 1e5, tol = .Machine$double.eps^0.25)
K |
number of co-primary endpoints |
n |
optional: sample size |
delta |
expected effect size (length |
Sigma |
known covariance matrix (dimension |
SD |
known standard deviations (length |
rho |
known correlations (length |
sig.level |
significance level (Type I error probability) |
power |
optional: power of test (1 minus Type II error probability) |
n.max |
upper end of the interval to be search for |
tol |
the desired accuracy for |
The function can be used to either compute sample size or power for continuous multiple co-primary endpoints with known covariance where a multivariate normal distribution is assumed. The implementation is based on the formulas given in the references below.
The null hypothesis reads for
at least one where Tk is treatment k,
Ck is control k and K is the number of co-primary endpoints.
One has to specify either n or power, the other parameter is
determined. Moreover, either covariance matrix Sigma or standard
deviations SD and correlations rho must be given.
Object of class power.mpe.test, a list of arguments (including the
computed one) augemented with method and note elements.
The function first appeared in package mpe, which is now archived on CRAN.
Srinath Kolampally, Matthias Kohl [email protected]
Sugimoto, T. and Sozu, T. and Hamasaki, T. (2012). A convenient formula for sample size calculations in clinical trials with multiple co-primary continuous endpoints. Pharmaceut. Statist., 11: 118-128. doi:10.1002/pst.505
Sozu, T. and Sugimoto, T. and Hamasaki, T. and Evans, S.R. (2015). Sample Size Determination in Clinical Trials with Multiple Endpoints. Springer Briefs in Statistics, ISBN 978-3-319-22005-5.
## compute power power.mpe.known.var(K = 2, n = 20, delta = c(1,1), Sigma = diag(c(1,1))) ## compute sample size power.mpe.known.var(K = 2, delta = c(1,1), Sigma = diag(c(2,2)), power = 0.9, sig.level = 0.025) ## known covariance matrix Sigma <- matrix(c(1.440, 0.840, 1.296, 0.840, 0.840, 1.960, 0.168, 1.568, 1.296, 0.168, 1.440, 0.420, 0.840, 1.568, 0.420, 1.960), ncol = 4) ## compute power power.mpe.known.var(K = 4, n = 60, delta = c(0.5, 0.75, 0.5, 0.75), Sigma = Sigma) ## equivalent: known SDs and correlation rho power.mpe.known.var(K = 4, n = 60,delta = c(0.5, 0.75, 0.5, 0.75), SD = c(1.2, 1.4, 1.2, 1.4), rho = c(0.5, 0.9, 0.5, 0.1, 0.8, 0.25))## compute power power.mpe.known.var(K = 2, n = 20, delta = c(1,1), Sigma = diag(c(1,1))) ## compute sample size power.mpe.known.var(K = 2, delta = c(1,1), Sigma = diag(c(2,2)), power = 0.9, sig.level = 0.025) ## known covariance matrix Sigma <- matrix(c(1.440, 0.840, 1.296, 0.840, 0.840, 1.960, 0.168, 1.568, 1.296, 0.168, 1.440, 0.420, 0.840, 1.568, 0.420, 1.960), ncol = 4) ## compute power power.mpe.known.var(K = 4, n = 60, delta = c(0.5, 0.75, 0.5, 0.75), Sigma = Sigma) ## equivalent: known SDs and correlation rho power.mpe.known.var(K = 4, n = 60,delta = c(0.5, 0.75, 0.5, 0.75), SD = c(1.2, 1.4, 1.2, 1.4), rho = c(0.5, 0.9, 0.5, 0.1, 0.8, 0.25))
The function calculates either sample size or power for continuous multiple co-primary endpoints with unknown covariance.
power.mpe.unknown.var(K, n = NULL, delta = NULL, Sigma, SD, rho, sig.level = 0.05, power = NULL, M = 10000, n.min = NULL, n.max = NULL, tol = .Machine$double.eps^0.25, use.uniroot = TRUE)power.mpe.unknown.var(K, n = NULL, delta = NULL, Sigma, SD, rho, sig.level = 0.05, power = NULL, M = 10000, n.min = NULL, n.max = NULL, tol = .Machine$double.eps^0.25, use.uniroot = TRUE)
K |
number of co-primary endpoints |
n |
optional: sample size |
delta |
expected effect size (length |
Sigma |
unknown covariance matrix (dimension |
SD |
unknown standard deviations (length |
rho |
unknown correlations (length |
sig.level |
significance level (Type I error probability) |
power |
optional: power of test (1 minus Type II error probability) |
M |
Number of replications for the required simulations. |
n.min |
Starting point of search interval for sample size |
n.max |
End point of search interval for sample size, must be larger than |
tol |
the desired accuracy for |
use.uniroot |
Finds one root of one equation |
The function can be used to either compute sample size or power for continuous multiple co-primary endpoints with unknown covariance. The implementation is based on the formulas given in the references below.
The null hypothesis reads for
at least one where Tk is treatment k,
Ck is control k and K is the number of co-primary endpoints.
One has to specify either n or power, the other parameter is
determined. An approach to calculate sample size n, is to first call
power.mpe.known.var and use the result as n.min. The input for
n.max must be larger then n.min. Moreover, either covariance
matrix Sigma or standard deviations SD and correlations rho
must be given.
The sample size is calculated by simulating Wishart distributed random matrices, hence the results include a certain random variation.
Object of class power.mpe.test, a list of arguments (including the
computed one) augmented with method and note elements.
The function first appeared in package mpe, which is now archived on CRAN.
Srinath Kolampally, Matthias Kohl [email protected]
Sugimoto, T. and Sozu, T. and Hamasaki, T. (2012). A convenient formula for sample size calculations in clinical trials with multiple co-primary continuous endpoints. Pharmaceut. Statist., 11: 118-128. doi:10.1002/pst.505
Sozu, T. and Sugimoto, T. and Hamasaki, T. and Evans, S.R. (2015). Sample Size Determination in Clinical Trials with Multiple Endpoints. Springer Briefs in Statistics, ISBN 978-3-319-22005-5.
## compute power ## Not run: power.mpe.unknown.var(K = 2, n = 20, delta = c(1,1), Sigma = diag(c(1,1))) ## To compute sample size, first assume covariance as known power.mpe.known.var(K = 2, delta = c(1,1), Sigma = diag(c(2,2)), power = 0.9, sig.level = 0.025) ## The value of n, which is 51, is used as n.min and n.max must be larger ## then n.min so we try 60. power.mpe.unknown.var(K = 2, delta = c(1,1), Sigma = diag(c(2,2)), power = 0.9, sig.level = 0.025, n.min = 51, n.max = 60) ## More complex example with unknown covariance matrix assumed to be Sigma <- matrix(c(1.440, 0.840, 1.296, 0.840, 0.840, 1.960, 0.168, 1.568, 1.296, 0.168, 1.440, 0.420, 0.840, 1.568, 0.420, 1.960), ncol = 4) ## compute power power.mpe.unknown.var(K = 4, n = 90, delta = c(0.5, 0.75, 0.5, 0.75), Sigma = Sigma) ## equivalent: unknown SDs and correlation rho power.mpe.unknown.var(K = 4, n = 90, delta = c(0.5, 0.75, 0.5, 0.75), SD = c(1.2, 1.4, 1.2, 1.4), rho = c(0.5, 0.9, 0.5, 0.1, 0.8, 0.25)) ## End(Not run)## compute power ## Not run: power.mpe.unknown.var(K = 2, n = 20, delta = c(1,1), Sigma = diag(c(1,1))) ## To compute sample size, first assume covariance as known power.mpe.known.var(K = 2, delta = c(1,1), Sigma = diag(c(2,2)), power = 0.9, sig.level = 0.025) ## The value of n, which is 51, is used as n.min and n.max must be larger ## then n.min so we try 60. power.mpe.unknown.var(K = 2, delta = c(1,1), Sigma = diag(c(2,2)), power = 0.9, sig.level = 0.025, n.min = 51, n.max = 60) ## More complex example with unknown covariance matrix assumed to be Sigma <- matrix(c(1.440, 0.840, 1.296, 0.840, 0.840, 1.960, 0.168, 1.568, 1.296, 0.168, 1.440, 0.420, 0.840, 1.568, 0.420, 1.960), ncol = 4) ## compute power power.mpe.unknown.var(K = 4, n = 90, delta = c(0.5, 0.75, 0.5, 0.75), Sigma = Sigma) ## equivalent: unknown SDs and correlation rho power.mpe.unknown.var(K = 4, n = 90, delta = c(0.5, 0.75, 0.5, 0.75), SD = c(1.2, 1.4, 1.2, 1.4), rho = c(0.5, 0.9, 0.5, 0.1, 0.8, 0.25)) ## End(Not run)
Compute sample size or power for comparing two negative binomial rates.
power.nb.test(n = NULL, mu0, mu1, RR, duration = 1, theta, ssize.ratio = 1, sig.level = 0.05, power = NULL, alternative = c("two.sided", "one.sided"), approach = 3)power.nb.test(n = NULL, mu0, mu1, RR, duration = 1, theta, ssize.ratio = 1, sig.level = 0.05, power = NULL, alternative = c("two.sided", "one.sided"), approach = 3)
n |
Sample size for group 0 (control group). |
mu0 |
expected rate of events per time unit for group 0 |
mu1 |
expected rate of events per time unit for group 1 |
RR |
ratio of expected event rates: mu1/mu0 |
duration |
(average) treatment duration |
theta |
theta parameter of negative binomial distribution; see |
ssize.ratio |
ratio of sample sizes: n1/n where n1 is sample size of group 1 |
sig.level |
Significance level (Type I error probability) |
power |
Power of test (1 minus Type II error probability) |
alternative |
one- or two-sided test |
approach |
1, 2, or 3; see Zhu and Lakkis (2014). |
Exactly one of the parameters n and power must be passed as
NULL, and that parameter is determined from the other.
The computations are based on the formulas given in Zhu and Lakkis (2014).
Please be careful, as we are using a slightly different parametrization
(theta = 1/k).
Zhu and Lakkis (2014) based on their simulation studies recommend to use their approach 2 or 3.
Object of class "power.htest", a list of the arguments
(including the computed one) augmented with a note element.
Matthias Kohl [email protected]
H. Zhu and H. Lakkis (2014). Sample size calculation for comparing two negative binomial rates. Statistics in Medicine, 33:376-387.
## examples from Table I in Zhu and Lakkis (2014) ## theta = 1/k, RR = rr, mu0 = r0, duration = mu_t power.nb.test(mu0 = 0.8, RR = 0.85, theta = 1/0.4, duration = 0.75, power = 0.8, approach = 1) power.nb.test(mu0 = 0.8, RR = 0.85, theta = 1/0.4, duration = 0.75, power = 0.8, approach = 2) power.nb.test(mu0 = 0.8, RR = 0.85, theta = 1/0.4, duration = 0.75, power = 0.8, approach = 3) power.nb.test(mu0 = 1.4, RR = 1.15, theta = 1/1.5, duration = 0.75, power = 0.8, approach = 1) power.nb.test(mu0 = 1.4, RR = 1.15, theta = 1/1.5, duration = 0.75, power = 0.8, approach = 2) power.nb.test(mu0 = 1.4, RR = 1.15, theta = 1/1.5, duration = 0.75, power = 0.8, approach = 3) ## examples from Table II in Zhu and Lakkis (2014) - seem to be total sample sizes ## can reproduce the results with mu_t = 1.0 (not 0.7!) power.nb.test(mu0 = 2.0, RR = 0.5, theta = 1, duration = 1.0, ssize.ratio = 1, power = 0.8, approach = 1) power.nb.test(mu0 = 2.0, RR = 0.5, theta = 1, duration = 1.0, ssize.ratio = 1, power = 0.8, approach = 2) power.nb.test(mu0 = 2.0, RR = 0.5, theta = 1, duration = 1.0, ssize.ratio = 1, power = 0.8, approach = 3) power.nb.test(mu0 = 10.0, RR = 1.5, theta = 1/5, duration = 1.0, ssize.ratio = 3/2, power = 0.8, approach = 1) power.nb.test(mu0 = 10.0, RR = 1.5, theta = 1/5, duration = 1.0, ssize.ratio = 3/2, power = 0.8, approach = 2) power.nb.test(mu0 = 10.0, RR = 1.5, theta = 1/5, duration = 1.0, ssize.ratio = 3/2, power = 0.8, approach = 3) ## examples from Table III in Zhu and Lakkis (2014) power.nb.test(mu0 = 5.0, RR = 2.0, theta = 1/0.5, duration = 1, power = 0.8, approach = 1) power.nb.test(mu0 = 5.0, RR = 2.0, theta = 1/0.5, duration = 1, power = 0.8, approach = 2) power.nb.test(mu0 = 5.0, RR = 2.0, theta = 1/0.5, duration = 1, power = 0.8, approach = 3) ## examples from Table IV in Zhu and Lakkis (2014) power.nb.test(mu0 = 5.9/3, RR = 0.4, theta = 0.49, duration = 3, power = 0.9, approach = 1) power.nb.test(mu0 = 5.9/3, RR = 0.4, theta = 0.49, duration = 3, power = 0.9, approach = 2) power.nb.test(mu0 = 5.9/3, RR = 0.4, theta = 0.49, duration = 3, power = 0.9, approach = 3) power.nb.test(mu0 = 13/6, RR = 0.2, theta = 0.52, duration = 6, power = 0.9, approach = 1) power.nb.test(mu0 = 13/6, RR = 0.2, theta = 0.52, duration = 6, power = 0.9, approach = 2) power.nb.test(mu0 = 13/6, RR = 0.2, theta = 0.52, duration = 6, power = 0.9, approach = 3) ## see Section 5 of Zhu and Lakkis (2014) power.nb.test(mu0 = 0.66, RR = 0.8, theta = 1/0.8, duration = 0.9, power = 0.9)## examples from Table I in Zhu and Lakkis (2014) ## theta = 1/k, RR = rr, mu0 = r0, duration = mu_t power.nb.test(mu0 = 0.8, RR = 0.85, theta = 1/0.4, duration = 0.75, power = 0.8, approach = 1) power.nb.test(mu0 = 0.8, RR = 0.85, theta = 1/0.4, duration = 0.75, power = 0.8, approach = 2) power.nb.test(mu0 = 0.8, RR = 0.85, theta = 1/0.4, duration = 0.75, power = 0.8, approach = 3) power.nb.test(mu0 = 1.4, RR = 1.15, theta = 1/1.5, duration = 0.75, power = 0.8, approach = 1) power.nb.test(mu0 = 1.4, RR = 1.15, theta = 1/1.5, duration = 0.75, power = 0.8, approach = 2) power.nb.test(mu0 = 1.4, RR = 1.15, theta = 1/1.5, duration = 0.75, power = 0.8, approach = 3) ## examples from Table II in Zhu and Lakkis (2014) - seem to be total sample sizes ## can reproduce the results with mu_t = 1.0 (not 0.7!) power.nb.test(mu0 = 2.0, RR = 0.5, theta = 1, duration = 1.0, ssize.ratio = 1, power = 0.8, approach = 1) power.nb.test(mu0 = 2.0, RR = 0.5, theta = 1, duration = 1.0, ssize.ratio = 1, power = 0.8, approach = 2) power.nb.test(mu0 = 2.0, RR = 0.5, theta = 1, duration = 1.0, ssize.ratio = 1, power = 0.8, approach = 3) power.nb.test(mu0 = 10.0, RR = 1.5, theta = 1/5, duration = 1.0, ssize.ratio = 3/2, power = 0.8, approach = 1) power.nb.test(mu0 = 10.0, RR = 1.5, theta = 1/5, duration = 1.0, ssize.ratio = 3/2, power = 0.8, approach = 2) power.nb.test(mu0 = 10.0, RR = 1.5, theta = 1/5, duration = 1.0, ssize.ratio = 3/2, power = 0.8, approach = 3) ## examples from Table III in Zhu and Lakkis (2014) power.nb.test(mu0 = 5.0, RR = 2.0, theta = 1/0.5, duration = 1, power = 0.8, approach = 1) power.nb.test(mu0 = 5.0, RR = 2.0, theta = 1/0.5, duration = 1, power = 0.8, approach = 2) power.nb.test(mu0 = 5.0, RR = 2.0, theta = 1/0.5, duration = 1, power = 0.8, approach = 3) ## examples from Table IV in Zhu and Lakkis (2014) power.nb.test(mu0 = 5.9/3, RR = 0.4, theta = 0.49, duration = 3, power = 0.9, approach = 1) power.nb.test(mu0 = 5.9/3, RR = 0.4, theta = 0.49, duration = 3, power = 0.9, approach = 2) power.nb.test(mu0 = 5.9/3, RR = 0.4, theta = 0.49, duration = 3, power = 0.9, approach = 3) power.nb.test(mu0 = 13/6, RR = 0.2, theta = 0.52, duration = 6, power = 0.9, approach = 1) power.nb.test(mu0 = 13/6, RR = 0.2, theta = 0.52, duration = 6, power = 0.9, approach = 2) power.nb.test(mu0 = 13/6, RR = 0.2, theta = 0.52, duration = 6, power = 0.9, approach = 3) ## see Section 5 of Zhu and Lakkis (2014) power.nb.test(mu0 = 0.66, RR = 0.8, theta = 1/0.8, duration = 0.9, power = 0.9)
Compute the power of the one-sample test for proportions, or determine parameters to obtain a target power.
power.prop1.test(n = NULL, p1 = NULL, p0 = 0.5, sig.level = 0.05, power = NULL, alternative = c("two.sided", "less", "greater"), cont.corr = TRUE, tol = .Machine$double.eps^0.25)power.prop1.test(n = NULL, p1 = NULL, p0 = 0.5, sig.level = 0.05, power = NULL, alternative = c("two.sided", "less", "greater"), cont.corr = TRUE, tol = .Machine$double.eps^0.25)
n |
number of observations (per group) |
p1 |
expected probability |
p0 |
probability under the null hypothesis |
sig.level |
significance level (Type I error probability) |
power |
power of test (1 minus Type II error probability) |
alternative |
one- or two-sided test. Can be abbreviated. |
cont.corr |
use continuity correction |
tol |
numerical tolerance used in root finding, the default providing (at least) four significant digits. |
Exactly one of the parameters n, p1, power, and
sig.level must be passed as NULL, and that parameter is determined
from the others. Notice that sig.level has a non-NULL default
so NULL must be explicitly passed if you want it computed.
The computation is based on the asymptotic formulas provided in Section 2.5.1
of Fleiss et al. (2003). If cont.corr = TRUE a continuity correction
is applied, which may lead to better approximations of the finite-sample
values.
Object of class "power.htest", a list of the arguments
(including the computed one) augmented with method and
note elements.
The documentation was adapted from power.prop.test.
Matthias Kohl [email protected]
J.L. Fleiss, B. Levin and M.C. Paik (2003). Statistical Methods for Rates and Proportions. Wiley Series in Probability and Statistics.
power.prop1.test(p1 = 0.4, power = 0.8) power.prop1.test(p1 = 0.4, power = 0.8, cont.corr = FALSE) power.prop1.test(p1 = 0.6, power = 0.8) power.prop1.test(n = 204, power = 0.8) power.prop1.test(n = 204, p1 = 0.4, power = 0.8, sig.level = NULL) power.prop1.test(n = 194, p1 = 0.4, power = 0.8, sig.level = NULL, cont.corr = FALSE) power.prop1.test(p1 = 0.1, p0 = 0.3, power = 0.8, alternative = "less") power.prop1.test(p1 = 0.1, p0 = 0.3, power = 0.8, alternative = "less", cont.corr = FALSE) power.prop1.test(n = 31, p0 = 0.3, power = 0.8, alternative = "less") power.prop1.test(n = 31, p1 = 0.1, p0 = 0.3, power = 0.8, sig.level = NULL, alternative = "less") power.prop1.test(p1 = 0.5, p0 = 0.3, power = 0.8, alternative = "greater") power.prop1.test(p1 = 0.5, p0 = 0.3, power = 0.8, alternative = "greater", cont.corr = FALSE) power.prop1.test(n = 40, p0 = 0.3, power = 0.8, alternative = "greater") power.prop1.test(n = 40, p1 = 0.5, p0 = 0.3, power = 0.8, sig.level = NULL, alternative = "greater")power.prop1.test(p1 = 0.4, power = 0.8) power.prop1.test(p1 = 0.4, power = 0.8, cont.corr = FALSE) power.prop1.test(p1 = 0.6, power = 0.8) power.prop1.test(n = 204, power = 0.8) power.prop1.test(n = 204, p1 = 0.4, power = 0.8, sig.level = NULL) power.prop1.test(n = 194, p1 = 0.4, power = 0.8, sig.level = NULL, cont.corr = FALSE) power.prop1.test(p1 = 0.1, p0 = 0.3, power = 0.8, alternative = "less") power.prop1.test(p1 = 0.1, p0 = 0.3, power = 0.8, alternative = "less", cont.corr = FALSE) power.prop1.test(n = 31, p0 = 0.3, power = 0.8, alternative = "less") power.prop1.test(n = 31, p1 = 0.1, p0 = 0.3, power = 0.8, sig.level = NULL, alternative = "less") power.prop1.test(p1 = 0.5, p0 = 0.3, power = 0.8, alternative = "greater") power.prop1.test(p1 = 0.5, p0 = 0.3, power = 0.8, alternative = "greater", cont.corr = FALSE) power.prop1.test(n = 40, p0 = 0.3, power = 0.8, alternative = "greater") power.prop1.test(n = 40, p1 = 0.5, p0 = 0.3, power = 0.8, sig.level = NULL, alternative = "greater")
Compute the power of the two-sample Welch t test, or determine parameters to obtain a target power.
power.welch.t.test(n = NULL, delta = NULL, sd1 = 1, sd2 = 1, sig.level = 0.05, power = NULL, alternative = c("two.sided", "one.sided"), strict = FALSE, tol = .Machine$double.eps^0.25) power.welch.t.test2(n = NULL, delta = NULL, sd1 = 1, sd2 = 1, sig.level = 0.05, power = NULL, alternative = c("two.sided", "one.sided"), strict = FALSE, tol = .Machine$double.eps^0.25, abstol = .Machine$double.eps^0.5) ssize.welch.t.test(delta = 1, sd1 = 1, sd2 = 1, sig.level = 0.05, power = 0.8, alternative = c("two.sided", "one.sided"), strict = FALSE)power.welch.t.test(n = NULL, delta = NULL, sd1 = 1, sd2 = 1, sig.level = 0.05, power = NULL, alternative = c("two.sided", "one.sided"), strict = FALSE, tol = .Machine$double.eps^0.25) power.welch.t.test2(n = NULL, delta = NULL, sd1 = 1, sd2 = 1, sig.level = 0.05, power = NULL, alternative = c("two.sided", "one.sided"), strict = FALSE, tol = .Machine$double.eps^0.25, abstol = .Machine$double.eps^0.5) ssize.welch.t.test(delta = 1, sd1 = 1, sd2 = 1, sig.level = 0.05, power = 0.8, alternative = c("two.sided", "one.sided"), strict = FALSE)
n |
number of observations (per group) |
delta |
(expected) true difference in means |
sd1 |
(expected) standard deviation of group 1 |
sd2 |
(expected) standard deviation of group 2 |
sig.level |
significance level (Type I error probability) |
power |
power of test (1 minus Type II error probability) |
alternative |
one- or two-sided test. Can be abbreviated. |
strict |
use strict interpretation in two-sided case |
tol |
numerical tolerance used in root finding, the default providing (at least) four significant digits. |
abstol |
controls the absolute convergence tolerance of the
|
In case of power.welch.t.test and power.welch.t.test2, exactly
one of the parameters n, delta, power,
sd1, sd2 and sig.level must be passed as NULL,
and that parameter is determined from the others. Notice that the last three
have non-NULL defaults, so NULL must be explicitly passed if you want to
compute them.
In the case of the power.welch.t.test function, both groups are assumed
to be of equal size (i.e. balanced design), whereas the power.welch.t.test2
function allows for unequal sample sizes per group (i.e. unbalanced design).
For unequal variances, the balanced design does not give the minimum total
sample size (sum of the group sample sizes).
The Function ssize.welch.t.test returns the optimal group sample sizes
by minimizing the total sample size (sum of the group sample sizes) that achieves
the given power. The solution of function power.welch.t.test2 might not
represent the global optimum in terms of total sample size.
If strict = TRUE is used, the power will include the probability of
rejection in the opposite direction of the true effect, in the two-sided
case. Without this the power will be half the significance level if the
true difference is zero.
Object of class "power.htest", a list of the arguments
(including the computed one) augmented with method and
note elements.
The function and its documentation was adapted from power.t.test
implemented by Peter Dalgaard and based on previous work by Claus Ekstroem.
uniroot is used to solve the power equation for unknowns, so
you may see errors from it, notably about inability to bracket the
root when invalid arguments are given.
Matthias Kohl [email protected]
S.L. Jan and G. Shieh (2011). Optimal sample sizes for Welch's test under various allocation and cost considerations. Behav Res Methods, 43, 4:1014-22.
## identical results as power.t.test, since sd = sd1 = sd2 = 1 power.welch.t.test(n = 20, delta = 1) power.welch.t.test(power = .90, delta = 1) power.welch.t.test(power = .90, delta = 1, alternative = "one.sided") power.welch.t.test2(n = c(20, 20), delta = 1) power.welch.t.test2(power = .90, delta = 1) power.welch.t.test2(power = .90, delta = 1, alternative = "one.sided") ssize.welch.t.test(power = 0.90, delta = 1) power.welch.t.test2(n = c(25, 20), delta = 1) power.welch.t.test2(n = c(20, 25), delta = 1) power.welch.t.test(n = 23, delta = 1) power.welch.t.test2(n = c(23, 23), delta = 1) ssize.welch.t.test(power = .90, delta = 1, alternative = "one.sided") power.welch.t.test2(n = c(19, 17), delta = 1, alternative = "one.sided") power.welch.t.test2(n = c(17, 19), delta = 1, alternative = "one.sided") power.welch.t.test(n = 18, delta = 1, alternative = "one.sided") power.welch.t.test2(n = c(18, 18), delta = 1, alternative = "one.sided") ########################################################### ## see Table 1 of Jan and Shieh (2011) ## delta = 1, sig.level = 0.05, power = 0.9 ########################################################### ####################################### ## sd1 = 1/3, sd2 = 1 ####################################### power.welch.t.test(delta = 1, sd1 = 1/3, sd2 = 1, power = 0.9) power.welch.t.test2(delta = 1, sd1 = 1/3, sd2 = 1, power = 0.9) ssize.welch.t.test(delta = 1, sd1 = 1/3, sd2 = 1, power = 0.9) ## n1 = 14, n2 = 14 power.welch.t.test(n = 14, delta = 1, sd1 = 1/3, sd2 = 1, power = NULL) power.welch.t.test2(n = c(14, 14), delta = 1, sd1 = 1/3, sd2 = 1, power = NULL) ## n1 = 8, n2 = 14 power.welch.t.test2(n = c(8, 16), delta = 1, sd1 = 1/3, sd2 = 1, power = NULL) ## n1 = 6, n2 = 18 power.welch.t.test2(n = c(6, 18), delta = 1, sd1 = 1/3, sd2 = 1, power = NULL) ####################################### ## sd1 = 1/2, sd2 = 1 ####################################### power.welch.t.test(delta = 1, sd1 = 1/2, sd2 = 1, power = 0.9) power.welch.t.test2(delta = 1, sd1 = 1/2, sd2 = 1, power = 0.9) ssize.welch.t.test(delta = 1, sd1 = 1/2, sd2 = 1, power = 0.9) ## n1 = 15, n2 = 15 power.welch.t.test(n = 15, delta = 1, sd1 = 1/2, sd2 = 1, power = NULL) power.welch.t.test2(n = c(15, 15), delta = 1, sd1 = 1/2, sd2 = 1, power = NULL) ## n1 = 9, n2 = 18 power.welch.t.test2(n = c(9, 18), delta = 1, sd1 = 1/2, sd2 = 1, power = NULL) ## n1 = 7, n2 = 21 power.welch.t.test2(n = c(7, 21), delta = 1, sd1 = 1/2, sd2 = 1, power = NULL) ####################################### ## sd1 = 1, sd2 = 1 ####################################### power.t.test(delta = 1, sd = 1, power = 0.9) power.welch.t.test(delta = 1, sd1 = 1, sd2 = 1, power = 0.9) power.welch.t.test2(delta = 1, sd1 = 1, sd2 = 1, power = 0.9) ssize.welch.t.test(delta = 1, sd1 = 1, sd2 = 1, power = 0.9) ## n1 = 23, n2 = 23 power.t.test(n = 23, delta = 1, sd = 1, power = NULL) power.welch.t.test(n = 23, delta = 1, sd1 = 1, sd2 = 1, power = NULL) power.welch.t.test2(n = c(23, 23), delta = 1, sd1 = 1, sd2 = 1, power = NULL) ## n1 = 17, n2 = 34 power.welch.t.test2(n = c(17, 34), delta = 1, sd1 = 1, sd2 = 1, power = NULL) ## n1 = 16, n2 = 48 power.welch.t.test2(n = c(16, 48), delta = 1, sd1 = 1, sd2 = 1, power = NULL) ####################################### ## sd1 = 2, sd2 = 1 ####################################### power.welch.t.test(delta = 1, sd1 = 2, sd2 = 1, power = 0.9) power.welch.t.test2(delta = 1, sd1 = 2, sd2 = 1, power = 0.9) ssize.welch.t.test(delta = 1, sd1 = 2, sd2 = 1, power = 0.9) ## n1 = 54, n2 = 54 power.welch.t.test(n = 54, delta = 1, sd1 = 2, sd2 = 1, power = NULL) power.welch.t.test2(n = c(54, 54), delta = 1, sd1 = 2, sd2 = 1, power = NULL) ## n1 = 49, n2 = 98 power.welch.t.test2(n = c(49, 98), delta = 1, sd1 = 2, sd2 = 1, power = NULL) ## n1 = 48, n2 = 144 power.welch.t.test2(n = c(48, 144), delta = 1, sd1 = 2, sd2 = 1, power = NULL) ####################################### ## sd1 = 3, sd2 = 1 ####################################### power.welch.t.test(delta = 1, sd1 = 3, sd2 = 1, power = 0.9) power.welch.t.test2(delta = 1, sd1 = 3, sd2 = 1, power = 0.9) ssize.welch.t.test(delta = 1, sd1 = 3, sd2 = 1, power = 0.9) ## n1 = 107, n2 = 107 power.welch.t.test(n = 107, delta = 1, sd1 = 3, sd2 = 1, power = NULL) power.welch.t.test2(n = c(107, 107), delta = 1, sd1 = 3, sd2 = 1, power = NULL) ## n1 = 102, n2 = 204 power.welch.t.test2(n = c(102, 204), delta = 1, sd1 = 3, sd2 = 1, power = NULL) ## n1 = 100, n2 = 300 power.welch.t.test2(n = c(100, 300), delta = 1, sd1 = 3, sd2 = 1, power = NULL)## identical results as power.t.test, since sd = sd1 = sd2 = 1 power.welch.t.test(n = 20, delta = 1) power.welch.t.test(power = .90, delta = 1) power.welch.t.test(power = .90, delta = 1, alternative = "one.sided") power.welch.t.test2(n = c(20, 20), delta = 1) power.welch.t.test2(power = .90, delta = 1) power.welch.t.test2(power = .90, delta = 1, alternative = "one.sided") ssize.welch.t.test(power = 0.90, delta = 1) power.welch.t.test2(n = c(25, 20), delta = 1) power.welch.t.test2(n = c(20, 25), delta = 1) power.welch.t.test(n = 23, delta = 1) power.welch.t.test2(n = c(23, 23), delta = 1) ssize.welch.t.test(power = .90, delta = 1, alternative = "one.sided") power.welch.t.test2(n = c(19, 17), delta = 1, alternative = "one.sided") power.welch.t.test2(n = c(17, 19), delta = 1, alternative = "one.sided") power.welch.t.test(n = 18, delta = 1, alternative = "one.sided") power.welch.t.test2(n = c(18, 18), delta = 1, alternative = "one.sided") ########################################################### ## see Table 1 of Jan and Shieh (2011) ## delta = 1, sig.level = 0.05, power = 0.9 ########################################################### ####################################### ## sd1 = 1/3, sd2 = 1 ####################################### power.welch.t.test(delta = 1, sd1 = 1/3, sd2 = 1, power = 0.9) power.welch.t.test2(delta = 1, sd1 = 1/3, sd2 = 1, power = 0.9) ssize.welch.t.test(delta = 1, sd1 = 1/3, sd2 = 1, power = 0.9) ## n1 = 14, n2 = 14 power.welch.t.test(n = 14, delta = 1, sd1 = 1/3, sd2 = 1, power = NULL) power.welch.t.test2(n = c(14, 14), delta = 1, sd1 = 1/3, sd2 = 1, power = NULL) ## n1 = 8, n2 = 14 power.welch.t.test2(n = c(8, 16), delta = 1, sd1 = 1/3, sd2 = 1, power = NULL) ## n1 = 6, n2 = 18 power.welch.t.test2(n = c(6, 18), delta = 1, sd1 = 1/3, sd2 = 1, power = NULL) ####################################### ## sd1 = 1/2, sd2 = 1 ####################################### power.welch.t.test(delta = 1, sd1 = 1/2, sd2 = 1, power = 0.9) power.welch.t.test2(delta = 1, sd1 = 1/2, sd2 = 1, power = 0.9) ssize.welch.t.test(delta = 1, sd1 = 1/2, sd2 = 1, power = 0.9) ## n1 = 15, n2 = 15 power.welch.t.test(n = 15, delta = 1, sd1 = 1/2, sd2 = 1, power = NULL) power.welch.t.test2(n = c(15, 15), delta = 1, sd1 = 1/2, sd2 = 1, power = NULL) ## n1 = 9, n2 = 18 power.welch.t.test2(n = c(9, 18), delta = 1, sd1 = 1/2, sd2 = 1, power = NULL) ## n1 = 7, n2 = 21 power.welch.t.test2(n = c(7, 21), delta = 1, sd1 = 1/2, sd2 = 1, power = NULL) ####################################### ## sd1 = 1, sd2 = 1 ####################################### power.t.test(delta = 1, sd = 1, power = 0.9) power.welch.t.test(delta = 1, sd1 = 1, sd2 = 1, power = 0.9) power.welch.t.test2(delta = 1, sd1 = 1, sd2 = 1, power = 0.9) ssize.welch.t.test(delta = 1, sd1 = 1, sd2 = 1, power = 0.9) ## n1 = 23, n2 = 23 power.t.test(n = 23, delta = 1, sd = 1, power = NULL) power.welch.t.test(n = 23, delta = 1, sd1 = 1, sd2 = 1, power = NULL) power.welch.t.test2(n = c(23, 23), delta = 1, sd1 = 1, sd2 = 1, power = NULL) ## n1 = 17, n2 = 34 power.welch.t.test2(n = c(17, 34), delta = 1, sd1 = 1, sd2 = 1, power = NULL) ## n1 = 16, n2 = 48 power.welch.t.test2(n = c(16, 48), delta = 1, sd1 = 1, sd2 = 1, power = NULL) ####################################### ## sd1 = 2, sd2 = 1 ####################################### power.welch.t.test(delta = 1, sd1 = 2, sd2 = 1, power = 0.9) power.welch.t.test2(delta = 1, sd1 = 2, sd2 = 1, power = 0.9) ssize.welch.t.test(delta = 1, sd1 = 2, sd2 = 1, power = 0.9) ## n1 = 54, n2 = 54 power.welch.t.test(n = 54, delta = 1, sd1 = 2, sd2 = 1, power = NULL) power.welch.t.test2(n = c(54, 54), delta = 1, sd1 = 2, sd2 = 1, power = NULL) ## n1 = 49, n2 = 98 power.welch.t.test2(n = c(49, 98), delta = 1, sd1 = 2, sd2 = 1, power = NULL) ## n1 = 48, n2 = 144 power.welch.t.test2(n = c(48, 144), delta = 1, sd1 = 2, sd2 = 1, power = NULL) ####################################### ## sd1 = 3, sd2 = 1 ####################################### power.welch.t.test(delta = 1, sd1 = 3, sd2 = 1, power = 0.9) power.welch.t.test2(delta = 1, sd1 = 3, sd2 = 1, power = 0.9) ssize.welch.t.test(delta = 1, sd1 = 3, sd2 = 1, power = 0.9) ## n1 = 107, n2 = 107 power.welch.t.test(n = 107, delta = 1, sd1 = 3, sd2 = 1, power = NULL) power.welch.t.test2(n = c(107, 107), delta = 1, sd1 = 3, sd2 = 1, power = NULL) ## n1 = 102, n2 = 204 power.welch.t.test2(n = c(102, 204), delta = 1, sd1 = 3, sd2 = 1, power = NULL) ## n1 = 100, n2 = 300 power.welch.t.test2(n = c(100, 300), delta = 1, sd1 = 3, sd2 = 1, power = NULL)
Printing objects of class "power.mpe.test" by simple print methods.
## S3 method for class 'power.mpe.test' print(x, digits = getOption("digits"), ...)## S3 method for class 'power.mpe.test' print(x, digits = getOption("digits"), ...)
x |
object of class |
digits |
number of significant digits to be used. |
... |
further arguments to be passed to or from methods. |
The print method is based on the respective method
print.power.htest of package stats.
A power.mpe.test object is just a named list of numbers and
character strings, supplemented with method and note
elements. The method is displayed as a title, the note
as a footnote, and the remaining elements are given in an aligned
‘name = value’ format.
the argument x, invisibly, as for all print
methods.
The function first appeared in package mpe, which is now archived on CRAN.
Srinath Kolampally, Matthias Kohl [email protected]
print.power.htest, power.mpe.known.var,
power.mpe.unknown.var
(pkv <- power.mpe.known.var(K = 2, delta = c(1,1), Sigma = diag(c(2,2)), power = 0.9, sig.level = 0.025)) print(pkv, digits = 4) # using less digits than default print(pkv, digits = 12) # using more digits than default(pkv <- power.mpe.known.var(K = 2, delta = c(1,1), Sigma = diag(c(2,2)), power = 0.9, sig.level = 0.025)) print(pkv, digits = 4) # using less digits than default print(pkv, digits = 12) # using more digits than default
Produce qq-plot(s) of the given effect size and p values assuming a uniform distribution.
qqunif(x, ...) ## Default S3 method: qqunif(x, min = 0, max = 1, ...) ## S3 method for class 'sim.power.ttest' qqunif(x, color.line = "orange", shape = 19, size = 1, alpha = 1, ...) ## S3 method for class 'sim.power.wtest' qqunif(x, color.line = "orange", shape = 19, size = 1, alpha = 1, ...)qqunif(x, ...) ## Default S3 method: qqunif(x, min = 0, max = 1, ...) ## S3 method for class 'sim.power.ttest' qqunif(x, color.line = "orange", shape = 19, size = 1, alpha = 1, ...) ## S3 method for class 'sim.power.wtest' qqunif(x, color.line = "orange", shape = 19, size = 1, alpha = 1, ...)
x |
numeric vector or data (object). |
min |
single numeric, lower limit of the distribution. |
max |
single numeric, upper limit of the distribution. |
color.line |
color of the line indicating the uniform distribution. |
shape |
point shape. |
size |
point size. |
alpha |
bleding factor (default: no blending. |
... |
further arguments that may be passed through). |
The plot generates a ggplot2 object that is shown.
Missing values are handled by the ggplot2 functions.
Object of class gg and ggplot.
Matthias Kohl [email protected]
## default qqunif(runif(100)) ## visualization of empirical power and type-I-error res1 <- sim.power.t.test(nx = 5, rx = rnorm, rx.H0 = rnorm, ny = 10, ry = function(x) rnorm(x, mean = 3, sd = 3), ry.H0 = function(x) rnorm(x, sd = 3), methods = c("student", "welch")) qqunif(res1, alpha = 0.1) res2 <- sim.power.wilcox.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm) qqunif(res2) res3 <- sim.power.cond.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm) qqunif(res3)## default qqunif(runif(100)) ## visualization of empirical power and type-I-error res1 <- sim.power.t.test(nx = 5, rx = rnorm, rx.H0 = rnorm, ny = 10, ry = function(x) rnorm(x, mean = 3, sd = 3), ry.H0 = function(x) rnorm(x, sd = 3), methods = c("student", "welch")) qqunif(res1, alpha = 0.1) res2 <- sim.power.wilcox.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm) qqunif(res2) res3 <- sim.power.cond.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm) qqunif(res3)
Simulate the empirical power and type-I-error of the two-sample test that follows pre-testing of assumptions.
sim.power.cond.test(nx, rx = NULL, rx.H0 = NULL, ny, ry = NULL, ry.H0 = NULL, sig.level = 0.05, conf.int = FALSE, mu = 0, alternative = c("two.sided", "less", "greater"), iter = 10000, normality = TRUE, variance = TRUE, var.equal = FALSE, exact = TRUE)sim.power.cond.test(nx, rx = NULL, rx.H0 = NULL, ny, ry = NULL, ry.H0 = NULL, sig.level = 0.05, conf.int = FALSE, mu = 0, alternative = c("two.sided", "less", "greater"), iter = 10000, normality = TRUE, variance = TRUE, var.equal = FALSE, exact = TRUE)
nx |
single numeric, sample size of first group. |
rx |
function to simulate the values of first group (assuming H1). |
rx.H0 |
|
ny |
single numeric, sample size of second group. |
ry |
function to simulate the values of second group (assuming H1). |
ry.H0 |
|
alternative |
one- or two-sided test. Can be abbreviated. |
sig.level |
significance level (type I error probability) |
conf.int |
logical, shall confidence intervals be computed. Strongly increases computation time! |
mu |
true value of the location shift for the null hypothesis. |
iter |
single positive integer, number of interations of the simulations. |
normality |
logical, group-wise pre-testing of the assumption of normality
applying the Shapiro-Wilk normality test; see |
variance |
logical, pre-testing of the assumption of equal variances
applying Levene's test; see |
var.equal |
logical, only relevant if |
exact |
logical, compute an exact p-value; see |
If specified, functions rx and ry are used to simulate the data under the
alternative hypothesis H1. If specified, functions rx.H0 and ry.H0
simulte the data unter the null hypothesis H0.
At least one pair of functions has to be provided, either rx and ry
or rx.H0 and ry.H0.
For fast computations functions from package matrixTests are used.
If the test for normality is significant, the Wilcoxon-Mann-Whitney test is computed.
If the test for equal variances is significant, the Welch t-test is used - otherwise, Student't t-test.
If normality is FALSE, the Shapiro-Wilk normality test and
the Wilcoxon-Mann-Whitney test are skipped.
If variance is FALSE, Levene's test is skipped and Student's or
Welch t-test is computed depending on var.equal.
Object of class "sim.power.cond.test" with the results of the
conditional two-sample tests. A list with elements H1 (optional),
H0 (optional) and SetUp.
Matthias Kohl [email protected]
Levene, H. (1960). Robust tests for equality of variances. In Olkin, I., Hotelling, H. et al. (eds.). Contributions to Probability and Statistics: Essays in Honor of Harold Hotelling. Stanford University Press. pp. 278-292.
Mann, H and Withney, D (1947). On a test of whether one of two random variables is stochastically larger than the other. Annals of mathematical Statistics, 18, 50-60.
Shapiro, S. S., Wilk, M. B. (1965). An Analysis of Variance Test for Normality (Complete Samples). Biometrika, 52(3/4), 591-611.
Student (1908). The Probable Error of a Mean. Biometrika, 6(1): 1-25.
Welch, B. L. (1947). The generalization of “Student's” problem when several different population variances are involved. Biometrika, 34 (1-2): 28-35.
Wilcoxon, F (1945). Individual Comparisons by Ranking Methods. Biometrics Bulletin, 1, 80-83.
wilcox.test, LocationTests,
wilcoxon
## Equal variance, small sample size power.t.test(n = 6, power = 0.8) (res <- sim.power.cond.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm)) ## test of normality not significant? table(res$H1$normality) ## test of equal variance not significant? ## (only if test for normality was not significant) table(res$H1$variance, useNA = "ifany") table(res$H0$normality) table(res$H0$variance, useNA = "ifany") ## only power (res <- sim.power.cond.test(nx = 6, rx = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2))) ## only type I error (res <- sim.power.cond.test(nx = 6, rx.H0 = rnorm, ny = 6, ry.H0 = rnorm)) ## Not run: ## only test for equal variances ## assume equal variances -> Student's t-test (res <- sim.power.cond.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm, normality = FALSE, var.equal = TRUE)) ## test of normality not significant? table(res$H1$normality, useNA = "ifany") ## test of equal variance not significant? ## (only if test for normality was not significant) table(res$H1$variance, useNA = "ifany") table(res$H0$normality) table(res$H0$variance, useNA = "ifany") ## only test for equal variances ## assume unequal variances -> Welch's t-test (res <- sim.power.cond.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm, normality = FALSE, var.equal = FALSE)) ## test of normality not significant? table(res$H1$normality, useNA = "ifany") ## test of equal variance not significant? ## (only if test for normality was not significant) table(res$H1$variance, useNA = "ifany") table(res$H0$normality) table(res$H0$variance, useNA = "ifany") ## only test of normality ## assume equal variances -> Student's t-test (res <- sim.power.cond.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm, variance = FALSE, var.equal = TRUE)) ## test of normality not significant? table(res$H1$normality, useNA = "ifany") ## test of equal variance not significant? ## (only if test for normality was not significant) table(res$H1$variance, useNA = "ifany") table(res$H0$normality) table(res$H0$variance, useNA = "ifany") ## only test of normality ## assume unequal variances -> Welch's t-test (res <- sim.power.cond.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm, variance = FALSE, var.equal = FALSE)) ## test of normality not significant? table(res$H1$normality, useNA = "ifany") ## test of equal variance not significant? ## (only if test for normality was not significant) table(res$H1$variance, useNA = "ifany") table(res$H0$normality) table(res$H0$variance, useNA = "ifany") ## End(Not run)## Equal variance, small sample size power.t.test(n = 6, power = 0.8) (res <- sim.power.cond.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm)) ## test of normality not significant? table(res$H1$normality) ## test of equal variance not significant? ## (only if test for normality was not significant) table(res$H1$variance, useNA = "ifany") table(res$H0$normality) table(res$H0$variance, useNA = "ifany") ## only power (res <- sim.power.cond.test(nx = 6, rx = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2))) ## only type I error (res <- sim.power.cond.test(nx = 6, rx.H0 = rnorm, ny = 6, ry.H0 = rnorm)) ## Not run: ## only test for equal variances ## assume equal variances -> Student's t-test (res <- sim.power.cond.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm, normality = FALSE, var.equal = TRUE)) ## test of normality not significant? table(res$H1$normality, useNA = "ifany") ## test of equal variance not significant? ## (only if test for normality was not significant) table(res$H1$variance, useNA = "ifany") table(res$H0$normality) table(res$H0$variance, useNA = "ifany") ## only test for equal variances ## assume unequal variances -> Welch's t-test (res <- sim.power.cond.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm, normality = FALSE, var.equal = FALSE)) ## test of normality not significant? table(res$H1$normality, useNA = "ifany") ## test of equal variance not significant? ## (only if test for normality was not significant) table(res$H1$variance, useNA = "ifany") table(res$H0$normality) table(res$H0$variance, useNA = "ifany") ## only test of normality ## assume equal variances -> Student's t-test (res <- sim.power.cond.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm, variance = FALSE, var.equal = TRUE)) ## test of normality not significant? table(res$H1$normality, useNA = "ifany") ## test of equal variance not significant? ## (only if test for normality was not significant) table(res$H1$variance, useNA = "ifany") table(res$H0$normality) table(res$H0$variance, useNA = "ifany") ## only test of normality ## assume unequal variances -> Welch's t-test (res <- sim.power.cond.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm, variance = FALSE, var.equal = FALSE)) ## test of normality not significant? table(res$H1$normality, useNA = "ifany") ## test of equal variance not significant? ## (only if test for normality was not significant) table(res$H1$variance, useNA = "ifany") table(res$H0$normality) table(res$H0$variance, useNA = "ifany") ## End(Not run)
Simulate the empirical power and type-I-error of two-sample t-tests; i.e., Student (equal variances), Welch, Hsu, Xiao, permutation and boostrap t-tests.
sim.power.t.test(nx, rx = NULL, rx.H0 = NULL, ny, ry = NULL, ry.H0 = NULL, sig.level = 0.05, conf.int = FALSE, mu = 0, alternative = c("two.sided", "less", "greater"), methods = c("student", "welch", "hsu", "xiao"), iter = 10000, R = 9999, useCombn = FALSE, parallel = FALSE, cl = NULL)sim.power.t.test(nx, rx = NULL, rx.H0 = NULL, ny, ry = NULL, ry.H0 = NULL, sig.level = 0.05, conf.int = FALSE, mu = 0, alternative = c("two.sided", "less", "greater"), methods = c("student", "welch", "hsu", "xiao"), iter = 10000, R = 9999, useCombn = FALSE, parallel = FALSE, cl = NULL)
nx |
single numeric, sample size of first group. |
rx |
|
rx.H0 |
|
ny |
single numeric, sample size of second group. |
ry |
|
ry.H0 |
|
sig.level |
significance level (type I error probability) |
conf.int |
logical, shall confidence intervals be computed. Increases computation time! |
mu |
true value of the location shift for the null hypothesis. |
alternative |
one- or two-sided test. Can be abbreviated. |
methods |
character, one or more of the following names: |
iter |
single integer, number of interations of the simulations. |
R |
number of (Monte-Carlo) permutations for permutation and bootstrap tests. |
useCombn |
logical; use all combinations instead of permutations in the
case of the two-sample t-test; see |
parallel |
logical; use parellel computing via function |
cl |
a cluster object, created by package |
If specified, functions rx and ry are used to simulate the data under the
alternative hypothesis H1. If specified, functions rx.H0 and ry.H0
simulte the data unter the null hypothesis H0.
At least one pair of functions has to be provided, either rx and ry
or rx.H0 and ry.H0.
For fast computations functions from package matrixTests are used.
Object of class "sim.power.ttest" with the results of the selected t-tests
in the list elements Student, Welch, Hsu, Xiao,
Perm, Boot. In addition, the simulation setup is saved in
element SetUp.
Matthias Kohl [email protected]
J. Hedderich, L. Sachs. Angewandte Statistik: Methodensammlung mit R. Springer 2018.
Hsu, P. (1938). Contribution to the theory of “student's” t-test as applied to the problem of two samples. Statistical Research Memoirs 2: 1-24.
Student (1908). The Probable Error of a Mean. Biometrika, 6(1): 1-25.
Welch, B. L. (1947). The generalization of “Student's” problem when several different population variances are involved. Biometrika, 34 (1-2): 28-35.
Y. Xiao (2018). On the Solution of a Generalized Behrens-Fisher Problem. Far East Journal of Theoretical Statistics, 54 (1), 21-140.
## Equal variance, small sample size power.t.test(n = 5, delta = 2) power.welch.t.test(n = 5, delta = 2) power.hsu.t.test(n = 5, delta = 2) sim.power.t.test(nx = 5, rx = rnorm, rx.H0 = rnorm, ny = 5, ry = function(n) rnorm(n, mean = 2), ry.H0 = rnorm, methods = c("student", "welch")) ## only power sim.power.t.test(nx = 5, rx = rnorm, ny = 5, ry = function(n) rnorm(n, mean = 2), methods = c("student", "welch")) ## only type I error sim.power.t.test(nx = 5, rx.H0 = rnorm, ny = 5, ry.H0 = rnorm, methods = c("student", "welch")) ## Unequal variance, small sample size power.welch.t.test(n = 5, delta = 5, sd1 = 1, sd2 = 3) power.hsu.t.test(n = 5, delta = 5, sd1 = 1, sd2 = 3) sim.power.t.test(nx = 5, rx = rnorm, rx.H0 = rnorm, ny = 5, ry = function(n) rnorm(n, mean = 5, sd = 3), ry.H0 = function(n) rnorm(n, sd = 3), methods = c("student", "welch")) ## Unequal variance, small unequal sample sizes power.welch.t.test2(n = c(5, 10), delta = 3, sd1 = 1, sd2 = 3) sim.power.t.test(nx = 5, rx = rnorm, rx.H0 = rnorm, ny = 10, ry = function(n) rnorm(n, mean = 3, sd = 3), ry.H0 = function(n) rnorm(n, sd = 3), methods = c("student", "welch")) ## Not run: ## with parallel computing in case of xiao as well as perm. and boot. tests sim.power.t.test(nx = 5, rx = rnorm, rx.H0 = rnorm, ny = 10, ry = function(n) rnorm(n, mean = 3, sd = 3), ry.H0 = function(n) rnorm(n, sd = 3), iter = 1000, methods = c("student", "welch", "hsu", "xiao", "perm.student", "perm.welch", "boot.student", "boot.welch"), parallel = TRUE) ## End(Not run)## Equal variance, small sample size power.t.test(n = 5, delta = 2) power.welch.t.test(n = 5, delta = 2) power.hsu.t.test(n = 5, delta = 2) sim.power.t.test(nx = 5, rx = rnorm, rx.H0 = rnorm, ny = 5, ry = function(n) rnorm(n, mean = 2), ry.H0 = rnorm, methods = c("student", "welch")) ## only power sim.power.t.test(nx = 5, rx = rnorm, ny = 5, ry = function(n) rnorm(n, mean = 2), methods = c("student", "welch")) ## only type I error sim.power.t.test(nx = 5, rx.H0 = rnorm, ny = 5, ry.H0 = rnorm, methods = c("student", "welch")) ## Unequal variance, small sample size power.welch.t.test(n = 5, delta = 5, sd1 = 1, sd2 = 3) power.hsu.t.test(n = 5, delta = 5, sd1 = 1, sd2 = 3) sim.power.t.test(nx = 5, rx = rnorm, rx.H0 = rnorm, ny = 5, ry = function(n) rnorm(n, mean = 5, sd = 3), ry.H0 = function(n) rnorm(n, sd = 3), methods = c("student", "welch")) ## Unequal variance, small unequal sample sizes power.welch.t.test2(n = c(5, 10), delta = 3, sd1 = 1, sd2 = 3) sim.power.t.test(nx = 5, rx = rnorm, rx.H0 = rnorm, ny = 10, ry = function(n) rnorm(n, mean = 3, sd = 3), ry.H0 = function(n) rnorm(n, sd = 3), methods = c("student", "welch")) ## Not run: ## with parallel computing in case of xiao as well as perm. and boot. tests sim.power.t.test(nx = 5, rx = rnorm, rx.H0 = rnorm, ny = 10, ry = function(n) rnorm(n, mean = 3, sd = 3), ry.H0 = function(n) rnorm(n, sd = 3), iter = 1000, methods = c("student", "welch", "hsu", "xiao", "perm.student", "perm.welch", "boot.student", "boot.welch"), parallel = TRUE) ## End(Not run)
Simulate the empirical power and type-I-error of Wilcoxon-Mann-Whitney tests.
sim.power.wilcox.test(nx, rx = NULL, rx.H0 = NULL, ny, ry = NULL, ry.H0 = NULL, alternative = c("two.sided", "less", "greater"), sig.level = 0.05, conf.int = FALSE, approximate = FALSE, ties = FALSE, iter = 10000, nresample = 10000, parallel = "no", ncpus = 1L, cl = NULL)sim.power.wilcox.test(nx, rx = NULL, rx.H0 = NULL, ny, ry = NULL, ry.H0 = NULL, alternative = c("two.sided", "less", "greater"), sig.level = 0.05, conf.int = FALSE, approximate = FALSE, ties = FALSE, iter = 10000, nresample = 10000, parallel = "no", ncpus = 1L, cl = NULL)
nx |
single numeric, sample size of first group. |
rx |
function to simulate the values of first group (assuming H1). |
rx.H0 |
|
ny |
single numeric, sample size of second group. |
ry |
function to simulate the values of second group (assuming H1). |
ry.H0 |
|
alternative |
one- or two-sided test. Can be abbreviated. |
sig.level |
significance level (type I error probability) |
conf.int |
logical, shall confidence intervals be computed. Strongly increases computation time! |
approximate |
logical, shall an approximate test be computed;
see |
ties |
logical, indicating whether ties may occur. Increases computation time! |
iter |
single positive integer, number of interations of the simulations. |
nresample |
single positive integer, the number of Monte Carlo replicates
used for the computation of the approximative reference distribution;
see |
parallel |
a character, the type of parallel operation: either |
ncpus |
a single integer, the number of processes to be used in parallel operation.
Defaults to 1L; see |
cl |
an object inheriting from class |
If specified, functions rx and ry are used to simulate the data under the
alternative hypothesis H1. If specified, functions rx.H0 and ry.H0
simulte the data unter the null hypothesis H0.
At least one pair of functions has to be provided, either rx and ry
or rx.H0 and ry.H0.
For fast computations functions from package matrixTests and package
coin are used.
Object of class "sim.power.wtest" with the results of the
Wilcoxon-Mann-Whitney tests. A list elements Exact, Asymptotic
and Approximate. In addition, the simulation setup is saved in
element SetUp.
Matthias Kohl [email protected]
Mann, H and Withney, D (1947). On a test of whether one of two random variables is stochastically larger than the other. Annals of mathematical Statistics, 18, 50-60.
Wilcoxon, F (1945). Individual Comparisons by Ranking Methods. Biometrics Bulletin, 1, 80-83.
wilcox.test, LocationTests,
wilcoxon
## Equal variance, small sample size power.t.test(n = 6, power = 0.8) sim.ssize.wilcox.test(rx = rnorm, ry = function(x) rnorm(x, mean = 2), power = 0.8, n.min = 3, n.max = 10, step.size = 1) sim.power.wilcox.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm) ## only power sim.power.wilcox.test(nx = 6, rx = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2)) ## only type I error sim.power.wilcox.test(nx = 6, rx.H0 = rnorm, ny = 6, ry.H0 = rnorm)## Equal variance, small sample size power.t.test(n = 6, power = 0.8) sim.ssize.wilcox.test(rx = rnorm, ry = function(x) rnorm(x, mean = 2), power = 0.8, n.min = 3, n.max = 10, step.size = 1) sim.power.wilcox.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm) ## only power sim.power.wilcox.test(nx = 6, rx = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2)) ## only type I error sim.power.wilcox.test(nx = 6, rx.H0 = rnorm, ny = 6, ry.H0 = rnorm)
Use Monte-Carlo simulations to determine the empirical power and sample size of the Xiao t-test.
sim.power.xiao.t.test(n1, n2, delta = 1, sd1 = 1, sd2 = 2, mu = 0, sig.level = 0.05, alternative = c("two.sided", "less", "greater"), iter = 10000, parallel = FALSE, cl = NULL) sim.ssize.xiao.t.test(delta = 1, sd1 = 1, sd2 = 2, mu = 0, sig.level = 0.05, power = 0.8, alternative = c("two.sided", "one.sided"), n.delta = 5, step.size = 1, iter = 10000, parallel = FALSE, cl = NULL)sim.power.xiao.t.test(n1, n2, delta = 1, sd1 = 1, sd2 = 2, mu = 0, sig.level = 0.05, alternative = c("two.sided", "less", "greater"), iter = 10000, parallel = FALSE, cl = NULL) sim.ssize.xiao.t.test(delta = 1, sd1 = 1, sd2 = 2, mu = 0, sig.level = 0.05, power = 0.8, alternative = c("two.sided", "one.sided"), n.delta = 5, step.size = 1, iter = 10000, parallel = FALSE, cl = NULL)
n1 |
number of observations of group 1 |
n2 |
number of observations of group 2 |
delta |
(expected) true difference in means |
sd1 |
(expected) standard deviation of group 1 |
sd2 |
(expected) standard deviation of group 2 |
mu |
true value of the location shift for the null hypothesis. |
sig.level |
significance level (Type I error probability) |
power |
power of test (1 minus Type II error probability) |
alternative |
one- or two-sided test. |
iter |
single integer, number of interations of the simulations. |
parallel |
logical; use parellel computing via function |
cl |
a cluster object, created by package |
n.delta |
single integer, neighborhood around the values obtained by
|
step.size |
step size for the neighborhood around the values obtained
by |
To find the (optimal) group sample sizes a grid search in a (small) neighborhood around the (optimal) values of the Welch t-test is performed.
Object of class "power.htest", a list of the arguments
(including the computed one) augmented with method and
note elements.
Matthias Kohl [email protected]
Y. Xiao (2018). On the Solution of a Generalized Behrens-Fisher Problem. Far East Journal of Theoretical Statistics, 54 (1), 21-140.
ssize.welch.t.test, xiao.t.test
## Not run: ## delta = 0.5, sd1 = 0.5, sd2 = 1, power.welch.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) power.hsu.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) power.welch.t.test2(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) # 22+87 = 109 power.hsu.t.test2(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) # 44+59 = 103 power.hsu.t.test2(n = c(44, 59), delta = 0.5, sd1 = 0.5, sd2 = 1) sim.power.xiao.t.test(n1 = 22, n2 = 87, delta = 0.5, sd1 = 0.5, sd2 = 1, parallel = TRUE) ssize.welch.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) # 35+62 = 97 ssize.hsu.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) # 39+62 = 101 sim.ssize.xiao.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9, parallel = TRUE) ## End(Not run)## Not run: ## delta = 0.5, sd1 = 0.5, sd2 = 1, power.welch.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) power.hsu.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) power.welch.t.test2(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) # 22+87 = 109 power.hsu.t.test2(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) # 44+59 = 103 power.hsu.t.test2(n = c(44, 59), delta = 0.5, sd1 = 0.5, sd2 = 1) sim.power.xiao.t.test(n1 = 22, n2 = 87, delta = 0.5, sd1 = 0.5, sd2 = 1, parallel = TRUE) ssize.welch.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) # 35+62 = 97 ssize.hsu.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9) # 39+62 = 101 sim.ssize.xiao.t.test(delta = 0.5, sd1 = 0.5, sd2 = 1, power = 0.9, parallel = TRUE) ## End(Not run)
Simulate the empirical power of bootstrap t tests for computing the required sample size.
sim.ssize.boot.t.test(rx, ry = NULL, mu = 0, sig.level = 0.05, power = 0.8, type = c("two.sample", "one.sample", "paired"), alternative = c("two.sided", "less", "greater"), var.equal = FALSE, R = 9999, symmetric = TRUE, n.min = 10, n.max = 200, step.size = 10, iter = 10000, BREAK = TRUE, parallel = FALSE, cl = NULL)sim.ssize.boot.t.test(rx, ry = NULL, mu = 0, sig.level = 0.05, power = 0.8, type = c("two.sample", "one.sample", "paired"), alternative = c("two.sided", "less", "greater"), var.equal = FALSE, R = 9999, symmetric = TRUE, n.min = 10, n.max = 200, step.size = 10, iter = 10000, BREAK = TRUE, parallel = FALSE, cl = NULL)
rx |
function to simulate the values of x, respectively x-y in the paired case. |
ry |
function to simulate the values of y in the two-sample case |
mu |
true values of the location shift for the null hypothesis. |
sig.level |
significance level (Type I error probability) |
power |
two-sample, one-sample or paired test |
type |
one- or two-sided test. Can be abbreviated. |
alternative |
one- or two-sided test. Can be abbreviated. |
var.equal |
a logical variable indicating whether to treat the
two variances as being equal. If |
R |
number of (Monte-Carlo) permutations. |
symmetric |
a logical variable indicating whether to assume symmetry
in the two-sided test. If |
n.min |
integer, start value of grid search. |
n.max |
integer, stop value of grid search. |
step.size |
integer, step size used in the grid search. |
iter |
integer, number of interations of the simulations. |
BREAK |
logical, grid search stops when the emperical power is larger than the requested power. |
parallel |
logical; use parellel computing via function |
cl |
a cluster object, created by package |
Functions rx and ry are used to simulate the data.
We recommend a two steps procedure: In the first step, start with a wide grid
and find out in which range of sample size values the intended power will
be achieved. In the second step, the interval identified in the first step
is used to find the sample size that leads to the required power setting
step.size = 1 and BREAK = FALSE. This approach is applied
in the examples below.
Object of class "power.htest", a list of the arguments
(including the computed one) augmented with method and
note elements.
Matthias Kohl [email protected]
B. Efron, R.J. Tibshirani. An Introduction to the Bootstrap. Chapman and Hall/CRC 1993.
A. Janssen, T. Pauls. How do bootstrap and permutation tests work?. Ann. Statist., 31(3), 768-806.
## Not run: ############################################################################### ## two-sample ## iter = 100 to reduce computation time ############################################################################### rx <- function(n) rnorm(n, mean = 0, sd = 1) ry <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.boot.t.test(rx = rx, ry = ry, n.min = 50, n.max = 100, iter = 1000, var.equal = TRUE) sim.ssize.boot.t.test(rx = rx, ry = ry, n.min = 50, n.max = 100, iter = 1000, var.equal = TRUE, parallel = TRUE) sim.ssize.boot.t.test(rx = rx, ry = ry, n.min = 60, n.max = 65, step.size = 1, iter = 1000, BREAK = FALSE, var.equal = TRUE, parallel = TRUE) ## compared to power.t.test(delta = 0.5, power = 0.8) rx <- function(n) rnorm(n, mean = 0, sd = 1) ry <- function(n) rnorm(n, mean = 0.5, sd = 1.5) sim.ssize.boot.t.test(rx = rx, ry = ry, n.max = 100, iter = 1000, alternative = "less") sim.ssize.boot.t.test(rx = rx, ry = ry, n.min = 85, n.max = 90, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "less") ## compared to power.welch.t.test(delta = 0.5, sd = 1, sd2 = 1.5, power = 0.8, alternative = "one.sided") rx <- function(n) rnorm(n, mean = 0.5, sd = 1) ry <- function(n) rnorm(n, mean = 0, sd = 1) sim.ssize.boot.t.test(rx = rx, ry = ry, n.max = 100, iter = 1000, alternative = "greater") sim.ssize.boot.t.test(rx = rx, ry = ry, n.min = 50, n.max = 55, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "greater") ## compared to power.t.test(delta = 0.5, power = 0.8, alternative = "one.sided") rx <- function(n) rgamma(n, scale = 10, shape = 1) ry <- function(n) rgamma(n, scale = 15, shape = 1) sim.ssize.boot.t.test(rx = rx, ry = ry, n.max = 200, iter = 1000) sim.ssize.boot.t.test(rx = rx, ry = ry, n.min = 125, n.max = 135, step.size = 1, iter = 1000, BREAK = FALSE) ############################################################################### ## one-sample ## iter = 1000 to reduce computation time ############################################################################### rx <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.boot.t.test(rx = rx, mu = 0, type = "one.sample", n.max = 100, iter = 1000) sim.ssize.boot.t.test(rx = rx, mu = 0, type = "one.sample", n.min = 33, n.max = 38, step.size = 1, iter = 1000, BREAK = FALSE) ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample") sim.ssize.boot.t.test(rx = rx, mu = 0, type = "one.sample", n.max = 100, iter = 1000, alternative = "greater") sim.ssize.boot.t.test(rx = rx, mu = 0, type = "one.sample", n.min = 25, n.max = 30, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "greater") ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample", alternative = "one.sided") sim.ssize.boot.t.test(rx = rx, mu = 1, type = "one.sample", n.max = 100, iter = 1000, alternative = "less") sim.ssize.boot.t.test(rx = rx, mu = 1, type = "one.sample", n.min = 20, n.max = 30, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "less") ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample", alternative = "one.sided") rx <- function(n) rgamma(n, scale = 10, shape = 1) sim.ssize.boot.t.test(rx = rx, mu = 5, type = "one.sample", n.max = 200, iter = 1000) sim.ssize.boot.t.test(rx = rx, mu = 5, type = "one.sample", n.min = 40, n.max = 50, step.size = 1, iter = 1000, BREAK = FALSE) ############################################################################### ## paired ## identical to one-sample, requires random number generating function ## that simulates the difference x-y ## iter = 1000 to reduce computation time ############################################################################### rxy <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.boot.t.test(rx = rxy, mu = 0, type = "paired", n.max = 100, iter = 1000) sim.ssize.boot.t.test(rx = rxy, mu = 0, type = "paired", n.min = 33, n.max = 38, step.size = 1, iter = 1000, BREAK = FALSE) ## compared to power.t.test(delta = 0.5, power = 0.8, type = "paired") ## End(Not run)## Not run: ############################################################################### ## two-sample ## iter = 100 to reduce computation time ############################################################################### rx <- function(n) rnorm(n, mean = 0, sd = 1) ry <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.boot.t.test(rx = rx, ry = ry, n.min = 50, n.max = 100, iter = 1000, var.equal = TRUE) sim.ssize.boot.t.test(rx = rx, ry = ry, n.min = 50, n.max = 100, iter = 1000, var.equal = TRUE, parallel = TRUE) sim.ssize.boot.t.test(rx = rx, ry = ry, n.min = 60, n.max = 65, step.size = 1, iter = 1000, BREAK = FALSE, var.equal = TRUE, parallel = TRUE) ## compared to power.t.test(delta = 0.5, power = 0.8) rx <- function(n) rnorm(n, mean = 0, sd = 1) ry <- function(n) rnorm(n, mean = 0.5, sd = 1.5) sim.ssize.boot.t.test(rx = rx, ry = ry, n.max = 100, iter = 1000, alternative = "less") sim.ssize.boot.t.test(rx = rx, ry = ry, n.min = 85, n.max = 90, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "less") ## compared to power.welch.t.test(delta = 0.5, sd = 1, sd2 = 1.5, power = 0.8, alternative = "one.sided") rx <- function(n) rnorm(n, mean = 0.5, sd = 1) ry <- function(n) rnorm(n, mean = 0, sd = 1) sim.ssize.boot.t.test(rx = rx, ry = ry, n.max = 100, iter = 1000, alternative = "greater") sim.ssize.boot.t.test(rx = rx, ry = ry, n.min = 50, n.max = 55, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "greater") ## compared to power.t.test(delta = 0.5, power = 0.8, alternative = "one.sided") rx <- function(n) rgamma(n, scale = 10, shape = 1) ry <- function(n) rgamma(n, scale = 15, shape = 1) sim.ssize.boot.t.test(rx = rx, ry = ry, n.max = 200, iter = 1000) sim.ssize.boot.t.test(rx = rx, ry = ry, n.min = 125, n.max = 135, step.size = 1, iter = 1000, BREAK = FALSE) ############################################################################### ## one-sample ## iter = 1000 to reduce computation time ############################################################################### rx <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.boot.t.test(rx = rx, mu = 0, type = "one.sample", n.max = 100, iter = 1000) sim.ssize.boot.t.test(rx = rx, mu = 0, type = "one.sample", n.min = 33, n.max = 38, step.size = 1, iter = 1000, BREAK = FALSE) ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample") sim.ssize.boot.t.test(rx = rx, mu = 0, type = "one.sample", n.max = 100, iter = 1000, alternative = "greater") sim.ssize.boot.t.test(rx = rx, mu = 0, type = "one.sample", n.min = 25, n.max = 30, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "greater") ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample", alternative = "one.sided") sim.ssize.boot.t.test(rx = rx, mu = 1, type = "one.sample", n.max = 100, iter = 1000, alternative = "less") sim.ssize.boot.t.test(rx = rx, mu = 1, type = "one.sample", n.min = 20, n.max = 30, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "less") ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample", alternative = "one.sided") rx <- function(n) rgamma(n, scale = 10, shape = 1) sim.ssize.boot.t.test(rx = rx, mu = 5, type = "one.sample", n.max = 200, iter = 1000) sim.ssize.boot.t.test(rx = rx, mu = 5, type = "one.sample", n.min = 40, n.max = 50, step.size = 1, iter = 1000, BREAK = FALSE) ############################################################################### ## paired ## identical to one-sample, requires random number generating function ## that simulates the difference x-y ## iter = 1000 to reduce computation time ############################################################################### rxy <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.boot.t.test(rx = rxy, mu = 0, type = "paired", n.max = 100, iter = 1000) sim.ssize.boot.t.test(rx = rxy, mu = 0, type = "paired", n.min = 33, n.max = 38, step.size = 1, iter = 1000, BREAK = FALSE) ## compared to power.t.test(delta = 0.5, power = 0.8, type = "paired") ## End(Not run)
Simulate the empirical power of permutation t tests for computing the required sample size.
sim.ssize.perm.t.test(rx, ry = NULL, mu = 0, sig.level = 0.05, power = 0.8, type = c("two.sample", "one.sample", "paired"), alternative = c("two.sided", "less", "greater"), var.equal = FALSE, R = 9999, symmetric = TRUE, useCombn = FALSE, n.min = 10, n.max = 200, step.size = 10, iter = 10000, BREAK = TRUE, parallel = FALSE, cl = NULL)sim.ssize.perm.t.test(rx, ry = NULL, mu = 0, sig.level = 0.05, power = 0.8, type = c("two.sample", "one.sample", "paired"), alternative = c("two.sided", "less", "greater"), var.equal = FALSE, R = 9999, symmetric = TRUE, useCombn = FALSE, n.min = 10, n.max = 200, step.size = 10, iter = 10000, BREAK = TRUE, parallel = FALSE, cl = NULL)
rx |
function to simulate the values of x, respectively x-y in the paired case. |
ry |
function to simulate the values of y in the two-sample case |
mu |
true values of the location shift for the null hypothesis. |
sig.level |
significance level (Type I error probability) |
power |
two-sample, one-sample or paired test |
type |
one- or two-sided test. Can be abbreviated. |
alternative |
one- or two-sided test. Can be abbreviated. |
var.equal |
a logical variable indicating whether to treat the
two variances as being equal. If |
R |
number of (Monte-Carlo) permutations. |
symmetric |
a logical variable indicating whether to assume symmetry
in the two-sided test. If |
useCombn |
logical; use all combinations instead of permutations in the
case of the two-sample t-test; see |
n.min |
integer, start value of grid search. |
n.max |
integer, stop value of grid search. |
step.size |
integer, step size used in the grid search. |
iter |
integer, number of interations of the simulations. |
BREAK |
logical, grid search stops when the emperical power is larger than the requested power. |
parallel |
logical; use parellel computing via function |
cl |
a cluster object, created by package |
Functions rx and ry are used to simulate the data.
We recommend a two steps procedure: In the first step, start with a wide grid
and find out in which range of sample size values the intended power will
be achieved. In the second step, the interval identified in the first step
is used to find the sample size that leads to the required power setting
step.size = 1 and BREAK = FALSE. This approach is applied
in the examples below.
Object of class "power.htest", a list of the arguments
(including the computed one) augmented with method and
note elements.
Matthias Kohl [email protected]
B. Efron, R.J. Tibshirani. An Introduction to the Bootstrap. Chapman and Hall/CRC 1993.
A. Janssen (1997). Studentized permutation tests for non-i.i.d, hypotheses and the generalized Behrens-Fisher problem. Statistics and Probability Letters, 36, 9-21.
A. Janssen, T. Pauls. How do bootstrap and permutation tests work?. Ann. Statist., 31(3), 768-806.
E. Chung, J.P. Romano (2013). Exact and asymptotically robust permutation tests. The Annals of Statistics, 41(2), 484-507.
## Not run: ############################################################################### ## two-sample ## iter = 100 to reduce computation time ############################################################################### rx <- function(n) rnorm(n, mean = 0, sd = 1) ry <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.perm.t.test(rx = rx, ry = ry, n.min = 50, n.max = 100, iter = 1000, var.equal = TRUE) sim.ssize.perm.t.test(rx = rx, ry = ry, n.min = 50, n.max = 100, iter = 1000, var.equal = TRUE, parallel = TRUE) sim.ssize.perm.t.test(rx = rx, ry = ry, n.min = 60, n.max = 65, step.size = 1, iter = 1000, BREAK = FALSE, var.equal = TRUE, parallel = TRUE) ## compared to power.t.test(delta = 0.5, power = 0.8) rx <- function(n) rnorm(n, mean = 0, sd = 1) ry <- function(n) rnorm(n, mean = 0.5, sd = 1.5) sim.ssize.perm.t.test(rx = rx, ry = ry, n.max = 100, iter = 1000, alternative = "less") sim.ssize.perm.t.test(rx = rx, ry = ry, n.min = 85, n.max = 90, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "less") ## compared to power.welch.t.test(delta = 0.5, sd = 1, sd2 = 1.5, power = 0.8, alternative = "one.sided") rx <- function(n) rnorm(n, mean = 0.5, sd = 1) ry <- function(n) rnorm(n, mean = 0, sd = 1) sim.ssize.perm.t.test(rx = rx, ry = ry, n.max = 100, iter = 1000, alternative = "greater") sim.ssize.perm.t.test(rx = rx, ry = ry, n.min = 50, n.max = 55, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "greater") ## compared to power.t.test(delta = 0.5, power = 0.8, alternative = "one.sided") rx <- function(n) rgamma(n, scale = 10, shape = 1) ry <- function(n) rgamma(n, scale = 15, shape = 1) sim.ssize.perm.t.test(rx = rx, ry = ry, n.max = 200, iter = 1000) sim.ssize.perm.t.test(rx = rx, ry = ry, n.min = 125, n.max = 135, step.size = 1, iter = 1000, BREAK = FALSE) ############################################################################### ## one-sample ## iter = 1000 to reduce computation time ############################################################################### rx <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.perm.t.test(rx = rx, mu = 0, type = "one.sample", n.max = 100, iter = 1000) sim.ssize.perm.t.test(rx = rx, mu = 0, type = "one.sample", n.min = 33, n.max = 38, step.size = 1, iter = 1000, BREAK = FALSE) ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample") sim.ssize.perm.t.test(rx = rx, mu = 0, type = "one.sample", n.max = 100, iter = 1000, alternative = "greater") sim.ssize.perm.t.test(rx = rx, mu = 0, type = "one.sample", n.min = 25, n.max = 30, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "greater") ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample", alternative = "one.sided") sim.ssize.perm.t.test(rx = rx, mu = 1, type = "one.sample", n.max = 100, iter = 1000, alternative = "less") sim.ssize.perm.t.test(rx = rx, mu = 1, type = "one.sample", n.min = 20, n.max = 30, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "less") ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample", alternative = "one.sided") rx <- function(n) rgamma(n, scale = 10, shape = 1) sim.ssize.perm.t.test(rx = rx, mu = 5, type = "one.sample", n.max = 200, iter = 1000) sim.ssize.perm.t.test(rx = rx, mu = 5, type = "one.sample", n.min = 40, n.max = 50, step.size = 1, iter = 1000, BREAK = FALSE) ############################################################################### ## paired ## identical to one-sample, requires random number generating function ## that simulates the difference x-y ## iter = 1000 to reduce computation time ############################################################################### rxy <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.perm.t.test(rx = rxy, mu = 0, type = "paired", n.max = 100, iter = 1000) sim.ssize.perm.t.test(rx = rxy, mu = 0, type = "paired", n.min = 33, n.max = 38, step.size = 1, iter = 1000, BREAK = FALSE) ## compared to power.t.test(delta = 0.5, power = 0.8, type = "paired") ## End(Not run)## Not run: ############################################################################### ## two-sample ## iter = 100 to reduce computation time ############################################################################### rx <- function(n) rnorm(n, mean = 0, sd = 1) ry <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.perm.t.test(rx = rx, ry = ry, n.min = 50, n.max = 100, iter = 1000, var.equal = TRUE) sim.ssize.perm.t.test(rx = rx, ry = ry, n.min = 50, n.max = 100, iter = 1000, var.equal = TRUE, parallel = TRUE) sim.ssize.perm.t.test(rx = rx, ry = ry, n.min = 60, n.max = 65, step.size = 1, iter = 1000, BREAK = FALSE, var.equal = TRUE, parallel = TRUE) ## compared to power.t.test(delta = 0.5, power = 0.8) rx <- function(n) rnorm(n, mean = 0, sd = 1) ry <- function(n) rnorm(n, mean = 0.5, sd = 1.5) sim.ssize.perm.t.test(rx = rx, ry = ry, n.max = 100, iter = 1000, alternative = "less") sim.ssize.perm.t.test(rx = rx, ry = ry, n.min = 85, n.max = 90, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "less") ## compared to power.welch.t.test(delta = 0.5, sd = 1, sd2 = 1.5, power = 0.8, alternative = "one.sided") rx <- function(n) rnorm(n, mean = 0.5, sd = 1) ry <- function(n) rnorm(n, mean = 0, sd = 1) sim.ssize.perm.t.test(rx = rx, ry = ry, n.max = 100, iter = 1000, alternative = "greater") sim.ssize.perm.t.test(rx = rx, ry = ry, n.min = 50, n.max = 55, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "greater") ## compared to power.t.test(delta = 0.5, power = 0.8, alternative = "one.sided") rx <- function(n) rgamma(n, scale = 10, shape = 1) ry <- function(n) rgamma(n, scale = 15, shape = 1) sim.ssize.perm.t.test(rx = rx, ry = ry, n.max = 200, iter = 1000) sim.ssize.perm.t.test(rx = rx, ry = ry, n.min = 125, n.max = 135, step.size = 1, iter = 1000, BREAK = FALSE) ############################################################################### ## one-sample ## iter = 1000 to reduce computation time ############################################################################### rx <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.perm.t.test(rx = rx, mu = 0, type = "one.sample", n.max = 100, iter = 1000) sim.ssize.perm.t.test(rx = rx, mu = 0, type = "one.sample", n.min = 33, n.max = 38, step.size = 1, iter = 1000, BREAK = FALSE) ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample") sim.ssize.perm.t.test(rx = rx, mu = 0, type = "one.sample", n.max = 100, iter = 1000, alternative = "greater") sim.ssize.perm.t.test(rx = rx, mu = 0, type = "one.sample", n.min = 25, n.max = 30, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "greater") ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample", alternative = "one.sided") sim.ssize.perm.t.test(rx = rx, mu = 1, type = "one.sample", n.max = 100, iter = 1000, alternative = "less") sim.ssize.perm.t.test(rx = rx, mu = 1, type = "one.sample", n.min = 20, n.max = 30, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "less") ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample", alternative = "one.sided") rx <- function(n) rgamma(n, scale = 10, shape = 1) sim.ssize.perm.t.test(rx = rx, mu = 5, type = "one.sample", n.max = 200, iter = 1000) sim.ssize.perm.t.test(rx = rx, mu = 5, type = "one.sample", n.min = 40, n.max = 50, step.size = 1, iter = 1000, BREAK = FALSE) ############################################################################### ## paired ## identical to one-sample, requires random number generating function ## that simulates the difference x-y ## iter = 1000 to reduce computation time ############################################################################### rxy <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.perm.t.test(rx = rxy, mu = 0, type = "paired", n.max = 100, iter = 1000) sim.ssize.perm.t.test(rx = rxy, mu = 0, type = "paired", n.min = 33, n.max = 38, step.size = 1, iter = 1000, BREAK = FALSE) ## compared to power.t.test(delta = 0.5, power = 0.8, type = "paired") ## End(Not run)
Simulate the empirical power of Wilcoxon rank sum and signed rank tests for computing the required sample size.
sim.ssize.wilcox.test(rx, ry = NULL, mu = 0, sig.level = 0.05, power = 0.8, type = c("two.sample", "one.sample", "paired"), alternative = c("two.sided", "less", "greater"), n.min = 10, n.max = 200, step.size = 10, iter = 10000, BREAK = TRUE, exact = NA, correct = TRUE)sim.ssize.wilcox.test(rx, ry = NULL, mu = 0, sig.level = 0.05, power = 0.8, type = c("two.sample", "one.sample", "paired"), alternative = c("two.sided", "less", "greater"), n.min = 10, n.max = 200, step.size = 10, iter = 10000, BREAK = TRUE, exact = NA, correct = TRUE)
rx |
function to simulate the values of x, respectively x-y in the paired case. |
ry |
function to simulate the values of y in the two-sample case |
mu |
true values of the location shift for the null hypothesis. |
sig.level |
significance level (Type I error probability) |
power |
two-sample, one-sample or paired test |
type |
one- or two-sided test. Can be abbreviated. |
alternative |
one- or two-sided test. Can be abbreviated. |
n.min |
integer, start value of grid search. |
n.max |
integer, stop value of grid search. |
step.size |
integer, step size used in the grid search. |
iter |
integer, number of interations of the simulations. |
BREAK |
logical, grid search stops when the emperical power is larger than the requested power. |
exact |
logical or NA (default) indicator whether an exact p-value should
be computed (see Details at |
correct |
logical indicator whether continuity correction should be applied
in the cases where p-values are obtained using normal approximation. A single
value or logical vector with values for each observation; see
|
Functions rx and ry are used to simulate the data and
functions row_wilcoxon_twosample and row_wilcoxon_onesample of
package matrixTests are used to efficiently compute the p values of the
respective test.
We recommend a two steps procedure: In the first step, start with a wide grid
and find out in which range of sample size values the intended power will
be achieved. In the second step, the interval identified in the first step
is used to find the sample size that leads to the required power setting
step.size = 1 and BREAK = FALSE. This approach is applied
in the examples below.
Object of class "power.htest", a list of the arguments
(including the computed one) augmented with method and
note elements.
Matthias Kohl [email protected]
Wilcoxon, F (1945). Individual Comparisons by Ranking Methods. Biometrics Bulletin, 1, 80-83.
############################################################################### ## two-sample ## iter = 1000 to reduce check time ############################################################################### rx <- function(n) rnorm(n, mean = 0, sd = 1) ry <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.wilcox.test(rx = rx, ry = ry, n.max = 100, iter = 1000) sim.ssize.wilcox.test(rx = rx, ry = ry, n.min = 65, n.max = 70, step.size = 1, iter = 1000, BREAK = FALSE) ## compared to power.t.test(delta = 0.5, power = 0.8) rx <- function(n) rnorm(n, mean = 0, sd = 1) ry <- function(n) rnorm(n, mean = 0.5, sd = 1.5) sim.ssize.wilcox.test(rx = rx, ry = ry, n.max = 100, iter = 1000, alternative = "less") sim.ssize.wilcox.test(rx = rx, ry = ry, n.min = 85, n.max = 90, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "less") ## compared to power.welch.t.test(delta = 0.5, sd = 1, sd2 = 1.5, power = 0.8, alternative = "one.sided") rx <- function(n) rnorm(n, mean = 0.5, sd = 1) ry <- function(n) rnorm(n, mean = 0, sd = 1) sim.ssize.wilcox.test(rx = rx, ry = ry, n.max = 100, iter = 1000, alternative = "greater") sim.ssize.wilcox.test(rx = rx, ry = ry, n.min = 50, n.max = 55, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "greater") ## compared to power.t.test(delta = 0.5, power = 0.8, alternative = "one.sided") rx <- function(n) rgamma(n, scale = 10, shape = 1) ry <- function(n) rgamma(n, scale = 15, shape = 1) sim.ssize.wilcox.test(rx = rx, ry = ry, n.max = 200, iter = 1000) sim.ssize.wilcox.test(rx = rx, ry = ry, n.min = 125, n.max = 135, step.size = 1, iter = 1000, BREAK = FALSE) ############################################################################### ## one-sample ## iter = 1000 to reduce check time ############################################################################### rx <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.wilcox.test(rx = rx, mu = 0, type = "one.sample", n.max = 100, iter = 1000) sim.ssize.wilcox.test(rx = rx, mu = 0, type = "one.sample", n.min = 33, n.max = 38, step.size = 1, iter = 1000, BREAK = FALSE) ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample") sim.ssize.wilcox.test(rx = rx, mu = 0, type = "one.sample", n.max = 100, iter = 1000, alternative = "greater") sim.ssize.wilcox.test(rx = rx, mu = 0, type = "one.sample", n.min = 25, n.max = 30, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "greater") ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample", alternative = "one.sided") sim.ssize.wilcox.test(rx = rx, mu = 1, type = "one.sample", n.max = 100, iter = 1000, alternative = "less") sim.ssize.wilcox.test(rx = rx, mu = 1, type = "one.sample", n.min = 20, n.max = 30, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "less") ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample", alternative = "one.sided") rx <- function(n) rgamma(n, scale = 10, shape = 1) sim.ssize.wilcox.test(rx = rx, mu = 5, type = "one.sample", n.max = 200, iter = 1000) sim.ssize.wilcox.test(rx = rx, mu = 5, type = "one.sample", n.min = 40, n.max = 50, step.size = 1, iter = 1000, BREAK = FALSE) ############################################################################### ## paired ## identical to one-sample, requires random number generating function ## that simulates the difference x-y ## iter = 1000 to reduce check time ############################################################################### rxy <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.wilcox.test(rx = rxy, mu = 0, type = "paired", n.max = 100, iter = 1000) sim.ssize.wilcox.test(rx = rxy, mu = 0, type = "paired", n.min = 33, n.max = 38, step.size = 1, iter = 1000, BREAK = FALSE) ## compared to power.t.test(delta = 0.5, power = 0.8, type = "paired")############################################################################### ## two-sample ## iter = 1000 to reduce check time ############################################################################### rx <- function(n) rnorm(n, mean = 0, sd = 1) ry <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.wilcox.test(rx = rx, ry = ry, n.max = 100, iter = 1000) sim.ssize.wilcox.test(rx = rx, ry = ry, n.min = 65, n.max = 70, step.size = 1, iter = 1000, BREAK = FALSE) ## compared to power.t.test(delta = 0.5, power = 0.8) rx <- function(n) rnorm(n, mean = 0, sd = 1) ry <- function(n) rnorm(n, mean = 0.5, sd = 1.5) sim.ssize.wilcox.test(rx = rx, ry = ry, n.max = 100, iter = 1000, alternative = "less") sim.ssize.wilcox.test(rx = rx, ry = ry, n.min = 85, n.max = 90, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "less") ## compared to power.welch.t.test(delta = 0.5, sd = 1, sd2 = 1.5, power = 0.8, alternative = "one.sided") rx <- function(n) rnorm(n, mean = 0.5, sd = 1) ry <- function(n) rnorm(n, mean = 0, sd = 1) sim.ssize.wilcox.test(rx = rx, ry = ry, n.max = 100, iter = 1000, alternative = "greater") sim.ssize.wilcox.test(rx = rx, ry = ry, n.min = 50, n.max = 55, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "greater") ## compared to power.t.test(delta = 0.5, power = 0.8, alternative = "one.sided") rx <- function(n) rgamma(n, scale = 10, shape = 1) ry <- function(n) rgamma(n, scale = 15, shape = 1) sim.ssize.wilcox.test(rx = rx, ry = ry, n.max = 200, iter = 1000) sim.ssize.wilcox.test(rx = rx, ry = ry, n.min = 125, n.max = 135, step.size = 1, iter = 1000, BREAK = FALSE) ############################################################################### ## one-sample ## iter = 1000 to reduce check time ############################################################################### rx <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.wilcox.test(rx = rx, mu = 0, type = "one.sample", n.max = 100, iter = 1000) sim.ssize.wilcox.test(rx = rx, mu = 0, type = "one.sample", n.min = 33, n.max = 38, step.size = 1, iter = 1000, BREAK = FALSE) ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample") sim.ssize.wilcox.test(rx = rx, mu = 0, type = "one.sample", n.max = 100, iter = 1000, alternative = "greater") sim.ssize.wilcox.test(rx = rx, mu = 0, type = "one.sample", n.min = 25, n.max = 30, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "greater") ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample", alternative = "one.sided") sim.ssize.wilcox.test(rx = rx, mu = 1, type = "one.sample", n.max = 100, iter = 1000, alternative = "less") sim.ssize.wilcox.test(rx = rx, mu = 1, type = "one.sample", n.min = 20, n.max = 30, step.size = 1, iter = 1000, BREAK = FALSE, alternative = "less") ## compared to power.t.test(delta = 0.5, power = 0.8, type = "one.sample", alternative = "one.sided") rx <- function(n) rgamma(n, scale = 10, shape = 1) sim.ssize.wilcox.test(rx = rx, mu = 5, type = "one.sample", n.max = 200, iter = 1000) sim.ssize.wilcox.test(rx = rx, mu = 5, type = "one.sample", n.min = 40, n.max = 50, step.size = 1, iter = 1000, BREAK = FALSE) ############################################################################### ## paired ## identical to one-sample, requires random number generating function ## that simulates the difference x-y ## iter = 1000 to reduce check time ############################################################################### rxy <- function(n) rnorm(n, mean = 0.5, sd = 1) sim.ssize.wilcox.test(rx = rxy, mu = 0, type = "paired", n.max = 100, iter = 1000) sim.ssize.wilcox.test(rx = rxy, mu = 0, type = "paired", n.min = 33, n.max = 38, step.size = 1, iter = 1000, BREAK = FALSE) ## compared to power.t.test(delta = 0.5, power = 0.8, type = "paired")
Compute sample size, power, delta, or significance level of a diagnostic test for an expected AUC.
ssize.auc.ci(AUC = NULL, delta = NULL, n = NULL, sig.level = 0.05, power = NULL, prev = NULL, NMAX = 1e4)ssize.auc.ci(AUC = NULL, delta = NULL, n = NULL, sig.level = 0.05, power = NULL, prev = NULL, NMAX = 1e4)
AUC |
Expected AUC. |
n |
Total sample size (number of cases + number of controls). |
delta |
|
sig.level |
Significance level (Type I error probability) |
power |
Assurance probability of confidence interval (1 minus Type II error probability) |
prev |
Expected prevalence, if |
NMAX |
Maximum sample size considered. |
Exactly one of the parameters n, delta, sig.level,
and power must be passed as NULL, and that parameter is determined
from the others. Notice that sig.level has a non-NULL default
so NULL must be explicitly passed if you want to compute it.
The computations use the variance of the AUC derived by Hanley and McNeil (1982) and incorporate an additional assurance probability (power) as in the approach of Flahault et al. (2005); see also Kohl and Muench (2024).
Object of class "power.htest", a list of the arguments
(including the computed one) augmented with method and
note elements.
uniroot is used to solve the equations for unknowns, so
you may see errors from it, notably about inability to bracket the
root when invalid arguments are given.
Matthias Kohl [email protected]
A. Flahault, M. Cadilhac, and G. Thomas (2005). Sample size calculation should be performed for design accuracy in diagnostic test studies. Journal of Clinical Epidemiology, 58(8):859-862.
J.A. Hanley, B.J. McNeil (1982). The meaning and use of the area under a receiver operating characteristic (ROC) curve. Radiology, 143(1):29-36.
M. Kohl, F. Muench (2024). Statistik Teil 12: Fallzahlplanung. Die Perfusiologie, 2024(4):124-127.
## compute n ssize.auc.ci(AUC = 0.9, delta = 0.05, power = 0.8) ## compute delta ssize.auc.ci(AUC = 0.9, n = 254, power = 0.8) ## compute power ssize.auc.ci(AUC = 0.9, n = 254, delta = 0.05) ## compute sig.level ssize.auc.ci(AUC = 0.9, n = 254, delta = 0.05, power = 0.8, sig.level = NULL)## compute n ssize.auc.ci(AUC = 0.9, delta = 0.05, power = 0.8) ## compute delta ssize.auc.ci(AUC = 0.9, n = 254, power = 0.8) ## compute power ssize.auc.ci(AUC = 0.9, n = 254, delta = 0.05) ## compute sig.level ssize.auc.ci(AUC = 0.9, n = 254, delta = 0.05, power = 0.8, sig.level = NULL)
Calculate sample size for training set in developing classifiers using high dimensional data. The calculation is based on the probability of correct classification (PCC).
ssize.pcc(gamma, stdFC, prev = 0.5, nrFeatures, sigFeatures = 20, verbose = FALSE)ssize.pcc(gamma, stdFC, prev = 0.5, nrFeatures, sigFeatures = 20, verbose = FALSE)
gamma |
tolerance between PCC(infty) and PCC(n). |
stdFC |
expected standardized fold-change; that is, expected fold-change devided by within class standard deviation. |
prev |
expected prevalence. |
nrFeatures |
number of features (variables) considered. |
sigFeatures |
number of significatn features; default (20) should be sufficient for most if not all cases. |
verbose |
print intermediate results. |
The computations are based the algorithm provided in Section~4.2 of Dobbin and Simon (2007). Prevalence is incorporated by the simple rough approach given in Section~4.4 (ibid.).
The results for prevalence equal to $50%$ are identical to the numbers computed by https://brb.nci.nih.gov/brb/samplesize/samplesize4GE.html. For other prevalences the numbers differ and are larger for our implementation.
Object of class "power.htest", a list of the arguments
(including the computed one) augmented with method and
note elements.
optimize is used to solve equation (4.3) of Dobbin and Simon (2007),
so you may see errors from it.
Matthias Kohl [email protected]
K. Dobbin and R. Simon (2007). Sample size planning for developing classifiers using high-dimensional DNA microarray data. Biostatistics, 8(1):101-117.
K. Dobbin, Y. Zhao, R. Simon (2008). How Large a Training Set is Needed to Develop a Classifier for Microarray Data? Clin Cancer Res., 14(1):108-114.
## see Table 2 of Dobbin et al. (2008) g <- 0.1 fc <- 1.6 ssize.pcc(gamma = g, stdFC = fc, nrFeatures = 22000) ## see Table 3 of Dobbin et al. (2008) g <- 0.05 fc <- 1.1 ssize.pcc(gamma = g, stdFC = fc, nrFeatures = 22000)## see Table 2 of Dobbin et al. (2008) g <- 0.1 fc <- 1.6 ssize.pcc(gamma = g, stdFC = fc, nrFeatures = 22000) ## see Table 3 of Dobbin et al. (2008) g <- 0.05 fc <- 1.1 ssize.pcc(gamma = g, stdFC = fc, nrFeatures = 22000)
Compute the sample size for the two-sided confidence interval of a single proportion.
ssize.prop.ci(prop, width, conf.level = 0.95, method = "wald-cc") ssize.propCI(prop, width, conf.level = 0.95, method = "wald-cc")ssize.prop.ci(prop, width, conf.level = 0.95, method = "wald-cc") ssize.propCI(prop, width, conf.level = 0.95, method = "wald-cc")
prop |
expected proportion |
width |
width of the confidence interval |
conf.level |
confidence level |
method |
method used to compute the confidence interval; see Details. |
The computation is based on the asymptotic formulas provided in Section 2.5.2
of Fleiss et al. (2003). If method = "wald-cc" a continuity correction
is applied.
There are also methods for Jeffreys, Clopper-Pearson, Wilson and the
Agresti-Coull interval; see also binomCI.
Object of class "power.htest", a list of the arguments
(including the computed one) augmented with method and
note elements.
Matthias Kohl [email protected]
J.L. Fleiss, B. Levin and M.C. Paik (2003). Statistical Methods for Rates and Proportions. Wiley Series in Probability and Statistics.
W.W. Piegorsch (2004). Sample sizes for improved binomial confidence intervals. Computational Statistics & Data Analysis, 46, 309-316.
M. Thulin (2014). The cost of using exact confidence intervals for a binomial proportion. Electronic Journal of Statistics, 8(1), 817-840.
ssize.propCI(prop = 0.1, width = 0.1) ssize.propCI(prop = 0.3, width = 0.1) ssize.propCI(prop = 0.3, width = 0.1, method = "wald") ssize.propCI(prop = 0.3, width = 0.1, method = "jeffreys") ssize.propCI(prop = 0.3, width = 0.1, method = "clopper-pearson") ssize.propCI(prop = 0.3, width = 0.1, method = "wilson") ssize.propCI(prop = 0.3, width = 0.1, method = "agresti-coull")ssize.propCI(prop = 0.1, width = 0.1) ssize.propCI(prop = 0.3, width = 0.1) ssize.propCI(prop = 0.3, width = 0.1, method = "wald") ssize.propCI(prop = 0.3, width = 0.1, method = "jeffreys") ssize.propCI(prop = 0.3, width = 0.1, method = "clopper-pearson") ssize.propCI(prop = 0.3, width = 0.1, method = "wilson") ssize.propCI(prop = 0.3, width = 0.1, method = "agresti-coull")
Compute the sample size for reference range studies, or determine parameters for a given sample size; see Jennen-Steinmetz and Wellek (2005).
ssize.reference.range(n = NULL, delta = NULL, ref.prob = 0.95, conf.prob = NULL, alternative = c("two.sided", "one.sided"), method = "parametric", exact = TRUE, tol = .Machine$double.eps^0.5)ssize.reference.range(n = NULL, delta = NULL, ref.prob = 0.95, conf.prob = NULL, alternative = c("two.sided", "one.sided"), method = "parametric", exact = TRUE, tol = .Machine$double.eps^0.5)
n |
number of observations |
delta |
difference between empirical and target coverage of reference range |
ref.prob |
target coverage of reference range |
conf.prob |
confidence probability to acchieve given difference between empirical and target coverage |
alternative |
a character string specifying |
method |
either |
exact |
use exact or approximate method |
tol |
numerical tolerance used in root finding, the default providing (at least) eight significant digits. |
Exactly one of the parameters n, delta, ref.prob
and conf.prob must be passed as NULL, and that parameter is
determined from the others. In case of ref.prob NULL must be explicitly
passed if you want to compute it.
If method "parametric" a normal distribution is assumed for the investigated
quantity.
If method "nonparametric" an arbitrary continuous probability distribution
is assumed.
If exact = TRUE is used, the computations use the exact formulas (5)
and (9) of Jennen-Steinmetz and Wellek (2005).
If exact = FALSE is used, the computations use the approximate formulas (6)
and (10) of Jennen-Steinmetz and Wellek (2005).
Object of class "power.htest", a list of the arguments
(including the computed one) augmented with method and
note elements.
uniroot is used to solve the equations for unknowns, so
you may see errors from it, notably about inability to bracket the
root when invalid arguments are given.
Matthias Kohl [email protected]
C. Jennen-Steinmetz, S. Wellek (2005). A new approach to sample size calculation for reference interval studies. Statistics in Medicine 24:3199-3212.
## see Table 1 in Jennen-Steinmetz and Wellek (2005) ssize.reference.range(delta = 0.03, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = TRUE) ## 135 vs 125 (error in Table 1) ssize.reference.range(delta = 0.03, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.03, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.03, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.025, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = TRUE) ssize.reference.range(delta = 0.025, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.025, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.025, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.02, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = TRUE) ## 314 vs. 305 (error Table 1?) ssize.reference.range(delta = 0.02, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.02, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.02, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.015, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = TRUE) ssize.reference.range(delta = 0.015, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.015, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.015, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.01, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = TRUE) ssize.reference.range(delta = 0.01, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.01, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.01, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.015, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = TRUE) ssize.reference.range(delta = 0.015, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.015, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.015, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.0125, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = TRUE) ssize.reference.range(delta = 0.0125, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.0125, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.0125, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.01, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = TRUE) ssize.reference.range(delta = 0.01, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.01, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.01, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.0075, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = TRUE) ssize.reference.range(delta = 0.0075, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.0075, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.0075, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.005, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = TRUE) ssize.reference.range(delta = 0.005, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.005, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.005, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ## results are equivalent to one-sided reference range with coverage of ## 95 percent instead of 90 percent; for example ssize.reference.range(delta = 0.03, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = TRUE, alternative = "one.sided") ## 135 vs 125 (error in Table 1) ssize.reference.range(delta = 0.03, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = TRUE, alternative = "one.sided") ssize.reference.range(delta = 0.03, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = FALSE, alternative = "one.sided") ssize.reference.range(delta = 0.03, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = FALSE, alternative = "one.sided")## see Table 1 in Jennen-Steinmetz and Wellek (2005) ssize.reference.range(delta = 0.03, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = TRUE) ## 135 vs 125 (error in Table 1) ssize.reference.range(delta = 0.03, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.03, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.03, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.025, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = TRUE) ssize.reference.range(delta = 0.025, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.025, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.025, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.02, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = TRUE) ## 314 vs. 305 (error Table 1?) ssize.reference.range(delta = 0.02, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.02, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.02, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.015, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = TRUE) ssize.reference.range(delta = 0.015, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.015, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.015, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.01, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = TRUE) ssize.reference.range(delta = 0.01, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.01, ref.prob = 0.9, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.01, ref.prob = 0.9, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.015, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = TRUE) ssize.reference.range(delta = 0.015, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.015, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.015, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.0125, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = TRUE) ssize.reference.range(delta = 0.0125, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.0125, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.0125, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.01, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = TRUE) ssize.reference.range(delta = 0.01, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.01, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.01, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.0075, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = TRUE) ssize.reference.range(delta = 0.0075, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.0075, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.0075, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ssize.reference.range(delta = 0.005, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = TRUE) ssize.reference.range(delta = 0.005, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = TRUE) ssize.reference.range(delta = 0.005, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = FALSE) ssize.reference.range(delta = 0.005, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = FALSE) ## results are equivalent to one-sided reference range with coverage of ## 95 percent instead of 90 percent; for example ssize.reference.range(delta = 0.03, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = TRUE, alternative = "one.sided") ## 135 vs 125 (error in Table 1) ssize.reference.range(delta = 0.03, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = TRUE, alternative = "one.sided") ssize.reference.range(delta = 0.03, ref.prob = 0.95, conf.prob = 0.9, method = "parametric", exact = FALSE, alternative = "one.sided") ssize.reference.range(delta = 0.03, ref.prob = 0.95, conf.prob = 0.9, method = "nonparametric", exact = FALSE, alternative = "one.sided")
Produce volcano plot(s) for simulations of power and type-I-error of tests.
## S3 method for class 'sim.power.ttest' volcano(x, alpha = 1, shape = 19, hex = FALSE, bins = 50, ...) ## S3 method for class 'sim.power.wtest' volcano(x, alpha = 1, shape = 19, hex = FALSE, bins = 50, ...)## S3 method for class 'sim.power.ttest' volcano(x, alpha = 1, shape = 19, hex = FALSE, bins = 50, ...) ## S3 method for class 'sim.power.wtest' volcano(x, alpha = 1, shape = 19, hex = FALSE, bins = 50, ...)
x |
object of class |
alpha |
bleding factor (default: no blending. |
shape |
point shape used. |
hex |
logical, should hexagonal binning be used. |
bins |
number of bins used for hexagonal binning. |
... |
further arguments that may be passed through). |
The plot generates a ggplot2 object that is shown.
Missing values are handled by the ggplot2 functions.
Object of class gg and ggplot.
Matthias Kohl [email protected]
Wikipedia contributors, Volcano plot (statistics), Wikipedia, The Free Encyclopedia, https://en.wikipedia.org/w/index.php?title=Volcano_plot_(statistics)&oldid=900217316 (accessed December 25, 2019).
For more sophisticated and flexible volcano plots see for instance: Blighe K, Rana S, Lewis M (2019). EnhancedVolcano: Publication-ready volcano plots with enhanced colouring and labeling. R/Bioconductor package. https://github.com/kevinblighe/EnhancedVolcano.
res1 <- sim.power.t.test(nx = 5, rx = rnorm, rx.H0 = rnorm, ny = 10, ry = function(x) rnorm(x, mean = 3, sd = 3), ry.H0 = function(x) rnorm(x, sd = 3), methods = c("student", "welch")) volcano(res1) ## Not run: ## low number of iterations to reduce computation time res2 <- sim.power.wilcox.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm, iter = 100, conf.int = TRUE) volcano(res2) ## End(Not run)res1 <- sim.power.t.test(nx = 5, rx = rnorm, rx.H0 = rnorm, ny = 10, ry = function(x) rnorm(x, mean = 3, sd = 3), ry.H0 = function(x) rnorm(x, sd = 3), methods = c("student", "welch")) volcano(res1) ## Not run: ## low number of iterations to reduce computation time res2 <- sim.power.wilcox.test(nx = 6, rx = rnorm, rx.H0 = rnorm, ny = 6, ry = function(x) rnorm(x, mean = 2), ry.H0 = rnorm, iter = 100, conf.int = TRUE) volcano(res2) ## End(Not run)