Merge pull request #35 from eduardgrebe/master

Update frrcal and bump inctools version to 1.0.15 for release

Author: Eduard Grebe

inctools-julia/Incidence.jl

@@ -7,14 +7,41 @@ import Distributions
 import Statistics
 import DataFrames
 
-function prevalence(pos, n, de = 1) #, f = 1
+function prevalence(pos, n, de = 1; ci = false, α = 0.05) #, f = 1
     if n == 0
         @warn "n = 0: prevalence undefined. n set to 1."
         n = 1
     end
+    if pos < 5
+        @warn "Too few successes for the normal approximation to be valid"
+    end
+    if n - pos < 5
+        @warn "Sample size twoo small for the normal apprximation to be valid"
+    end
     p = pos/n
     σ = sqrt( (p * (1 - p)) / n ) * de #* sqrt(1 - f)
-    return p, σ
+    if !ci
+        return (p = p, σ = σ)
+    elseif ci
+        Fd =
+        if pos == 0
+            lb = 0
+            ub = 1-(α/2)^(1/n)
+        elseif pos == n
+            lb = (α/2)^(1/n)
+            ub = 1
+        else
+            lb = ( 1+(n-pos+1)/(pos * Distributions.quantile(Distributions.FDist(2*pos, 2*(n-pos+1)), α/2)) ) ^ (-1)
+            ub = ( 1+(n-pos)/((pos + 1) * Distributions.quantile(Distributions.FDist(2*(pos + 1), 2*(n-pos)), 1-α/2)) ) ^ (-1)
+        end
+
+        # Normal approximation
+        # pprime = (pos + 2) / (n + 4)
+        # sigmaprime = sqrt( (pprime * (1 - pprime)) / (n + 4) )
+        #lb = max(pprime - 1.96 * sigmaprime,0)
+        # ub = min(pprime + 1.96 * sigmaprime,1)
+        return (p = p, σ = σ, ci = [lb, ub])
+    end
 end
 
 function kassanjee(prev::Float64, prevR::Float64, mdri::Float64, frr::Float64, T::Float64)

inctools/DESCRIPTION

@@ -2,10 +2,10 @@ Encoding: UTF-8
 Package: inctools
 Type: Package
 Title: Incidence Estimation Tools
-Version: 1.0.14.9005
-Date: 2019-08-28
+Version: 1.0.15
+Date: 2019-11-05
 Authors@R: c(
-  person("Eduard","Grebe",email="[email protected]",role=c("cre","aut")),
+  person("Eduard","Grebe",email="[email protected]",role=c("cre","aut")),
   person("Alex","Welte", email="[email protected]", role="aut"),
   person("Avery","McIntosh",email="[email protected]",role="aut"),
   person("Petra","Bäumler",email="[email protected]",role="aut"),
@@ -43,7 +43,8 @@ Imports:
     pracma,
     foreach,
     parallel,
-    doParallel
+    doParallel,
+    binom
 Suggests:
     survey,
     knitr,

inctools/R/calibration_frr.R

@@ -29,6 +29,7 @@
 #' @param recency_vars Variables to be used in determining recency outcomes
 #' @param recency_params Vector of numeric parameters (e.g. thresholds) for determining recency according to the relevant rule
 #' @param alpha Confidence level, default=0.05.
+#' @param method Method for computing confidence interval on binomial probability (passed to binom::binom.confint). Default is Clopper-Pearson 'exact' method. Accepted values: `c("exact", "ac", "asymptotic", "wilson", "prop.test", "bayes", "logit", "cloglog", "probit")`.
 #' @param debug Enable debugging mode (browser)
 #' @details The package contains long form documentation in the form of vignettes that cover the use of the main fucntions. Use browseVignettes(package="inctools") to access them.
 #'
@@ -51,6 +52,7 @@
 #'        recency_rule = "independent_thresholds",
 #'        recency_vars = c("Result","VL"),
 #'        recency_params = c(10,0,1000,1),
+#'        method = "exact",
 #'        alpha = 0.05)
 #' @export
 frrcal <- function(data = NULL,
@@ -61,6 +63,7 @@ frrcal <- function(data = NULL,
                    recency_vars = NULL,
                    recency_params = NULL,
                    alpha = 0.05,
+                   method = "exact",
                    debug = FALSE) {
 
   if (debug) {browser()}
@@ -101,6 +104,13 @@ frrcal <- function(data = NULL,
   if (is.null(time_var)) {stop("Error: No time variable provided.")}
   if (is.null(time_var)) {stop("Error: No recency variables provided variable provided.")}
 
+  if (is.null(method)) {stop("Error: Confidence interval method must be specified.")}
+  
+  if (length(method) != 1) {stop("Error: Exactly one confidence interval method must be specified.")}
+  
+  if ( !(method %in% c("exact", "ac", "asymptotic", "wilson", "prop.test", "bayes", "logit", "cloglog", "probit"))) {
+    stop("Confidence interval method must be one of 'exact', 'ac', 'asymptotic', 'wilson', 'prop.test', 'bayes', 'logit', 'cloglog', 'probit'. See help of binom::binom.test() for further details.")
+    }
 
   # check that subject id, time and recency variables exist
   variables <- colnames(data)
@@ -143,13 +153,13 @@ frrcal <- function(data = NULL,
     }
   }
   subjectdata <- subjectdata[!is.na(subjectdata$sid),]
-  nr <- ceiling(sum(subjectdata$recent))
+  nr <- as.integer(ceiling(sum(subjectdata$recent)))
   n <- nrow(subjectdata)
   p <- nr / n
   sigma <- sqrt( (p * (1 - p)) / n )
-  binom_test <- stats::binom.test(nr, n, p = 0, conf.level = 1 - alpha)
-  FRR <- data.frame(p,  sigma, binom_test$conf.int[1], binom_test$conf.int[2], alpha, nr, n, nrow(data))
-  colnames(FRR) <- c("FRRest", "SE", "LB", "UB", "alpha", "n_recent", "n_subjects", "n_observations")
-  rownames(FRR) <- ""
+  binom_ci <- binom::binom.confint(nr, n, conf.level = 1 - alpha, methods = method)
+  FRR <- tibble::tibble(FRRest = p, SE = sigma, LB = binom_ci$lower, UB = binom_ci$upper, 
+                        alpha = alpha, n_recent = nr, n_subjects = n, n_observations = nrow(data), 
+                        ci_method = method)
   return(FRR)
 }

inctools/man/frrcal.Rd

@@ -155,6 +155,7 @@ frr_untreated <- frrcal(data = untreated_all,
                         recency_rule = "independent_thresholds",
                         recency_vars = c("LAg_ODn","viral_load"),
                         recency_params = c(1.5,0,75,1),
+                        method = "exact",
                         debug = TRUE)
 print(frr_untreated)
 frr_treated <- frrcal(data = treated_all,