implemented analytical CV estimation

DESCRIPTION

@@ -1,15 +1,14 @@
 Package: ObsCovgTools
 Type: Package
 Title: Evaluate Fishery Observer Coverage for Bycatch Estimation
-Version: 3.2.1-2
-Date: 2022-09-13
+Version: 3.3.0
+Date: 2022-09-14
 URL: https://kacurtis.github.io/ObsCovgTools, https://github.com/kacurtis/ObsCovgTools
 Depends: R (>= 3.4.0)
 Imports: 
 	dplyr,
 	magrittr,
 	rlang,
-	Runuran,
 	tibble
 Suggests: 
   shiny

NAMESPACE

@@ -1,10 +1,9 @@
 # Generated by roxygen2: do not edit by hand
 
-export(plot_cv_obscov)
+export(plot_cv)
 export(plot_probposobs)
 export(plot_uclnegobs)
 export(probnzeros)
 export(run_shiny)
-export(sim_cv_obscov)
 importFrom(magrittr,"%>%")
 importFrom(rlang,.data)

NEWS.md

@@ -1,3 +1,8 @@
+# ObsCovgTools 3.3.0
+
+* Implemented an analytical solution for CV corresponding to observer coverage, combining sim_cv_obscov.r and plot_cv_obscov.r into one function analogous to those for the other two management objectives. Elimination of Monte Carlo simulations provides major efficiency gains and increases precision of result. 
+
+
 # ObsCovgTools 3.2.1-2
 
 * Added a `NEWS.md` file to track changes to the package.

R/plot_cv.r

@@ -0,0 +1,108 @@
+#' Plot bycatch estimation CV vs. observer coverage
+#' 
+#' \code{plot_cv} plots projected bycatch estimation CVs vs observer 
+#' coverage, and returns minimum observer coverage needed to achieve 
+#' user-specified target CV and percentile. 
+#'   
+#' @param te an integer greater than 1. Total effort in fishery (e.g., trips 
+#'   or sets).
+#' @param bpue a positive number. Bycatch per unit effort.
+#' @param d a number greater than or equal to 1. Dispersion index. The dispersion 
+#'   index corresponds to the variance-to-mean ratio of effort-unit-level bycatch, 
+#'   so \code{d = 1} corresponds to Poisson-distributed bycatch, and \code{d > 1} 
+#'   to overdispersed bycatch.
+#' @param targetcv a non-negative number less than 1. Target CV (as a proportion). 
+#'   If set to 0, no corresponding minimum observer coverage will be 
+#'   highlighted or returned.
+#' @param silent logical. If \code{TRUE}, print output to terminal is suppressed.
+#' @param showplot logical. If \code{FALSE}, plotting is suppressed.
+#' @param ... additional arguments for compatibility with Shiny. 
+#' 
+#' @details
+#' \strong{Caveat:} \code{plot_cv} assumes that (1) observer coverage is 
+#' representative, (2) bycatch (\code{bpue}) is in terms of individuals (not 
+#' weight) per unit effort, and (3) the specified dispersion index reflects 
+#' the highest level of any hierarchical variance (e.g., using dispersion index 
+#' at trip level if greater than that at set level). Violating these assumptions 
+#' will likely result in negatively biased projections of bycatch estimation CV 
+#' for a given level of observer coverage. More conservative projections can be 
+#' obtained by using a higher dispersion index \code{d}. Users may want to 
+#' explore uncertainty in dispersion index and in bycatch per unit effort by 
+#' varying those inputs.
+#'   
+#' @return If \code{targetcv} is non-zero, a list with one component:
+#'   \item{targetoc}{minimum observer coverage in terms of percentage.} 
+#' @return Returned invisibly. 
+#'   
+#' @export 
+plot_cv <- function(te, bpue, d = 2, targetcv = 0.3, 
+                           showplot = TRUE, silent = FALSE, ...) {
+
+  # check input values
+  if ((ceiling(te) != floor(te)) || te<2) stop("te must be a positive integer > 1")
+  if (bpue<=0) stop("bpue must be > 0")
+  if (d<1) stop("d must be >= 1")
+  if(targetcv<0 || targetcv>=1) stop("targetcv must be >= 0 and < 1")
+
+  # get shiny flag if specified or set to FALSE
+  myArgs <- match.call()
+  if (!("as.shiny" %in% names(myArgs))) as.shiny <- FALSE
+  else as.shiny <- myArgs$as.shiny
+
+  # get CV of total bycatch estimates for range of observer coverage
+  if (te<1000) { oc <- 1:te 
+  } else { oc <- round(seq(0.001,1,0.001)*te) }
+  df <- data.frame(nobs = oc, pobs = oc/te)
+  fpc <- (te-df$nobs)/(te-1)   # finite population correction, Cochran 1973
+  s2 <- d*bpue   # variance of NB or Poisson distribution
+  df$cv <- sqrt((fpc * s2)/df$nobs) / bpue
+
+  # get minimum required observer coverage if targetcv provided
+  if (targetcv) {
+    itarget <- min(which(df$cv <= targetcv))
+    targetoc <- 100*ceiling_dec(df$pobs[itarget], 3)
+  }
+
+  # plot 
+  if (showplot) {
+    with(df, graphics::plot(100*pobs, cv, type="l", lty=1, lwd=2,
+                            xlim=c(0,100), ylim=c(0,1), 
+                            xaxs="i", yaxs="i", xaxp=c(0,100,10), yaxp=c(0,1,10),
+                            xlab="Observer Coverage (%)", ylab="CV of Bycatch Estimate",
+                            main="CV of Bycatch Estimate vs Observer Coverage"))
+    with(df, graphics::lines(100*pobs, cv))
+    graphics::abline(h=1, v=100)
+    # add minimum required observer coverage
+    if (targetcv) {
+      graphics::abline(h=targetcv, col=2, lwd=2, lty=2)
+      graphics::par(xpd=TRUE)
+      graphics::points(targetoc, targetcv, pch=8, col=2, cex=1.5, lwd=2)
+      graphics::par(xpd=FALSE)
+      legpos <- ifelse(any(df$pobs > 0.7 & df$cv > 0.5), 
+                       "bottomleft", "topright")
+      graphics::legend(legpos, lty=c(2,0), pch=c(NA,8), col=c(2,2), lwd=c(2,2), 
+                       pt.cex=1.5, y.intersp=1.1, legend=c("target CV", "min coverage"))
+    }
+  }
+
+  # print recommended minimum observer coverage
+  if (targetcv)  {
+    if (!is.na(targetoc)) {
+      rec <- paste0("Minimum observer coverage to achieve CV \u2264 ", targetcv, 
+                     " is ", format(targetoc, nsmall=1), "%.\n")
+    } else {
+      rec <- paste0("Simulated observer coverage levels do not include range corresponding to ",
+                 "minimum observer coverage to achieve CV \u2264 ", targetcv, ".\n")
+    }
+  } else { rec <- "" }
+  if (!as.shiny & !silent) 
+    cat(paste0(rec, "Please review the caveats in the documentation.\n"))
+
+  # return recommended minimum observer coverage
+  if (as.shiny) {
+    return(invisible(list(rec = rec)))
+  } else { 
+    if (targetcv)
+      return(invisible(list(targetoc = targetoc)))
+  }
+}

R/plot_probposobs.r

@@ -3,7 +3,9 @@
 #' \code{plot_probposobs} plots (1) probability of observing at least one bycatch
 #'   event vs observer coverage and (2) probability of any bycatch occurring in 
 #'   total fishery effort, given total fishery effort, bycatch per unit effort, 
-#'   and dispersion index. 
+#'   and dispersion index. The function returns returns minimum observer coverage 
+#'   needed to achieve user-specified probability of observing bycatch if it 
+#'   occurs.
 #'   
 #' @param te an integer greater than 1. Total effort in fishery (e.g., trips 
 #'   or sets).
@@ -63,13 +65,14 @@ plot_probposobs <- function(te, bpue, d = 2, targetppos = 95, showplot = TRUE,
   if (!("as.shiny" %in% names(myArgs))) as.shiny <- FALSE
   else as.shiny <- myArgs$as.shiny
 
-  # percent probability of positive observed bycatch
+  # get probability of positive observed bycatch for range of observer coverage
   if (te<1000) { oc <- 1:te 
   } else { oc <- round(seq(0.001,1,0.001)*te) }
   df <- data.frame(nobs = oc, pobs = oc/te)
   df$pp <- 1-probnzeros(df$nobs, bpue, d)   # probability of positive observed bycatch
   ppt <- utils::tail(df$pp,1)   # probability of positive bycatch in total effort
   df$ppc <- df$pp/ppt
+  # get minimum required observer coverage if targetppos provided
   if (targetppos) {
     itarget <- min(which(df$ppc >= targetppos/100))
     targetoc <- 100*ceiling_dec(df$pobs[itarget], 3)
@@ -80,7 +83,8 @@ plot_probposobs <- function(te, bpue, d = 2, targetppos = 95, showplot = TRUE,
     opar <- graphics::par(no.readonly = TRUE)
     graphics::par(xpd=TRUE)
     graphics::plot(100*df$pobs, 100*(df$ppc), type="l", lty=1, lwd=2,
-                   xlim=c(0,100), ylim=c(0,100), xaxs="i", yaxs="i", xaxp=c(0,100,10), yaxp=c(0,100,10),
+                   xlim=c(0,100), ylim=c(0,100), 
+                   xaxs="i", yaxs="i", xaxp=c(0,100,10), yaxp=c(0,100,10),
                    xlab="Observer Coverage (%)", ylab="Probability of Positive Bycatch (%)",
                    main="Probability of Positive Bycatch")
     graphics::lines(x=c(0,100),y=rep(100*ppt,2),lwd=3, lty=3)

R/plot_uclnegobs.r

@@ -2,7 +2,10 @@
 #' 
 #' \code{plot_uclnegobs} plots upper confidence limit of total bycatch vs 
 #'   observer coverage when no bycatch is observed, given total fishery effort, 
-#'   dispersion index, and confidence level.
+#'   dispersion index, and confidence level. The function returns (1) minimum 
+#'   observer coverage needed to fall within user-specified upper confidence 
+#'   limit for bycatch when none was observed, and/or (2) the upper confidence 
+#'   limit for bycatch given specified observer coverage and no observed bycatch. 
 #'   
 #' @param te an integer greater than 1. Total effort in fishery (e.g., trips 
 #'   or sets).
@@ -91,7 +94,7 @@ plot_uclnegobs <- function(te, d = 2, cl = 95, targetucl = 0, fixedoc = 0,
   df$ucl <- df$fpc * te * solveucl(a=a, d=dv[2], n=df$nobs)
   ucl.dh <- df$fpc * te * solveucl(a=a, d=dv[3], n=df$nobs)
 
-  # identify target observer coverage if target ucl specified
+  # get minimum required observer coverage if target ucl specified
   if (targetucl) {
     itarget <- min(which(df$ucl <= targetucl))
     targetoc <- 100*ceiling_dec(df$pobs[itarget], 3)