Better support metabolites in `sawtooth`

Also, add more example code to make it easier to use
`one_box`, `sawtooth` and `pfm_degradation`.

However, for generating sawtooth curves with specified
output times, overlaying the output of several calls to
`mkinpredict` is probably still easier.

1 files changed, 31 insertions, 4 deletions

diff --git a/R/pfm_degradation.R b/R/pfm_degradation.R
index 832a797..63a6e19 100644
--- a/R/pfm_degradation.R
+++ b/R/pfm_degradation.R
@@ -9,21 +9,48 @@
 #' @param years For how many years should the degradation be predicted?
 #' @param step_days What step size in days should the output have?
 #' @param times The output times
+#' @return A data frame containing the output times and the concentrations
+#'   assuming initial concentrations of 1 for the parent and zero for
+#'   metabolites, if any.
 #' @export
 #' @author Johannes Ranke
 #' @examples
-#' head(pfm_degradation("SFO", DT50 = 10))
+#' # Simple example of an SFO decline curve
+#' sfo_out <- pfm_degradation("SFO", DT50 = 10)
+#' head(sfo_out)
+#'
+#' # Fictive example with a metabolite where we first generate an SFO-SFO model
+#' sfo_sfo <- mkinmod(
+#'   parent = mkinsub("SFO", to = "metabolite"),
+#'   metabolite = mkinsub("SFO"))
+#'
+#' sfo_sfo_out <- pfm_degradation(sfo_sfo,
+#'   parms = c(k_parent = 0.1, f_parent_to_metabolite = 0.5, k_metabolite = 0.02))
+#'
+#' plot(
+#'   sfo_sfo_out[, "time"],
+#'   sfo_sfo_out[, "parent"], type = "l",
+#'   xlab = "Time", ylab = "Relative concentration",
+#'   xlim = c(0, 100))
+#' lines(
+#'   sfo_sfo_out[, "time"],
+#'   sfo_sfo_out[, "metabolite"], lty = 2)
+#'
 pfm_degradation <- function(model = "SFO",
   DT50 = 1000, parms = c(k_parent = log(2)/DT50),
   years = 1, step_days = 1,
   times = seq(0, years * 365, by = step_days))
 {
-  if (model %in% c("SFO", "FOMC", "DFOP", "HS", "IORE")) {
-    model <- mkinmod(parent = list(type = model))
+  if (!inherits(model, "mkinmod")) {
+    if (model[1] %in% c("SFO", "FOMC", "DFOP", "HS", "IORE")) {
+      model <- mkinmod(parent = list(type = model))
+    } else {
+      stop("Please specify the model with a suitable name or an mkinmod object")
+    }
   }
   initial_state = c(1, rep(0, length(model$diffs) - 1))
   names(initial_state) <- names(model$diffs)
-  time_course <- mkinpredict(model, odeparms = parms, 
+  time_course <- mkinpredict(model, odeparms = parms,
     odeini = initial_state,
     outtimes = times,
     solution_type = ifelse(length(model$spec) == 1, "analytical", "deSolve"))