Add the logistic model

NEWS.md

@@ -1,4 +1,4 @@
-# mkin 0.9.47.6 (2019-01-31)
+# mkin 0.9.48.1 (2019-02-21)
 
 - Add the function 'logLik.mkinfit' which makes it possible to calculate an AIC for mkinfit objects
 
@@ -16,6 +16,8 @@
 
 - Add the function 'CAKE_export' to facilitate cross-checking of results
 
+- Implement the logistic model (only tested for parent fits)
+
 # mkin 0.9.47.5 (2018-09-14)
 
 - Make the two-component error model stop in cases where it is inadequate to avoid nls crashes on windows

R/endpoints.R

@@ -159,6 +159,14 @@ endpoints <- function(fit) {
       ep$distimes[obs_var, c(paste("DT50", obs_var, "b1", sep = "_"))] = DT50_b1
       ep$distimes[obs_var, c(paste("DT50", obs_var, "b2", sep = "_"))] = DT50_b2
     }
+    if (type == "logistic") {
+      # FOCUS kinetics (2014) p. 67
+      kmax = parms.all["kmax"]
+      k0 = parms.all["k0"]
+      r = parms.all["r"]
+      DT50 = (1/r) * log(1 - ((kmax/k0) * (1 - 2^(r/kmax))))
+      DT90 = (1/r) * log(1 - ((kmax/k0) * (1 - 10^(r/kmax))))
+    }
     ep$distimes[obs_var, c("DT50", "DT90")] = c(DT50, DT90)
   }
   return(ep)

R/logistic.solution.R

@@ -0,0 +1,4 @@
+logistic.solution <- function(t, parent.0, kmax, k0, r)
+{
+	parent = parent.0 * (kmax / (kmax - k0 + k0 * exp (r * t))) ^(kmax/r)
+}

R/mkinfit.R

@@ -1,4 +1,4 @@
-# Copyright (C) 2010-2018 Johannes Ranke
+# Copyright (C) 2010-2019 Johannes Ranke
 # Portions of this code are copyright (C) 2013 Eurofins Regulatory AG
 # Contact: jranke@uni-bremen.de
 # The summary function is an adapted and extended version of summary.modFit
@@ -48,7 +48,7 @@ mkinfit <- function(mkinmod, observed,
 {
   # Check mkinmod and generate a model for the variable whith the highest value
   # if a suitable string is given
-  parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE")
+  parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE", "logistic")
   if (class(mkinmod) != "mkinmod") {
     presumed_parent_name = observed[which.max(observed$value), "name"]
     if (mkinmod[[1]] %in% parent_models_available) {
@@ -153,6 +153,9 @@ mkinfit <- function(mkinmod, observed,
     if (parmname == "k2") parms.ini[parmname] = 0.01
     if (parmname == "tb") parms.ini[parmname] = 5
     if (parmname == "g") parms.ini[parmname] = 0.5
+    if (parmname == "kmax") parms.ini[parmname] = 0.1
+    if (parmname == "k0") parms.ini[parmname] = 0.0001
+    if (parmname == "r") parms.ini[parmname] = 0.2
   }
   # Default values for formation fractions in case they are present
   for (box in mod_vars) {
@@ -376,7 +379,7 @@ mkinfit <- function(mkinmod, observed,
     lower[index_k] <- 0
     index_k__iore <- grep("^k__iore_", names(lower))
     lower[index_k__iore] <- 0
-    other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2", "tb"), names(lower))
+    other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2", "tb", "r"), names(lower))
     lower[other_rate_parms] <- 0
   }
 

R/mkinmod.R

@@ -1,4 +1,4 @@
-# Copyright (C) 2010-2015 Johannes Ranke {{{
+# Copyright (C) 2010-2015,2019 Johannes Ranke {{{
 # Contact: jranke@uni-bremen.de
 
 # This file is part of the R package mkin
@@ -42,8 +42,8 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb
   parms <- vector()
   # }}}
 
-  # Do not return a coefficient matrix mat when FOMC, IORE, DFOP or HS is used for the parent {{{
-  if(spec[[1]]$type %in% c("FOMC", "IORE", "DFOP", "HS")) {
+  # Do not return a coefficient matrix mat when FOMC, IORE, DFOP, HS or logistic is used for the parent {{{
+  if(spec[[1]]$type %in% c("FOMC", "IORE", "DFOP", "HS", "logistic")) {
     mat = FALSE
   } else mat = TRUE
   #}}}
@@ -57,10 +57,10 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb
     # Check the type component of the compartment specification {{{
     if(is.null(spec[[varname]]$type)) stop(
       "Every part of the model specification must be a list containing a type component")
-    if(!spec[[varname]]$type %in% c("SFO", "FOMC", "IORE", "DFOP", "HS", "SFORB")) stop(
-      "Available types are SFO, FOMC, IORE, DFOP, HS and SFORB only")
-    if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS") & match(varname, obs_vars) != 1) {
-        stop(paste("Types FOMC, DFOP and HS are only implemented for the first compartment,",
+    if(!spec[[varname]]$type %in% c("SFO", "FOMC", "IORE", "DFOP", "HS", "SFORB", "logistic")) stop(
+      "Available types are SFO, FOMC, IORE, DFOP, HS, SFORB and logistic only")
+    if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS", "logistic") & match(varname, obs_vars) != 1) {
+        stop(paste("Types FOMC, DFOP, HS and logistic are only implemented for the first compartment,",
                    "which is assumed to be the source compartment"))
     }
     #}}}
@@ -71,6 +71,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb
       IORE = varname,
       DFOP = varname,
       HS = varname,
+      logistic = varname,
       SFORB = paste(varname, c("free", "bound"), sep = "_")
     )
     map[[varname]] <- new_boxes
@@ -141,6 +142,11 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb
       decline_term <- paste(HS_decline, "*", box_1)
       parms <- c(parms, "k1", "k2", "tb")
     } #}}}
+    if(spec[[varname]]$type == "logistic") { # {{{ Add logistic decline term
+      # From p. 67 of the FOCUS kinetics report (2014)
+      decline_term <- paste("(k0 * kmax)/(k0 + (kmax - k0) * exp(-r * time)) *", box_1)
+      parms <- c(parms, "kmax", "k0", "r")
+    } #}}}
     # Add origin decline term to box 1 (usually the only box, unless type is SFORB)#{{{
     diffs[[box_1]] <- paste(diffs[[box_1]], "-", decline_term)#}}}
     if(spec[[varname]]$type == "SFORB") { # {{{ Add SFORB reversible binding terms

R/mkinpredict.R

@@ -1,4 +1,4 @@
-# Copyright (C) 2010-2016,2018 Johannes Ranke
+# Copyright (C) 2010-2016,2018,2019 Johannes Ranke
 # Some lines in this code are copyright (C) 2013 Eurofins Regulatory AG
 # Contact: jranke@uni-bremen.de
 
@@ -83,7 +83,11 @@ mkinpredict.mkinmod <- function(x,
           evalparse(parent.name),
           evalparse(paste("k", parent.name, "bound", sep="_")),
           evalparse(paste("k", sub("free", "bound", parent.name), "free", sep="_")),
-          evalparse(paste("k", parent.name, "sink", sep="_")))
+          evalparse(paste("k", parent.name, "sink", sep="_"))),
+      logistic = logistic.solution(outtimes,
+          evalparse(parent.name),
+          evalparse("kmax"), evalparse("k0"),
+          evalparse("r"))
     )
     out <- data.frame(outtimes, o)
     names(out) <- c("time", sub("_free", "", parent.name))

R/transform_odeparms.R

@@ -1,4 +1,4 @@
-# Copyright (C) 2010-2014 Johannes Ranke
+# Copyright (C) 2010-2014,2019 Johannes Ranke
 # Contact: jranke@uni-bremen.de
 
 # This file is part of the R package mkin
@@ -71,8 +71,9 @@ transform_odeparms <- function(parms, mkinmod,
   }
 
   # Transform also FOMC parameters alpha and beta, DFOP and HS rates k1 and k2
-  # and HS parameter tb if transformation of rates is requested
-  for (pname in c("alpha", "beta", "k1", "k2", "tb")) {
+  # and HS parameter tb as well as logistic model parameters kmax, k0 and r if
+  # transformation of rates is requested
+  for (pname in c("alpha", "beta", "k1", "k2", "tb", "kmax", "k0", "r")) {
     if (!is.na(parms[pname])) {
       if (transform_rates) {
         transparms[paste0("log_", pname)] <- log(parms[pname])
@@ -151,8 +152,8 @@ backtransform_odeparms <- function(transparms, mkinmod,
     }
   }
 
-  # Transform parameters also for FOMC, DFOP and HS models
-  for (pname in c("alpha", "beta", "k1", "k2", "tb")) {
+  # Transform parameters also for FOMC, DFOP, HS and logistic models
+  for (pname in c("alpha", "beta", "k1", "k2", "tb", "kmax", "k0", "r")) {
     if (transform_rates) {
       pname_trans = paste0("log_", pname)
       if (!is.na(transparms[pname_trans])) {

_pkgdown.yml

@@ -61,6 +61,7 @@ reference:
       - SFORB.solution
       - HS.solution
       - IORE.solution
+      - logistic.solution
   - title: Generate synthetic datasets
     contents:
       - add_err