Population correction

DESCRIPTION

@@ -36,7 +36,6 @@ Suggests:
   parallel,
   data.table
 Encoding: UTF-8
-RoxygenNote: 7.3.3
 RdMacros: Rdpack
 Depends: 
     R (>= 3.5.0)
@@ -45,3 +44,4 @@ BugReports: https://github.com/SymbolixAU/collision/issues
 LazyData: true
 VignetteBuilder: knitr
 Roxygen: list(markdown = TRUE)
+Config/roxygen2/version: 8.0.0

NEWS.md

@@ -2,6 +2,7 @@
 
 # collision 1.1 (in development)
 
+* [issue24](https://github.com/SymbolixAU/collision/issues/24) added population correction as an optional parameter in `n_collision` function.
 * [issue46a](https://github.com/SymbolixAU/collision/issues/46a) added spatial correction to account for turbine clustering when turbine spacing is smaller than the effective survey width/area.
 * [issue40](https://github.com/SymbolixAU/collision/issues/40) improved handling of height of observed flux window to align with the max height of the turbine.
 * [issue41](https://github.com/SymbolixAU/collision/issues/41) bug fix to ensure the proportion of time a turbine is operational is accounted for in the `p_collision_dynamic` function.

R/n_collisions.R

@@ -4,18 +4,22 @@
 #' avoidance rate of static components (tower + static presented area)
 #' @param avoidance_rate_dynamic numeric; Number between 0 and 1 representing the
 #' avoidance rate of dynamic components (moving blade edge)
-#' @param n_flights numeric; Flights though the turbine per unit time. Calculated from 
+#' @param n_flights numeric; Flights though the turbine per unit time. Calculated from
 #'                [turbine_flights_year()] or similar.
-#' @param p_coll_static numeric; the probability of collision with static 
-#'                      turbine components if an interaction occurs and 
-#'                      avoidance is zero. Calculated from 
+#' @param p_coll_static numeric; the probability of collision with static
+#'                      turbine components if an interaction occurs and
+#'                      avoidance is zero. Calculated from
 #'                      [prob_collision_static()]
-#' @param p_coll_dynamic numeric; the probability of collision with dynamic 
-#'                       blade edge turbine components if an interaction occurs 
-#'                       and avoidance is zero. Calculated from 
+#' @param p_coll_dynamic numeric; the probability of collision with dynamic
+#'                       blade edge turbine components if an interaction occurs
+#'                       and avoidance is zero. Calculated from
 #'                       [prob_collision_dynamic()]
-#' 
-#' @return numeric; number of collisions per unit time. Time interval 
+#' @param population numeric or NULL; the regional population size available to
+#'                    replace collided birds. When `NULL` (default) the function
+#'                    returns the standard infinite-sink estimate. When supplied, a
+#'                    finite-population correction is applied.
+#'
+#' @return numeric; number of collisions per unit time. Time interval
 #'    is the same as referenced by `n_flights` input.
 #'
 #' @examples
@@ -27,25 +31,42 @@
 #'   p_coll_dynamic = 0.03
 #' )
 #'
+#' # With finite-population correction
+#' n_collision(
+#'   avoidance_rate_static = 0.99,
+#'   avoidance_rate_dynamic = 0.90,
+#'   n_flights = 100,
+#'   p_coll_static = 0.05,
+#'   p_coll_dynamic = 0.03,
+#'   population = 500
+#' )
+#'
 #' @export
-n_collision <- function(
-    avoidance_rate_static,
-    avoidance_rate_dynamic,
-    n_flights,
-    p_coll_static,
-    p_coll_dynamic) {
-
+n_collision <- function(avoidance_rate_static,
+                        avoidance_rate_dynamic,
+                        n_flights,
+                        p_coll_static,
+                        p_coll_dynamic,
+                        population = NULL) {
   check_num_bounds(avoidance_rate_static, min = 0, max = 1)
   check_num_bounds(avoidance_rate_dynamic, min = 0, max = 1)
   check_num_bounds(n_flights, min = 0)
   check_num_bounds(p_coll_static, min = 0, max = 1)
   check_num_bounds(p_coll_dynamic, min = 0, max = 1)
-
-
-  return(
-    n_flights *(
-      (1 - avoidance_rate_static) * p_coll_static +
-        (1 - avoidance_rate_dynamic) * p_coll_dynamic
-    )
+
+  n_risk <- n_flights * (
+    (1 - avoidance_rate_static) * p_coll_static +
+      (1 - avoidance_rate_dynamic) * p_coll_dynamic
   )
+
+  if (is.null(population)) {
+    return(n_risk)
+  }
+
+  check_num_bounds(population, min = 0)
+
+  S <- n_risk / n_flights
+  m_ind <- n_flights / population
+
+  return(population * (1 - (1 - S)^m_ind))
 }

R/populationcorrection.R

@@ -0,0 +1,12 @@
+#' Optional population correction
+#' 
+#' The default calculation predicts the number of flight movements that will
+#' end in collision. This has an inherent assumption that the number of birds
+#' generating the flights will be automatically replenished following for
+#' collision. For small local populations this will be overly conservative.  
+#' This function converts the number of movements ending in collision to
+#' the number of birds, by accounting for the population at risk
+#' 
+#' 
+#' 
+#' @export

inst/tinytest/test_n_collision.R

@@ -121,3 +121,86 @@ expect_error(
   "variable out of bounds"
 )
 
+# --- population (finite-population correction) ---
+
+# NULL population returns the standard result
+expect_equal(
+  n_collision(
+    avoidance_rate_static = 0.99,
+    avoidance_rate_dynamic = 0.90,
+    n_flights = 100,
+    p_coll_static = 0.05,
+    p_coll_dynamic = 0.03,
+    population = NULL
+  ),
+  0.35,
+  tol = 1e-6
+)
+
+# finite population gives a corrected count
+expect_equal(
+  n_collision(
+    avoidance_rate_static = 0.99,
+    avoidance_rate_dynamic = 0.90,
+    n_flights = 100,
+    p_coll_static = 0.05,
+    p_coll_dynamic = 0.03,
+    population = 500
+  ),
+  500 * (1 - (1 - 0.0035)^0.2),
+  tol = 1e-6
+)
+
+# vectorised avoidance_dynamic works with population
+expect_equal(
+  n_collision(
+    avoidance_rate_static = 0.99,
+    avoidance_rate_dynamic = c(0.90, 0.95),
+    n_flights = 100,
+    p_coll_static = 0.05,
+    p_coll_dynamic = 0.03,
+    population = 500
+  ),
+  c(500 * (1 - (1 - 0.0035)^0.2), 500 * (1 - (1 - 0.002)^0.2)),
+  tol = 1e-6
+)
+
+# very large population converges to the infinite sink result
+expect_equal(
+  n_collision(
+    avoidance_rate_static = 0.99,
+    avoidance_rate_dynamic = 0.90,
+    n_flights = 100,
+    p_coll_static = 0.05,
+    p_coll_dynamic = 0.03,
+    population = 1e8
+  ),
+  0.35,
+  tol = 1e-2
+)
+
+# negative population is rejected
+expect_error(
+  n_collision(
+    avoidance_rate_static = 0.99,
+    avoidance_rate_dynamic = 0.90,
+    n_flights = 100,
+    p_coll_static = 0.05,
+    p_coll_dynamic = 0.03,
+    population = -100
+  ),
+  "variable out of bounds"
+)
+
+# non-numeric population is rejected
+expect_error(
+  n_collision(
+    avoidance_rate_static = 0.99,
+    avoidance_rate_dynamic = 0.90,
+    n_flights = 100,
+    p_coll_static = 0.05,
+    p_coll_dynamic = 0.03,
+    population = "500"
+  ),
+  "Numeric input expected"
+)