Merge pull request #48 from NIEHS/sm

Multi-assay support added

DESCRIPTION

@@ -21,6 +21,7 @@ Imports:
     stats,
     tibble,
     tidyr,
+    tidyselect,
     truncnorm,
     utils
 Suggests: 

R/GeoTox.R

@@ -24,7 +24,7 @@ GeoTox <- function() {
         exposure = list(
           expos_mean  = "mean",
           expos_sd    = "sd",
-          expos_label = "casn"
+          expos_label = "CASRN"
         ),
         internal_dose = list(
           time    = 1,
@@ -43,7 +43,7 @@ GeoTox <- function() {
 print.GeoTox <- function(x, ...) {
 
   names_simulated <- c("age", "IR", "obesity", "C_ext", "C_ss")
-  names_computed <- c("D_int", "C_invitro", "resp")
+  names_computed <- c("D_int", "C_invitro", "resp", "sensitivity")
   names_other <- setdiff(names(x),
                          c(names_simulated, names_computed))
 
@@ -84,6 +84,23 @@ print.GeoTox <- function(x, ...) {
   info_computed <- info_computed[info_computed$Class != "", , drop = FALSE]
 
   cat("GeoTox object\n")
+  if (is.null(x$hill_params)) {
+    n_assays <- 0
+    n_chems <- 0
+  } else {
+    if ("assay" %in% names(x$hill_params)) {
+      n_assays <- length(unique(x$hill_params$assay))
+    } else {
+      n_assays <- 1
+    }
+    if ("chem" %in% names(x$hill_params)) {
+      n_chems <- length(unique(x$hill_params$chem))
+    } else {
+      n_chems <- 1
+    }
+  }
+  cat("Assays: ", n_assays, "\n", sep = "")
+  cat("Chemicals: ", n_chems, "\n", sep = "")
   if (nrow(info_simulated) > 0) {
     n_regions <- length(x[[info_simulated$Name[1]]])
   } else if (nrow(info_computed) > 0) {
@@ -121,15 +138,26 @@ plot.GeoTox <- function(x,
                         type = c("resp", "hill", "exposure", "sensitivity"),
                         ...) {
   type <- match.arg(type)
-  switch(type,
-         resp = plot_resp(x$resp,
-                          region_boundary = x$boundaries$region,
-                          group_boundary  = x$boundaries$group,
-                          ...),
-         hill = plot_hill(x$hill_params),
-         exposure = plot_exposure(x$exposure,
-                                  region_boundary = x$boundaries$region,
-                                  group_boundary  = x$boundaries$group,
-                                  ...),
-         sensitivity = plot_sensitivity(x, ...))
+  if (type == "resp") {
+    plot_resp(x$resp,
+              region_boundary = x$boundaries$region,
+              group_boundary  = x$boundaries$group,
+              ...)
+  } else if (type == "hill") {
+    plot_hill(x$hill_params,
+              ...)
+  } else if (type == "exposure") {
+    dots = list(...)
+    chem_label = dots$chem_label %||% x$par$exposure$expos_label
+    ncol = dots$ncol %||% 2
+    plot_exposure(x$exposure,
+                  region_boundary = x$boundaries$region,
+                  group_boundary  = x$boundaries$group,
+                  chem_label      = chem_label,
+                  ncol            = ncol)
+  } else if (type == "sensitivity") {
+    plot_sensitivity(x, ...)
+  } else {
+    stop("Invalid type.")
+  }
 }

R/calc_concentration_response.R

@@ -19,19 +19,46 @@ calc_concentration_response <- function(C_invitro,
                                         tp_b_mult = 1.5,
                                         fixed = FALSE) {
 
-  if (inherits(C_invitro, "matrix")) {
-    .calc_concentration_response(C_invitro, hill_params, tp_b_mult, fixed)
+  if (!any(c("matrix", "list") %in% class(C_invitro))) {
+    stop("C_invitro must be a matrix or list")
+  }
+  if (!is.list(C_invitro)) C_invitro <- list(C_invitro)
+
+  # Split hill_params by assay
+  if ("assay" %in% names(hill_params)) {
+    hill_params <- split(hill_params, ~assay)
   } else {
-    mapply(
-      .calc_concentration_response,
-      C_invitro = C_invitro,
-      hill_params = list(hill_params),
-      tp_b_mult = tp_b_mult,
-      fixed = fixed,
-      SIMPLIFY = FALSE
-    )
+    hill_params <- list(hill_params)
   }
-
+
+  # Calculate response for each assay
+  lapply(C_invitro, \(C_invitro_i) {
+    lapply(hill_params, \(hill_params_j) {
+      if (ncol(C_invitro_i) == 1 & nrow(hill_params_j) == 1) {
+        .calc_concentration_response(C_invitro_i, hill_params_j, tp_b_mult, fixed)
+      } else {
+        if (!"chem" %in% names(hill_params_j)) {
+          stop("'hill_params' must contain a 'chem' column", call. = FALSE)
+        }
+        chems <- hill_params_j$chem
+        if (!all(chems %in% colnames(C_invitro_i))) {
+          stop("'hill_params' chemicals missing in 'C_invitro'", call. = FALSE)
+        }
+        C_invitro_i <- C_invitro_i[, chems, drop = FALSE]
+        res <- .calc_concentration_response(C_invitro_i,
+                                            hill_params_j,
+                                            tp_b_mult,
+                                            fixed) |> 
+          dplyr::mutate(sample = dplyr::row_number(), .before = 1)
+        if ("assay" %in% names(hill_params_j)) {
+          res <- res |> 
+            dplyr::mutate(assay = hill_params_j$assay[[1]], .before = 1)
+        }
+        res
+      }
+    }) |> 
+      dplyr::bind_rows()
+  })
 }
 
 .calc_concentration_response <- function(
@@ -42,32 +69,46 @@ calc_concentration_response <- function(C_invitro,
 
   # TODO value of b not consistent
   # grep "tp.sim <-" ~/github/GeoToxMIE/*.R
-  tp <- t(sapply(1:nrow(C_invitro), function(x) {
+  tp <- lapply(1:nrow(C_invitro), function(x) {
     truncnorm::rtruncnorm(
       1,
       a    = 0,
       b    = hill_params$resp_max * tp_b_mult,
       mean = hill_params$tp,
       sd   = if (fixed) 0 else hill_params$tp.sd
     )
-  }))
+  }) |> unlist()
+  tp <- matrix(tp, nrow = nrow(C_invitro), byrow = TRUE)
+  # Replace NAs with 0 (can occur when tp and tp.sd are 0)
+  tp[is.na(tp)] <- 0
 
-  logAC50 <- t(sapply(1:nrow(C_invitro), function(x) {
+  logAC50 <- lapply(1:nrow(C_invitro), function(x) {
     truncnorm::rtruncnorm(
       1,
       a    = hill_params$logc_min - 2.0001,
       b    = hill_params$logc_max + 0.5001,
       mean = hill_params$logAC50,
       sd   = if (fixed) 0 else hill_params$logAC50.sd
     )
-  }))
+  }) |> unlist()
+  logAC50 <- matrix(logAC50, nrow = nrow(C_invitro), byrow = TRUE)
 
   GCA.eff <- IA.eff <- GCA.HQ.10 <- IA.HQ.10 <- rep(NA, nrow(C_invitro))
   for (i in 1:nrow(C_invitro)) {
 
     C_i <- C_invitro[i, ]
     tp_i <- tp[i, ]
     AC50_i <- 10^logAC50[i, ]
+
+    if (all(C_i == 0)) {
+      GCA.eff[i] <- IA.eff[i] <- GCA.HQ.10[i] <- IA.HQ.10[i] <- NA
+      next
+    } else {
+      idx <- which(C_i > 0)
+      C_i <- C_i[idx]
+      tp_i <- tp_i[idx]
+      AC50_i <- AC50_i[idx]
+    }
 
     mixture.result <- stats::optimize(
       obj_GCA, interval = interval, Ci = C_i, tp = tp_i, AC50 = AC50_i

R/calculate_response.R

@@ -14,18 +14,6 @@
 #'
 #' @seealso [calc_internal_dose], [calc_invitro_concentration],
 #' [calc_concentration_response]
-#' @examples
-#' # For information about geo_tox_data, see:
-#' # vignette("package_data", package = "GeoTox")
-#' 
-#' geoTox <- GeoTox() |> 
-#'   simulate_population(age = split(geo_tox_data$age, ~FIPS)[1:5],
-#'                       obesity = geo_tox_data$obesity[1:5, ],
-#'                       exposure = split(geo_tox_data$exposure, ~FIPS)[1:5],
-#'                       simulated_css = geo_tox_data$simulated_css,
-#'                       n = 10) |> 
-#'   set_hill_params(fit_hill(split(geo_tox_data$dose_response, ~casn))) |>
-#'   calculate_response()
 calculate_response <- function(x, ...) {
 
   # Update parameters

R/data.R

@@ -1,4 +1,4 @@
-#' Geo Tox data
+#' GeoTox data
 #'
 #' Sample data for use in vignettes and function examples. See the
 #' "package_data" vignette for details on how this data was gathered.