jranke / mkin
1
utils::globalVariables(c("name", "value_mean"))
2

3
#' Calculate the minimum error to assume in order to pass the variance test
4
#'
5
#' This function finds the smallest relative error still resulting in passing
6
#' the chi-squared test as defined in the FOCUS kinetics report from 2006.
7
#'
8
#' This function is used internally by \code{\link{summary.mkinfit}}.
9
#'
10
#' @param fit an object of class \code{\link{mkinfit}}.
11
#' @param alpha The confidence level chosen for the chi-squared test.
12
#' @importFrom stats qchisq aggregate
13
#' @return A dataframe with the following components: \item{err.min}{The
14
#' relative error, expressed as a fraction.} \item{n.optim}{The number of
15
#' optimised parameters attributed to the data series.} \item{df}{The number of
16
#' remaining degrees of freedom for the chi2 error level calculations.  Note
17
#' that mean values are used for the chi2 statistic and therefore every time
18
#' point with observed values in the series only counts one time.} The
19
#' dataframe has one row for the total dataset and one further row for each
20
#' observed state variable in the model.
21
#' @references FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence
22
#' and Degradation Kinetics from Environmental Fate Studies on Pesticides in EU
23
#' Registration} Report of the FOCUS Work Group on Degradation Kinetics, EC
24
#' Document Reference Sanco/10058/2005 version 2.0, 434 pp,
25
#' \url{http://esdac.jrc.ec.europa.eu/projects/degradation-kinetics}
26
#' @keywords manip
27
#' @examples
28
#'
29
#' SFO_SFO = mkinmod(parent = mkinsub("SFO", to = "m1"),
30
#'                   m1 = mkinsub("SFO"),
31
#'                   use_of_ff = "max")
32
#'
33
#' fit_FOCUS_D = mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE)
34
#' round(mkinerrmin(fit_FOCUS_D), 4)
35
#' \dontrun{
36
#'   fit_FOCUS_E = mkinfit(SFO_SFO, FOCUS_2006_E, quiet = TRUE)
37
#'   round(mkinerrmin(fit_FOCUS_E), 4)
38
#' }
39
#'
40
#' @export
41
mkinerrmin <- function(fit, alpha = 0.05)
42
{
43 1
  parms.optim <- fit$par
44

45 1
  kinerrmin <- function(errdata, n.parms) {
46 1
    means.mean <- mean(errdata$observed, na.rm = TRUE)
47 1
    df = nrow(errdata) - n.parms
48

49 1
    err.min <- sqrt((1 / qchisq(1 - alpha, df)) *
50 1
               sum((errdata$observed - errdata$predicted)^2)/(means.mean^2))
51

52 1
    return(list(err.min = err.min, n.optim = n.parms, df = df))
53
  }
54

55 1
  errdata <- aggregate(cbind(observed, predicted) ~ time + variable, data = fit$data, mean, na.rm=TRUE)
56 1
  errdata <- errdata[order(errdata$time, errdata$variable), ]
57

58
  # Remove values at time zero for variables whose value for state.ini is fixed,
59
  # as these will not have any effect in the optimization and should therefore not
60
  # be counted as degrees of freedom.
61 1
  fixed_initials = gsub("_0$", "", rownames(subset(fit$fixed, type = "state")))
62 1
  errdata <- subset(errdata, !(time == 0 & variable %in% fixed_initials))
63

64 1
  n.optim.overall <- length(parms.optim) - length(fit$errparms)
65

66 1
  errmin.overall <- kinerrmin(errdata, n.optim.overall)
67 1
  errmin <- data.frame(err.min = errmin.overall$err.min,
68 1
    n.optim = errmin.overall$n.optim, df = errmin.overall$df)
69 1
  rownames(errmin) <- "All data"
70

71
  # The degrees of freedom are counted according to FOCUS kinetics (2011, p. 164)
72 1
  for (obs_var in fit$obs_vars)
73
  {
74 1
    errdata.var <- subset(errdata, variable == obs_var)
75

76
    # Check if initial value is optimised
77 1
    n.initials.optim <- length(grep(paste(obs_var, ".*", "_0", sep=""), names(parms.optim)))
78

79
    # Rate constants and IORE exponents are attributed to the source variable
80 1
    n.k.optim <- length(grep(paste("^k", obs_var, sep="_"), names(parms.optim)))
81 1
    n.k.optim <- n.k.optim + length(grep(paste("^log_k", obs_var, sep="_"),
82 1
                                         names(parms.optim)))
83 1
    n.k__iore.optim <- length(grep(paste("^k__iore", obs_var, sep="_"), names(parms.optim)))
84 1
    n.k__iore.optim <- n.k__iore.optim + length(grep(paste("^log_k__iore",
85 1
          obs_var, sep = "_"), names(parms.optim)))
86

87 1
    n.N.optim <- length(grep(paste("^N", obs_var, sep="_"), names(parms.optim)))
88

89 1
    n.ff.optim <- 0
90
    # Formation fractions are attributed to the target variable, so look
91
    # for source compartments with formation fractions
92 1
    for (source_var in fit$obs_vars) {
93 1
      n.ff.source = length(grep(paste("^f", source_var, sep = "_"),
94 1
                                 names(parms.optim)))
95 1
      n.paths.source = length(fit$mkinmod$spec[[source_var]]$to)
96 1
      for (target_var in fit$mkinmod$spec[[source_var]]$to) {
97 1
        if (obs_var == target_var) {
98 1
          n.ff.optim <- n.ff.optim + n.ff.source/n.paths.source
99
        }
100
      }
101
    }
102

103 1
    n.optim <- sum(n.initials.optim, n.k.optim, n.k__iore.optim, n.N.optim, n.ff.optim)
104

105
    # FOMC, DFOP and HS parameters are only counted if we are looking at the
106
    # first variable in the model which is always the source variable
107 1
    if (obs_var == fit$obs_vars[[1]]) {
108 1
      special_parms = c("alpha", "log_alpha", "beta", "log_beta",
109 1
                        "k1", "log_k1", "k2", "log_k2",
110 1
                        "g", "g_ilr", "g_qlogis", "tb", "log_tb")
111 1
      n.optim <- n.optim + length(intersect(special_parms, names(parms.optim)))
112
    }
113

114
    # Calculate and add a line to the dataframe holding the results
115 1
    errmin.tmp <- kinerrmin(errdata.var, n.optim)
116 1
    errmin[obs_var, c("err.min", "n.optim", "df")] <- errmin.tmp
117
  }
118

119 1
  return(errmin)
120
}

Read our documentation on viewing source code .

Loading