R/confIntProp.R In userfriendlyscience: Quantitative Analysis Made Accessible

Documented in confIntProp

```#' Confidence intervals for proportions, vectorized over all arguments
#'
#' This function simply computes confidence intervals for proportions.
#'
#' uses \code{\link{pbeta}}, with some lines of code taken from the
#' \code{\link{binom.test}} source. Specifically, the count for the low
#' category is specified as first 'shape argument' to \code{\link{pbeta}}, and
#' the total count (either the sum of the count for the low category and the
#' count for the high category, or the total number of cases if
#' \code{compareHiToLo} is \code{FALSE}) minus the count for the low category
#' as the second 'shape argument'.
#'
#' @param x The number of 'successes', i.e. the number of events, observations,
#' or cases that one is interested in.
#' @param n The total number of cases or observatons.
#' @param conf.level The confidence level.
#' @return The confidence interval bounds in a twodimensional matrix, with the
#' first column containing the lower bound and the second column containing the
#' upper bound.
#'
#' Maintainer: Gjalt-Jorn Peters <gjalt-jorn@@userfriendlyscience.com>
#' function for which this was written.}
#' @keywords univar htest
#' @examples
#'
#'   ### Simple case
#'   confIntProp(84, 200);
#'
#'   ### Using vectors
#'   confIntProp(c(2,3), c(10, 20), conf.level=c(.90, .95, .99));
#'
#' @export confIntProp
confIntProp <- function(x, n, conf.level = .95) {

originalN <- n;
originalX <- x;

if (length(x) > length(n)) {
if (length(n) == 1) {
n <- rep(n, length(x));
} else {
stop("Length of 'x' is ", length(x),
", but length of 'n' is ", length(n), ".");
}
} else if (length(x) < length(n)) {
if (length(x) == 1) {
x <- rep(x, length(n));
} else {
stop("Length of 'x' is ", length(x),
", but length of 'n' is ", length(n), ".");
}
}

p.L <- function(x, n, alpha) {
if (x == 0)
0
else qbeta(alpha, x, n - x + 1)
}
p.U <- function(x, n, alpha) {
if (x == n)
1
else qbeta(1 - alpha, x + 1, n - x)
}
### END binom.test code

### Vectorized both for x/n pairs and for confidence levels;
### first process for each confidence level
res <- sapply(conf.level, function(confLevel) {
### Then for x/n pairs for this given confidence level
return(sapply(seq_along(originalX), function(i) {
### Then for n
return(sapply(seq_along(originalN), function(j) {
return(c(x[i] / n[j],
p.L(x=x[i],
n=n[j],
alpha=(1 - confLevel)/2),
p.U(x=x[i],
n=n[j],
alpha=(1 - confLevel)/2)));
}));
}));
});

if (ncol(res) > 1) {
res <- lapply(1:ncol(res), function(confI) {
rslt <- matrix(res[, confI], ncol=3, byrow=TRUE);
rownames(rslt) <- paste0(rslt[, 1], ", ",
rep(100*conf.level[confI], each=nrow(rslt)), "%");
rslt <- rslt[, -1];

if (is.null(ncol(rslt))) {
names(rslt) <- c('ci.lo', 'ci.hi');
} else {
colnames(rslt) <- c('ci.lo', 'ci.hi');
}
return(rslt);
});
res <- do.call(rbind, res);
if (is.null(row.names(res))) {
row.names(res) <- paste0(x/n, ", ", 100*conf.level, "%");
}
} else {
res <- matrix(res, ncol=3, byrow=TRUE);
rownames(res) <- paste0(res[, 1], ", ",
rep(100*conf.level, each=nrow(res)), "%");
res <- res[, -1, drop=FALSE];
colnames(res) <- c('ci.lo', 'ci.hi');
}
return(res);
}
```

Try the userfriendlyscience package in your browser

Any scripts or data that you put into this service are public.

userfriendlyscience documentation built on May 2, 2019, 1:09 p.m.