Skip to content

Commit

Permalink
DEV: Fix log_one_minus_phi
Browse files Browse the repository at this point in the history
- Moved to utils.R
- Delete duplicate definition in detection.R
- Fix namespacing
  • Loading branch information
fredjaya committed Jun 28, 2024
1 parent f459cb5 commit dcc028b
Show file tree
Hide file tree
Showing 3 changed files with 16 additions and 30 deletions.
18 changes: 1 addition & 17 deletions R/detection.R
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
}
}

14 changes: 1 addition & 13 deletions R/fisher_information.R
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
14 changes: 14 additions & 0 deletions R/util.R
Original file line number Diff line number Diff line change
Expand Up @@ -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){
Expand Down

0 comments on commit dcc028b

Please sign in to comment.