jared
38cdaa8ea7
R-CMD-check / R CMD check (push) Failing after 3s
- Add real R CMD check workflow via Gitea Actions (rocker/r-ver:4.4) - Remove Jenkinsfile and demo workflow - Modernize DESCRIPTION: Authors@R, R >= 4.1.0, URL/BugReports fields, move tidycensus from Imports to Suggests, testthat edition 3 - Fix agresti_coull_interval: correct implementation, export it - Convert get_fips/get_stabbr to use requireNamespace for tidycensus - Remove civilytics:: self-reference in rnh()
93 lines
3.1 KiB
R
93 lines
3.1 KiB
R
|
|
|
|
# https://github.com/cran/binom/blob/master/R/binom.confint.R
|
|
# Consider importing and crediting this code ^^
|
|
# https://towardsdatascience.com/five-confidence-intervals-for-proportions-that-you-should-know-about-7ff5484c024f
|
|
# https://andrewpwheeler.com/2020/11/30/confidence-intervals-around-proportions/
|
|
#' Get a simple Clopper Pearson interval
|
|
#'
|
|
#' @param num number of successes
|
|
#' @param den number of trials
|
|
#' @param conf.level default 0.95, set the confidence interval to return
|
|
#'
|
|
#' @return three values forming the upper and lower bounds of the confidence region and the true value
|
|
#' @export
|
|
clopper_pearson <- function(num, den, conf.level = 0.95) {
|
|
# Same results as binom.test in base R
|
|
quant <- (1 - conf.level) / 2
|
|
low <- qbeta(quant, num, den-num+1)
|
|
hi <- qbeta(1-quant, num+1, den-num)
|
|
obs <- num/den
|
|
return(c("low" = low, "observed" = obs,"high" = hi))
|
|
}
|
|
|
|
|
|
# z_gap_test_v <- Vectorize(z_gap_test,
|
|
# SIMPLIFY = TRUE)# we only want to return a scalar
|
|
|
|
|
|
#z_gap_test(a_prop = 0.051, a_count = 2000, b_prop = 0.11, b_count = 100)
|
|
|
|
|
|
#' Calculate a univariate z score by comparing to a population
|
|
#'
|
|
#' @param unit_prop proportion for the group we are comparing
|
|
#' @param global_prop the global proportion
|
|
#' @param unit_denom the population size for the group we are comparing
|
|
#'
|
|
#' @return a z-score
|
|
#' @export
|
|
#'
|
|
#' @examples
|
|
#' z_univariate(unit_prop = 0.13, global_prop = 0.11, unit_denom = 2500)
|
|
z_univariate <- function(unit_prop, global_prop, unit_denom) {
|
|
num <- unit_prop - global_prop
|
|
denom <- sqrt(
|
|
(global_prop * (1-global_prop))/unit_denom
|
|
)
|
|
z = num / denom
|
|
return(z)
|
|
|
|
}
|
|
|
|
#' Calculate a Wald interval
|
|
#'
|
|
#' @param x the numerator, number of times the event occurs
|
|
#' @param n the denominator, the number of trials
|
|
#' @param conf.level default 0.95, set the confidence interval to return
|
|
#'
|
|
#' @return two values forming the upper and lower bounds of the confidence region
|
|
#' @export
|
|
#'
|
|
#' @examples
|
|
#' waldInterval(x = 20, n =40) #this will return 0.345 and 0.655
|
|
waldInterval <- function(x, n, conf.level = 0.95){
|
|
p <- x/n
|
|
sd <- sqrt(p*((1-p)/n))
|
|
z <- qnorm(c( (1 - conf.level)/2, 1 - (1-conf.level)/2)) #returns the value of thresholds at which conf.level has to be cut at. for 95% CI, this is -1.96 and +1.96
|
|
ci <- p + z*sd
|
|
names(ci) <- c('lwr', 'upr')
|
|
return(ci)
|
|
}
|
|
|
|
#' Calculate the Agresti-Coull interval
|
|
#'
|
|
#' @param num number of successes
|
|
#' @param den number of trials
|
|
#' @param conf.level default 0.95, confidence level for the interval
|
|
#'
|
|
#' @return three values forming the lower bound, observed proportion, and upper bound
|
|
#' @export
|
|
#'
|
|
#' @examples
|
|
#' agresti_coull_interval(20, 40)
|
|
#' agresti_coull_interval(2, 100, conf.level = 0.99)
|
|
agresti_coull_interval <- function(num, den, conf.level = 0.95) {
|
|
z <- qnorm(1 - (1 - conf.level) / 2)
|
|
n_tilde <- den + z^2
|
|
p_tilde <- (num + z^2 / 2) / n_tilde
|
|
margin <- z * sqrt(p_tilde * (1 - p_tilde) / n_tilde)
|
|
obs <- num / den
|
|
return(c("low" = p_tilde - margin, "observed" = obs, "high" = p_tilde + margin))
|
|
}
|