jranke / mkin
1
#' Add normally distributed errors to simulated kinetic degradation data
2
#' 
3
#' Normally distributed errors are added to data predicted for a specific
4
#' degradation model using \code{\link{mkinpredict}}. The variance of the error
5
#' may depend on the predicted value and is specified as a standard deviation.
6
#' 
7
#' @param prediction A prediction from a kinetic model as produced by
8
#'   \code{\link{mkinpredict}}.
9
#' @param sdfunc A function taking the predicted value as its only argument and
10
#'   returning a standard deviation that should be used for generating the
11
#'   random error terms for this value.
12
#' @param secondary The names of state variables that should have an initial
13
#'   value of zero
14
#' @param n The number of datasets to be generated.
15
#' @param LOD The limit of detection (LOD). Values that are below the LOD after
16
#'   adding the random error will be set to NA.
17
#' @param reps The number of replicates to be generated within the datasets.
18
#' @param digits The number of digits to which the values will be rounded.
19
#' @param seed The seed used for the generation of random numbers. If NA, the
20
#'   seed is not set.
21
#' @importFrom stats rnorm
22
#' @return A list of datasets compatible with \code{\link{mmkin}}, i.e. the
23
#'   components of the list are datasets compatible with \code{\link{mkinfit}}.
24
#' @author Johannes Ranke
25
#' @references Ranke J and Lehmann R (2015) To t-test or not to t-test, that is
26
#' the question. XV Symposium on Pesticide Chemistry 2-4 September 2015,
27
#' Piacenza, Italy
28
#' https://jrwb.de/posters/piacenza_2015.pdf
29
#' @examples
30
#' 
31
#' # The kinetic model
32
#' m_SFO_SFO <- mkinmod(parent = mkinsub("SFO", "M1"),
33
#'                      M1 = mkinsub("SFO"), use_of_ff = "max")
34
#' 
35
#' # Generate a prediction for a specific set of parameters
36
#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
37
#' 
38
#' # This is the prediction used for the "Type 2 datasets" on the Piacenza poster
39
#' # from 2015
40
#' d_SFO_SFO <- mkinpredict(m_SFO_SFO,
41
#'                          c(k_parent = 0.1, f_parent_to_M1 = 0.5,
42
#'                            k_M1 = log(2)/1000),
43
#'                          c(parent = 100, M1 = 0),
44
#'                          sampling_times)
45
#' 
46
#' # Add an error term with a constant (independent of the value) standard deviation
47
#' # of 10, and generate three datasets
48
#' d_SFO_SFO_err <- add_err(d_SFO_SFO, function(x) 10, n = 3, seed = 123456789 )
49
#' 
50
#' # Name the datasets for nicer plotting
51
#' names(d_SFO_SFO_err) <- paste("Dataset", 1:3)
52
#' 
53
#' # Name the model in the list of models (with only one member in this case) for
54
#' # nicer plotting later on.  Be quiet and use only one core not to offend CRAN
55
#' # checks
56
#' \dontrun{
57
#' f_SFO_SFO <- mmkin(list("SFO-SFO" = m_SFO_SFO),
58
#'                    d_SFO_SFO_err, cores = 1,
59
#'                    quiet = TRUE)
60
#' 
61
#' plot(f_SFO_SFO)
62
#' 
63
#' # We would like to inspect the fit for dataset 3 more closely
64
#' # Using double brackets makes the returned object an mkinfit object
65
#' # instead of a list of mkinfit objects, so plot.mkinfit is used
66
#' plot(f_SFO_SFO[[3]], show_residuals = TRUE)
67
#' 
68
#' # If we use single brackets, we should give two indices (model and dataset),
69
#' # and plot.mmkin is used
70
#' plot(f_SFO_SFO[1, 3])
71
#' }
72
#' 
73
#' @export
74
add_err <- function(prediction, sdfunc, secondary = c("M1", "M2"),
75
  n = 1000, LOD = 0.1, reps = 2, digits = 1, seed = NA)
76
{
77 3
  if (!is.na(seed)) set.seed(seed)
78

79 3
  prediction <- as.data.frame(prediction)
80

81
  # The output of mkinpredict is in wide format
82 3
  d_long = mkin_wide_to_long(prediction, time = "time")
83

84
  # Set up the list to be returned
85 3
  d_return = list()
86

87
  # Generate datasets one by one in a loop
88 3
  for (i in 1:n) {
89 3
    d_rep = data.frame(lapply(d_long, rep, each = reps))
90 3
    d_rep$value = rnorm(length(d_rep$value), d_rep$value, sdfunc(d_rep$value))
91

92 3
    d_rep[d_rep$time == 0 & d_rep$name %in% secondary, "value"] <- 0
93

94
    # Set values below the LOD to NA
95 3
    d_NA <- transform(d_rep, value = ifelse(value < LOD, NA, value))
96

97
    # Round the values for convenience
98 3
    d_NA$value <- round(d_NA$value, digits)
99

100 3
    d_return[[i]] <- d_NA
101
  }
102

103 3
  return(d_return)
104
}

Read our documentation on viewing source code .

Loading