Skip to content
Snippets Groups Projects

dissect ecxsys()

Merged dissect ecxsys()
refactor/dissect-ecxsys into v2.0
Merged Sebastian Henz requested to merge refactor/dissect-ecxsys into v2.0
Compare
and
Show latest version
3 files
+ 73
− 84
Compare changes
  • Side-by-side
  • Inline
Files
3
R/ecxsys.R
+ 72
− 61
@@ -122,6 +122,9 @@ ecxsys <- function(concentration,
conc_shift <- 2 # Powers of ten to shift the control downwards from the
# second lowest concentration. This is required to approximate 0 because
# of the logarithmic axis.
if (is.unsorted(concentration)) {
stop("The values must be sorted by increasing concentration.")
}
if (any(concentration < 0)) {
stop("Concentrations must be >= 0")
} else if (min(concentration) > 0) {
@@ -129,9 +132,6 @@ ecxsys <- function(concentration,
} else {
min_conc <- 10 ^ floor(log10(concentration[2]) - conc_shift)
}
if (is.unsorted(concentration)) {
stop("The values must be sorted by increasing concentration.")
}
# scale observed effect -----------------------------------------------
@@ -141,53 +141,18 @@ ecxsys <- function(concentration,
effect_tox_env_observed <- effect_tox_env_observed / effect_max
# prepare adjusted control concentration ------------------------------
conc_adjust_factor <- 10^-5
output$conc_adjust_factor <- conc_adjust_factor
# traditional simple model (LL.5) -------------------------------------
conc_with_control_shifted <- c(min_conc, concentration[-1])
effect_tox_observed_averaged <- moving_weighted_mean(effect_tox_observed)
effect_tox_env_observed_averaged <- moving_weighted_mean(effect_tox_env_observed)
temp <- approx(
log10(conc_with_control_shifted),
effect_tox_observed_averaged,
n = 10
)
conc_interpolated <- 10^temp$x
effect_tox_observed_interpolated_simple_model <- temp$y
option_warn_original <- getOption("warn")
effect_tox_mod_simple <- drc::drm(
effect_tox_observed_interpolated_simple_model ~ conc_interpolated,
fct = drc::LL.5(fixed = c(
NA, 0, effect_tox_observed_averaged[1], NA, NA
))
)
options(warn = option_warn_original)
effect_tox_simple <- predict(effect_tox_mod_simple, data.frame(concentration))
output$effect_tox_mod_simple <- effect_tox_mod_simple
output$effect_tox_simple <- effect_tox_simple * effect_max
simple_tox <- fit_simple_model(min_conc, concentration,
effect_tox_observed,
original_options)
output$effect_tox_mod_simple <- simple_tox$effect_simple_mod
output$effect_tox_simple <- simple_tox$effect_simple * effect_max
if (with_env) {
effect_tox_env_observed_interpolated_simple_model <- approx(
log10(conc_with_control_shifted),
effect_tox_env_observed_averaged,
xout = temp$x
)$y
option_warn_original <- getOption("warn")
effect_tox_env_mod_simple <- drc::drm(
effect_tox_env_observed_interpolated_simple_model ~ conc_interpolated,
fct = drc::LL.5(fixed = c(
NA, 0, effect_tox_env_observed_averaged[1], NA, NA
))
)
options(warn = option_warn_original)
effect_tox_env_simple <- predict(
effect_tox_env_mod_simple,
data.frame(concentration)
)
output$effect_tox_env_mod_simple <- effect_tox_env_mod_simple
output$effect_tox_env_simple <- effect_tox_env_simple * effect_max
simple_tox_env <- fit_simple_model(min_conc, concentration,
effect_tox_env_observed,
original_options)
output$effect_tox_env_mod_simple <- simple_tox_env$effect_simple_mod
output$effect_tox_env_simple <- simple_tox_env$effect_simple * effect_max
}
@@ -269,7 +234,7 @@ ecxsys <- function(concentration,
sys_stress_tox[hormesis_index:length(sys_stress_tox)] <- 0
# Add sys_stress_tox to the output before it is fitted:
output$sys_stress_tox <- sys_stress_tox[keep]
sys_stress_tox_mod <- fit_sys(stress_tox, sys_stress_tox)
sys_stress_tox_mod <- fit_sys(stress_tox, sys_stress_tox, original_options)
output$sys_stress_tox_mod <- sys_stress_tox_mod
sys_stress_tox <- unname(
predict(sys_stress_tox_mod, data.frame(stress_tox))
@@ -302,7 +267,8 @@ ecxsys <- function(concentration,
sys_stress_tox_env[hormesis_index:length(sys_stress_tox_env)] <- 0
# Add sys_stress_tox_env to the output before it is fitted:
output$sys_stress_tox_env <- sys_stress_tox_env[keep]
sys_stress_tox_env_mod <- fit_sys(stress_tox, sys_stress_tox_env)
sys_stress_tox_env_mod <- fit_sys(stress_tox, sys_stress_tox_env,
original_options)
output$sys_stress_tox_env_mod <- sys_stress_tox_env_mod
sys_stress_tox_env <- unname(
predict(sys_stress_tox_env_mod, data.frame(stress_tox))
@@ -322,6 +288,8 @@ ecxsys <- function(concentration,
# 0 but because it's a log axis I have to make the values just really
# small. The concentrations in the gap won't be used for plotting later.
n_smooth <- 1000 # number of points to approximate the curves
conc_adjust_factor <- 10^-5
output$conc_adjust_factor <- conc_adjust_factor
concentration_smooth <- 10 ^ seq(
log10(min_conc * conc_adjust_factor),
log10(max(concentration)),
@@ -339,12 +307,11 @@ ecxsys <- function(concentration,
reset_options <- function(original_options) {
# Reset all the options which have changed.
# You may ask why I don't just cache and reset the options like this:
# original_options <- options()
# options(original_options)
# The reason is that when I do this and ecxsys generates warnings then those
# warnings don't show up in the console. I don't know why, but resetting
# only the options which have changed alleviates that problem.
# You may ask why I don't just reset the options using
# options(original_options). The reason is that when I do this and ecxsys
# generates warnings then those warnings don't show up in the console. I
# don't know why, but resetting only the options which have changed
# alleviates that problem.
changed <- list()
for (n in names(original_options)) {
@@ -357,14 +324,11 @@ reset_options <- function(original_options) {
}
fit_sys <- function(stress_tox, sys) {
fit_sys <- function(stress_tox, sys, original_options) {
result_name <- paste0(deparse(substitute(sys)), "_mod")
# There is no other way to suppress that one error message
# from optim except by changing the options temporarily.
warn_error_original <- list(
warn = getOption("warn"),
show.error.messages = getOption("show.error.messages")
)
warn_error_original <- original_options[c("warn", "show.error.messages")]
options(show.error.messages = FALSE)
tryCatch(
drc::drm(sys ~ stress_tox, fct = drc::W1.3()),
@@ -387,3 +351,50 @@ fit_sys <- function(stress_tox, sys) {
)
}
moving_weighted_mean <- function(x) {
# This is used to smooth out points which are jumping up and down.
count <- rep(1, length(x))
x_diff <- diff(x) * -1
while (any(x_diff < 0, na.rm = TRUE)) {
i <- which(x_diff < 0)[1]
j <- i + 1
x[i] <- weighted.mean(x[i:j], count[i:j])
x <- x[-j]
count[i] <- count[i] + count[j]
count <- count[-j]
x_diff <- diff(x) * -1
}
rep(x, count)
}
fit_simple_model <- function(min_conc,
concentration,
effect_observed,
original_options) {
# The traditional log-logistic model, here using the five-parameter
# log-logistic function drc::LL.5.
conc_with_control_shifted <- c(min_conc, concentration[-1])
effect_observed_averaged <- moving_weighted_mean(effect_observed)
interpolated <- approx(
log10(conc_with_control_shifted),
effect_observed_averaged,
n = 10
)
conc_interpolated <- 10^interpolated$x
effect_interpolated <- interpolated$y
on.exit(options(original_options["warn"]))
effect_simple_mod <- drc::drm(
effect_interpolated ~ conc_interpolated,
fct = drc::LL.5(fixed = c(NA, 0, effect_observed_averaged[1], NA, NA))
)
effect_simple <- predict(
effect_simple_mod,
data.frame(conc_interpolated = concentration)
)
list(
effect_simple_mod = effect_simple_mod,
effect_simple = effect_simple
)
}