Pull sm branch to SET branch (#8)

* Added conc resp vignette

* Removed tcpl comparison data from vignette

* Adjust age order to match GeoToxMIE script

NAMESPACE

@@ -3,6 +3,7 @@
 export(calc_concentration_response)
 export(calc_internal_dose)
 export(calc_invitro_concentration)
+export(fit_hill)
 export(simulate_age)
 export(simulate_exposure)
 export(simulate_inhalation_rate)

R/fit_hill.R

@@ -2,8 +2,10 @@
 #'
 #' @param log10_conc base-10 log scale concentration
 #' @param resp response
+#' @param fixed_slope if TRUE, slope is fixed at 1
 #'
 #' @return fit and other stats
+#' @export
 fit_hill <- function(log10_conc, resp, fixed_slope = TRUE) {
 
   # Compute initial values

vignettes/dev-conc-resp.Rmd

@@ -0,0 +1,114 @@
+---
+title: "dev-conc-resp"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{dev-conc-resp}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+```{r setup}
+library(GeoToxPackage)
+options(tidyverse.quiet = TRUE)
+library(tidyverse)
+```
+
+### Load data
+
+```{r load_data}
+# Load dose response data
+load("~/dev/GeoTox/data/LTEA_HepaRG_CYP1A1_up 41 chems for Kyle 220131.RData")
+ice_data <- cdat; rm(cdat)
+```
+
+### Format data
+
+```{r format_data}
+# Split dose response data by chemical
+ice_conc_resp <- split(as.data.frame(ice_data), ~casn)
+```
+
+### 2-parameter Hill fit
+
+```{r fits}
+fits <- lapply(ice_conc_resp, function(df) {
+  suppressWarnings(
+    fit_hill(df$logc, df$resp)
+  )
+})
+```
+
+### Extract fit parameters
+
+```{r fit_params}
+fit_params <- do.call(
+  rbind,
+  lapply(fits, function(fit) {
+    as_tibble(t(unlist(fit))) %>%
+      rename(
+        tp         = par.tp,
+        tp.sd      = sds.tp,
+        logAC50    = par.logAC50,
+        logAC50.sd = sds.logAC50
+      ) %>%
+      select(
+        tp, tp.sd, logAC50, logAC50.sd,
+        logc_min, logc_max, resp_min, resp_max, AIC
+      )
+  })
+)
+```
+
+### Replace NA sd with mean
+
+```{r replace_na}
+# TODO is this a good idea? See plots below
+fit_params <- fit_params %>%
+  mutate(
+    tp.na      = is.na(tp.sd),      # for plot below
+    logAC50.na = is.na(logAC50.sd), # for plot below
+    tp.sd      = if_else(is.na(tp.sd), tp, tp.sd),
+    logAC50.sd = if_else(is.na(logAC50.sd), logAC50, logAC50.sd)
+  )
+```
+
+```{r}
+xylim <- range(with(fit_params, c(tp, tp.sd)), na.rm = T)
+ggplot(fit_params, aes(tp, tp.sd)) +
+  geom_abline(linetype = 3) +
+  geom_point(aes(color = tp.na), show.legend = FALSE) +
+  coord_cartesian(xlim = xylim, ylim = xylim)
+```
+
+```{r}
+xylim <- range(with(fit_params, c(logAC50, logAC50.sd)), na.rm = T)
+ggplot(fit_params, aes(logAC50, logAC50.sd)) +
+  geom_abline(linetype = 3) +
+  geom_point(aes(color = logAC50.na), show.legend = FALSE) +
+  coord_cartesian(xlim = xylim, ylim = xylim)
+```
+
+### Plot fits
+
+```{r plot_fits}
+log10_x <- seq(-3, 3, length.out = 100)
+
+y <- as.matrix(
+  apply(fit_params[, c("tp", "logAC50")], 1, function(par) {
+    par["tp"] / (1 + 10^(par["logAC50"] - log10_x))
+  })
+)
+colnames(y) <- names(ice_conc_resp)
+y <- as_tibble(y) %>% mutate(x = 10^log10_x, .before = 1)
+
+ggplot(y %>% pivot_longer(!x), aes(x, value, color = name)) +
+  geom_line(show.legend = FALSE) +
+  scale_x_log10(labels = scales::label_math(10^.x, format = log10))
+```

vignettes/dev-sensitivity.Rmd

@@ -100,7 +100,16 @@ age <- age.data %>%
 ### Simulate age
 
 ```{r sim_age}
-simulated_age <- lapply(split(age, ~FIPS), simulate_age, n = MC_iter)
+age_split <- split(age, ~FIPS)
+
+# Adjust order due to changing FIPS 46102 to 46113
+# TODO used to match GeoToxMIE output, but does this adjustment line up with
+#      other input data, e.g. C_ss?
+idx <- 1:length(age_split)
+idx[2379:2384] <- c(2384, 2379:2383)
+age_split <- age_split[idx]
+
+simulated_age <- lapply(age_split, simulate_age, n = MC_iter)
 
 # Replace sampled age values with median when age is fixed
 if (step != 1) {