Merge pull request #348 from USEPA/343-mc5-ow-and-loec

343/345 add err for multiple ow methods, update mc5 and loec.coff method to ow_loec.coff and include_loec.coff methods

Author: cthunes

Diff for: R/mc5.R

@@ -78,10 +78,10 @@ mc5 <- function(ae, wr = FALSE) {
 
   ## Load cutoff methods
   ms <- tcplMthdLoad(lvl = 5L, id = ae, type = "mc")
-  if ('loec.coff' %in% ms$mthd) {
-    # using the loec method
+  if (any(c("ow_loec.coff", "include_loec.coff") %in% ms$mthd)) {
+    # using a loec method
     loec.mthd = TRUE
-    ms <- ms[!mthd=='loec.coff']
+    ms <- ms[!mthd=='include_loec.coff']
   }
 
   #special case where osd needs to be overwritten
@@ -93,6 +93,9 @@ mc5 <- function(ae, wr = FALSE) {
 
   ## Extract methods that need to overwrite hitc and hit_val
   ms_overwrite <- ms[grepl("ow_",mthd),]
+  if (nrow(ms_overwrite) > 1) {
+    stop(paste0("Only one level 5 hit-call override method may be assigned concurrently. Currently assigned: ", paste0(ms_overwrite$mthd, collapse = ", ")))
+  }
   ## Extract methods that don't overwrite
   ms <- ms[!grepl("ow_",mthd),]
 
@@ -124,7 +127,6 @@ mc5 <- function(ae, wr = FALSE) {
 
     ## Complete the loec calculations
     if (loec.mthd) {
-      
       all_resp_gt_conc <-function(resp) {
         # all resp > coff
         return(as.integer(all(abs(resp) > cutoff))) # All responses must be greater than coff
@@ -137,13 +139,18 @@ mc5 <- function(ae, wr = FALSE) {
       mc3[is.infinite(loec), loec := NA] #convert Inf to NA
       mc3[, loec_hitc := max(loec_coff), by = spid] # is there a loec? used for hitc
       mc3 <- mc3[dat, mult='first', on='spid', nomatch=0L]
-      
-      dat <- dat[mc3[,c("spid","loec", "loec_hitc")],on = "spid"]
-      
-    } else {
-      # if we're using v3 schema and not loec method we want to tcplfit2
-      dat <- tcplHit2(dat, coff = cutoff)
+      mc3 <- mc3[,c("m4id","aeid","coff","loec")] |> melt(measure.vars = c("loec"), variable.name = "hit_param", value.name = "hit_val")
+    }
+    
+    # if we're using v3 schema we want to tcplfit2
+    dat <- tcplHit2(dat, coff = cutoff)
+    
+    if (loec.mthd) {
+      exprs <- lapply(mthd_funcs["include_loec.coff"], do.call, args = list())
+      fenv <- environment()
+      invisible(rapply(exprs, eval, envir = fenv))
     }
+    
   } else {
     # Legacy fitting
 
@@ -381,30 +388,21 @@ mc5 <- function(ae, wr = FALSE) {
     dat <- dat[ , .SD, .SDcols = outcols]
   }
 
-  
-  # apply overwrite methods
+  # apply overwrite method
   if (nrow(ms_overwrite) > 0) {
     exprs <- lapply(mthd_funcs[ms_overwrite$mthd], do.call, args = list())
     fenv <- environment()
     invisible(rapply(exprs, eval, envir = fenv))
   }
 
-  # apply loec.coff 
-  if (loec.mthd) {
-    exprs <- lapply(mthd_funcs[c("loec.coff")], do.call, args = list())
-    fenv <- environment()
-    invisible(rapply(exprs, eval, envir = fenv))
-  }
-  
 
   ttime <- round(difftime(Sys.time(), stime, units = "sec"), 2)
   ttime <- paste(unclass(ttime), units(ttime))
   cat("Processed L5 AEID", ae, " (", nrow(dat), 
       " rows; ", ttime, ")\n", sep = "")
 
   res <- TRUE
-  
-  
+
 
   ## Load into mc5 table -- else return results
   if (wr) {

Diff for: R/mc5_mthds.R

@@ -106,17 +106,20 @@
 #'   \item{maxmed20pct}{Add a cutoff value of 20 percent of the maximum of all endpoint maximal 
 #'   average response values (max_med).}
 #'   \item{coff_2.32}{Add a cutoff value of 2.32.}
-#'   \item{loec.coff}{Identify the lowest observed effective concentration (loec) where the values 
-#'   of all responses are outside the cutoff band (i.e. abs(resp) > cutoff). If loec exists, assume 
-#'   hit call = 1, fitc = 100, model_type = 1. Winning model is not selected based on curve fits 
-#'   and therefore additional potency estimates are not derived.}
 #'   \item{ow_bidirectional_loss}{Multiply winning model hitcall (hitc) by -1 for models fit in the 
 #'   positive analysis direction. Typically used for endpoints where only negative responses are 
 #'   biologically relevant.}
 #'   \item{ow_bidirectional_gain}{Multiply winning model hitcall (hitc) by -1 for models fit in the 
 #'   negative analysis direction. Typically used for endpoints where only positive responses are 
 #'   biologically relevant.}
 #'   \item{osd_coff_bmr}{Overwrite the osd value so that bmr equals cutoff.}
+#'   \item{ow_loec.coff}{Identify the lowest observed effective concentration (loec) where the values 
+#'   of all responses are outside the cutoff band (i.e. abs(resp) > cutoff). loec is stored alongside 
+#'   winning model and potency estimates. If loec exists, assume hit call = 1, fitc = 100, 
+#'   model_type = 1, and if not, assume hit call = 0.}
+#'   \item{include_loec.coff}{Identify the lowest observed effective concentration (loec) where the 
+#'   values of all responses are outside the cutoff band (i.e. abs(resp) > cutoff). loec is 
+#'   stored alongside winning model and potency estimates.}
 #'  }
 #' }
 #'
@@ -354,12 +357,24 @@ mc5_mthds <- function(ae) {
 
     },
 
-    loec.coff = function() {
+    ow_loec.coff = function() {
 
-      e1 <- bquote(dat[, c("modl", "fitc", "model_type", "hitc") := list("loec", 100L, 1, loec_hitc)])
-      e2 <- bquote(dat <- dat |> melt(measure.vars = c("loec"), variable.name = "hit_param", value.name = "hit_val"))
-      e3 <- bquote(dat <- dat[,c("m4id", "aeid", "modl", "hitc", "fitc", "coff", "model_type", "hit_param", "hit_val")])
-      list(e1, e2, e3)
+      # get all endpoint sample m4ids where the loec param is not na
+      e1 <- bquote(loec.m4ids <- dat[(hit_param == "loec") & !is.na(hit_val), unique(m4id)])
+      # set hitcall and hitc param to 1 if found in m4id list and 0 if not
+      e2 <- bquote(dat[hit_param == "hitcall", hit_val:=ifelse(m4id %in% loec.m4ids, 1, 0)])
+      e3 <- bquote(dat[, hitc:=ifelse(m4id %in% loec.m4ids, 1, 0)])
+      # update modl to loec, fitc to 100, model_type to 1
+      e4 <- bquote(dat[, c("modl", "fitc", "model_type") := list("loec", 100L, 1)])
+      list(e1, e2, e3, e4)
+      
+    },
+    
+    include_loec.coff = function() {
+      
+      e1 <- bquote(dat <- rbind(dat, mc3, fill = TRUE) |> arrange(m4id))
+      e2 <- bquote(dat <- dat |> group_by(m4id) |> tidyr::fill(modl,hitc,fitc,model_type) |> as.data.table())
+      list(e1, e2)
 
     }