-
Notifications
You must be signed in to change notification settings - Fork 31
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
metacont(): bug fix for exact SMD method if total sample size is abov…
…e 343 (#42)
- Loading branch information
Showing
8 changed files
with
312 additions
and
153 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,79 +1,168 @@ | ||
| ciAgrestiCoull <- function(event, n, level = 0.95) { | ||
|
|
||
| if (level <= 0 | level >= 1) | ||
| stop("no valid level for confidence interval") | ||
|
|
||
| ciAgrestiCoull <- function(event, n, level) { | ||
| chknumeric(event, 0) | ||
| chknumeric(n, 0, zero = TRUE) | ||
| chklevel(level) | ||
| ## | ||
| if (length(event) == 1 & length(n) > 1) | ||
| event <- rep(event, length(n)) | ||
| else if (length(event) > 1 & length(n) == 1) | ||
| n <- rep(n, length(event)) | ||
| else | ||
| if (length(event) != length(n)) | ||
| stop("Arguments 'event' and 'n' must be of same length.", call. = FALSE) | ||
| ## | ||
| if (any(event > n, na.rm = TRUE)) | ||
| stop("Number of events must be smaller equal sample size.", call. = FALSE) | ||
| ## | ||
| z <- qnorm(1 - (1 - level) / 2) | ||
| ## | ||
| n <- n + z^2 | ||
| p <- 1 / n * (event + 0.5 * z^2) | ||
|
|
||
| lower <- p - z * sqrt(1 / n * p * (1 - p)) | ||
| upper <- p + z * sqrt(1 / n * p * (1 - p)) | ||
| prop <- 1 / n * (event + 0.5 * z^2) | ||
| ## | ||
| lower <- prop - z * sqrt(1 / n * prop * (1 - prop)) | ||
| upper <- prop + z * sqrt(1 / n * prop * (1 - prop)) | ||
| ## | ||
| lower[lower < 0] <- 0 | ||
| upper[upper > 1] <- 1 | ||
|
|
||
| list(n = n, p = p, lower = lower, upper = upper, level = level) | ||
| ## | ||
| list(event = event, n = n, | ||
| prop = prop, lower = lower, upper = upper, | ||
| statistic = NA, p = NA, level = level, | ||
| df = NA, null.effect = NA) | ||
| } | ||
|
|
||
|
|
||
| ciClopperPearson <- function(event, n, level, null.effect) { | ||
| chknumeric(event, 0) | ||
| chknumeric(n, 0, zero = TRUE) | ||
| chklevel(level) | ||
| chknumeric(null.effect, min = 0, max = 1, length = 1) | ||
| ## | ||
| if (length(event) == 1 & length(n) > 1) | ||
| event <- rep(event, length(n)) | ||
| else if (length(event) > 1 & length(n) == 1) | ||
| n <- rep(n, length(event)) | ||
| else | ||
| if (length(event) != length(n)) | ||
| stop("Arguments 'event' and 'n' must be of same length.", call. = FALSE) | ||
| ## | ||
| if (any(event > n, na.rm = TRUE)) | ||
| stop("Number of events must be smaller equal sample size.", call. = FALSE) | ||
| ## | ||
| k <- length(event) | ||
| lower <- upper <- statistic <- pval <- rep(NA, k) | ||
| ## | ||
| for (i in seq_len(k)) { | ||
| if (!is.na(event[i] & !is.na(n[i]))) { | ||
| cint <- binom.test(event[i], n[i], conf.level = level, | ||
| p = if (!is.na(null.effect)) null.effect else 0.5) | ||
| ## | ||
| lower[i] <- cint$conf.int[[1]] | ||
| upper[i] <- cint$conf.int[[2]] | ||
| if (!is.na(null.effect)) | ||
| pval[i] <- cint$p.value | ||
| } | ||
| else { | ||
| lower[i] <- NA | ||
| upper[i] <- NA | ||
| pval[i] <- NA | ||
| } | ||
| } | ||
| ## | ||
| list(event = event, n = n, | ||
| prop = event / n, lower = lower, upper = upper, | ||
| statistic = statistic, p = pval, level = level, | ||
| df = NA, null.effect = null.effect) | ||
| } | ||
| ciSimpleAsymptotic <- function(event, n, level = 0.95, correct = FALSE) { | ||
|
|
||
| ## Reference: | ||
|
|
||
|
|
||
| ciSimpleAsymptotic <- function(event, n, level, correct = FALSE) { | ||
| ## | ||
| ## Newcombe RG. Two-sided confidence intervals for the single | ||
| ## proportion: Comparison of seven methods. | ||
| ## Stat Med 1998, Apr 30;17(8): 857-72 | ||
|
|
||
| if (level <= 0 | level >= 1) | ||
| stop("no valid level for confidence interval") | ||
|
|
||
| p <- event / n | ||
| q <- 1 - p | ||
| ## | ||
| chknumeric(event, 0) | ||
| chknumeric(n, 0, zero = TRUE) | ||
| chklevel(level) | ||
| chklogical(correct) | ||
| ## | ||
| if (length(event) == 1 & length(n) > 1) | ||
| event <- rep(event, length(n)) | ||
| else if (length(event) > 1 & length(n) == 1) | ||
| n <- rep(n, length(event)) | ||
| else | ||
| if (length(event) != length(n)) | ||
| stop("Arguments 'event' and 'n' must be of same length.", call. = FALSE) | ||
| ## | ||
| if (any(event > n, na.rm = TRUE)) | ||
| stop("Number of events must be smaller equal sample size.", call. = FALSE) | ||
| ## | ||
| prop <- event / n | ||
| z <- qnorm(1 - (1 - level) / 2) | ||
|
|
||
| ## | ||
| if (!correct) { | ||
| lower <- p - z * sqrt(p * q / n) | ||
| upper <- p + z * sqrt(p * q / n) | ||
| lower <- prop - z * sqrt(prop * (1 - prop) / n) | ||
| upper <- prop + z * sqrt(prop * (1 - prop) / n) | ||
| } | ||
| else { | ||
| lower <- p - (z * sqrt(p * q / n) + 1 / (2 * n)) | ||
| upper <- p + (z * sqrt(p * q / n) + 1 / (2 * n)) | ||
| lower <- prop - (z * sqrt(prop * (1 - prop) / n) + 1 / (2 * n)) | ||
| upper <- prop + (z * sqrt(prop * (1 - prop) / n) + 1 / (2 * n)) | ||
| } | ||
| ## | ||
| lower[lower < 0] <- 0 | ||
| upper[upper > 1] <- 1 | ||
|
|
||
| list(n = n, p = p, lower = lower, upper = upper, level = level) | ||
| ## | ||
| list(event = event, n = n, | ||
| prop = prop, lower = lower, upper = upper, | ||
| statistic = NA, p = NA, level = level, | ||
| df = NA, null.effect = NA) | ||
| } | ||
| ciWilsonScore <- function(event, n, level = 0.95, correct = FALSE) { | ||
|
|
||
| ## Reference: | ||
|
|
||
|
|
||
| ciWilsonScore <- function(event, n, level, correct = FALSE) { | ||
| ## | ||
| ## Newcombe RG. Two-sided confidence intervals for the single | ||
| ## proportion: Comparison of seven methods. | ||
| ## Stat Med 1998, Apr 30;17(8): 857-72 | ||
|
|
||
| if (level <= 0 | level >= 1) | ||
| stop("no valid level for confidence interval") | ||
|
|
||
| p <- event / n | ||
| q <- 1 - p | ||
| ## | ||
| chknumeric(event, 0) | ||
| chknumeric(n, 0, zero = TRUE) | ||
| chklevel(level) | ||
| ## | ||
| if (length(event) == 1 & length(n) > 1) | ||
| event <- rep(event, length(n)) | ||
| else if (length(event) > 1 & length(n) == 1) | ||
| n <- rep(n, length(event)) | ||
| else | ||
| if (length(event) != length(n)) | ||
| stop("Arguments 'event' and 'n' must be of same length.", call. = FALSE) | ||
| ## | ||
| if (any(event > n, na.rm = TRUE)) | ||
| stop("Number of events must be smaller equal sample size.", call. = FALSE) | ||
| ## | ||
| prop <- event / n | ||
| z <- qnorm(1 - (1 - level) / 2) | ||
|
|
||
| ## | ||
| if (!correct) { | ||
| lower <- (2 * n * p + z^2 - | ||
| z * sqrt(z^2 + 4 * n * p * q)) / (2 * (n + z^2)) | ||
| upper <- (2 * n * p + z^2 + | ||
| z * sqrt(z^2 + 4 * n * p * q)) / (2 * (n + z^2)) | ||
| lower <- (2 * n * prop + z^2 - | ||
| z * sqrt(z^2 + 4 * n * prop * (1 - prop))) / (2 * (n + z^2)) | ||
| upper <- (2 * n * prop + z^2 + | ||
| z * sqrt(z^2 + 4 * n * prop * (1 - prop))) / (2 * (n + z^2)) | ||
| } | ||
| else { | ||
| lower <- (2 * n * p + z^2 - 1 - | ||
| z * sqrt(z^2 - 2 - 1/n + 4 * p * (n * q + 1))) / | ||
| lower <- (2 * n * prop + z^2 - 1 - | ||
| z * sqrt(z^2 - 2 - 1/n + 4 * prop * (n * (1 - prop) + 1))) / | ||
| (2 * (n + z^2)) | ||
| lower[lower < 0] <- 0 | ||
| upper <- (2 * n * p + z^2 + 1 + | ||
| z * sqrt(z^2 + 2 - 1/n + 4 * p * (n * q - 1))) / | ||
| upper <- (2 * n * prop + z^2 + 1 + | ||
| z * sqrt(z^2 + 2 - 1/n + 4 * prop * (n * (1 - prop) - 1))) / | ||
| (2 * (n + z^2)) | ||
| upper[upper > 1] <- 1 | ||
| } | ||
|
|
||
| list(n = n, p = p, lower = lower, upper = upper, level = level) | ||
| ## | ||
| list(event = event, n = n, | ||
| prop = prop, lower = lower, upper = upper, | ||
| statistic = NA, p = NA, level = level, | ||
| df = NA, null.effect = NA) | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.