summaryrefslogblamecommitdiff
path: root/CakeErrMin.R
blob: ff49ad5d5ddfa06110d3eb7f350856a4b3c7e389 (plain) (tree) 
5b3daf3





ae64167


5b3daf3




















































































1
2
3
4
5

6
7

8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91





                                                                                

                                                                                    



















































































                                                                                                             
## -----------------------------------------------------------------------------
## Chi squared errmin function.
## -----------------------------------------------------------------------------
# Some of the CAKE R modules are based on mkin.
#
# Modifications developed by Tessella for Syngenta, Copyright (C) 2011-2020 Syngenta
# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091
#
# The CAKE R modules are free software: you can
# redistribute them and/or modify them 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/> 

if(getRversion() >= '2.15.1') utils::globalVariables(c("name"))

CakeErrMin <- function(fit, alpha = 0.05)
{
  parms.optim <- fit$par
  kinerrmin <- function(errdata, n.parms) {
    means.mean <- mean(errdata$value_mean, na.rm=TRUE)
    df = length(errdata$value_mean) - n.parms
  
    err.min <- sqrt((1 / qchisq(1 - alpha, df)) *
                sum((errdata$value_mean - errdata$value_pred)^2)/(means.mean^2))
    
    return(list(err.min = err.min, n.optim = n.parms, df = df))
  }

  means <- aggregate(value ~ time + name, data = fit$observed, mean, na.rm=TRUE)
  errdata <- merge(means, fit$predicted, by = c("time", "name"), 
    suffixes = c("_mean", "_pred"))
  errdata <- errdata[order(errdata$time, errdata$name), ]
  
  # Remove values at time zero for variables whose value for state.ini is fixed,
  # as these will not have any effect in the optimization and should therefore not
  # be counted as degrees of freedom.
  fixed_initials = gsub("_0$", "", rownames(subset(fit$fixed, type = "state")))
  errdata <- subset(errdata, !(time == 0 & name %in% fixed_initials))
  
  n.optim.overall <- length(parms.optim)

  errmin.overall <- kinerrmin(errdata, n.optim.overall)
  errmin <- data.frame(err.min = errmin.overall$err.min,
    n.optim = errmin.overall$n.optim, df = errmin.overall$df)
  rownames(errmin) <- "All data"

  for (obs_var in fit$obs_vars)
  {
    errdata.var <- subset(errdata, name == obs_var)

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

    # Rate constants and DFOP parameters are attributed to the source variable
    n.k.optim <- length(grep(paste("^k", obs_var, sep="_"), names(parms.optim)))
    n.k1.dfop.optim <- length(grep(paste("^k1", obs_var, sep="_"), names(parms.optim)))
    n.k2.dfop.optim <- length(grep(paste("^k2", obs_var, sep="_"), names(parms.optim)))
    n.g.dfop.optim <- length(grep(paste("^g", obs_var, sep="_"), names(parms.optim)))

    # Formation fractions are attributed to the target variable
    n.ff.optim <- length(grep(paste("^f", ".*", obs_var, "$", sep=""), names(parms.optim)))

    n.optim <- n.k.optim + n.k1.dfop.optim + n.k2.dfop.optim + n.g.dfop.optim + n.initials.optim + n.ff.optim

    # FOMC, HS and IORE parameters are only counted if we are looking at the
    # first variable in the model which is always the source variable
    if (obs_var == fit$obs_vars[[1]]) {
      if ("alpha" %in% names(parms.optim)) n.optim <- n.optim + 1
      if ("beta" %in% names(parms.optim)) n.optim <- n.optim + 1
      if ("k1" %in% names(parms.optim)) n.optim <- n.optim + 1
      if ("k2" %in% names(parms.optim)) n.optim <- n.optim + 1
      if ("tb" %in% names(parms.optim)) n.optim <- n.optim + 1
      if ("N" %in% names(parms.optim)) n.optim <- n.optim + 1
    }

    # Calculate and add a line to the results
    errmin.tmp <- kinerrmin(errdata.var, n.optim)
    errmin[obs_var, c("err.min", "n.optim", "df")] <- errmin.tmp
  }

  return(errmin)
}

Contact - Imprint