[med-svn] [r-cran-epir] 02/05: New upstream version 0.9-91
Andreas Tille
tille at debian.org
Sat Nov 11 13:56:25 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-epir.
commit bac6ed14f898e70298d183d87ad8c5abd1e638d8
Author: Andreas Tille <tille at debian.org>
Date: Sat Nov 11 14:50:32 2017 +0100
New upstream version 0.9-91
---
DESCRIPTION | 10 +-
MD5 | 22 +--
R/epi.ccsize.R | 498 ++++++++++++++++++++++++++------------------------
R/epi.cohortsize.R | 134 ++++++++------
R/epi.cp.R | 49 ++---
R/epi.meansize.R | 97 ++++++----
R/epi.propsize.R | 153 ++++++++++------
man/epi.ccsize.Rd | 57 +++---
man/epi.cohortsize.Rd | 40 ++--
man/epi.cp.Rd | 12 +-
man/epi.meansize.Rd | 16 +-
man/epi.propsize.Rd | 42 +++--
12 files changed, 632 insertions(+), 498 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 22e2f15..66d7cb3 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
Package: epiR
-Version: 0.9-87
-Date: 2017-06-30
+Version: 0.9-91
+Date: 2017-10-26
Title: Tools for the Analysis of Epidemiological Data
-Author: Mark Stevenson <mark.stevenson1 at unimelb.edu.au> with contributions from Telmo Nunes, Cord Heuer, Jonathon Marshall, Javier Sanchez, Ron Thornton, Jeno Reiczigel, Jim Robison-Cox, Paola Sebastiani, Peter Solymos, Kazuki Yoshida, Geoff Jones, Sarah Pirikahu, Simon Firestone and Ryan Kyle.
+Author: Mark Stevenson <mark.stevenson1 at unimelb.edu.au> with contributions from Telmo Nunes, Cord Heuer, Jonathon Marshall, Javier Sanchez, Ron Thornton, Jeno Reiczigel, Jim Robison-Cox, Paola Sebastiani, Peter Solymos, Kazuki Yoshida, Geoff Jones, Sarah Pirikahu, Simon Firestone, Ryan Kyle and Johann Popp.
Maintainer: Mark Stevenson <mark.stevenson1 at unimelb.edu.au>
Description: Tools for the analysis of epidemiological data. Contains functions for directly and indirectly adjusting measures of disease frequency, quantifying measures of association on the basis of single or multiple strata of count data presented in a contingency table, and computing confidence intervals around incidence risk and incidence rate estimates. Miscellaneous functions for use in meta-analysis, diagnostic test interpretation, and sample size calculations.
Depends: R (>= 3.0.0), survival
@@ -11,6 +11,6 @@ Suggests: MASS (>= 3.1-20)
License: GPL (>= 2)
URL: http://fvas.unimelb.edu.au/veam
NeedsCompilation: no
-Packaged: 2017-06-30 07:54:55 UTC; Mark
+Packaged: 2017-10-26 04:12:40 UTC; marks1
Repository: CRAN
-Date/Publication: 2017-06-30 22:15:10 UTC
+Date/Publication: 2017-10-27 13:07:08 UTC
diff --git a/MD5 b/MD5
index 8cca1cd..cd3df70 100644
--- a/MD5
+++ b/MD5
@@ -1,4 +1,4 @@
-3153e775b66b924b25a378df61848059 *DESCRIPTION
+3fc0a923a43866a0192811642e0aeb94 *DESCRIPTION
7fef6704a714f9fc47a477f6b751aaab *NAMESPACE
36914e5e5904e60bd5bd33af6cb05a45 *R/epi.2by2.R
90797cfa53b823b298af4c32eb3435bf *R/epi.RtoBUGS.R
@@ -7,14 +7,14 @@ a9a4deaaec9cd9435529455159385225 *R/epi.about.R
88fc1fad29654a3b8cff28d58bd332d6 *R/epi.betabuster.R
89e1ade6b7a25729a3e7db833abd2739 *R/epi.bohning.R
b8e64c719988731209aa1fbc7fdcc5aa *R/epi.ccc.R
-fdcbfe2ac48ef09a6028cb3ab86293e2 *R/epi.ccsize.R
+dad0d24004b6f34d5e4563c0452e1835 *R/epi.ccsize.R
52927695dc6e0d0002c6ab70a12ac926 *R/epi.cluster1size.R
aed994d4e585e092089481c32bc12835 *R/epi.cluster2size.R
c3ff42af84bb01e927867d54b4ac4334 *R/epi.clustersize.R
-e52bded5499df4fecafc409235593c79 *R/epi.cohortsize.R
+e88d3240c668b29504f247bf2c824ebd *R/epi.cohortsize.R
8e08de8792494503161a98b7e38d9ef9 *R/epi.conf.R
d535d0aaae250f95c38e00ccbba0b7ca *R/epi.convgrid.R
-4eb1a96f7427ea81382d37754d717f49 *R/epi.cp.R
+6fcf38898074b7da8d20c1b83a76f730 *R/epi.cp.R
0bce821540346b1e58e640d7334c043b *R/epi.cpresids.R
e3d7b38e88eb943c976edb0f38350174 *R/epi.descriptives.R
7ca563951950b4476e001ce3dba64043 *R/epi.detectsize.R
@@ -33,7 +33,7 @@ e935ed017215d39b2225f61b23de3bfd *R/epi.herdtest.R
c6d1bbc1733a57b2ed6973bbc9a2689a *R/epi.iv.R
0231853e687a5153f9e9bb186bfd5fed *R/epi.kappa.R
dd5a6b971d1652170e84fa8a709215a9 *R/epi.ltd.R
-20acf4c5b5334be8d1d66a8c83bcdace *R/epi.meansize.R
+52fd3efc9736ce101747cf779233ae84 *R/epi.meansize.R
27df85ca4570d0077ec68ff552739d13 *R/epi.mh.R
064a3d44cfc040bc957d37d9dd6ac64b *R/epi.nomogram.R
fda5a6132e84888d9ae9ec0e3054aee6 *R/epi.noninfb.R
@@ -44,7 +44,7 @@ b1c0bdf50b56a4b9df09ac89f2de0b98 *R/epi.pooled.R
3a35fac64c1bcda32c3dcabe7b7d14b2 *R/epi.popsize.R
9984eb80ee76e03e982e401c2c3a2490 *R/epi.prcc.R
3f91e8c64f3dda591ed7f2d98dc5d39c *R/epi.prev.R
-aeafef81811418812f4d895466a1e480 *R/epi.propsize.R
+fbe3c8b6aebc5541f474d11443e427d4 *R/epi.propsize.R
b149e431ad851995742d356767b1a1ab *R/epi.simplesize.R
4e719a2b9c4ba572f0ee3c3d70d54190 *R/epi.smd.R
32b646ecabba4110b2566b44d1879d49 *R/epi.stratasize.R
@@ -64,14 +64,14 @@ c0207f2c0a218b3c799cfd6e5b0df220 *man/epi.about.Rd
a3f56b6074f71a9f6a5356d298d20a01 *man/epi.betabuster.Rd
7e06c95f246ea6a933847b73c2c0288a *man/epi.bohning.Rd
3ba388722aeff534ebc8661ac89051f0 *man/epi.ccc.Rd
-7142895a2342a283b2c732428536f340 *man/epi.ccsize.Rd
+1606caafb53acc3ec969752708f4702c *man/epi.ccsize.Rd
2758ba80c5bde26752e8b5de29e03161 *man/epi.cluster1size.Rd
b86f88dbb6fdfca830a2afbc9a764f30 *man/epi.cluster2size.Rd
572fe905ca4759a69445fda1de10d48e *man/epi.clustersize.Rd
-a83f74d10ec9b3bf0f9cb140d8f6fcba *man/epi.cohortsize.Rd
+f702487d058e34fa75a4173710ab44c4 *man/epi.cohortsize.Rd
ddc1befc74ee119935f03d93a3301e04 *man/epi.conf.Rd
07e702fd65b9cc6710d27c133b94f89f *man/epi.convgrid.Rd
-b8c620a7bf19cdff691279488646ab37 *man/epi.cp.Rd
+3650dbdffa8726d7eb0859f10c4b7e57 *man/epi.cp.Rd
11dc595f8d992cd9b43c0c1eb83879e7 *man/epi.cpresids.Rd
4227b839d52e6052284da81817ec836a *man/epi.descriptives.Rd
9d9251304690d373ee190fb27626904c *man/epi.detectsize.Rd
@@ -92,7 +92,7 @@ f5eb9a339ceb32787522e1b10cfd2693 *man/epi.interaction.Rd
ede50c148a4b3ab71a74733a0ac953db *man/epi.iv.Rd
866584c8d65a7cf56e3e3a85b0250096 *man/epi.kappa.Rd
85eb0f47054239387557a52af03a674f *man/epi.ltd.Rd
-28437d1309ead1e42782b863a71b9b6f *man/epi.meansize.Rd
+9078d0569ac7fd0f05b04103c7a83088 *man/epi.meansize.Rd
867e8fc4d3f310a9f044d02309d6a5b4 *man/epi.mh.Rd
22c11e69793964a066401af0d99903dd *man/epi.nomogram.Rd
f10943cdc98bd05521d53f5218844743 *man/epi.noninfb.Rd
@@ -103,7 +103,7 @@ cc7f75a85351f5542916421fa4b7e648 *man/epi.offset.Rd
94d5f2ced08ba7556c21e6360d5f54b8 *man/epi.popsize.Rd
eb4fa4e5bf5996b7f46b9bb9e67f57e0 *man/epi.prcc.Rd
1f01142f1724a5c1c88162aa96b8befe *man/epi.prev.Rd
-246145a7c8c19f7ff7298ee636627f61 *man/epi.propsize.Rd
+0cc4da187b5d7bf97eda273ed0d622b3 *man/epi.propsize.Rd
8e0668a9c1a7376e32521204873b8659 *man/epi.simplesize.Rd
b6389f885f3fdf71f58c115334f230d4 *man/epi.smd.Rd
96e8daea0c1e36beb684ebbd85b7f2c6 *man/epi.stratasize.Rd
diff --git a/R/epi.ccsize.R b/R/epi.ccsize.R
index 974d31c..c901880 100644
--- a/R/epi.ccsize.R
+++ b/R/epi.ccsize.R
@@ -1,258 +1,274 @@
"epi.ccsize" <- function(OR, p0, n, power, r = 1, rho = 0, design = 1, sided.test = 2, conf.level = 0.95, method = "unmatched", fleiss = FALSE) {
- # p0: proportion of controls exposed.
- # rho: correlation between case and control exposures in matched pairs (defaults to 0).
-
- # https://www2.ccrb.cuhk.edu.hk/stat/epistudies/cc2.htm
+ # p0: proportion of controls exposed.
+ # rho: correlation between case and control exposures in matched pairs (defaults to 0).
- alpha.new <- (1 - conf.level) / sided.test
- z.alpha <- qnorm(1 - alpha.new, mean = 0, sd = 1)
-
-# ------------------------------------------------------------------------------------------------------
- # Unmatched case-control - sample size:
- if(method == "unmatched" & !is.na(OR) & is.na(n) & !is.na(power)){
+ # https://www2.ccrb.cuhk.edu.hk/stat/epistudies/cc2.htm
- beta <- 1 - power
- z.beta <- qnorm(p = beta, lower.tail = FALSE)
-
- # Sample size without Fleiss correction:
- q0 <- 1 - p0
-
- p1 <- (p0 * OR) / (1 + p0 * (OR - 1))
- q1 <- 1 - p1
-
- pbar <- (p1 + (r * p0)) / (r + 1)
- qbar <- 1 - pbar
+ alpha.new <- (1 - conf.level) / sided.test
+ z.alpha <- qnorm(1 - alpha.new, mean = 0, sd = 1)
- Np1 <- z.alpha * sqrt((r + 1) * pbar * qbar)
- Np2 <- z.beta * sqrt((r * p1 * q1) + (p0 * q0))
- Dp1 <- r * (p1 - p0)^2
- n. <- (Np1 + Np2)^2 / Dp1
-
- if(fleiss == TRUE){
- d <- 1 + ((2 * (r + 1)) / (n. * r * abs(p0 - p1)))
- n1 <- (n. / 4) * (1 + sqrt(d))^2
+ if (!is.na(OR) & !is.na(n) & !is.na(power)){
+ stop("Error: at least one of OR, n and power must be NA.")
}
- if(fleiss == FALSE){
- n1 <- n.
+ # ------------------------------------------------------------------------------------------------------
+ # Unmatched case-control - sample size:
+ if(method == "unmatched" & !is.na(OR) & is.na(n) & !is.na(power)){
+
+ beta <- 1 - power
+ z.beta <- qnorm(p = beta, lower.tail = FALSE)
+
+ # Sample size without Fleiss correction - from Woodward's spreadsheet:
+ Pc <- (p0 / (r + 1)) * ((r * OR / (1 + (OR - 1) * p0)) + 1)
+ t1 <- (r + 1) * (1 + (OR - 1) * p0)^2
+ t2 <- r * p0^2 * (p0 - 1)^2 * (OR - 1)^2
+ t3 <- z.alpha * sqrt((r + 1) * Pc * (1 - Pc))
+ t4 <- OR * p0 * (1 - p0)
+ t5 <- 1 + (OR - 1) * p0
+ t6 <- t4 / (t5^2)
+ t7 <- r * p0 * (1 - p0)
+ t8 <- z.beta * sqrt(t6 + t7)
+ n. <- (t1 / t2) * (t3 + t8)^2
+
+ # q0 <- 1 - p0
+ # p1 <- (p0 * OR) / (1 + p0 * (OR - 1))
+ # q1 <- 1 - p1
+
+ # pbar <- (p1 + (r * p0)) / (r + 1)
+ # qbar <- 1 - pbar
+
+ # Np1 <- z.alpha * sqrt((r + 1) * pbar * qbar)
+ # Np2 <- z.beta * sqrt((r * p1 * q1) + (p0 * q0))
+ # Dp1 <- r * (p1 - p0)^2
+ # n. <- (Np1 + Np2)^2 / Dp1
+
+ if(fleiss == TRUE){
+ d <- 1 + ((2 * (r + 1)) / (n. * r * abs(p0 - p1)))
+ n <- (n. / 4) * (1 + sqrt(d))^2
+ }
+
+ if(fleiss == FALSE){
+ n <- n.
+ }
+
+ # Account for the design effect:
+ n1 <- n / (1 + r)
+ n1 <- ceiling(n1 * design)
+ n2 <- ceiling(r * n1)
+
+ rval <- list(n.total = n1 + n2, n.case = n1, n.control = n2, power = power, OR = OR)
}
- # Account for the design effect:
- n1 <- ceiling(n1 * design)
- n2 <- ceiling(r * n1)
+ # Unmatched case-control - power:
+ else
+ if(method == "unmatched" & !is.na(OR) & !is.na(n) & is.na(power)){
+
+ # From Woodward's spreadsheet:
+ Pc <- (p0 / (r + 1)) * ((r * OR / (1 + (OR - 1) * p0)) + 1)
+ M <- abs((OR - 1) * (p0 - 1) / (1 + (OR - 1) * p0))
+ termn1 <- M * p0 * sqrt(n * r) / sqrt(r + 1)
+ termn2 <- z.alpha * sqrt((r + 1) * Pc *(1 - Pc))
+ termd1 <- OR * p0 * (1 - p0) / (1 + (OR - 1) * p0)^2
+ termd2 <- r * p0 * (1 - p0)
+ z.beta <- (termn1 - termn2) / sqrt(termd1 + termd2)
- rval <- list(n.total = n1 + n2, n.case = n1, n.control = n2, power = power, OR = OR)
+ # q0 <- 1 - p0
+ # p1 <- (p0 * OR) / (1 + p0 * (OR - 1))
+ # q1 <- 1 - p1
+
+ # pbar <- (p1 + (r * p0)) / (r + 1)
+ # qbar <- 1 - pbar
+
+ n1 <- n / (1 + r)
+
+ # z.beta <- sqrt((n1 * r * (p1 - p0)^2) / (pbar * qbar * (r + 1))) - z.alpha
+ power <- pnorm(z.beta, mean = 0, sd = 1)
+
+ # Account for the design effect:
+ n1 <- ceiling(n1 * design)
+ n2 <- n - n1
+ # n2 <- ceiling(r * n1)
+
+ rval <- list(n.total = n1 + n2, n.case = n1, n.control = n2, power = power, OR = OR)
+ }
+
+ # Unmatched case-control - effect:
+ else
+ if(method == "unmatched" & is.na(OR) & !is.na(n) & !is.na(power)){
+
+ n1 <- n / (r + 1)
+ n1 <- ceiling(n1 * design)
+ n2 <- n - n1
+ # n2 <- ceiling(r * n1)
+
+ Pfun <- function(OR, power, p0, r, n, design, z.alpha){
+ q0 <- 1 - p0
+ p1 <- (p0 * OR) / (1 + p0 * (OR - 1))
+ q1 <- 1 - p1
+
+ pbar <- (p1 + (r * p0)) / (r + 1)
+ qbar <- 1 - pbar
+
+ # Account for the design effect:
+ n1 <- n / (1 + r)
+ n1 <- ceiling(n1 * design)
+ n2 <- n - n1
+ # n2 <- ceiling(r * n1)
+
+ z.beta <- sqrt((n1 * r * (p1 - p0)^2) / (pbar * qbar * (r + 1))) - z.alpha
+
+ # Take the calculated value of the power and subtract the power entered by the user:
+ pnorm(z.beta, mean = 0, sd = 1) - power
+ }
+
+ # Find the value of OR that matches the power entered by the user:
+ OR.up <- uniroot(Pfun, power = power, p0 = p0, r = r, n = n, design = design, z.alpha = z.alpha, interval = c(1,1E06))$root
+ OR.lo <- uniroot(Pfun, power = power, p0 = p0, r = r, n = n, design = design, z.alpha = z.alpha, interval = c(0.0001,1))$root
+
+ # x = seq(from = 0.01, to = 100, by = 0.01)
+ # y = Pfun(x, power = 0.8, p0 = 0.15, r = 1, n = 150, design = 1, z.alpha = 1.96)
+ # windows(); plot(x, y, xlim = c(0,5)); abline(h = 0, lty = 2)
+ # Two possible values for OR meet the conditions of Pfun. So hence we set the lower bound of the search interval to 1.
+
+ rval <- list(n.total = n1 + n2, n.case = n1, n.control = n2, power = power, OR = c(OR.lo, OR.up))
+ }
+
+
+ # ------------------------------------------------------------------------------------------------------
+ # Matched case-control - sample size:
+
+ if(method == "matched" & !is.na(OR) & is.na(n) & !is.na(power)){
+
+ beta <- 1 - power
+ z.beta <- qnorm(p = beta, lower.tail = FALSE)
+
+ odds0 = p0 / (1 - p0)
+ odds1 = odds0 * OR
+ p1 = odds1 / (1 + odds1)
+ delta.p = p1 - p0
+
+ psi <- OR
+ pq <- p1 * (1 - p1) * p0 * (1 - p0)
+ p0.p <- (p1 * p0 + rho * sqrt(pq)) / p1
+ p0.n <- (p0 * (1 - p1) - rho * sqrt(pq)) / (1 - p1)
+
+ tm <- ee.psi <- ee.one <- nu.psi <- nu.one <- rep(NA, r)
+ for(m in 1:r){
+ tm[m] <- p1 * choose(r, m - 1) * ((p0.p)^(m - 1)) * ((1 - p0.p)^(r - m + 1)) + (1 - p1) * choose(r,m) * ((p0.n)^m) * ((1 - p0.n))^(r - m)
+ ee.psi[m] <- m * tm[m] * psi / (m * psi + r - m + 1)
+ ee.one[m] <- m * tm[m] / (r + 1)
+ nu.psi[m] <- m * tm[m] * psi * (r - m + 1) / ((m * psi + r - m + 1)^2)
+ nu.one[m] <- m * tm[m] * (r - m + 1) / ((r + 1)^2)
+ }
+
+ ee.psi <- sum(ee.psi)
+ ee.one <- sum(ee.one)
+ nu.psi <- sum(nu.psi)
+ nu.one <- sum(nu.one)
+
+ n1 <- ceiling(((z.beta * sqrt(nu.psi) + z.alpha * sqrt(nu.one))^2) / ((ee.one - ee.psi)^2))
+ n2 <- ceiling(r * n1)
+
+ rval <- list(n.total = n1 + n2, n.case = n1, n.control = n2, power = power, OR = OR)
}
-
- # Unmatched case-control - power:
- else
- if(method == "unmatched" & !is.na(OR) & !is.na(n) & is.na(power)){
- q0 <- 1 - p0
-
- p1 <- (p0 * OR) / (1 + p0 * (OR - 1))
- q1 <- 1 - p1
-
- pbar <- (p1 + (r * p0)) / (r + 1)
- qbar <- 1 - pbar
-
- n1 <- n / (r + 1) * r
- z.beta <- sqrt((n1 * r * (p1 - p0)^2) / (pbar * qbar * (r + 1))) - z.alpha
- power <- pnorm(z.beta, mean = 0, sd = 1)
-
- # Account for the design effect:
- n1 <- ceiling(n1 * design)
- n2 <- ceiling(r * n1)
-
- rval <- list(n.total = n1 + n2, n.case = n1, n.control = n2, power = power, OR = OR)
-}
-
- # Unmatched case-control - effect:
- else
- if(method == "unmatched" & is.na(OR) & !is.na(n) & !is.na(power)){
-
- n1 <- n / (r + 1) * r
- n1 <- ceiling(n1 * design)
- n2 <- ceiling(r * n1)
-
- Pfun <- function(OR, power, p0, r, n, design, z.alpha){
- q0 <- 1 - p0
- p1 <- (p0 * OR) / (1 + p0 * (OR - 1))
- q1 <- 1 - p1
-
- pbar <- (p1 + (r * p0)) / (r + 1)
- qbar <- 1 - pbar
-
- # Account for the design effect:
- n1 <- n / (r + 1) * r
- n1 <- ceiling(n1 * design)
- n2 <- ceiling(r * n1)
-
- z.beta <- sqrt((n1 * r * (p1 - p0)^2) / (pbar * qbar * (r + 1))) - z.alpha
-
- # Take the calculated value of the power and subtract the power entered by the user:
- pnorm(z.beta, mean = 0, sd = 1) - power
- }
-
- # Find the value of OR that matches the power entered by the user:
- OR.up <- uniroot(Pfun, power = power, p0 = p0, r = r, n = n, design = design, z.alpha = z.alpha, interval = c(1,1E06))$root
- OR.lo <- uniroot(Pfun, power = power, p0 = p0, r = r, n = n, design = design, z.alpha = z.alpha, interval = c(0.0001,1))$root
-
- # x = seq(from = 0.01, to = 100, by = 0.01)
- # y = Pfun(x, power = 0.8, p0 = 0.15, r = 1, n = 150, design = 1, z.alpha = 1.96)
- # windows(); plot(x, y, xlim = c(0,5)); abline(h = 0, lty = 2)
- # Two possible values for OR meet the conditions of Pfun. So hence we set the lower bound of the search interval to 1.
+ # Matched case-control - power:
+ else
+ if(method == "matched" & !is.na(OR) & !is.na(n) & is.na(power)){
- rval <- list(n.total = n1 + n2, n.case = n1, n.control = n2, power = power, OR = c(OR.lo, OR.up))
- }
-
-
-# ------------------------------------------------------------------------------------------------------
- # Matched case-control - sample size:
+ odds0 = p0 / (1 - p0)
+ odds1 = odds0 * OR
+ p1 = odds1 / (1 + odds1)
+ delta.p = p1 - p0
+
+ psi <- OR
+ beta <- 1 - power
+ z.beta <- qnorm(p = beta,lower.tail = FALSE)
+ pq <- p1 * (1 - p1) * p0 * (1 - p0)
+ p0.p <- (p1 * p0 + rho * sqrt(pq)) / p1
+ p0.n <- (p0 * (1 - p1) - rho * sqrt(pq)) / (1 - p1)
+
+ tm <- ee.psi <- ee.one <- nu.psi <- nu.one <- rep(NA, r)
+ for(m in 1:r){
+ tm[m] <- p1 * choose(r, m - 1) * ((p0.p)^(m - 1)) * ((1 - p0.p)^(r - m + 1)) + (1 - p1) * choose(r,m) * ((p0.n)^m) * ((1 - p0.n))^(r - m)
+ ee.psi[m] <- m * tm[m] * psi / (m * psi + r - m + 1)
+ ee.one[m] <- m * tm[m] / (r + 1)
+ nu.psi[m] <- m * tm[m] * psi * (r - m + 1) / ((m * psi + r - m + 1)^2)
+ nu.one[m] <- m * tm[m] * (r - m + 1) / ((r + 1)^2)
+ }
+
+ ee.psi <- sum(ee.psi)
+ ee.one <- sum(ee.one)
+ nu.psi <- sum(nu.psi)
+ nu.one <- sum(nu.one)
+
+ # Account for the design effect:
+ n1 <- n / (1 + r)
+ n1 <- ceiling(n1 * design)
+ n2 <- n - n1
+ # n2 <- ceiling(r * n1)
- if(method == "matched" & !is.na(OR) & is.na(n) & !is.na(power)){
-
- beta <- 1 - power
- z.beta <- qnorm(p = beta, lower.tail = FALSE)
-
- odds0 = p0 / (1 - p0)
- odds1 = odds0 * OR
- p1 = odds1 / (1 + odds1)
- delta.p = p1 - p0
-
- psi <- OR
- pq <- p1 * (1 - p1) * p0 * (1 - p0)
- p0.p <- (p1 * p0 + rho * sqrt(pq)) / p1
- p0.n <- (p0 * (1 - p1) - rho * sqrt(pq)) / (1 - p1)
-
- tm <- ee.psi <- ee.one <- nu.psi <- nu.one <- rep(NA, r)
- for(m in 1:r){
- tm[m] <- p1 * choose(r, m - 1) * ((p0.p)^(m - 1)) * ((1 - p0.p)^(r - m + 1)) + (1 - p1) * choose(r,m) * ((p0.n)^m) * ((1 - p0.n))^(r - m)
- ee.psi[m] <- m * tm[m] * psi / (m * psi + r - m + 1)
- ee.one[m] <- m * tm[m] / (r + 1)
- nu.psi[m] <- m * tm[m] * psi * (r - m + 1) / ((m * psi + r - m + 1)^2)
- nu.one[m] <- m * tm[m] * (r - m + 1) / ((r + 1)^2)
- }
-
- ee.psi <- sum(ee.psi)
- ee.one <- sum(ee.one)
- nu.psi <- sum(nu.psi)
- nu.one <- sum(nu.one)
-
- n.treat <- ceiling(((z.beta * sqrt(nu.psi) + z.alpha * sqrt(nu.one))^2) / ((ee.one - ee.psi)^2))
-
- # Account for the design effect:
- n.treat <- n.treat * design
- n.control <- n.treat * r
-
- n.crude <- ceiling(n.treat + n.control)
- n.treat <- ceiling(n.treat)
- n.control <- ceiling(n.control)
- n.total <- n.treat + n.control
-
- rval <- list(n.total = n.total, n.case = n.treat, n.control = n.control, power = power, OR = OR)
- }
-
- # Matched case-control - power:
- else
- if(method == "matched" & !is.na(OR) & !is.na(n) & is.na(power)){
-
- odds0 = p0 / (1 - p0)
- odds1 = odds0 * OR
- p1 = odds1 / (1 + odds1)
- delta.p = p1 - p0
-
- psi <- OR
- beta <- 1 - power
- z.beta <- qnorm(p = beta,lower.tail = FALSE)
- pq <- p1 * (1 - p1) * p0 * (1 - p0)
- p0.p <- (p1 * p0 + rho * sqrt(pq)) / p1
- p0.n <- (p0 * (1 - p1) - rho * sqrt(pq)) / (1 - p1)
-
- tm <- ee.psi <- ee.one <- nu.psi <- nu.one <- rep(NA, r)
- for(m in 1:r){
- tm[m] <- p1 * choose(r, m - 1) * ((p0.p)^(m - 1)) * ((1 - p0.p)^(r - m + 1)) + (1 - p1) * choose(r,m) * ((p0.n)^m) * ((1 - p0.n))^(r - m)
- ee.psi[m] <- m * tm[m] * psi / (m * psi + r - m + 1)
- ee.one[m] <- m * tm[m] / (r + 1)
- nu.psi[m] <- m * tm[m] * psi * (r - m + 1) / ((m * psi + r - m + 1)^2)
- nu.one[m] <- m * tm[m] * (r - m + 1) / ((r + 1)^2)
- }
-
- ee.psi <- sum(ee.psi)
- ee.one <- sum(ee.one)
- nu.psi <- sum(nu.psi)
- nu.one <- sum(nu.one)
-
- # Account for the design effect:
- n <- n / design
- n.treat <- ceiling(n / (r + 1)) * r
- n.control <- ceiling(n / (r + 1)) * 1
-
- z.beta <- (sqrt(n.treat * ((ee.one - ee.psi)^2)) - z.alpha * sqrt(nu.one)) / sqrt(nu.psi)
- power <- 1 - pnorm(q = z.beta, lower.tail = FALSE)
-
- rval <- list(n.total = n.treat + n.control, n.case = n.treat, n.control = n.control, power = power, OR = OR)
- }
-
- # Matched case-control - effect:
- else
- if(method == "matched" & is.na(OR) & !is.na(n) & !is.na(power)){
-
- n1 <- n / (r + 1) * r
- n1 <- ceiling(n1 * design)
- n2 <- ceiling(r * n1)
-
- Pfun <- function(OR, power, p0, r, n, design, z.alpha){
- odds0 = p0 / (1 - p0)
- odds1 = odds0 * OR
- p1 = odds1 / (1 + odds1)
- delta.p = p1 - p0
-
- psi <- OR
- beta <- 1 - power
- z.beta <- qnorm(p = beta, lower.tail = FALSE)
- pq <- p1 * (1 - p1) * p0 * (1 - p0)
- p0.p <- (p1 * p0 + rho * sqrt(pq)) / p1
- p0.n <- (p0 * (1 - p1) - rho * sqrt(pq)) / (1 - p1)
-
- tm <- ee.psi <- ee.one <- nu.psi <- nu.one <- rep(NA, r)
- for(m in 1:r){
- tm[m] <- p1 * choose(r, m - 1) * ((p0.p)^(m - 1)) * ((1 - p0.p)^(r - m + 1)) + (1 - p1) * choose(r,m) * ((p0.n)^m) * ((1 - p0.n))^(r - m)
- ee.psi[m] <- m * tm[m] * psi / (m * psi + r - m + 1)
- ee.one[m] <- m * tm[m] / (r + 1)
- nu.psi[m] <- m * tm[m] * psi * (r - m + 1) / ((m * psi + r - m + 1)^2)
- nu.one[m] <- m * tm[m] * (r - m + 1) / ((r + 1)^2)
- }
-
- ee.psi <- sum(ee.psi)
- ee.one <- sum(ee.one)
- nu.psi <- sum(nu.psi)
- nu.one <- sum(nu.one)
-
- # Account for the design effect:
- n <- n / design
- n1 <- ceiling(n / (r + 1)) * r
- n2 <- ceiling(n / (r + 1)) * 1
-
- z.beta <- (sqrt(n1 * ((ee.one - ee.psi)^2)) - z.alpha * sqrt(nu.one)) / sqrt(nu.psi)
-
- # Take the calculated value of the power and subtract the power entered by the user:
- pnorm(z.beta, mean = 0, sd = 1) - power
- }
-
- # Find the value of OR that matches the power entered by the user:
- OR.up <- uniroot(Pfun, power = power, p0 = p0, r = r, n = n, design = design, z.alpha = z.alpha, interval = c(1,1E06))$root
- OR.lo <- uniroot(Pfun, power = power, p0 = p0, r = r, n = n, design = design, z.alpha = z.alpha, interval = c(0.0001,1))$root
-
- # x = seq(from = 0.01, to = 100, by = 0.01)
- # y = Pfun(x, power = 0.8, p0 = 0.15, r = 1, n = 150, design = 1, z.alpha = 1.96)
- # windows(); plot(x, y, xlim = c(0,5)); abline(h = 0, lty = 2)
- # Two possible values for OR meet the conditions of Pfun. So hence we set the lower bound of the search interval to 1.
+ z.beta <- (sqrt(n1 * ((ee.one - ee.psi)^2)) - z.alpha * sqrt(nu.one)) / sqrt(nu.psi)
+ power <- 1 - pnorm(q = z.beta, lower.tail = FALSE)
+
+ rval <- list(n.total = n1 + n2, n.case = n1, n.control = n2, power = power, OR = OR)
+ }
+
+ # Matched case-control - effect:
+ else
+ if(method == "matched" & is.na(OR) & !is.na(n) & !is.na(power)){
- rval <- list(n.total = n1 + n2, n.case = n1, n.control = n2, power = power, OR = c(OR.lo, OR.up))
- }
+ n1 <- n / (1 + r)
+ n1 <- ceiling(n1 * design)
+ n2 <- n - n1
+ # n2 <- ceiling(r * n1)
+
+ Pfun <- function(OR, power, p0, r, n, design, z.alpha){
+ odds0 = p0 / (1 - p0)
+ odds1 = odds0 * OR
+ p1 = odds1 / (1 + odds1)
+ delta.p = p1 - p0
+
+ psi <- OR
+ beta <- 1 - power
+ z.beta <- qnorm(p = beta, lower.tail = FALSE)
+ pq <- p1 * (1 - p1) * p0 * (1 - p0)
+ p0.p <- (p1 * p0 + rho * sqrt(pq)) / p1
+ p0.n <- (p0 * (1 - p1) - rho * sqrt(pq)) / (1 - p1)
+
+ tm <- ee.psi <- ee.one <- nu.psi <- nu.one <- rep(NA, r)
+ for(m in 1:r){
+ tm[m] <- p1 * choose(r, m - 1) * ((p0.p)^(m - 1)) * ((1 - p0.p)^(r - m + 1)) + (1 - p1) * choose(r,m) * ((p0.n)^m) * ((1 - p0.n))^(r - m)
+ ee.psi[m] <- m * tm[m] * psi / (m * psi + r - m + 1)
+ ee.one[m] <- m * tm[m] / (r + 1)
+ nu.psi[m] <- m * tm[m] * psi * (r - m + 1) / ((m * psi + r - m + 1)^2)
+ nu.one[m] <- m * tm[m] * (r - m + 1) / ((r + 1)^2)
+ }
+
+ ee.psi <- sum(ee.psi)
+ ee.one <- sum(ee.one)
+ nu.psi <- sum(nu.psi)
+ nu.one <- sum(nu.one)
- # ------------------------------------------------------------------------------------------------------
- rval
+ z.beta <- (sqrt(n1 * ((ee.one - ee.psi)^2)) - z.alpha * sqrt(nu.one)) / sqrt(nu.psi)
+
+ # Take the calculated value of the power and subtract the power entered by the user:
+ pnorm(z.beta, mean = 0, sd = 1) - power
+ }
+
+ # Find the value of OR that matches the power entered by the user:
+ OR.up <- uniroot(Pfun, power = power, p0 = p0, r = r, n = n, design = design, z.alpha = z.alpha, interval = c(1,1E06))$root
+ OR.lo <- uniroot(Pfun, power = power, p0 = p0, r = r, n = n, design = design, z.alpha = z.alpha, interval = c(0.0001,1))$root
+
+ # x = seq(from = 0.01, to = 100, by = 0.01)
+ # y = Pfun(x, power = 0.8, p0 = 0.15, r = 1, n = 150, design = 1, z.alpha = 1.96)
+ # windows(); plot(x, y, xlim = c(0,5)); abline(h = 0, lty = 2)
+ # Two possible values for OR meet the conditions of Pfun. So hence we set the lower bound of the search interval to 1.
+
+ rval <- list(n.total = n1 + n2, n.case = n1, n.control = n2, power = power, OR = c(OR.lo, OR.up))
+ }
+
+ # ------------------------------------------------------------------------------------------------------
+ rval
}
diff --git a/R/epi.cohortsize.R b/R/epi.cohortsize.R
index a5adaa8..4490048 100644
--- a/R/epi.cohortsize.R
+++ b/R/epi.cohortsize.R
@@ -1,68 +1,84 @@
"epi.cohortsize" <- function(exposed, unexposed, n, power, r = 1, design = 1, sided.test = 2, conf.level = 0.95) {
- alpha.new <- (1 - conf.level) / sided.test
- z.alpha <- qnorm(1 - alpha.new, mean = 0, sd = 1)
-
- if(!is.na(exposed) & !is.na(unexposed) & is.na(n) & !is.na(power)){
- # Sample size estimate. From Woodward p 405:
- z.beta <- qnorm(power, mean = 0, sd = 1)
- lambda <- exposed / unexposed
- pi <- unexposed
- pc <- (pi * ((r * lambda) + 1)) / (r + 1)
- p1 <- (r + 1) / (r * (lambda - 1)^2 * pi^2)
- p2 <- z.alpha * sqrt((r + 1) * pc * (1 - pc))
- p3 <- z.beta * sqrt((lambda * pi * (1 - (lambda * pi))) + (r * pi * (1 - pi)))
- n <- p1 * (p2 + p3)^2
+ alpha.new <- (1 - conf.level) / sided.test
+ z.alpha <- qnorm(1 - alpha.new, mean = 0, sd = 1)
- # Account for the design effect:
- n <- n * design
+ if (!is.na(exposed) & !is.na(n) & !is.na(power)){
+ stop("Error: at least one of exposed, n and power must be NA.")
+ }
- n.crude <- ceiling(n)
- n.exposed <- ceiling(n / (r + 1)) * r
- n.unexposed <- ceiling(n / (r + 1)) * 1
- n.total <- n.exposed + n.unexposed
- rval <- list(n.crude = n.crude, n.total = n.total, n.exposed = n.exposed, n.unexposed = n.unexposed)
+ # Sample size:
+ if(!is.na(exposed) & !is.na(unexposed) & is.na(n) & !is.na(power)){
+ # Sample size estimate. From Woodward p 405:
+ z.beta <- qnorm(power, mean = 0, sd = 1)
+ lambda <- exposed / unexposed
+ pi <- unexposed
+ pc <- (pi * ((r * lambda) + 1)) / (r + 1)
+ p1 <- (r + 1) / (r * (lambda - 1)^2 * pi^2)
+ p2 <- z.alpha * sqrt((r + 1) * pc * (1 - pc))
+ p3 <- z.beta * sqrt((lambda * pi * (1 - (lambda * pi))) + (r * pi * (1 - pi)))
+ n <- p1 * (p2 + p3)^2
+
+ # Account for the design effect:
+ n <- n * design
+
+ n.crude <- ceiling(n)
+ n.exposed <- ceiling(n / (r + 1)) * r
+ n.unexposed <- ceiling(n / (r + 1)) * 1
+ n.total <- n.exposed + n.unexposed
+
+ rval <- list(n.total = n.total, n.exposed = n.exposed, n.unexposed = n.unexposed, power = power, lambda = lambda)
}
-
- else
- if(!is.na(exposed) & !is.na(unexposed) & !is.na(n) & is.na(power)){
- # Study power. From Woodward p 409:
- lambda <- exposed / unexposed
- pi <- unexposed
- pc <- (pi * ((r * lambda) + 1)) / (r + 1)
-
- # Account for the design effect:
- n <- n / design
- t1 <- ifelse(lambda >= 1,
- (pi * (lambda - 1) * sqrt(n * r)),
- (pi * (1 - lambda) * sqrt(n * r)))
-
- t2 <- z.alpha * (r + 1) * sqrt(pc * (1 - pc))
- t3 <- (r + 1) * (lambda * pi * (1 - lambda * pi) + r * pi * (1 - pi))
- z.beta <- (t1 - t2) / sqrt(t3)
- power <- pnorm(z.beta, mean = 0, sd = 1)
- rval <- list(power = power)
- }
-
- else
- if(is.na(exposed) & !is.na(unexposed) & !is.na(n) & !is.na(power)){
- # Risk ratio to be detected - requires a value for unexposed. From Woodward p 409:
- z.beta <- qnorm(power, mean = 0, sd = 1)
- pi <- unexposed
+ # Power:
+ else
+ if(!is.na(exposed) & !is.na(unexposed) & !is.na(n) & is.na(power)){
+ # Study power. From Woodward p 409:
+ lambda <- exposed / unexposed
+ pi <- unexposed
+ pc <- (pi * ((r * lambda) + 1)) / (r + 1)
+
+ # Account for the design effect:
+ n <- n / design
+ n.exposed <- ceiling(n / (r + 1)) * r
+ n.unexposed <- ceiling(n / (r + 1)) * 1
+ n.total <- n.exposed + n.unexposed
+
+ t1 <- ifelse(lambda >= 1,
+ (pi * (lambda - 1) * sqrt(n * r)),
+ (pi * (1 - lambda) * sqrt(n * r)))
+
+ t2 <- z.alpha * (r + 1) * sqrt(pc * (1 - pc))
+ t3 <- (r + 1) * (lambda * pi * (1 - lambda * pi) + r * pi * (1 - pi))
+ z.beta <- (t1 - t2) / sqrt(t3)
+ power <- pnorm(z.beta, mean = 0, sd = 1)
+
+ rval <- list(n.total = n.total, n.exposed = n.exposed, n.unexposed = n.unexposed, power = power, lambda = lambda)
+ }
- # Account for the design effect:
- n <- n / design
+ # Lambda:
+ else
+ if(is.na(exposed) & !is.na(unexposed) & !is.na(n) & !is.na(power)){
+ # Risk ratio to be detected - requires a value for unexposed. From Woodward p 409:
+ z.beta <- qnorm(power, mean = 0, sd = 1)
+ pi <- unexposed
+
+ # Account for the design effect:
+ n <- n / design
+ n.exposed <- ceiling(n / (r + 1)) * r
+ n.unexposed <- ceiling(n / (r + 1)) * 1
+ n.total <- n.exposed + n.unexposed
+
+ Y <- r * n * pi^2
+ Z <- (r + 1) * pi * (z.alpha + z.beta)^2
+ a <- Y + (pi * Z)
+ b <- (2 * Y) + Z
+ c <- Y - (r * (1 - pi) * Z)
+ lambda.pos <- (1 / (2 * a)) * (b + sqrt(b^2 - 4 * a * c))
+ lambda.neg <- (1 / (2 * a)) * (b - sqrt(b^2 - 4 * a * c))
+
+ rval <- list(n.total = n.total, n.exposed = n.exposed, n.unexposed = n.unexposed, power = power, lambda = sort(c(lambda.neg, lambda.pos)))
+ }
- Y <- r * n * pi^2
- Z <- (r + 1) * pi * (z.alpha + z.beta)^2
- a <- Y + (pi * Z)
- b <- (2 * Y) + Z
- c <- Y - (r * (1 - pi) * Z)
- lambda.pos <- (1 / (2 * a)) * (b + sqrt(b^2 - 4 * a * c))
- lambda.neg <- (1 / (2 * a)) * (b - sqrt(b^2 - 4 * a * c))
- rval <- list(lambda = sort(c(lambda.neg, lambda.pos)))
- }
-
- rval
+ rval
}
\ No newline at end of file
diff --git a/R/epi.cp.R b/R/epi.cp.R
index 4446bcc..b8a78ee 100644
--- a/R/epi.cp.R
+++ b/R/epi.cp.R
@@ -1,23 +1,30 @@
epi.cp <- function(dat){
- obs <- as.data.frame(cbind(id = 1:nrow(dat), cp = rep(0, times = nrow(dat))))
- nvar <- dim(dat)[2]
- cp <- unique(dat)
- cp <- cbind(id = 1:nrow(cp),cp)
-
- for(i in 1:nrow(cp)){
- tmp <- rep(0, times = nrow(dat))
- for(j in 1:nvar){
- tmp. <- as.numeric(dat[,j] == cp[i,(j+1)])
- tmp <- cbind(tmp, tmp.)
- }
- tmp <- apply(tmp, MARGIN = 1, FUN = sum)
- id <- tmp == nvar
- obs$cp[id] <- cp$id[i]
- }
- n <- hist(obs$cp, breaks = seq(from = 0, to = nrow(cp), by = 1), plot = FALSE)
- n <- n$counts
- end <- nvar + 1
- cov.pattern <- as.data.frame(cbind(id = cp[,1], n, cp[,2:end]))
- rval <- list(cov.pattern = cov.pattern, id = obs$cp)
- rval
+
+ # Re-write of function following email from Johann Popp 11 October 2017.
+ ndat <- data.frame(id = 1:nrow(dat), dat)
+
+ # Add an indicator variable for covariate patterns:
+ ndat$indi <- apply(X = ndat[,ncol(ndat):2], MARGIN = 1, FUN = function(x) as.factor(paste(x, collapse = "")))
+
+ # Order the data according to the indicator variable:
+ ndat <- ndat[order(ndat$indi),]
+
+ # Create a variable that indicates all the cases of each covariate pattern:
+ cp.id <- tapply(ndat$id, ndat$indi, function(x) paste(x, collapse = ","))
+
+ # Create a data frame of unique covariate patterns:
+ cp <- unique(ndat[,2:ncol(ndat)])
+ n <- as.numeric(unlist(lapply(strsplit(cp.id, ","), length)))
+
+ id <- tapply(ndat$id, ndat$indi, function(x) (x)[1])
+ lookup <- data.frame(id = 1:length(n), indi = row.names(id))
+
+ cov.pattern <- data.frame(id = 1:length(n), n, cp[,-ncol(cp)])
+ rownames(cov.pattern) <- rownames(cp)
+
+ # Create a vector with the covariate pattern for each case:
+ id <- lookup$id[match(ndat$indi, lookup$indi)]
+ # id <- as.numeric(unlist(lapply(strsplit(cp.id, ","), function(x) rep(min(as.numeric(unlist(x))), length(x)))))[order(ndat$id)]
+
+ list(cov.pattern = cov.pattern, id = id)
}
diff --git a/R/epi.meansize.R b/R/epi.meansize.R
index 21aec35..d94f1c4 100644
--- a/R/epi.meansize.R
+++ b/R/epi.meansize.R
@@ -1,48 +1,67 @@
"epi.meansize" <- function(treat, control, n, sigma, power, r = 1, design = 1, sided.test = 2, conf.level = 0.95) {
- alpha.new <- (1 - conf.level) / sided.test
- z.alpha <- qnorm(1 - alpha.new, mean = 0, sd = 1)
-
- if(!is.na(treat) & !is.na(control) & is.na(n) & !is.na(sigma) & !is.na(power)){
- # Sample size. From Woodward p 398:
- z.beta <- qnorm(power, mean = 0, sd = 1)
- delta <- abs(treat - control)
- n <- ((r + 1)^2 * (z.alpha + z.beta)^2 * sigma^2) / (delta^2 * r)
+ alpha.new <- (1 - conf.level) / sided.test
+ z.alpha <- qnorm(1 - alpha.new, mean = 0, sd = 1)
- # Account for the design effect:
- n <- n * design
-
- n.crude <- ceiling(n)
- n.treat <- ceiling(n / (r + 1)) * r
- n.control <- ceiling(n / (r + 1)) * 1
- n.total <- n.treat + n.control
- rval <- list(n.crude = n.crude, n.total = n.total, n.treat = n.treat, n.control = n.control)
+ if (!is.na(treat) & !is.na(control) & !is.na(n) & !is.na(power)){
+ stop("Error: at least one of treat, n and power must be NA.")
}
-
- else
- if(!is.na(treat) & !is.na(control) & !is.na(n) & !is.na(sigma) & is.na(power)){
- # Study power. From Woodward p 401:
- delta <- abs(treat - control)
-
- # Account for the design effect:
- n <- n / design
- z.beta <- ((delta * sqrt(n * r)) / ((r + 1) * sigma)) - z.alpha
- power <- pnorm(z.beta, mean = 0, sd = 1)
- rval <- list(power = power)
+ # Sample size:
+ if(!is.na(treat) & !is.na(control) & is.na(n) & !is.na(sigma) & !is.na(power)){
+ # Sample size. From Woodward p 398:
+ z.beta <- qnorm(power, mean = 0, sd = 1)
+ delta <- abs(treat - control)
+ n <- ((r + 1)^2 * (z.alpha + z.beta)^2 * sigma^2) / (delta^2 * r)
+
+ # Account for the design effect:
+ n <- n * design
+
+ n.crude <- ceiling(n)
+ n.treat <- ceiling(n / (r + 1)) * r
+ n.control <- ceiling(n / (r + 1)) * 1
+ n.total <- n.treat + n.control
+
+ rval <- list(n.total = n.total, n.treat = n.treat, n.control = n.control, power = power, delta = delta)
}
-
- else
- if(is.na(treat) & is.na(control) & !is.na(n) & !is.na(sigma) & !is.na(power)){
- # Maximum detectable difference. From Woodward p 401:
- z.beta <- qnorm(power, mean = 0, sd = 1)
- # Account for the design effect:
- n <- n / design
+ # Power:
+ else
+ if(!is.na(treat) & !is.na(control) & !is.na(n) & !is.na(sigma) & is.na(power)){
+ # Study power. From Woodward p 401:
+ delta <- abs(treat - control)
+
+ # Account for the design effect:
+ n <- n / design
+
+ n.crude <- ceiling(n)
+ n.treat <- ceiling(n / (r + 1)) * r
+ n.control <- ceiling(n / (r + 1)) * 1
+ n.total <- n.treat + n.control
+
+ z.beta <- ((delta * sqrt(n * r)) / ((r + 1) * sigma)) - z.alpha
+ power <- pnorm(z.beta, mean = 0, sd = 1)
+
+ rval <- list(n.total = n.total, n.treat = n.treat, n.control = n.control, power = power, delta = delta)
+ }
- delta <- ((r + 1) * (z.alpha + z.beta) * sigma) / (sqrt(n * r))
- rval <- list(delta = delta)
- }
-
- rval
+ # Delta:
+ else
+ if(is.na(treat) & is.na(control) & !is.na(n) & !is.na(sigma) & !is.na(power)){
+ # Maximum detectable difference. From Woodward p 401:
+ z.beta <- qnorm(power, mean = 0, sd = 1)
+
+ # Account for the design effect:
+ n <- n / design
+
+ n.crude <- ceiling(n)
+ n.treat <- ceiling(n / (r + 1)) * r
+ n.control <- ceiling(n / (r + 1)) * 1
+ n.total <- n.treat + n.control
+
+ delta <- ((r + 1) * (z.alpha + z.beta) * sigma) / (sqrt(n * r))
+ rval <- list(n.total = n.total, n.treat = n.treat, n.control = n.control, power = power, delta = delta)
+ }
+
+ rval
}
diff --git a/R/epi.propsize.R b/R/epi.propsize.R
index a959cde..9fb642e 100644
--- a/R/epi.propsize.R
+++ b/R/epi.propsize.R
@@ -1,63 +1,106 @@
"epi.propsize" <- function(treat, control, n, power, r = 1, design = 1, sided.test = 2, conf.level = 0.95) {
- alpha.new <- (1 - conf.level) / sided.test
- z.alpha <- qnorm(1 - alpha.new, mean = 0, sd = 1)
-
- if (!is.na(treat) & !is.na(control) & is.na(n) & !is.na(power)) {
- # Sample size.
- z.beta <- qnorm(power, mean = 0, sd = 1)
- # delta <- abs(treat - control)
- # n <- (1/delta^2) * ((z.alpha * sqrt(treat * (1 - treat))) + (z.beta * sqrt(control * (1 - control))))^2
-
- # From Woodward's spreadsheet. Changed 130814:
- lambda <- treat / control
- Pc <- control * (r * lambda + 1) / (r + 1)
- T1 <- (r + 1) / (r * (lambda - 1)^2 * control^2)
- T2 <- (r + 1) * Pc *(1 - Pc)
- T3 <- lambda * control * (1 - lambda * control) + r * control * (1 - control)
- n <- T1 * (z.alpha * sqrt(T2) + z.beta * sqrt(T3))^2
-
- # Account for the design effect:
- n <- n * design
-
- # n.total <- 2 * ceiling(0.5 * n)
- # rval <- list(n.total = n.total)
- n.crude <- ceiling(n)
- n.treat <- ceiling(n / (r + 1)) * r
- n.control <- ceiling(n / (r + 1)) * 1
- n.total <- n.treat + n.control
- rval <- list(n.crude = n.crude, n.total = n.total, n.treat = n.treat, n.control = n.control)
- }
+ alpha.new <- (1 - conf.level) / sided.test
+ z.alpha <- qnorm(1 - alpha.new, mean = 0, sd = 1)
- else
- if (!is.na(treat) & !is.na(control) & !is.na(n) & is.na(power)) {
- # Power.
- # Account for the design effect:
- n <- n / design
-
- # From Woodward's spreadsheet. Changed 130814:
- lambda <- control / treat
- Pc <- treat * (r * lambda + 1) / (r + 1)
- T1 <- ifelse(lambda >= 1, treat * (lambda - 1) * sqrt(n * r), treat * (1 - lambda) * sqrt(n * r))
- T2 <- z.alpha * (r + 1) * sqrt(Pc * (1 - Pc))
- T3 <- (r + 1) * (lambda * treat * (1 - lambda * treat) + r * treat * (1 - treat))
- z.beta <- (T1 - T2) / sqrt(T3)
- # z.beta <- ((delta * sqrt(n)) - (z.alpha * sqrt(treat * (1 - treat))))/(sqrt(control * (1 - control)))
- power <- pnorm(z.beta, mean = 0, sd = 1)
- rval <- list(power = power)
+ if (!is.na(treat) & !is.na(control) & !is.na(n) & !is.na(power)){
+ stop("Error: at least one of treat, n and power must be NA.")
}
- else
- if (!is.na(treat) & !is.na(control) & !is.na(n) & !is.na(power)) {
- # Maximum detectable difference.
- z.beta <- qnorm(power, mean = 0, sd = 1)
-
- # Account for the design effect:
- n <- n / design
-
- delta <- 1/sqrt(n) * ((z.alpha * sqrt(treat * (1 - treat))) + (z.beta * sqrt(control * (1 - control))))
- rval <- list(delta = delta)
+ # Sample size.
+ if (!is.na(treat) & !is.na(control) & is.na(n) & !is.na(power)) {
+
+ z.beta <- qnorm(power, mean = 0, sd = 1)
+ # delta <- abs(treat - control)
+ # n <- (1/delta^2) * ((z.alpha * sqrt(treat * (1 - treat))) + (z.beta * sqrt(control * (1 - control))))^2
+
+ # From Woodward's spreadsheet. Changed 130814:
+ lambda <- treat / control
+ Pc <- control * (r * lambda + 1) / (r + 1)
+ T1 <- (r + 1) / (r * (lambda - 1)^2 * control^2)
+ T2 <- (r + 1) * Pc *(1 - Pc)
+ T3 <- lambda * control * (1 - lambda * control) + r * control * (1 - control)
+ n <- T1 * (z.alpha * sqrt(T2) + z.beta * sqrt(T3))^2
+
+ # Account for the design effect:
+ n1 <- n / (r + 1)
+ n1 <- ceiling(n1 * design)
+ n2 <- ceiling(r * n1)
+
+ rval <- list(n.total = n1 + n2, n.treat = n1, n.control = n2, power = power, lambda = lambda)
+
}
- rval
+ # Power.
+ else
+ if (!is.na(treat) & !is.na(control) & !is.na(n) & is.na(power)) {
+
+ # Account for the design effect:
+ n1 <- n / (r + 1)
+ n1 <- ceiling(n1 * design)
+ n2 <- ceiling(r * n1)
+
+ # From Woodward's spreadsheet. Changed 130814:
+ lambda <- control / treat
+ Pc <- treat * (r * lambda + 1) / (r + 1)
+ T1 <- ifelse(lambda >= 1, treat * (lambda - 1) * sqrt(n * r), treat * (1 - lambda) * sqrt(n * r))
+ T2 <- z.alpha * (r + 1) * sqrt(Pc * (1 - Pc))
+ T3 <- (r + 1) * (lambda * treat * (1 - lambda * treat) + r * treat * (1 - treat))
+ z.beta <- (T1 - T2) / sqrt(T3)
+ # z.beta <- ((delta * sqrt(n)) - (z.alpha * sqrt(treat * (1 - treat))))/(sqrt(control * (1 - control)))
+ power <- pnorm(z.beta, mean = 0, sd = 1)
+
+ rval <- list(n.total = n1 + n2, n.treat = n1, n.control = n2, power = power, lambda = 1 / lambda)
+ }
+
+ # Lambda:
+ else
+ if (is.na(treat) & !is.na(control) & !is.na(n) & !is.na(power)) {
+
+ z.beta <- qnorm(power, mean = 0, sd = 1)
+
+ # Account for the design effect:
+ n1 <- n / (r + 1)
+ n1 <- ceiling(n1 * design)
+ n2 <- ceiling(r * n1)
+
+ # delta <- 1/sqrt(n) * ((z.alpha * sqrt(treat * (1 - treat))) + (z.beta * sqrt(control * (1 - control))))
+
+ # Here we use the formulae for study power (from Woodward's spreadsheet) and then solve for treat
+ # (which then allows us to calculate lambda).
+ # Note lambda defined as control / treat (hence we take the inverse of lambda for reporting):
+
+ # Where lambda > 1:
+ Pfun <- function(treat, control, n, r, z.alpha){
+ lambda <- control / treat
+ Pc <- treat * (r * lambda + 1) / (r + 1)
+ T1 <- treat * (lambda - 1) * sqrt(n * r)
+ T2 <- z.alpha * (r + 1) * sqrt(Pc * (1 - Pc))
+ T3 <- (r + 1) * (lambda * treat * (1 - lambda * treat) + r * treat * (1 - treat))
+ z.beta <- (T1 - T2) / sqrt(T3)
+
+ # Take the calculated value of the power and subtract the power entered by the user:
+ pnorm(z.beta, mean = 0, sd = 1) - power
+ }
+ treatu <- uniroot(Pfun, control = control, n = n, r = r, z.alpha = z.alpha, interval = c(1E-6,1))$root
+
+ # Where lambda < 1:
+ Pfun <- function(treat, control, n, r, z.alpha){
+ lambda <- control / treat
+ Pc <- treat * (r * lambda + 1) / (r + 1)
+ T1 <- treat * (1 - lambda) * sqrt(n * r)
+ T2 <- z.alpha * (r + 1) * sqrt(Pc * (1 - Pc))
+ T3 <- (r + 1) * (lambda * treat * (1 - lambda * treat) + r * treat * (1 - treat))
+ z.beta <- (T1 - T2) / sqrt(T3)
+
+ # Take the calculated value of the power and subtract the power entered by the user:
+ pnorm(z.beta, mean = 0, sd = 1) - power
+ }
+ treatl <- uniroot(Pfun, control = control, n = n, r = r, z.alpha = z.alpha, interval = c(1E-6,1))$root
+
+ rval <- list(n.total = n1 + n2, n.treat = n1, n.control = n2, power = power, lambda = sort(c(treatu / control, treatl / control)))
+
+ }
+
+ rval
}
diff --git a/man/epi.ccsize.Rd b/man/epi.ccsize.Rd
index 03e935c..5d8d732 100644
--- a/man/epi.ccsize.Rd
+++ b/man/epi.ccsize.Rd
@@ -3,11 +3,11 @@
\alias{epi.ccsize}
\title{
-Sample size or power for an unmatched or matched case-control study
+Sample size, power or minimum detectable odds ratio for an unmatched or matched case-control study
}
\description{
-Computes the sample size or power for an unmatched or matched case-control study.
+Calculates the sample size, power or minimum detectable odds ratio for an unmatched or matched case-control study.
}
\usage{
@@ -126,8 +126,8 @@ epi.ccsize(OR = 2.0, p0 = 0.20, n = NA, power = 0.90, r = 3, rho = 0,
design = 1, sided.test = 2, conf.level = 0.95, method = "unmatched",
fleiss = FALSE)
-## A total of 600 subjects need to be enrolled in the study: 150 cases and
-## 450 controls.
+## A total of 620 subjects need to be enrolled in the study: 155 cases and
+## 465 controls.
## An alternative is to conduct a matched case-control study rather than the
## unmatched design outlined above. One case will be matched to one control
@@ -146,31 +146,42 @@ epi.ccsize(OR = 2.0, p0 = 0.20, n = NA, power = 0.90, r = 1, rho = 0.01,
## EXAMPLE 5:
## Code to reproduce the isograph shown in Figure 2 in Dupont (1988):
r <- 1
-p0 <- seq(from = 0.05, to = 0.95, length = 50)
+p0 = seq(from = 0.05, to = 0.95, length = 50)
OR <- seq(from = 1.05, to = 6, length = 100)
-vals <- expand.grid(p0, OR)
-names(vals) <- c("p0", "OR")
-
-vals$n.case <- NA; vals$n.control <- NA; vals$n.total <- NA
-
-for(i in 1:nrow(vals)){
- trval <- epi.ccsize(OR = vals$OR[i], p0 = vals$p0[i], power = 0.80,
- n = NA, rho = 0, design = 1, sided.test = 2, conf.level = 0.95,
- method = "matched")
- vals$n.case[i] <- trval$n.case
- vals$n.control[i] <- trval$n.control
- vals$n.total[i] <- trval$n.total
-}
+dat <- expand.grid(p0 = p0, OR = OR)
+dat$n.total <- NA
+
+for(i in 1:nrow(dat)){
+ dat$n.total[i] <- epi.ccsize(OR = dat$OR[i], p0 = dat$p0[i], n = NA,
+ power = 0.80, r = 1, rho = 0, design = 1, sided.test = 2,
+ conf.level = 0.95, method = "unmatched", fleiss = FALSE)$n.total
+}
-grid.n <- matrix(vals$n.case, nrow = length(p0))
-lev <- c(22:30,32,34,36,40,45,50,55,60,70,80,90,100,125,150,175,
+grid.n <- matrix(dat$n.total, nrow = length(p0))
+breaks <- c(22:30,32,34,36,40,45,50,55,60,70,80,90,100,125,150,175,
200,300,500,1000)
par(mar = c(5,5,0,5), bty = "n")
-contour(x = p0, y = OR, z = log10(grid.n), add = FALSE, levels = log10(lev),
- labels = lev, xlim = c(0,1), ylim = c(1,6), las = 1, method = "flattest",
+contour(x = p0, y = OR, z = log10(grid.n), add = FALSE, levels = log10(breaks),
+ labels = breaks, xlim = c(0,1), ylim = c(1,6), las = 1, method = "flattest",
xlab = 'Proportion of controls exposed', ylab = "Minimum OR to detect")
+\dontrun{
+## The same plot using ggplot2:
+
+library(ggplot2); library(directlabels)
+
+p <- ggplot(data = dat, aes(x = p0, y = OR, z = n.total)) +
+ geom_contour(aes(colour = ..level..), breaks = breaks) +
+ xlab("Proportion of controls exposed") +
+ ylab("Minimum OR to detect") +
+ xlim(0,1) +
+ ylim(1,6)
+
+print(direct.label(p, list("far.from.others.borders", "calc.boxes",
+ "enlarge.box", hjust = 1, vjust = 1, box.color = NA,
+ fill = "transparent", "draw.rects")))
+}
## EXAMPLE 6:
## From page 1164 of Dupont (1988). A matched case control study is to be
@@ -194,7 +205,7 @@ epi.ccsize(OR = 3.0, p0 = 0.60, n = NA, power = 0.80, r = 3, rho = 0.2,
fleiss = FALSE)
## A total of 204 subjects need to be enrolled in the study: 51 cases and
-## 153 controls (204 subjects in total).
+## 153 controls.
}
diff --git a/man/epi.cohortsize.Rd b/man/epi.cohortsize.Rd
index 6040a2c..c83b155 100644
--- a/man/epi.cohortsize.Rd
+++ b/man/epi.cohortsize.Rd
@@ -3,11 +3,11 @@
\alias{epi.cohortsize}
\title{
-Sample size, power or minimum detectable effect for a cohort study
+Sample size, power or minimum detectable risk ratio for a cohort study
}
\description{
-Computes sample size, power or minimum detectable effect for a cohort study.
+Calculates the sample size, power or minimum detectable risk ratio for a cohort study.
}
\usage{
@@ -31,12 +31,12 @@ The methodology in this function follows the approach described in Chapter 8 of
}
\value{
-A list containing one or more of the following:
- \item{n.crude}{the crude estimated total number of subjects required for the specified level of confidence and power.}
- \item{n.total}{the total estimated number of subjects required for the specified level of confidence and power, respecting the requirement for \code{r} times as many individuals in the treatment group compared with the control group.}
- \item{delta}{the minimum detectable difference given the specified level of confidence and power.}
- \item{lambda}{the minimum detectable risk ratio >1 and the maximum detectable risk ratio <1.}
- \item{power}{the power of the study given the number of study subjects, the expected risk ratio and level of confidence.}
+A list containing the following:
+ \item{n.total}{the total number of subjects required for the specified level of confidence and power, respecting the requirement for \code{r} times as many individuals in the treatment group compared with the control group.}
+ \item{n.treat}{the total number of subjects in the treatment group for the specified level of confidence and power, respecting the requirement for \code{r} times as many individuals in the treatment group compared with the control group.}
+ \item{n.control}{the total number of subjects in the control group for the specified level of confidence and power, respecting the requirement for \code{r} times as many individuals in the treatment group compared with the control group.}
+ \item{power}{the power of the study given the number of study subjects, the expected effect size and level of confidence.}
+ \item{lambda}{the outcome proportion in the exposed group divided by the outcome proportion in the unexposed group (a risk ratio).}
}
\references{
@@ -47,6 +47,10 @@ Woodward M (2005). Epidemiology Study Design and Data Analysis. Chapman & Hall/C
\note{
The power of a study is its ability to demonstrate the presence of an association, given that an association actually exists.
+
+Values need to be entered for \code{unexposed}, \code{n}, and \code{power} to return a value for \code{lambda}. In this situation, the lower value of lambda represents the maximum detectable risk ratio that is less than 1; the upper value of lambda represents the minimum detectable risk ratio greater than 1.
+
+
}
\examples{
@@ -61,9 +65,8 @@ The power of a study is its ability to demonstrate the presence of an associatio
## 100,000 per year. Assuming equal numbers of smokers and non-smokers are
## sampled, how many men should be sampled overall?
-exposed = 1.4 * (5 * 413)/100000
-unexposed = (5 * 413)/100000
-epi.cohortsize(exposed = exposed, unexposed = unexposed, n = NA, power = 0.90,
+e1 = 1.4 * (5 * 413)/100000; e0 = (5 * 413)/100000
+epi.cohortsize(exposed = e1, unexposed = e0, n = NA, power = 0.90,
r = 1, design = 1, sided.test = 1, conf.level = 0.95)
## A total of 12,130 men need to be sampled (6065 smokers and 6065 non-smokers).
@@ -73,8 +76,8 @@ epi.cohortsize(exposed = exposed, unexposed = unexposed, n = NA, power = 0.90,
## Say, for example, we are only able to enrol 5000 subjects into the study
## described above. What is the minimum and maximum detectable risk ratio?
-unexposed = (5 * 413)/100000
-epi.cohortsize(exposed = NA, unexposed = unexposed, n = 5000, power = 0.90,
+e0 = (5 * 413)/100000
+epi.cohortsize(exposed = NA, unexposed = e0, n = 5000, power = 0.90,
r = 1, design = 1, sided.test = 1, conf.level = 0.95)
## The minimum detectable risk ratio >1 is 1.65. The maximum detectable
@@ -97,13 +100,14 @@ epi.cohortsize(exposed = 0.30, unexposed = 0.15, n = NA, power = 0.80,
## Assume now that this study is going to be carried out using animals from a
## number of herds. What is the required sample size when you account for the
-## observation that response to treatment is likely to cluster across herds.
+## observation that response to treatment is likely to cluster within herds.
## For the exercise, assume that the intra-cluster correlation coefficient
-## (the rate of homogeneity, rho) is 0.05 and the average number of cows per
-## herd is 30. Calculate the design effect, given
-## rho = (design - 1) / (nbar - 1), where nbar equals the average number of
-## individuals per cluster:
+## (the rate of homogeneity, rho) for this treatment is 0.05 and the
+## average number of cows sampled per herd will be 30.
+
+## Calculate the design effect, given rho = (design - 1) / (nbar - 1),
+## where nbar equals the average number of individuals per cluster:
design <- 0.05 * (30 - 1) + 1
epi.cohortsize(exposed = 0.30, unexposed = 0.15, n = NA, power = 0.80,
diff --git a/man/epi.cp.Rd b/man/epi.cp.Rd
index f0518aa..f48a44c 100644
--- a/man/epi.cp.Rd
+++ b/man/epi.cp.Rd
@@ -15,17 +15,21 @@ epi.cp(dat)
}
\arguments{
- \item{dat}{an \emph{i} row by \emph{j} column data frame where the \emph{i} rows represent individual observations and the \emph{m} columns represent covariates.}
+ \item{dat}{an \emph{i} row by \emph{j} column data frame where the \emph{i} rows represent individual observations and the \emph{m} columns represent a set of \emph{m} covariates.}
}
\details{
-A covariate pattern is a unique combination of values of predictor variables. For example, if a model contains two dichotomous predictors, there will be four covariate patterns possible: \code{(1,1)}, \code{(1,0)}, \code{(0,1)}, and \code{(0,0)}. This function extracts the \emph{n} unique covariate patterns from a data set comprised of \emph{i} observations, labelling them from 1 to \emph{n}. A vector of length \emph{m} is also returned, listing the covariate pattern identifier for each [...]
+This function provides a summary of the covariate patterns (unique combinations of values of covariates) present in a data. This function extracts the \emph{k} unique covariate patterns in a data set comprised of \emph{i} observations, labelling them from 1 to \emph{k}. The frequency of occurrence of each covariate pattern is listed. A vector of length \emph{i} is also returned, listing the 1:\emph{k} covariate pattern identifier for each observation.
}
\value{
A list containing the following:
- \item{cov.pattern}{a data frame with columns: \code{id} the unique covariate patterns, \code{n} the number of occasions each of the listed covariate pattern appears in the data, and the unique covariate combinations.}
- \item{id}{a vector listing the covariate pattern identifier for each observation.}
+ \item{cov.pattern}{a data frame with columns: \code{id} the unique covariate pattern identifier (labelled 1 to \emph{k}), \code{n} the number of occasions each of the listed covariate pattern appears in the data, and the unique covariate combinations.}
+ \item{id}{a vector of length \emph{i} listing the 1:\emph{k} covariate pattern identifier for each observation.}
+}
+
+\author{
+Thanks to Johann Popp for providing code and suggestions to enhance the speed of this function.
}
\references{
diff --git a/man/epi.meansize.Rd b/man/epi.meansize.Rd
index 3ad70a7..4c7d0f5 100644
--- a/man/epi.meansize.Rd
+++ b/man/epi.meansize.Rd
@@ -7,7 +7,7 @@ Sample size, power and minimum detectable difference when comparing means
}
\description{
-Computes the sample size, power or minimum detectable difference when comparing means.
+Calculates the sample size, power or minimum detectable difference when comparing means.
}
\usage{
@@ -32,11 +32,12 @@ The methodology in this function follows the approach described in Chapter 8 of
}
\value{
-A list containing one or more of the following:
- \item{n.crude}{the crude estimated total number of subjects required for the specified level of confidence and power.}
- \item{n.total}{the total estimated number of subjects required for the specified level of confidence and power, respecting the requirement for \code{r} times as many individuals in the treatment group compared with the control group.}
+A list containing the following:
+ \item{n.total}{the total number of subjects required for the specified level of confidence and power, respecting the requirement for \code{r} times as many individuals in the treatment group compared with the control group.}
+ \item{n.treat}{the total number of subjects in the treatment group for the specified level of confidence and power, respecting the requirement for \code{r} times as many individuals in the treatment group compared with the control group.}
+ \item{n.control}{the total number of subjects in the control group for the specified level of confidence and power, respecting the requirement for \code{r} times as many individuals in the treatment group compared with the control group.}
+ \item{power}{the power of the study given the number of study subjects, the expected effect size and level of confidence.}
\item{delta}{the minimum detectable difference given the specified level of confidence and power.}
- \item{power}{the power of the study given the number of study subjects, the expected effect size and level of confidence.}
}
\references{
@@ -89,9 +90,8 @@ epi.meansize(treat = 5, control = 4.5, n = NA, sigma = 1.4, power = 0.95,
epi.meansize(treat = 8, control = 4, n = NA, sigma = 4, power = 0.90,
r = 4, design = 1, sided.test = 2, conf.level = 0.95)
-## The estimated sample size is 66. We round this up to the nearest multiple
-## of 5, to 70. We allocate 70/5 = 14 women to the placebo group and four
-## times as many (56) to the iron treatment group.
+## The estimated sample size is 70. We allocate 70/5 = 14 women to the
+## placebo group and four times as many (56) to the iron treatment group.
}
diff --git a/man/epi.propsize.Rd b/man/epi.propsize.Rd
index 9cb74a6..7ee9d5d 100644
--- a/man/epi.propsize.Rd
+++ b/man/epi.propsize.Rd
@@ -3,11 +3,11 @@
\alias{epi.propsize}
\title{
-Sample size, power and minimum detectable difference when comparing proportions
+Sample size, power and minimum detectable risk ratio when comparing proportions
}
\description{
-Computes the sample size, power or minimum detectable difference when comparing proportions.
+Calculates the sample size, power or minimum detectable risk ratio when comparing proportions.
}
\usage{
@@ -16,9 +16,9 @@ epi.propsize(treat, control, n, power, r = 1, design = 1,
}
\arguments{
- \item{treat}{the expected value for the treatment group (see below).}
- \item{control}{the expected value for the control group (see below).}
- \item{n}{scalar, defining the total number of subjects in the study (i.e. the number in the treatment and control group).}
+ \item{treat}{the expected proportion for the treatment group (see below).}
+ \item{control}{the expected proportion for the control group (see below).}
+ \item{n}{scalar, defining the total number of subjects in the study (i.e. the number in the treatment plus the number in the control group).}
\item{power}{scalar, the required study power.}
\item{r}{scalar, the number in the treatment group divided by the number in the control group.}
\item{design}{scalar, the estimated design effect.}
@@ -27,7 +27,7 @@ epi.propsize(treat, control, n, power, r = 1, design = 1,
}
\details{
-The methodology in this function follows closely the approach described in Chapter 8 of Woodward (2005).
+The methodology in this function follows the approach described in Chapter 8 of Woodward (2005).
With this function it is assumed that one of the two proportions is known and we want to test the null hypothesis that the second proportion is equal to the first. Users are referred to the \code{\link{epi.cohortsize}} function which relates to the two-sample problem where neither proportion is known (or assumed, at least).
@@ -35,11 +35,12 @@ Because there is much more uncertainty in the two sample problem where neither p
}
\value{
-A list containing one or more of the following:
- \item{n.crude}{the crude estimated total number of subjects required for the specified level of confidence and power.}
- \item{n.total}{the total estimated number of subjects required for the specified level of confidence and power, respecting the requirement for \code{r} times as many individuals in the treatment group compared with the control group.}
- \item{delta}{the minimum detectable difference given the specified level of confidence and power.}
- \item{power}{the power of the study given the number of study subjects, the expected effect size and level of confidence.}
+A list containing the following:
+ \item{n.total}{the total number of subjects required for the specified level of confidence and power, respecting the requirement for \code{r} times as many individuals in the treatment group compared with the control group.}
+ \item{n.treat}{the total number of subjects in the treatment group for the specified level of confidence and power, respecting the requirement for \code{r} times as many individuals in the treatment group compared with the control group.}
+ \item{n.control}{the total number of subjects in the control group for the specified level of confidence and power, respecting the requirement for \code{r} times as many individuals in the treatment group compared with the control group.}
+ \item{power}{the power of the study given the number of study subjects, the expected effect size and level of confidence.}
+ \item{lambda}{the proportion in the treatment group divided by the proportion in the control group (a risk ratio).}
}
\references{
@@ -53,7 +54,7 @@ Woodward M (2005). Epidemiology Study Design and Data Analysis. Chapman & Hall/C
\note{
The power of a study is its ability to demonstrate the presence of an association, given that an association actually exists.
-Values need to be entered for \code{control}, \code{n}, and \code{power} to return a value for \code{delta}.
+Values need to be entered for \code{control}, \code{n}, and \code{power} to return a value for \code{lambda}. In this situation, the lower value of lambda represents the maximum detectable risk ratio that is less than 1; the upper value of lambda represents the minimum detectable risk ratio greater than 1.
}
\examples{
@@ -68,8 +69,21 @@ Values need to be entered for \code{control}, \code{n}, and \code{power} to retu
epi.propsize(treat = 0.30, control = 0.32, n = NA, power = 0.90,
r = 1, design = 1, sided.test = 1, conf.level = 0.95)
-## ## A total of 18,315 men should be sampled: 9158 in the treatment group and
-## 9158 in the control group.
+## ## A total of 18,316 men should be sampled: 9158 in the treatment group and
+## 9158 in the control group. The risk ratio (that is, the prevalence of
+## smoking in males post government initiative divided by the prevalence of
+## smoking in males pre inititative is 0.94.
+
+
+## EXAMPLE 2:
+## If we sample only 10,000 men (5000 in the treatment group and 5000 in the
+## control group) what is the maximum detectable risk ratio that is less
+## than 1?
+
+epi.propsize(treat = NA, control = 0.32, n = 10000, power = 0.90,
+ r = 1, design = 1, sided.test = 1, conf.level = 0.95)
+
+## If we sample only 10,000 men the maximum detectable risk ratio will be 0.91.
}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/r-cran-epir.git
More information about the debian-med-commit
mailing list