diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/endpoints.R | 28 | ||||
| -rw-r--r-- | R/mkinerrmin.R | 8 | ||||
| -rw-r--r-- | R/mkinmod.R | 26 | ||||
| -rw-r--r-- | R/mkinparplot.R | 8 | ||||
| -rw-r--r-- | R/mkinpredict.R | 4 |
5 files changed, 56 insertions, 18 deletions
diff --git a/R/endpoints.R b/R/endpoints.R index 676831a0..ae9aa3ec 100644 --- a/R/endpoints.R +++ b/R/endpoints.R @@ -1,11 +1,11 @@ endpoints <- function(fit) {
# Calculate dissipation times DT50 and DT90 and formation
# fractions as well as SFORB eigenvalues from optimised parameters
- # Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from HS and DFOP,
- # as well as from Eigenvalues b1 and b2 of any SFORB models
+ # Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from
+ # HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models
ep <- list()
obs_vars <- fit$obs_vars
- parms.all <- fit$bparms.ode
+ parms.all <- c(fit$bparms.optim, fit$bparms.fixed)
ep$ff <- vector()
ep$SFORB <- vector()
ep$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)),
@@ -49,6 +49,28 @@ endpoints <- function(fit) { DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
ep$distimes[obs_var, c("DT50back")] = DT50_back
}
+ if (type == "IORE") {
+ k_names = grep(paste("k.iore", obs_var, sep="_"), names(parms.all), value=TRUE)
+ k_tot = sum(parms.all[k_names])
+ # From the NAFTA kinetics guidance, p. 5
+ n = parms.all[paste("N", obs_var, sep = "_")]
+ k = k_tot
+ # Use the initial concentration of the parent compound
+ source_name = fit$mkinmod$map[[1]][[1]]
+ c0 = parms.all[paste(source_name, "0", sep = "_")]
+ alpha = 1 / (n - 1)
+ beta = (c0^(1 - n))/(k * (n - 1))
+ DT50 = beta * (2^(1/alpha) - 1)
+ DT90 = beta * (10^(1/alpha) - 1)
+ DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
+ ep$distimes[obs_var, c("DT50back")] = DT50_back
+ if (fit$mkinmod$use_of_ff == "min") {
+ for (k_name in k_names)
+ {
+ ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot
+ }
+ }
+ }
if (type == "DFOP") {
k1 = parms.all["k1"]
k2 = parms.all["k2"]
diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R index 671bcaab..d89a2f91 100644 --- a/R/mkinerrmin.R +++ b/R/mkinerrmin.R @@ -55,13 +55,15 @@ mkinerrmin <- function(fit, alpha = 0.05) # Check if initial value is optimised
n.initials.optim <- length(grep(paste(obs_var, ".*", "_0", sep=""), names(parms.optim)))
- # Rate constants are attributed to the source variable
+ # Rate constants and IORE exponents are attributed to the source variable
n.k.optim <- length(grep(paste("^k", obs_var, sep="_"), names(parms.optim)))
+ n.k.iore.optim <- length(grep(paste("^k.iore", obs_var, sep="_"), names(parms.optim)))
+ n.N.optim <- length(grep(paste("^N", 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.initials.optim + n.ff.optim
+ n.optim <- sum(n.initials.optim, n.k.optim, n.k.iore.optim, n.N.optim, n.ff.optim)
# FOMC, DFOP and HS parameters are only counted if we are looking at the
# first variable in the model which is always the source variable
@@ -74,7 +76,7 @@ mkinerrmin <- function(fit, alpha = 0.05) if ("tb" %in% names(parms.optim)) n.optim <- n.optim + 1
}
- # Calculate and add a line to the results
+ # Calculate and add a line to the dataframe holding the results
errmin.tmp <- kinerrmin(errdata.var, n.optim)
errmin[obs_var, c("err.min", "n.optim", "df")] <- errmin.tmp
}
diff --git a/R/mkinmod.R b/R/mkinmod.R index 321887fc..8b86303f 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -1,4 +1,4 @@ -# Copyright (C) 2010-2013 Johannes Ranke {{{
+# Copyright (C) 2010-2014 Johannes Ranke {{{
# Contact: jranke@uni-bremen.de
# This file is part of the R package mkin
@@ -33,12 +33,10 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) stop("The use of formation fractions 'use_of_ff' can only be 'min' or 'max'")
# The returned model will be a list of character vectors, containing {{{
- # differential equations, parameter names and a mapping from model variables
- # to observed variables. If possible, a matrix representation of the
- # differential equations is included
+ # differential equations (if supported), parameter names and a mapping from
+ # model variables to observed variables. If possible, a matrix representation
+ # of the differential equations is included
parms <- vector()
- diffs <- vector()
- map <- list()
# }}}
# Do not return a coefficient matrix mat when FOMC, IORE, DFOP or HS is used for the parent {{{
@@ -49,6 +47,8 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) # Establish a list of differential equations as well as a map from observed {{{
# compartments to differential equations
+ diffs <- vector()
+ map <- list()
for (varname in obs_vars)
{
# Check the type component of the compartment specification {{{
@@ -59,7 +59,8 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS") & match(varname, obs_vars) != 1) {
stop(paste("Types FOMC, DFOP and HS are only implemented for the first compartment,",
"which is assumed to be the source compartment"))
- } #}}}
+ }
+ #}}}
# New (sub)compartments (boxes) needed for the model type {{{
new_boxes <- switch(spec[[varname]]$type,
SFO = varname,
@@ -88,7 +89,11 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) if(spec[[varname]]$type %in% c("SFO", "IORE", "SFORB")) { # {{{ Add decline term
if (use_of_ff == "min") { # Minimum use of formation fractions
if(spec[[varname]]$sink) {
- # If sink is required, add first-order sink term
+ if(spec[[varname]]$type == "IORE" && length(spec[[varname]]$to) > 0) {
+ stop("Transformation reactions from compounds modelled with IORE\n",
+ "are only supported with formation fractions (use_of_ff = 'max')")
+ }
+ # If sink is required, add first-order/IORE sink term
k_compound_sink <- paste("k", box_1, "sink", sep = "_")
if(spec[[varname]]$type == "IORE") {
k_compound_sink <- paste("k.iore", box_1, "sink", sep = "_")
@@ -173,9 +178,10 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) for (target in to) {
target_box <- switch(spec[[target]]$type,
SFO = target,
+ IORE = target,
SFORB = paste(target, "free", sep = "_"))
- if (use_of_ff == "min" && spec[[varname]]$type %in% c("SFO",
- "SFORB")) {
+ if (use_of_ff == "min" && spec[[varname]]$type %in% c("SFO", "SFORB"))
+ {
k_from_to <- paste("k", origin_box, target_box, sep = "_")
parms <- c(parms, k_from_to)
diffs[[origin_box]] <- paste(diffs[[origin_box]], "-",
diff --git a/R/mkinparplot.R b/R/mkinparplot.R index 683a8b24..4abecab9 100644 --- a/R/mkinparplot.R +++ b/R/mkinparplot.R @@ -19,11 +19,13 @@ mkinparplot <- function(object) { state.optim = rownames(subset(object$start, type == "state")) deparms.optim = rownames(subset(object$start, type == "deparm")) fractions.optim = grep("^f_", deparms.optim, value = TRUE) + N.optim = grep("^N_", deparms.optim, value = TRUE) if ("g" %in% deparms.optim) fractions.optim <- c("g", fractions.optim) - rates.optim.unsorted = setdiff(deparms.optim, fractions.optim) + rates.optim.unsorted = setdiff(deparms.optim, union(fractions.optim, N.optim)) rates.optim <- rownames(object$start[rates.optim.unsorted, ]) n.plot <- c(state.optim = length(state.optim), rates.optim = length(rates.optim), + N.optim = length(N.optim), fractions.optim = length(fractions.optim)) n.plot <- n.plot[n.plot > 0] @@ -42,6 +44,8 @@ mkinparplot <- function(object) { na.rm = TRUE, finite = TRUE), rates.optim = range(c(0, unlist(values)), na.rm = TRUE, finite = TRUE), + N.optim = range(c(0, 1, unlist(values)), + na.rm = TRUE, finite = TRUE), fractions.optim = range(c(0, 1, unlist(values)), na.rm = TRUE, finite = TRUE)) stripchart(values["Estimate", ][length(parnames):1], @@ -49,7 +53,7 @@ mkinparplot <- function(object) { ylim = c(0.5, length(get(type)) + 0.5), yaxt = "n") if (type %in% c("rates.optim", "fractions.optim")) abline(v = 0, lty = 2) - if (type %in% c("fractions.optim")) abline(v = 1, lty = 2) + if (type %in% c("N.optim", "fractions.optim")) abline(v = 1, lty = 2) position <- ifelse(values["Estimate", ] < mean(xlim), "right", "left") text(ifelse(position == "left", min(xlim), max(xlim)), length(parnames):1, parnames, diff --git a/R/mkinpredict.R b/R/mkinpredict.R index 43b5d2f3..da013d50 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -109,6 +109,10 @@ mkinpredict <- function(mkinmod, odeparms, odeini, rtol = rtol, ... ) + if (sum(is.na(out)) > 0) { + stop("Differential equations were not integrated for all output times because\n", + "NaN values occurred in output from ode()") + } } if (map_output) { # Output transformation for models with unobserved compartments like SFORB |