From dcc028b1c6d5d11118a046ac86b727308326d9c2 Mon Sep 17 00:00:00 2001
From: Fred Jaya <36436914+fredjaya@users.noreply.github.com>
Date: Fri, 28 Jun 2024 12:35:41 +1000
Subject: [PATCH] DEV: Fix log_one_minus_phi

- Moved to utils.R
- Delete duplicate definition in detection.R
- Fix namespacing
---
 R/detection.R          | 18 +-----------------
 R/fisher_information.R | 14 +-------------
 R/util.R               | 14 ++++++++++++++
 3 files changed, 16 insertions(+), 30 deletions(-)

diff --git a/R/detection.R b/R/detection.R
index 749cbaa..4507d82 100644
--- a/R/detection.R
+++ b/R/detection.R
@@ -98,7 +98,7 @@ detection_errors_cluster <- function(pool_size, pool_number, cluster_number,
       #stats::dnorm(x, mean = mu, sd = sigma) * (1 - (1 - sensitivity - specificity) * (1-invlink(x))^pool_size - sensitivity)^pool_number
       #equivalent to the above, but more numerically stable
       exp(stats::dnorm(x, mean = mu, sd = sigma, log = TRUE) +
-            log_1_minus_phi(invlink(x), pool_size, sensitivity, specificity) * pool_number)
+            log_one_minus_phi(invlink(x), pool_size, sensitivity, specificity) * pool_number)
     }
     
     typeII <- stats::integrate(f, -Inf, Inf)$value ^ cluster_number
@@ -215,19 +215,3 @@ detection_errors_cluster <- function(pool_size, pool_number, cluster_number,
 #   }
 #   list(typeI = typeI, typeII = typeII)
 # }
-
-
-# Helper function giving log of 1 minus the pool positivity probability for a
-# given prevalence (p), pool size (s), sensitivity (spec), and specificity
-# (spec). Not to be exported
-log_one_minus_phi <- function(p, s, sens, spec){
-  # log(1 - sens - (1 - spec - sens) * (1 - p)^s)
-  if(p %in% 0:1){
-    log1p(- (1 - (1 - p)^s) * sens - (1 - p)^s * (1 - spec))
-  }else if(sens == 1){
-    log1p(-p) * s + log(spec)
-  }else{
-    log(expm1(log1p(-p) * s) * (sens - 1) + exp(log1p(-p) * s) * spec)
-  }
-}
-
diff --git a/R/fisher_information.R b/R/fisher_information.R
index 8ce867d..1f05a32 100644
--- a/R/fisher_information.R
+++ b/R/fisher_information.R
@@ -149,18 +149,6 @@ fi_pool_cluster <- function(pool_size,
     }
   }
   
-  log_one_minus_phi <- function(p){
-    # log(1 - varphi - (1 - psi - varphi) * (1 - p)^s)
-    if(p %in% 0:1){
-      log1p(- (1 - (1 - p)^s) * varphi - (1 - p)^s * (1 - psi))
-    }else if(varphi == 1){
-      log1p(-p) * s + log(psi)
-    }else{
-      log(expm1(log1p(-p) * s) * (varphi - 1) + exp(log1p(-p) * s) * psi)
-    }
-  }
-  
-
   if (form %in% c("discrete", "beta")) {
     # In these cases Fisher information matrix is calculated directly in terms
     # of theta and rho and checks on sums of likelihoods and likelihood
@@ -429,7 +417,7 @@ fi_pool_cluster <- function(pool_size,
         }else{
           out[j] <- stats::dnorm(zj, mean = mu, sd = sigma, log = TRUE) +
             sum(log(phi(pj)) * y) +
-            sum(log_one_minus_phi(pj) * (N - y))
+            sum(log_one_minus_phi(pj, s, varphi, psi) * (N - y))
         }
       }
       out <- exp(out)
diff --git a/R/util.R b/R/util.R
index 5a20b4a..4b3674f 100644
--- a/R/util.R
+++ b/R/util.R
@@ -8,6 +8,20 @@ cloglog_inv <- function(x){
   - expm1(-exp(x))
 }
 
+# Helper function giving log of 1 minus the pool positivity probability for a
+# given prevalence (p), pool size (s), sensitivity (spec), and specificity
+# (spec). Not to be exported
+log_one_minus_phi <- function(p, s, sens, spec){
+  # log(1 - sens - (1 - spec - sens) * (1 - p)^s)
+  if(p %in% 0:1){
+    log1p(- (1 - (1 - p)^s) * sens - (1 - p)^s * (1 - spec))
+  }else if(sens == 1){
+    log1p(-p) * s + log(spec)
+  }else{
+    log(expm1(log1p(-p) * s) * (sens - 1) + exp(log1p(-p) * s) * spec)
+  }
+}
+
 # Calculating parameters for link-normal style distributions
 
 mu_sigma_linknorm <- function(.mean, .var,link,invlink){