added the possibility to calculate a variable IC value

OpenStats/R/Backend_History.R

@@ -146,6 +146,7 @@ eval_entry_V1_2 <- function(entry, DataModelState,
     DoseResponse = {
       res <- dose_response_V1_2$new(
         DataModelState$df,
+        entry[["IC [%]"]],
         entry[["Log transform x-axis"]],
         entry[["Log transform y-axis"]],
         entry[["Column containing the names"]],

OpenStats/R/Backend_LC50.R

@@ -79,7 +79,7 @@ false_discovery_rate <- function(res) {
 }
 env_lc_V1_2$false_discovery_rate <- false_discovery_rate
 
-check_fit <- function(model, min_conc, max_conc,
+check_fit <- function(model, ic_percentage, min_conc, max_conc,
                       min_abs, max_abs, substance_name, unit) {
   if (model$fit$convergence != TRUE) {
     return(errorClass$new(paste(
@@ -93,57 +93,68 @@ check_fit <- function(model, min_conc, max_conc,
   e <- coefficients(model)[4] # IC50
   RSE <- summary(model)$rseMat[1] # residual standard error estimated
   Response_lowestdose_predicted <- predict(
-    model, data.frame(concentration = min_conc),
+    model, data.frame(conc = min_conc),
     se.fit = FALSE
   )[1]
   Response_highestdose_predicted <- predict(
-    model, data.frame(concentration = max_conc),
+    model, data.frame(conc = max_conc),
     se.fit = FALSE
   )[1]
   Response_difference <- abs(
     Response_lowestdose_predicted - Response_highestdose_predicted
   )
   HillCoefficient <- b
-  IC50_relative <- e
+  ed_res <- drc::ED(
+    model,
+    respLev = ic_percentage,
+    interval = "delta",
+    level = 0.95,
+    type = "relative",
+    display = FALSE
+  )
+  IC_relative <- ed_res[1, 1]
   Problems <- ""
   if (Response_difference < 0.25) {
     Problems <- paste(Problems,
-      "Response Difference lower than 25%",
-      collapse = " , "
+      "Response Difference lower than 25%", collapse = " , "
     )
-  } else if (IC50_relative > max_conc) {
+  } else if (IC_relative > max_conc) {
     Problems <- paste(Problems,
-      "IC50 larger than highest measured concentration",
-      collapse = " , "
+      "IC larger than highest measured concentration", collapse = " , "
     )
-  } else if (IC50_relative < min_conc) {
+  } else if (IC_relative < min_conc) {
     Problems <- paste(Problems,
-      "IC50 lower than lowest measured concentration",
-      collapse = " , "
+      "IC lower than lowest measured concentration", collapse = " , "
     )
   }
-  confidence_interval <- confint(model, parm = c("e"), level = 0.95)
-  IC50_relative_lower <- confidence_interval[1]
-  IC50_relative_higher <- confidence_interval[2]
+
+  IC_relative_lower <- ed_res[1, 3]
+  IC_relative_higher <- ed_res[1, 4]
   p_value <- drc::noEffect(model)[3]
   Response_lowestdose_predicted <- env_lc_V1_2$shapenumber(Response_lowestdose_predicted)
   Response_highestdose_predicted <- env_lc_V1_2$shapenumber(Response_highestdose_predicted)
   HillCoefficient <- env_lc_V1_2$shapenumber(HillCoefficient)
-  IC50_relative <- env_lc_V1_2$shapenumber(IC50_relative)
-  IC50_relative_lower <- env_lc_V1_2$shapenumber(IC50_relative_lower)
-  IC50_relative_higher <- env_lc_V1_2$shapenumber(IC50_relative_higher)
-  pIC50 <- env_lc_V1_2$shapenumber(-log10(IC50_relative))
+  IC_relative <- env_lc_V1_2$shapenumber(IC_relative)
+  IC_relative_lower <- env_lc_V1_2$shapenumber(IC_relative_lower)
+  IC_relative_higher <- env_lc_V1_2$shapenumber(IC_relative_higher)
+  pIC <- env_lc_V1_2$shapenumber(-log10(IC_relative))
   p_value <- env_lc_V1_2$shapenumber(p_value)
   outvar <- data.frame(
     name = substance_name,
     Response_lowestdose_predicted = Response_lowestdose_predicted,
     Response_highestdose_predicted = Response_highestdose_predicted,
     HillCoefficient = HillCoefficient,
     asymptote_one = c, asymptote_two = d,
-    IC50_relative = IC50_relative, IC50_relative_lower = IC50_relative_lower,
-    IC50_relative_higher = IC50_relative_higher, unit = unit, pIC50 = pIC50,
+    IC_relative = IC_relative, IC_relative_lower = IC_relative_lower, IC_relative_higher = IC_relative_higher,
+    unit = unit, pIC = pIC,
     RSE = RSE, p_value = p_value, Problems = Problems
   )
+  names(outvar)[7:9] <- c(
+    paste0("IC_", ic_percentage, "_relative"),
+    paste0("IC_", ic_percentage, "_relative_lower"),
+    paste0("IC_", ic_percentage, "_relative_higher")
+  )
+
   return(outvar)
 }
 env_lc_V1_2$check_fit <- check_fit
@@ -169,7 +180,7 @@ drawplot_only_raw_data <- function(df, abs_col, conc_col, title, unit) {
 env_lc_V1_2$drawplot_only_raw_data <- drawplot_only_raw_data
 
 drawplot <- function(df, abs_col, conc_col, unit, model, valid_points, title,
-                     IC50_relative, IC50_relative_lower, IC50_relative_higher,
+                     IC_relative, IC_relative_lower, IC_relative_higher,
                      islog_x, islog_y) {
 
   conc <- function() stop("Should never be called") # Please R CMD check
@@ -188,8 +199,8 @@ drawplot <- function(df, abs_col, conc_col, unit, model, valid_points, title,
   max_conc <- max(df[, conc_col]) +
     0.1 * (max(df[, conc_col]) - min(df[, conc_col]))
   min_conc <- min(df[, conc_col]) - 0.1 * min(df[, conc_col])
-  xmin <- IC50_relative_lower
-  xmax <- IC50_relative_higher
+  xmin <- IC_relative_lower
+  xmax <- IC_relative_higher
   if (!is.na(xmin) && !is.na(xmax)) {
     ymin <- min(df[, abs_col])
     ymax <- max(df[, abs_col])
@@ -228,7 +239,7 @@ drawplot <- function(df, abs_col, conc_col, unit, model, valid_points, title,
 }
 env_lc_V1_2$drawplot <- drawplot
 
-ic50_internal <- function(df, abs, conc,
+ic_internal <- function(df, ic_percentage, abs, conc,
                           title, islog_x, islog_y, unit) {
   model <- drc::drm(abs ~ conc,
     data = df, fct = drc::LL.4(),
@@ -242,25 +253,25 @@ ic50_internal <- function(df, abs, conc,
     fct = drc::LL.4(), robust = "mean",
   )
   res <- env_lc_V1_2$check_fit(
-    model, min(df[, conc]),
+    model, ic_percentage, min(df[, conc]),
     max(df[, conc]), min(df[, abs]), max(df[, abs]), title, unit
   )
   p <- env_lc_V1_2$drawplot(
-    df, abs, conc, unit, model, valid_points, title, res$IC50_relative,
-    res$IC50_relative_lower, res$IC50_relative_higher,
+    df, abs, conc, unit, model, valid_points, title,
+    res[[7]], res[[8]], res[[9]],
     islog_x, islog_y
   )
   return(list(res, p))
 }
-env_lc_V1_2$ic50_internal <- ic50_internal
+env_lc_V1_2$ic_internal <- ic_internal
 
 check_dr_df <- function(df, abs_col,
                         conc_col, substance_name_col) {
   if (!is.character(df[, substance_name_col]) &&
     !is.factor(df[, substance_name_col])) {
     return(errorClass$new("The substance names are not character"))
   }
-  substances <- unique(df[, substance_name_col]) # TODO: is this even possible?
+  substances <- unique(df[, substance_name_col]) # is this even possible?
   if (length(substances) < 1) {
     return(errorClass$new("The data for compounds seems to be missing"))
   }
@@ -279,13 +290,13 @@ transform_conc_dr <- function(conc_col) {
     ))
   }
   if (!is.numeric(temp_conc)) {
-    return(errorClass$new("The concentration data is not numerical")) # TODO: is this even possible?
+    return(errorClass$new("The concentration data is not numerical")) # is this even possible?
   }
   return(temp_conc)
 }
 env_lc_V1_2$transform_conc_dr <- transform_conc_dr
 
-ic50 <- function(df, abs_col, conc_col,
+ic <- function(df, ic_percentage, abs_col, conc_col,
                  substance_name_col, unit_col,
                  islog_x, islog_y) {
   # Checks
@@ -314,8 +325,8 @@ ic50 <- function(df, abs_col, conc_col,
 
     m <- tryCatch(
       {
-        m <- env_lc_V1_2$ic50_internal(
-          df_temp,
+        m <- env_lc_V1_2$ic_internal(
+          df_temp, ic_percentage,
           "abs", "conc",
           substances[i],
           islog_x, islog_y, unit
@@ -335,4 +346,4 @@ ic50 <- function(df, abs_col, conc_col,
   }
   return(res)
 }
-env_lc_V1_2$ic50 <- ic50
+env_lc_V1_2$ic <- ic

OpenStats/R/Backend_V1_2_Engine.R

@@ -1325,12 +1325,14 @@ dose_response_V1_2 <- R6::R6Class(
     com = NULL,
     res_df = NULL,
     res_p = NULL,
+    ic_percentage = NULL,
 
-    initialize = function(df,
+    initialize = function(df, ic_percentage,
                           is_xlog, is_ylog,
                           substance_names, unit_names,
                           formula, com = communicator_V1_2) {
       self$df <- df
+      self$ic_percentage <- ic_percentage
       self$df[, substance_names] <- self$df[, substance_names] |> as.character()
       self$is_xlog <- is_xlog
       self$is_ylog <- is_ylog
@@ -1347,7 +1349,7 @@ dose_response_V1_2 <- R6::R6Class(
         expr = {
           promise_history_entry <- self$create_history(new_name)
           ResultsState$bgp$start(
-            fun = function(df, formula, substance_names, unit_names, is_xlog, is_ylog) {
+            fun = function(df, ic_percentage, formula, substance_names, unit_names, is_xlog, is_ylog) {
               f <- as.character(formula)
               dep <- f[2]
               indep <- f[3]
@@ -1358,9 +1360,8 @@ dose_response_V1_2 <- R6::R6Class(
                 OpenStats:::env_check_ast_V1_2$check_ast(str2lang(indep), colnames(df))
                 OpenStats:::env_check_ast_V1_2$check_ast(str2lang(dep), colnames(df))
 
-                res <- OpenStats:::env_lc_V1_2$ic50(
-                  df, dep,
-                  indep, substance_names, unit_names,
+                res <- OpenStats:::env_lc_V1_2$ic(
+                  df, ic_percentage, dep, indep, substance_names, unit_names,
                   is_xlog, is_ylog
                 )
                 if (inherits(res, "errorClass")) {
@@ -1393,7 +1394,7 @@ dose_response_V1_2 <- R6::R6Class(
                 current_page = 1L)
             },
             args = list(
-              df = self$df, formula = self$formula, substance_names = self$substance_names,
+              df = self$df, ic_percentage = self$ic_percentage, formula = self$formula, substance_names = self$substance_names,
               unit_names = self$unit_names, is_xlog = self$is_xlog, is_ylog = self$is_ylog
             ),
             promise_result_name = new_name,
@@ -1411,6 +1412,7 @@ dose_response_V1_2 <- R6::R6Class(
     create_history = function(new_name) {
       list(
         type = "DoseResponse",
+        "IC [%]" = self$ic_percentage,
         "Column containing the names" = self$substance_names,
         "Column containing the units" = self$unit_names,
         "Log transform x-axis" = self$is_xlog,

OpenStats/R/Server_DoseResponse.R

@@ -80,7 +80,7 @@ DoseResponseServer <- function(id, DataModelState, ResultsState) {
 
     run_dr <- function(df, new_name) {
       dr <- get_dose_response()$new(
-        df, input$xTransform, input$yTransform,
+        df, input$ic_percentage, input$xTransform, input$yTransform,
         input$substanceNames, input$unitNames, DataModelState$formula
       )
       dr$eval(ResultsState, new_name)

OpenStats/R/UI_DoseResponse.R

@@ -3,6 +3,9 @@ DoseResponseSidebarUI <- function(id) {
     "Dose Response analysis",
     shiny::uiOutput(shiny::NS(id, "substanceNamesUI")),
     shiny::uiOutput(shiny::NS(id, "unitNamesUI")),
+    shiny::sliderInput(shiny::NS(id, "ic_percentage"), "Percentage if IC",
+      min = 1, max = 99, value = 50
+    ),
     shiny::uiOutput(shiny::NS(id, "DoseResponseUI"))
   )
 }

OpenStats/inst/test_data/2dose_response_history.json

@@ -17,6 +17,7 @@
   },
   {
     "type": "DoseResponse",
+    "IC [%]": 50,
     "Column containing the names": "names",
     "Column containing the units": "units",
     "Log transform x-axis": false,
@@ -26,6 +27,7 @@
   },
   {
     "type": "DoseResponse",
+    "IC [%]": 50,
     "Column containing the names": "names",
     "Column containing the units": "units",
     "Log transform x-axis": false,
@@ -35,6 +37,7 @@
   },
   {
     "type": "DoseResponse",
+    "IC [%]": 50,
     "Column containing the names": "names",
     "Column containing the units": "units",
     "Log transform x-axis": false,

OpenStats/inst/tinytest/test_Backend_DoseResponse.R

@@ -40,7 +40,7 @@ test_check_fit <- function(ic50_true) {
   abs <- "abs"
   title = "Bla"
   res <- OpenStats:::check_fit(
-    model, min(df[, conc]),
+    model, 50, min(df[, conc]),
     max(df[, conc]), min(df[, abs]), max(df[, abs]), title, "M"
   )
   expect_true(is.data.frame(res))
@@ -50,13 +50,13 @@ test_check_fit(10)
 # Test ic50 internal
 test_ic50_internal <- function(ic50_true) {
   data <- simulate("A", 7, ic50_true)
-  res <- OpenStats:::ic50_internal(data, "abs", "conc", "substance", FALSE, FALSE, "M")
+  res <- OpenStats:::ic_internal(data, 50, "abs", "conc", "substance", FALSE, FALSE, "M")
   res_df <- res[[1]]
   tol_percentage <- 0.1
   rel_error <- function(a, b) {
     abs(b - a) / b
   }
-  expect_true(rel_error(ic50_true, res_df$IC50_relative) < tol_percentage)
+  expect_true(rel_error(ic50_true, res_df$IC_50_relative) < tol_percentage)
 }
 test_ic50_internal(0.6)
 test_ic50_internal(5)
@@ -187,12 +187,12 @@ test_drawplot <- function() {
   abs <- "abs"
   title = "Bla"
   res <- OpenStats:::check_fit(
-    model, min(df[, conc]),
+    model, 50, min(df[, conc]),
     max(df[, conc]), min(df[, abs]), max(df[, abs]), title, "M"
   )
   p <- OpenStats:::drawplot(
-    df, abs, conc, "M", model, valid_points, title, res$IC50_relative,
-    res$IC50_relative_lower, res$IC50_relative_higher,
+    df, abs, conc, "M", model, valid_points, title, res$IC_50_relative,
+    res$IC_50_relative_lower, res$IC_50_relative_higher,
     FALSE, FALSE
   )
   layers <- p$layers
@@ -274,7 +274,7 @@ test_ic50 <- function() {
     names = c("A", "A", "A", "A", "A"),
     unit = "M"
   )
-  result <- OpenStats:::ic50(data, "abs", "conc", "names", "unit", FALSE, FALSE)
+  result <- OpenStats:::ic(data, 50, "abs", "conc", "names", "unit", FALSE, FALSE)
   checks[[1]] <- expect_true(is.list(result))
   checks[[2]] <- expect_true(is.data.frame(result[[1]][[1]]))
   expect_true(all(unlist(checks)))

OpenStats/inst/tinytest/test_Backend_TestEngine.R

@@ -718,6 +718,7 @@ test_dose_response <- function(in_background) {
 
   dr <- OpenStats:::dose_response_V1_2$new(
     df = df,
+    ic_percentage = 50,
     is_xlog = FALSE,
     is_ylog = FALSE,
     substance_names = "name",

OpenStats/inst/tinytest/test_Server_DoseResponse.R

@@ -2,8 +2,8 @@ df <- read.csv(system.file("/test_data/DoseResponse.csv", package = "OpenStats")
 
 # Expected
 # =================================================================
-expected <- OpenStats:::ic50(
-  df, "abs", "conc",
+expected <- OpenStats:::ic(
+  df, 50, "abs", "conc",
   "names", "units",
   FALSE, FALSE
 )
@@ -37,6 +37,7 @@ test_dose_response <- function(app, srv) {
   shiny::testServer(srv, {
     DataModelState$df      <- df
     DataModelState$formula <- new("LinearFormula", formula = abs ~ conc)
+    session$setInputs(`DOSERESPONSE-ic_percentage` = 50)
     session$setInputs(`DOSERESPONSE-substanceNames` = "names")
     session$setInputs(`DOSERESPONSE-unitNames` = "units")
     session$setInputs(`DOSERESPONSE-yTransform` = FALSE)