Convert main vignette to Rmd/html, add_err(), fixes

2 files changed, 56 insertions, 7 deletions

diff --git a/R/add_err.R b/R/add_err.R
new file mode 100644
index 00000000..0995e634
--- /dev/null
+++ b/R/add_err.R
@@ -0,0 +1,48 @@
+# Copyright (C) 2015-2016 Johannes Ranke
+# Contact: jranke@uni-bremen.de
+
+# This file is part of the R package mkin
+
+# mkin is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+
+# You should have received a copy of the GNU General Public License along with
+# this program. If not, see <http://www.gnu.org/licenses/>
+
+add_err = function(prediction, sdfunc, 
+                   n = 1000, LOD = 0.1, reps = 2, 
+                   digits = 1, seed = NA)
+{
+  if (!is.na(seed)) set.seed(seed)
+
+  # The output of mkinpredict is in wide format
+  d_long = mkin_wide_to_long(prediction, time = "time")
+
+  # Set up the list to be returned
+  d_return = list()
+
+  # Generate datasets one by one in a loop
+  for (i in 1:n) {
+    d_rep = data.frame(lapply(d_long, rep, each = 2))
+    d_rep$value = rnorm(length(d_rep$value), d_rep$value, sdfunc(d_rep$value))
+         
+    d_rep[d_rep$time == 0 & d_rep$name %in% c("M1", "M2"), "value"] <- 0
+
+    # Set values below the LOD to NA
+    d_NA <- transform(d_rep, value = ifelse(value < LOD, NA, value))
+
+    # Round the values for convenience
+    d_NA$value <- round(d_NA$value, digits)
+
+    d_return[[i]] <- d_NA
+  }
+
+  return(d_return)
+}
diff --git a/R/plot.mmkin.R b/R/plot.mmkin.R
index 02bf4d6e..7b54be4b 100644
--- a/R/plot.mmkin.R
+++ b/R/plot.mmkin.R
@@ -20,12 +20,17 @@ plot.mmkin <- function(x, main = "auto", legends = 1, errmin_var = "All data", e
                        cex = 0.7, rel.height.middle = 0.9, ...) {
   n.m <- nrow(x)
   n.d <- ncol(x)
+
+  # We can handle either a row (different models, same dataset) 
+  # or a column (same model, different datasets)
   if (n.m > 1 & n.d > 1) stop("Please select fits either for one model or for one dataset")
   if (n.m == 1 & n.d == 1) loop_over = "none"
   if (n.m > 1) loop_over <- "models"
   if (n.d > 1) loop_over <- "datasets"
   n.fits <- length(x)
 
+  # Set the main plot titles from the names of the models or the datasets
+  # Will be integer indexes if no other names are present in the mmkin object
   if (main == "auto") {
     main = switch(loop_over,
                   none = paste(rownames(x), colnames(x)),
@@ -34,11 +39,13 @@ plot.mmkin <- function(x, main = "auto", legends = 1, errmin_var = "All data", e
   }
 
   oldpar <- par(no.readonly = TRUE)
+
+  # Set relative plot heights, so the first and the last plot are the norm
+  # and the middle plots (if n.fits >2) are smaller by rel.height.middle
   rel.heights <- if (n.fits > 2) c(1, rep(rel.height.middle, n.fits - 2), 1)
                  else rep(1, n.fits)
   layout(matrix(1:(2 * n.fits), n.fits, 2, byrow = TRUE), heights = rel.heights)
 
-  #par(mfrow = c(n.fits, 2))
   par(mar = c(3.0, 4.1, 4.1, 2.1)) # Reduce bottom margin by 2.1 - hides x axis legend
   par(cex = cex)
 
@@ -48,12 +55,6 @@ plot.mmkin <- function(x, main = "auto", legends = 1, errmin_var = "All data", e
       # reduced plot height, therefore we need rel.height.middle in the layout
       par(mar = c(3.0, 4.1, 2.1, 2.1))
     }
-    if (i.fit == n.fits) {
-      # Reduce top margin by 2 after the first plot as we have no main title, 
-      # plot height remains about constant
-      par(mar = c(5.1, 4.1, 2.1, 2.1))
-
-    }
     fit <- x[[i.fit]]
     plot(fit, legend = legends == i.fit, ...)