aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2020-12-07 10:50:29 +0100
committerJohannes Ranke <jranke@uni-bremen.de>2020-12-07 10:50:29 +0100
commit22633a64469cec54cdfbeb74ff640409e57651ba (patch)
tree96ebb0576b639788490ea928d9c4cb03249586c5 /R
parentd3531ffa0d1f190920174105facf799dd6c3a851 (diff)
Possibility to use saemix transformations in some cases
Namely for SFO, DFOP, SFO-SFO and DFOP-SFO
Diffstat (limited to 'R')
-rw-r--r--R/saem.R172
1 files changed, 127 insertions, 45 deletions
diff --git a/R/saem.R b/R/saem.R
index 3376d3c6..e5634ac7 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -17,6 +17,13 @@ utils::globalVariables(c("predicted", "std"))
#' [mkinmod] model to different datasets
#' @param verbose Should we print information about created objects of
#' type [saemix::SaemixModel] and [saemix::SaemixData]?
+#' @param transformations Per default, all parameter transformations are done
+#' in mkin. If this argument is set to 'saemix', parameter transformations
+#' are done in 'saemix' for the supported cases. Currently this is only
+#' supported in cases where the initial concentration of the parent is not fixed,
+#' SFO or DFOP is used for the parent and there is either no metabolite or one.
+#' @param solution_type Possibility to specify the solution type in case the
+#' automatic choice is not desired
#' @param quiet Should we suppress the messages saemix prints at the beginning
#' and the end of the optimisation process?
#' @param cores The number of cores to be used for multicore processing using
@@ -90,12 +97,15 @@ saem <- function(object, control, ...) UseMethod("saem")
#' @rdname saem
#' @export
saem.mmkin <- function(object,
+ transformations = c("mkin", "saemix"),
+ solution_type = "auto",
control = list(displayProgress = FALSE, print = FALSE,
save = FALSE, save.graphs = FALSE),
cores = 1,
verbose = FALSE, suppressPlot = TRUE, quiet = FALSE, ...)
{
- m_saemix <- saemix_model(object, cores = cores, verbose = verbose, ...)
+ m_saemix <- saemix_model(object, cores = cores, verbose = verbose,
+ solution_type = solution_type, transformations = transformations, ...)
d_saemix <- saemix_data(object, verbose = verbose)
if (suppressPlot) {
# We suppress the log-likelihood curve that saemix currently
@@ -113,10 +123,15 @@ saem.mmkin <- function(object,
}
transparms_optim <- f_saemix@results@fixed.effects
names(transparms_optim) <- f_saemix@results@name.fixed
- bparms_optim <- backtransform_odeparms(transparms_optim,
- object[[1]]$mkinmod,
- object[[1]]$transform_rates,
- object[[1]]$transform_fractions)
+
+ if (transformations == "mkin") {
+ bparms_optim <- backtransform_odeparms(transparms_optim,
+ object[[1]]$mkinmod,
+ object[[1]]$transform_rates,
+ object[[1]]$transform_fractions)
+ } else {
+ bparms_optim <- transparms_optim
+ }
return_data <- nlme_data(object)
@@ -136,6 +151,7 @@ saem.mmkin <- function(object,
mkinmod = object[[1]]$mkinmod,
mmkin = object,
solution_type = object[[1]]$solution_type,
+ transformations = transformations,
transform_rates = object[[1]]$transform_rates,
transform_fractions = object[[1]]$transform_fractions,
so = f_saemix,
@@ -188,13 +204,20 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
#' @rdname saem
#' @return An [saemix::SaemixModel] object.
#' @export
-saemix_model <- function(object, cores = 1, verbose = FALSE, ...) {
+saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"),
+ cores = 1, verbose = FALSE, ...)
+{
if (nrow(object) > 1) stop("Only row objects allowed")
mkin_model <- object[[1]]$mkinmod
- solution_type <- object[[1]]$solution_type
degparms_optim <- mean_degparms(object)
+ if (transformations == "saemix") {
+ degparms_optim <- backtransform_odeparms(degparms_optim,
+ object[[1]]$mkinmod,
+ object[[1]]$transform_rates,
+ object[[1]]$transform_fractions)
+ }
degparms_fixed <- object[[1]]$bparms.fixed
# Transformations are done in the degradation function
@@ -249,8 +272,15 @@ saemix_model <- function(object, cores = 1, verbose = FALSE, ...) {
}
} else {
if (parent_type == "SFO") {
- model_function <- function(psi, id, xidep) {
- psi[id, 1] * exp( - exp(psi[id, 2]) * xidep[, "time"])
+ if (transformations == "mkin") {
+ model_function <- function(psi, id, xidep) {
+ psi[id, 1] * exp( - exp(psi[id, 2]) * xidep[, "time"])
+ }
+ } else {
+ model_function <- function(psi, id, xidep) {
+ psi[id, 1] * exp( - psi[id, 2] * xidep[, "time"])
+ }
+ transform.par = c(0, 1)
}
}
if (parent_type == "FOMC") {
@@ -259,11 +289,21 @@ saemix_model <- function(object, cores = 1, verbose = FALSE, ...) {
}
}
if (parent_type == "DFOP") {
- model_function <- function(psi, id, xidep) {
- g <- plogis(psi[id, 4])
- t <- xidep[, "time"]
- psi[id, 1] * (g * exp(- exp(psi[id, 2]) * t) +
- (1 - g) * exp(- exp(psi[id, 3]) * t))
+ if (transformations == "mkin") {
+ model_function <- function(psi, id, xidep) {
+ g <- plogis(psi[id, 4])
+ t <- xidep[, "time"]
+ psi[id, 1] * (g * exp(- exp(psi[id, 2]) * t) +
+ (1 - g) * exp(- exp(psi[id, 3]) * t))
+ }
+ } else {
+ model_function <- function(psi, id, xidep) {
+ g <- psi[id, 4]
+ t <- xidep[, "time"]
+ psi[id, 1] * (g * exp(- psi[id, 2] * t) +
+ (1 - g) * exp(- psi[id, 3] * t))
+ }
+ transform.par = c(0, 1, 1, 3)
}
}
if (parent_type == "HS") {
@@ -282,53 +322,95 @@ saemix_model <- function(object, cores = 1, verbose = FALSE, ...) {
# Parent with one metabolite
# Parameter names used in the model functions are as in
# https://nbviewer.jupyter.org/urls/jrwb.de/nb/Symbolic%20ODE%20solutions%20for%20mkin.ipynb
- if (length(mkin_model$spec) == 2) {
- types <- unname(sapply(mkin_model$spec, function(x) x$type))
+ types <- unname(sapply(mkin_model$spec, function(x) x$type))
+ if (length(mkin_model$spec) == 2 &! "SFORB" %in% types ) {
# Initial value for the metabolite (n20) must be fixed
if (names(odeini_fixed) == names(mkin_model$spec)[2]) {
n20 <- odeini_fixed
parent_name <- names(mkin_model$spec)[1]
if (identical(types, c("SFO", "SFO"))) {
- model_function <- function(psi, id, xidep) {
- t <- xidep[, "time"]
- n10 <- psi[id, 1]
- k1 <- exp(psi[id, 2])
- k2 <- exp(psi[id, 3])
- f12 <- plogis(psi[id, 4])
- ifelse(xidep[, "name"] == parent_name,
- n10 * exp(- k1 * t),
- (((k2 - k1) * n20 - f12 * k1 * n10) * exp(- k2 * t)) / (k2 - k1) +
- (f12 * k1 * n10 * exp(- k1 * t)) / (k2 - k1)
- )
+ if (transformations == "mkin") {
+ model_function <- function(psi, id, xidep) {
+ t <- xidep[, "time"]
+ n10 <- psi[id, 1]
+ k1 <- exp(psi[id, 2])
+ k2 <- exp(psi[id, 3])
+ f12 <- plogis(psi[id, 4])
+ ifelse(xidep[, "name"] == parent_name,
+ n10 * exp(- k1 * t),
+ (((k2 - k1) * n20 - f12 * k1 * n10) * exp(- k2 * t)) / (k2 - k1) +
+ (f12 * k1 * n10 * exp(- k1 * t)) / (k2 - k1)
+ )
+ }
+ } else {
+ model_function <- function(psi, id, xidep) {
+ t <- xidep[, "time"]
+ n10 <- psi[id, 1]
+ k1 <- psi[id, 2]
+ k2 <- psi[id, 3]
+ f12 <- psi[id, 4]
+ ifelse(xidep[, "name"] == parent_name,
+ n10 * exp(- k1 * t),
+ (((k2 - k1) * n20 - f12 * k1 * n10) * exp(- k2 * t)) / (k2 - k1) +
+ (f12 * k1 * n10 * exp(- k1 * t)) / (k2 - k1)
+ )
+ }
+ transform.par = c(0, 1, 1, 3)
}
}
if (identical(types, c("DFOP", "SFO"))) {
- model_function <- function(psi, id, xidep) {
- t <- xidep[, "time"]
- n10 <- psi[id, 1]
- k2 <- exp(psi[id, 2])
- f12 <- plogis(psi[id, 3])
- l1 <- exp(psi[id, 4])
- l2 <- exp(psi[id, 5])
- g <- plogis(psi[id, 6])
- ifelse(xidep[, "name"] == parent_name,
- n10 * (g * exp(- l1 * t) + (1 - g) * exp(- l2 * t)),
- ((f12 * g - f12) * l2 * n10 * exp(- l2 * t)) / (l2 - k2) -
- (f12 * g * l1 * n10 * exp(- l1 * t)) / (l1 - k2) +
- ((((l1 - k2) * l2 - k2 * l1 + k2^2) * n20 +
- ((f12 * l1 + (f12 * g - f12) * k2) * l2 -
- f12 * g * k2 * l1) * n10) * exp( - k2 * t)) /
- ((l1 - k2) * l2 - k2 * l1 + k2^2)
- )
+ if (transformations == "mkin") {
+ model_function <- function(psi, id, xidep) {
+ t <- xidep[, "time"]
+ n10 <- psi[id, 1]
+ k2 <- exp(psi[id, 2])
+ f12 <- plogis(psi[id, 3])
+ l1 <- exp(psi[id, 4])
+ l2 <- exp(psi[id, 5])
+ g <- plogis(psi[id, 6])
+ ifelse(xidep[, "name"] == parent_name,
+ n10 * (g * exp(- l1 * t) + (1 - g) * exp(- l2 * t)),
+ ((f12 * g - f12) * l2 * n10 * exp(- l2 * t)) / (l2 - k2) -
+ (f12 * g * l1 * n10 * exp(- l1 * t)) / (l1 - k2) +
+ ((((l1 - k2) * l2 - k2 * l1 + k2^2) * n20 +
+ ((f12 * l1 + (f12 * g - f12) * k2) * l2 -
+ f12 * g * k2 * l1) * n10) * exp( - k2 * t)) /
+ ((l1 - k2) * l2 - k2 * l1 + k2^2)
+ )
+ }
+ } else {
+ model_function <- function(psi, id, xidep) {
+ t <- xidep[, "time"]
+ n10 <- psi[id, 1]
+ k2 <- psi[id, 2]
+ f12 <- psi[id, 3]
+ l1 <- psi[id, 4]
+ l2 <- psi[id, 5]
+ g <- psi[id, 6]
+ ifelse(xidep[, "name"] == parent_name,
+ n10 * (g * exp(- l1 * t) + (1 - g) * exp(- l2 * t)),
+ ((f12 * g - f12) * l2 * n10 * exp(- l2 * t)) / (l2 - k2) -
+ (f12 * g * l1 * n10 * exp(- l1 * t)) / (l1 - k2) +
+ ((((l1 - k2) * l2 - k2 * l1 + k2^2) * n20 +
+ ((f12 * l1 + (f12 * g - f12) * k2) * l2 -
+ f12 * g * k2 * l1) * n10) * exp( - k2 * t)) /
+ ((l1 - k2) * l2 - k2 * l1 + k2^2)
+ )
+ }
+ transform.par = c(0, 1, 3, 1, 1, 3)
}
}
}
}
}
- if (is.function(model_function)) {
+ if (is.function(model_function) & solution_type == "auto") {
solution_type = "analytical saemix"
} else {
+
+ if (solution_type == "auto")
+ solution_type <- object[[1]]$solution_type
+
model_function <- function(psi, id, xidep) {
uid <- unique(id)

Contact - Imprint