aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
Diffstat (limited to 'R')
-rw-r--r--R/DFOP.solution.R5
-rw-r--r--R/FOMC.solution.R4
-rw-r--r--R/HS.solution.R6
-rw-r--r--R/SFO.solution.R4
-rw-r--r--R/SFORB.solution.R9
-rw-r--r--R/endpoints.R108
-rw-r--r--R/gmkin.R10
-rw-r--r--R/ilr.R51
-rw-r--r--R/mkin_long_to_wide.R30
-rw-r--r--R/mkin_wide_to_long.R34
-rw-r--r--R/mkinerrmin.R85
-rw-r--r--R/mkinfit.R482
-rw-r--r--R/mkinmod.R242
-rw-r--r--R/mkinparplot.R62
-rw-r--r--R/mkinplot.R4
-rw-r--r--R/mkinpredict.R121
-rw-r--r--R/mkinresplot.R53
-rw-r--r--R/plot.mkinfit.R98
-rw-r--r--R/transform_odeparms.R101
19 files changed, 1 insertions, 1508 deletions
diff --git a/R/DFOP.solution.R b/R/DFOP.solution.R
deleted file mode 100644
index 8531cfe..0000000
--- a/R/DFOP.solution.R
+++ /dev/null
@@ -1,5 +0,0 @@
-DFOP.solution <- function(t, parent.0, k1, k2, g)
-{
- parent = g * parent.0 * exp(-k1 * t) +
- (1 - g) * parent.0 * exp(-k2 * t)
-}
diff --git a/R/FOMC.solution.R b/R/FOMC.solution.R
deleted file mode 100644
index 8bd13d6..0000000
--- a/R/FOMC.solution.R
+++ /dev/null
@@ -1,4 +0,0 @@
-FOMC.solution <- function(t, parent.0, alpha, beta)
-{
- parent = parent.0 / (t/beta + 1)^alpha
-}
diff --git a/R/HS.solution.R b/R/HS.solution.R
deleted file mode 100644
index 4651a6a..0000000
--- a/R/HS.solution.R
+++ /dev/null
@@ -1,6 +0,0 @@
-HS.solution <- function(t, parent.0, k1, k2, tb)
-{
- parent = ifelse(t <= tb,
- parent.0 * exp(-k1 * t),
- parent.0 * exp(-k1 * tb) * exp(-k2 * (t - tb)))
-}
diff --git a/R/SFO.solution.R b/R/SFO.solution.R
deleted file mode 100644
index 3a376e4..0000000
--- a/R/SFO.solution.R
+++ /dev/null
@@ -1,4 +0,0 @@
-SFO.solution <- function(t, parent.0, k)
-{
- parent = parent.0 * exp(-k * t)
-}
diff --git a/R/SFORB.solution.R b/R/SFORB.solution.R
deleted file mode 100644
index 45a4533..0000000
--- a/R/SFORB.solution.R
+++ /dev/null
@@ -1,9 +0,0 @@
-SFORB.solution = function(t, parent.0, k_12, k_21, k_1output) {
- sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21)
- b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
- b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp
-
- parent = parent.0 *
- (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +
- ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t))
-}
diff --git a/R/endpoints.R b/R/endpoints.R
deleted file mode 100644
index c9a6a51..0000000
--- a/R/endpoints.R
+++ /dev/null
@@ -1,108 +0,0 @@
-endpoints <- function(fit, pseudoDT50 = FALSE) {
- # Calculate dissipation times DT50 and DT90 and, if necessary, formation
- # fractions and SFORB eigenvalues from optimised parameters
- ep <- list()
- obs_vars <- fit$obs_vars
- parms.all <- fit$bparms.ode
- ep$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)),
- DT90 = rep(NA, length(obs_vars)),
- row.names = obs_vars)
- ep$ff <- vector()
- ep$SFORB <- vector()
- for (obs_var in obs_vars) {
- type = names(fit$mkinmod$map[[obs_var]])[1]
-
- # Get formation fractions if directly fitted, and calculate remaining fraction to sink
- f_names = grep(paste("f", obs_var, sep = "_"), names(parms.all), value=TRUE)
- f_values = parms.all[f_names]
- f_to_sink = 1 - sum(f_values)
- names(f_to_sink) = ifelse(type == "SFORB",
- paste(obs_var, "free", "sink", sep = "_"),
- paste(obs_var, "sink", sep = "_"))
- for (f_name in f_names) {
- ep$ff[[sub("f_", "", sub("_to_", "_", f_name))]] = f_values[[f_name]]
- }
- ep$ff = append(ep$ff, f_to_sink)
-
- # Get the rest
- if (type == "SFO") {
- k_names = grep(paste("k", obs_var, sep="_"), names(parms.all), value=TRUE)
- k_tot = sum(parms.all[k_names])
- DT50 = log(2)/k_tot
- DT90 = log(10)/k_tot
- for (k_name in k_names)
- {
- ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot
- }
- }
- if (type == "FOMC") {
- alpha = parms.all["alpha"]
- beta = parms.all["beta"]
- DT50 = beta * (2^(1/alpha) - 1)
- DT90 = beta * (10^(1/alpha) - 1)
- }
- if (type == "DFOP") {
- k1 = parms.all["k1"]
- k2 = parms.all["k2"]
- g = parms.all["g"]
- f <- function(t, x) {
- fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t)
- (fraction - (1 - x/100))^2
- }
- DTmax <- 1000
- DT50.o <- optimize(f, c(0.001, DTmax), x=50)$minimum
- DT50 = ifelse(DTmax - DT50.o < 0.1, NA, DT50.o)
- DT90.o <- optimize(f, c(0.001, DTmax), x=90)$minimum
- DT90 = ifelse(DTmax - DT90.o < 0.1, NA, DT90.o)
- }
- if (type == "HS") {
- k1 = parms.all["k1"]
- k2 = parms.all["k2"]
- tb = parms.all["tb"]
- DTx <- function(x) {
- DTx.a <- (log(100/(100 - x)))/k1
- DTx.b <- tb + (log(100/(100 - x)) - k1 * tb)/k2
- if (DTx.a < tb) DTx <- DTx.a
- else DTx <- DTx.b
- return(DTx)
- }
- DT50 <- DTx(50)
- DT90 <- DTx(90)
- }
- if (type == "SFORB") {
- # FOCUS kinetics (2006), p. 60 f
- k_out_names = grep(paste("k", obs_var, "free", sep="_"), names(parms.all), value=TRUE)
- k_out_names = setdiff(k_out_names, paste("k", obs_var, "free", "bound", sep="_"))
- k_1output = sum(parms.all[k_out_names])
- k_12 = parms.all[paste("k", obs_var, "free", "bound", sep="_")]
- k_21 = parms.all[paste("k", obs_var, "bound", "free", sep="_")]
-
- sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21)
- b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
- b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp
-
- SFORB_fraction = function(t) {
- ((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +
- ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)
- }
- f_50 <- function(t) (SFORB_fraction(t) - 0.5)^2
- max_DT <- 1000
- DT50.o <- optimize(f_50, c(0.01, max_DT))$minimum
- if (abs(DT50.o - max_DT) < 0.01) DT50 = NA else DT50 = DT50.o
- f_90 <- function(t) (SFORB_fraction(t) - 0.1)^2
- DT90.o <- optimize(f_90, c(0.01, max_DT))$minimum
- if (abs(DT90.o - max_DT) < 0.01) DT90 = NA else DT90 = DT90.o
-
- for (k_out_name in k_out_names)
- {
- ep$ff[[sub("k_", "", k_out_name)]] = parms.all[[k_out_name]] / k_1output
- }
-
- # Return the eigenvalues for comparison with DFOP rate constants
- ep$SFORB[[paste(obs_var, "b1", sep="_")]] = b1
- ep$SFORB[[paste(obs_var, "b2", sep="_")]] = b2
- }
- ep$distimes[obs_var, ] = c(DT50, DT90)
- }
- return(ep)
-}
diff --git a/R/gmkin.R b/R/gmkin.R
index 912d31c..31e8cf4 100644
--- a/R/gmkin.R
+++ b/R/gmkin.R
@@ -1,11 +1,3 @@
gmkin <- function() {
- if (require(gWidgetsWWW2)) load_app(system.file("GUI/gmkin.R", package = "mkin"))
- else {
- message(
- "\nYou need to install gWidgetsWWW2 in order to run the mkin GUI gmkin.\n",
- "This package is not currently on CRAN but can be installed from github\n",
- "using the devtools package:\n",
- "library(devtools)\n",
- "install_github('gWidgetsWWW2', 'jverzani')")
- }
+ load_app(system.file("GUI/gmkin.R", package = "gmkin"))
}
diff --git a/R/ilr.R b/R/ilr.R
deleted file mode 100644
index 389653e..0000000
--- a/R/ilr.R
+++ /dev/null
@@ -1,51 +0,0 @@
-# $Id$
-
-# Copyright (C) 2012 René Lehmann, 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/>
-
-ilr <- function(x) {
- z <- vector()
- for (i in 1:(length(x) - 1)) {
- z[i] <- sqrt(i/(i+1)) * log((prod(x[1:i]))^(1/i) / x[i+1])
- }
- return(z)
-}
-
-invilr<-function(x) {
- D <- length(x) + 1
- z <- c(x, 0)
- y <- rep(0, D)
- s <- sqrt(1:D*2:(D+1))
- q <- z/s
- y[1] <- sum(q[1:D])
- for (i in 2:D) {
- y[i] <- sum(q[i:D]) - sqrt((i-1)/i) * z[i-1]
- }
- z <- vector()
- for (i in 1:D) {
- z[i] <- exp(y[i])/sum(exp(y))
- }
-
- # Work around a numerical problem with NaN values returned by the above
- # Only works if there is only one NaN value: replace it with 1
- # if the sum of the other components is < 1e-10
- if (sum(is.na(z)) == 1 && sum(z[!is.na(z)]) < 1e-10)
- z = ifelse(is.na(z), 1, z)
-
- return(z)
-}
diff --git a/R/mkin_long_to_wide.R b/R/mkin_long_to_wide.R
deleted file mode 100644
index 87e1afd..0000000
--- a/R/mkin_long_to_wide.R
+++ /dev/null
@@ -1,30 +0,0 @@
-# $Id$
-
-# Copyright (C) 2010-2011 Johannes Ranke
-# Contact: mkin-devel@lists.berlios.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/>
-
-mkin_long_to_wide <- function(long_data, time = "time", outtime = "time")
-{
- colnames <- unique(long_data$name)
- wide_data <- data.frame(time = subset(long_data, name == colnames[1], time))
- names(wide_data) <- outtime
- for (var in colnames) {
- wide_data[var] <- subset(long_data, name == var, value)
- }
- return(wide_data)
-}
diff --git a/R/mkin_wide_to_long.R b/R/mkin_wide_to_long.R
deleted file mode 100644
index f1814fc..0000000
--- a/R/mkin_wide_to_long.R
+++ /dev/null
@@ -1,34 +0,0 @@
-# $Id$
-
-# Copyright (C) 2010-2013 Johannes Ranke
-# Contact: mkin-devel@lists.berlios.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/>
-if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "value"))
-
-mkin_wide_to_long <- function(wide_data, time = "t")
-{
- colnames <- names(wide_data)
- if (!(time %in% colnames)) stop("The data in wide format have to contain a variable named ", time, ".")
- vars <- subset(colnames, colnames != time)
- n <- length(colnames) - 1
- long_data <- data.frame(
- name = rep(vars, each = length(wide_data[[time]])),
- time = as.numeric(rep(wide_data[[time]], n)),
- value = as.numeric(unlist(wide_data[vars])),
- row.names = NULL)
- return(long_data)
-}
diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R
deleted file mode 100644
index 074cc56..0000000
--- a/R/mkinerrmin.R
+++ /dev/null
@@ -1,85 +0,0 @@
-# $Id$
-
-# Copyright (C) 2010-2013 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/>
-if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "value_mean"))
-
-mkinerrmin <- 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), ]
-
- # Any value that is set to exactly zero is not really an observed value
- # Remove those at time 0 - those are caused by the FOCUS recommendation
- # to set metabolites occurring at time 0 to 0
- errdata <- subset(errdata, !(time == 0 & value_mean == 0))
-
- 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 are attributed to the source variable
- n.k.optim <- length(grep(paste("^k", 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
-
- # 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
- 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 ("g" %in% names(parms.optim)) n.optim <- n.optim + 1
- if ("tb" %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)
-}
diff --git a/R/mkinfit.R b/R/mkinfit.R
deleted file mode 100644
index f44c67d..0000000
--- a/R/mkinfit.R
+++ /dev/null
@@ -1,482 +0,0 @@
-# Copyright (C) 2010-2014 Johannes Ranke
-# Portions of this code are copyright (C) 2013 Eurofins Regulatory AG
-# Contact: jranke@uni-bremen.de
-# The summary function is an adapted and extended version of summary.modFit
-# from the FME package, v 1.1 by Soetart and Petzoldt, which was in turn
-# inspired by summary.nls.lm
-
-# 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/>
-if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "value"))
-
-mkinfit <- function(mkinmod, observed,
- parms.ini = "auto",
- state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)),
- fixed_parms = NULL,
- fixed_initials = names(mkinmod$diffs)[-1],
- solution_type = "auto",
- method.ode = "lsoda",
- method.modFit = "Marq",
- control.modFit = list(),
- transform_rates = TRUE,
- transform_fractions = TRUE,
- plot = FALSE, quiet = FALSE,
- err = NULL, weight = "none", scaleVar = FALSE,
- atol = 1e-8, rtol = 1e-10, n.outtimes = 100,
- reweight.method = NULL,
- reweight.tol = 1e-8, reweight.max.iter = 10,
- trace_parms = FALSE,
- ...)
-{
- # Get the names of the state variables in the model
- mod_vars <- names(mkinmod$diffs)
-
- # Get the names of observed variables
- obs_vars = names(mkinmod$spec)
-
- # Subset observed data with names of observed data in the model
- observed <- subset(observed, name %in% obs_vars)
-
- # Define starting values for parameters where not specified by the user
- if (parms.ini[[1]] == "auto") parms.ini = vector()
-
- # Prevent inital parameter specifications that are not in the model
- wrongpar.names <- setdiff(names(parms.ini), mkinmod$parms)
- if (length(wrongpar.names) > 0) {
- stop("Initial parameter(s) ", paste(wrongpar.names, collapse = ", "),
- " not used in the model")
- }
-
- k_salt = 0
- defaultpar.names <- setdiff(mkinmod$parms, names(parms.ini))
- for (parmname in defaultpar.names) {
- # Default values for rate constants, depending on the parameterisation
- if (substr(parmname, 1, 2) == "k_") {
- parms.ini[parmname] = 0.1 + k_salt
- k_salt = k_salt + 1e-4
- }
- # Default values for rate constants for reversible binding
- if (grepl("free_bound$", parmname)) parms.ini[parmname] = 0.1
- if (grepl("bound_free$", parmname)) parms.ini[parmname] = 0.02
- # Default values for formation fractions
- if (substr(parmname, 1, 2) == "f_") parms.ini[parmname] = 0.2
- # Default values for the FOMC, DFOP and HS models
- if (parmname == "alpha") parms.ini[parmname] = 1
- if (parmname == "beta") parms.ini[parmname] = 10
- if (parmname == "k1") parms.ini[parmname] = 0.1
- if (parmname == "k2") parms.ini[parmname] = 0.01
- if (parmname == "tb") parms.ini[parmname] = 5
- if (parmname == "g") parms.ini[parmname] = 0.5
- }
-
- # Name the inital state variable values if they are not named yet
- if(is.null(names(state.ini))) names(state.ini) <- mod_vars
-
- # Transform initial parameter values for fitting
- transparms.ini <- transform_odeparms(parms.ini, mod_vars,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
-
- # Parameters to be optimised:
- # Kinetic parameters in parms.ini whose names are not in fixed_parms
- parms.fixed <- transparms.ini[fixed_parms]
- parms.optim <- transparms.ini[setdiff(names(transparms.ini), fixed_parms)]
-
- # Inital state variables in state.ini whose names are not in fixed_initials
- state.ini.fixed <- state.ini[fixed_initials]
- state.ini.optim <- state.ini[setdiff(names(state.ini), fixed_initials)]
-
- # Preserve names of state variables before renaming initial state variable
- # parameters
- state.ini.optim.boxnames <- names(state.ini.optim)
- state.ini.fixed.boxnames <- names(state.ini.fixed)
- if(length(state.ini.optim) > 0) {
- names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_")
- }
- if(length(state.ini.fixed) > 0) {
- names(state.ini.fixed) <- paste(names(state.ini.fixed), "0", sep="_")
- }
-
- # Decide if the solution of the model can be based on a simple analytical
- # formula, the spectral decomposition of the matrix (fundamental system)
- # or a numeric ode solver from the deSolve package
- if (!solution_type %in% c("auto", "analytical", "eigen", "deSolve"))
- stop("solution_type must be auto, analytical, eigen or de Solve")
- if (solution_type == "analytical" && length(mkinmod$spec) > 1)
- stop("Analytical solution not implemented for models with metabolites.")
- if (solution_type == "eigen" && !is.matrix(mkinmod$coefmat))
- stop("Eigenvalue based solution not possible, coefficient matrix not present.")
- if (solution_type == "auto") {
- if (length(mkinmod$spec) == 1) {
- solution_type = "analytical"
- } else {
- if (is.matrix(mkinmod$coefmat)) {
- solution_type = "eigen"
- if (max(observed$value, na.rm = TRUE) < 0.1) {
- stop("The combination of small observed values (all < 0.1) and solution_type = eigen is error-prone")
- }
- } else {
- solution_type = "deSolve"
- }
- }
- }
-
- cost.old <- 1e100 # The first model cost should be smaller than this value
- calls <- 0 # Counter for number of model solutions
- out_predicted <- NA
- # Define the model cost function
- cost <- function(P)
- {
- assign("calls", calls+1, inherits=TRUE) # Increase the model solution counter
-
- # Trace parameter values if requested
- if(trace_parms) cat(P, "\n")
-
- # Time points at which observed data are available
- # Make sure we include time 0, so initial values for state variables are for time 0
- outtimes = sort(unique(c(observed$time, seq(min(observed$time),
- max(observed$time),
- length.out = n.outtimes))))
-
- if(length(state.ini.optim) > 0) {
- odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed)
- names(odeini) <- c(state.ini.optim.boxnames, state.ini.fixed.boxnames)
- } else {
- odeini <- state.ini.fixed
- names(odeini) <- state.ini.fixed.boxnames
- }
-
- odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed)
-
- parms <- backtransform_odeparms(odeparms, mod_vars,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
-
- # Solve the system with current transformed parameter values
- out <- mkinpredict(mkinmod, parms, odeini, outtimes,
- solution_type = solution_type,
- method.ode = method.ode,
- atol = atol, rtol = rtol, ...)
-
- assign("out_predicted", out, inherits=TRUE)
-
- mC <- modCost(out, observed, y = "value",
- err = err, weight = weight, scaleVar = scaleVar)
-
- # Report and/or plot if the model is improved
- if (mC$model < cost.old) {
- if(!quiet) cat("Model cost at call ", calls, ": ", mC$model, "\n")
-
- # Plot the data and current model output if requested
- if(plot) {
- outtimes_plot = seq(min(observed$time), max(observed$time), length.out=100)
-
- out_plot <- mkinpredict(mkinmod, parms, odeini, outtimes_plot,
- solution_type = solution_type,
- method.ode = method.ode,
- atol = atol, rtol = rtol, ...)
-
- plot(0, type="n",
- xlim = range(observed$time), ylim = range(observed$value, na.rm=TRUE),
- xlab = "Time", ylab = "Observed")
- col_obs <- pch_obs <- 1:length(obs_vars)
- lty_obs <- rep(1, length(obs_vars))
- names(col_obs) <- names(pch_obs) <- names(lty_obs) <- obs_vars
- for (obs_var in obs_vars) {
- points(subset(observed, name == obs_var, c(time, value)),
- pch = pch_obs[obs_var], col = col_obs[obs_var])
- }
- matlines(out_plot$time, out_plot[-1], col = col_obs, lty = lty_obs)
- legend("topright", inset=c(0.05, 0.05), legend=obs_vars,
- col=col_obs, pch=pch_obs, lty=1:length(pch_obs))
- }
-
- assign("cost.old", mC$model, inherits=TRUE)
- }
- return(mC)
- }
-
- lower <- rep(-Inf, length(c(state.ini.optim, parms.optim)))
- names(lower) <- c(names(state.ini.optim), names(parms.optim))
- if (!transform_rates) {
- index_k <- grep("^k_", names(lower))
- lower[index_k] <- 0
- other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2"), names(lower))
- lower[other_rate_parms] <- 0
- }
-
- fit <- modFit(cost, c(state.ini.optim, parms.optim),
- method = method.modFit, control = control.modFit, lower = lower, ...)
-
- # Reiterate the fit until convergence of the variance components (IRLS)
- # if requested by the user
- weight.ini = weight
- if (!is.null(err)) weight.ini = "manual"
-
- if (!is.null(reweight.method)) {
- if (reweight.method != "obs") stop("Only reweighting method 'obs' is implemented")
- if(!quiet) {
- cat("IRLS based on variance estimates for each observed variable\n")
- }
- if (!quiet) {
- cat("Initial variance estimates are:\n")
- print(signif(fit$var_ms_unweighted, 8))
- }
- reweight.diff = 1
- n.iter <- 0
- if (!is.null(err)) observed$err.ini <- observed[[err]]
- err = "err.irls"
- while (reweight.diff > reweight.tol & n.iter < reweight.max.iter) {
- n.iter <- n.iter + 1
- sigma.old <- sqrt(fit$var_ms_unweighted)
- observed[err] <- sqrt(fit$var_ms_unweighted)[as.character(observed$name)]
- fit <- modFit(cost, fit$par, method = method.modFit,
- control = control.modFit, lower = lower, ...)
- reweight.diff = sum((sqrt(fit$var_ms_unweighted) - sigma.old)^2)
- if (!quiet) {
- cat("Iteration", n.iter, "yields variance estimates:\n")
- print(signif(fit$var_ms_unweighted, 8))
- cat("Sum of squared differences to last variance estimates:",
- signif(reweight.diff, 2), "\n")
- }
- }
- }
-
- # We need to return some more data for summary and plotting
- fit$solution_type <- solution_type
- fit$transform_rates <- transform_rates
- fit$transform_fractions <- transform_fractions
-
- # We also need the model for summary and plotting
- fit$mkinmod <- mkinmod
-
- # We need data and predictions for summary and plotting
- fit$observed <- observed
- fit$obs_vars <- obs_vars
- fit$predicted <- mkin_wide_to_long(out_predicted, time = "time")
-
- # Collect initial parameter values in two dataframes
- fit$start <- data.frame(value = c(state.ini.optim,
- backtransform_odeparms(parms.optim, mod_vars,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)))
- fit$start$type = c(rep("state", length(state.ini.optim)),
- rep("deparm", length(parms.optim)))
- fit$start$transformed = c(state.ini.optim, parms.optim)
-
- fit$fixed <- data.frame(value = c(state.ini.fixed,
- backtransform_odeparms(parms.fixed, mod_vars,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)))
- fit$fixed$type = c(rep("state", length(state.ini.fixed)),
- rep("deparm", length(parms.fixed)))
-
- bparms.optim = backtransform_odeparms(fit$par, mod_vars,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
- bparms.fixed = backtransform_odeparms(c(state.ini.fixed, parms.fixed),
- mod_vars,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
- bparms.all = c(bparms.optim, bparms.fixed)
-
- # Collect observed, predicted, residuals and weighting
- data <- merge(fit$observed, fit$predicted, by = c("time", "name"))
- data$name <- ordered(data$name, levels = obs_vars)
- data <- data[order(data$name, data$time), ]
-
- fit$data <- data.frame(time = data$time,
- variable = data$name,
- observed = data$value.x,
- predicted = data$value.y)
- fit$data$residual <- fit$data$observed - fit$data$predicted
- if (!is.null(data$err.ini)) fit$data$err.ini <- data$err.ini
- if (!is.null(err)) fit$data[[err]] <- data[[err]]
-
- fit$atol <- atol
- fit$rtol <- rtol
- fit$weight.ini <- weight.ini
- fit$reweight.method <- reweight.method
- fit$reweight.tol <- reweight.tol
- fit$reweight.max.iter <- reweight.max.iter
-
- # Return all backtransformed parameters for summary
- fit$bparms.optim <- bparms.optim
- fit$bparms.fixed <- bparms.fixed
- fit$bparms.ode <- bparms.all[mkinmod$parms] # Return ode parameters for further fitting
- fit$date <- date()
-
- class(fit) <- c("mkinfit", "modFit")
- return(fit)
-}
-
-summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...) {
- param <- object$par
- pnames <- names(param)
- p <- length(param)
- mod_vars <- names(object$mkinmod$diffs)
- covar <- try(solve(0.5*object$hessian), silent = TRUE) # unscaled covariance
- rdf <- object$df.residual
- resvar <- object$ssr / rdf
- if (!is.numeric(covar)) {
- covar <- NULL
- se <- lci <- uci <- tval <- pval1 <- pval2 <- rep(NA, p)
- } else {
- rownames(covar) <- colnames(covar) <- pnames
- se <- sqrt(diag(covar) * resvar)
- lci <- param + qt(alpha/2, rdf) * se
- uci <- param + qt(1-alpha/2, rdf) * se
- tval <- param/se
- pval1 <- 2 * pt(abs(tval), rdf, lower.tail = FALSE)
- pval2 <- pt(abs(tval), rdf, lower.tail = FALSE)
- }
-
- names(se) <- pnames
- modVariance <- object$ssr / length(object$residuals)
-
- param <- cbind(param, se, lci, uci, tval, pval1, pval2)
- dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "Lower", "Upper",
- "t value", "Pr(>|t|)", "Pr(>t)"))
-
- blci <- buci <- numeric()
- # Only transform boundaries of CI for one parameter at a time
- for (pname in pnames) {
- par.lower <- par.upper <- object$par
- par.lower[pname] <- param[pname, "Lower"]
- par.upper[pname] <- param[pname, "Upper"]
- blci[pname] <- backtransform_odeparms(par.lower, mod_vars,
- object$transform_rates, object$transform_fractions)[pname]
- buci[pname] <- backtransform_odeparms(par.upper, mod_vars,
- object$transform_rates, object$transform_fractions)[pname]
- }
- bparam <- cbind(object$bparms.optim, blci, buci)
- dimnames(bparam) <- list(pnames, c("Estimate", "Lower", "Upper"))
-
- ans <- list(
- version = as.character(packageVersion("mkin")),
- Rversion = paste(R.version$major, R.version$minor, sep="."),
- date.fit = object$date,
- date.summary = date(),
- solution_type = object$solution_type,
- use_of_ff = object$mkinmod$use_of_ff,
- weight.ini = object$weight.ini,
- reweight.method = object$reweight.method,
- residuals = object$residuals,
- residualVariance = resvar,
- sigma = sqrt(resvar),
- modVariance = modVariance,
- df = c(p, rdf),
- cov.unscaled = covar,
- cov.scaled = covar * resvar,
- info = object$info,
- niter = object$iterations,
- stopmess = message,
- par = param,
- bpar = bparam)
-
- ans$diffs <- object$mkinmod$diffs
- if(data) ans$data <- object$data
- ans$start <- object$start
-
- ans$fixed <- object$fixed
-
- ans$errmin <- mkinerrmin(object, alpha = 0.05)
-
- ans$bparms.ode <- object$bparms.ode
- ep <- endpoints(object)
- if (length(ep$ff) != 0)
- ans$ff <- ep$ff
- if(distimes) ans$distimes <- ep$distimes
- if(length(ep$SFORB) != 0) ans$SFORB <- ep$SFORB
- class(ans) <- c("summary.mkinfit", "summary.modFit")
- return(ans)
-}
-
-# Expanded from print.summary.modFit
-print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), ...) {
- cat("mkin version: ", x$version, "\n")
- cat("R version: ", x$Rversion, "\n")
- cat("Date of fit: ", x$date.fit, "\n")
- cat("Date of summary:", x$date.summary, "\n")
-
- cat("\nEquations:\n")
- print(noquote(as.character(x[["diffs"]])))
- df <- x$df
- rdf <- df[2]
-
- cat("\nMethod used for solution of differential equation system:\n")
- cat(x$solution_type, "\n")
-
- cat("\nWeighting:", x$weight.ini)
- if(!is.null(x$reweight.method)) cat(" then iterative reweighting method",
- x$reweight.method)
- cat("\n")
-
- cat("\nStarting values for optimised parameters:\n")
- print(x$start)
-
- cat("\nFixed parameter values:\n")
- if(length(x$fixed$value) == 0) cat("None\n")
- else print(x$fixed)
-
- cat("\nOptimised, transformed parameters:\n")
- print(signif(x$par, digits = digits))
-
- cat("\nBacktransformed parameters:\n")
- print(signif(x$bpar, digits = digits))
-
- cat("\nResidual standard error:",
- format(signif(x$sigma, digits)), "on", rdf, "degrees of freedom\n")
-
- cat("\nChi2 error levels in percent:\n")
- x$errmin$err.min <- 100 * x$errmin$err.min
- print(x$errmin, digits=digits,...)
-
- printdistimes <- !is.null(x$distimes)
- if(printdistimes){
- cat("\nEstimated disappearance times:\n")
- print(x$distimes, digits=digits,...)
- }
-
- printff <- !is.null(x$ff)
- if(printff & x$use_of_ff == "min"){
- cat("\nEstimated formation fractions:\n")
- print(data.frame(ff = x$ff), digits=digits,...)
- }
-
- printSFORB <- !is.null(x$SFORB)
- if(printSFORB){
- cat("\nEstimated Eigenvalues of SFORB model(s):\n")
- print(x$SFORB, digits=digits,...)
- }
-
- cat("\nParameter correlation:\n")
- if (!is.null(x$cov.unscaled)){
- Corr <- cov2cor(x$cov.unscaled)
- rownames(Corr) <- colnames(Corr) <- rownames(x$par)
- print(Corr, digits = digits, ...)
- } else {
- cat("Could not estimate covariance matrix; singular system:\n")
- }
-
- printdata <- !is.null(x$data)
- if (printdata){
- cat("\nData:\n")
- print(format(x$data, digits = digits, ...), row.names = FALSE)
- }
-
- invisible(x)
-}
-# vim: set ts=2 sw=2 expandtab:
diff --git a/R/mkinmod.R b/R/mkinmod.R
deleted file mode 100644
index 2bc6ae3..0000000
--- a/R/mkinmod.R
+++ /dev/null
@@ -1,242 +0,0 @@
-# Copyright (C) 2010-2013 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/> }}}
-
-mkinmod <- function(..., use_of_ff = "min", speclist = NULL)
-{
- if (is.null(speclist)) spec <- list(...)
- else spec <- speclist
- obs_vars <- names(spec)
-
- # Check if any of the names of the observed variables contains any other
- for (obs_var in obs_vars) {
- if (length(grep(obs_var, obs_vars)) > 1) stop("Sorry, variable names can not contain each other")
- if (grepl("_to_", obs_var)) stop("Sorry, names of observed variables can not contain _to_")
- }
-
- if (!use_of_ff %in% c("min", "max"))
- 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
- parms <- vector()
- diffs <- vector()
- map <- list()
- # }}}
-
- # Give a warning when a model with time dependent degradation uses formation {{{
- # fractions
- if(spec[[1]]$type %in% c("FOMC", "DFOP", "HS")) {
- mat = FALSE
- } else mat = TRUE
- #}}}
-
- # Establish a list of differential equations as well as a map from observed {{{
- # compartments to differential equations
- for (varname in obs_vars)
- {
- # Check the type component of the compartment specification {{{
- if(is.null(spec[[varname]]$type)) stop(
- "Every part of the model specification must be a list containing a type component")
- if(!spec[[varname]]$type %in% c("SFO", "FOMC", "DFOP", "HS", "SFORB")) stop(
- "Available types are SFO, FOMC, DFOP, HS and SFORB only")
- 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,
- FOMC = varname,
- DFOP = varname,
- HS = varname,
- SFORB = paste(varname, c("free", "bound"), sep="_")
- )
- map[[varname]] <- new_boxes
- names(map[[varname]]) <- rep(spec[[varname]]$type, length(new_boxes)) #}}}
- # Start a new differential equation for each new box {{{
- new_diffs <- paste("d_", new_boxes, " =", sep="")
- names(new_diffs) <- new_boxes
- diffs <- c(diffs, new_diffs) #}}}
- } #}}}
-
- # Create content of differential equations and build parameter list {{{
- for (varname in obs_vars)
- {
- # Get the name of the box(es) we are working on for the decline term(s)
- box_1 = map[[varname]][[1]] # This is the only box unless type is SFORB
- # Turn on sink if this is not explicitly excluded by the user by
- # specifying sink=FALSE
- if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE
- if(spec[[varname]]$type %in% c("SFO", "SFORB")) { # {{{ Add SFO or SFORB 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
- k_compound_sink <- paste("k", box_1, "sink", sep="_")
- parms <- c(parms, k_compound_sink)
- decline_term <- paste(k_compound_sink, "*", box_1)
- } else { # otherwise no decline term needed here
- decline_term = "0"
- }
- } else {
- k_compound <- paste("k", box_1, sep="_")
- parms <- c(parms, k_compound)
- decline_term <- paste(k_compound, "*", box_1)
- }
- } #}}}
- if(spec[[varname]]$type == "FOMC") { # {{{ Add FOMC decline term
- # From p. 53 of the FOCUS kinetics report
- decline_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", box_1)
- parms <- c(parms, "alpha", "beta")
- } #}}}
- if(spec[[varname]]$type == "DFOP") { # {{{ Add DFOP decline term
- # From p. 57 of the FOCUS kinetics report
- decline_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", box_1)
- parms <- c(parms, "k1", "k2", "g")
- } #}}}
- if(spec[[varname]]$type == "HS") { # {{{ Add HS decline term
- # From p. 55 of the FOCUS kinetics report
- decline_term <- paste("ifelse(time <= tb, k1, k2)", "*", box_1)
- parms <- c(parms, "k1", "k2", "tb")
- } #}}}
- # Add origin decline term to box 1 (usually the only box, unless type is SFORB)#{{{
- diffs[[box_1]] <- paste(diffs[[box_1]], "-", decline_term)#}}}
- if(spec[[varname]]$type == "SFORB") { # {{{ Add SFORB reversible binding terms
- box_2 = map[[varname]][[2]]
- if (use_of_ff == "min") { # Minimum use of formation fractions
- k_free_bound <- paste("k", varname, "free", "bound", sep="_")
- k_bound_free <- paste("k", varname, "bound", "free", sep="_")
- parms <- c(parms, k_free_bound, k_bound_free)
- reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+",
- k_bound_free, "*", box_2)
- reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-",
- k_bound_free, "*", box_2)
- } else { # Use formation fractions also for the free compartment
- stop("The maximum use of formation fractions is not supported for SFORB models")
- # The problems were: Calculation of dissipation times did not work in this case
- # and the coefficient matrix is not generated correctly by the code present
- # in this file in this case
- f_free_bound <- paste("f", varname, "free", "bound", sep="_")
- k_bound_free <- paste("k", varname, "bound", "free", sep="_")
- parms <- c(parms, f_free_bound, k_bound_free)
- reversible_binding_term_1 <- paste("+", k_bound_free, "*", box_2)
- reversible_binding_term_2 <- paste("+", f_free_bound, "*", k_compound, "*", box_1, "-",
- k_bound_free, "*", box_2)
- }
- diffs[[box_1]] <- paste(diffs[[box_1]], reversible_binding_term_1)
- diffs[[box_2]] <- paste(diffs[[box_2]], reversible_binding_term_2)
- } #}}}
-
- # Transfer between compartments#{{{
- to <- spec[[varname]]$to
- if(!is.null(to)) {
- # Name of box from which transfer takes place
- origin_box <- box_1
-
- # Add transfer terms to listed compartments
- for (target in to) {
- target_box <- switch(spec[[target]]$type,
- SFO = target,
- SFORB = paste(target, "free", sep="_"))
- 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]], "-",
- k_from_to, "*", origin_box)
- diffs[[target_box]] <- paste(diffs[[target_box]], "+",
- k_from_to, "*", origin_box)
- } else {
- if (!spec[[varname]]$sink) {
- stop("Turning off the sink when using formation fractions is not supported")
- }
- fraction_to_target = paste("f", origin_box, "to", target, sep="_")
- parms <- c(parms, fraction_to_target)
- diffs[[target_box]] <- paste(diffs[[target_box]], "+",
- fraction_to_target, "*", decline_term)
- }
- }
- } #}}}
- } #}}}
-
- model <- list(diffs = diffs, parms = parms, map = map, spec = spec, use_of_ff = use_of_ff)
-
- # Create coefficient matrix if appropriate#{{{
- if (mat) {
- boxes <- names(diffs)
- n <- length(boxes)
- m <- matrix(nrow=n, ncol=n, dimnames=list(boxes, boxes))
-
- if (use_of_ff == "min") { # {{{ Minimum use of formation fractions
- for (from in boxes) {
- for (to in boxes) {
- if (from == to) { # diagonal elements
- k.candidate = paste("k", from, c(boxes, "sink"), sep="_")
- k.candidate = sub("free.*bound", "free_bound", k.candidate)
- k.candidate = sub("bound.*free", "bound_free", k.candidate)
- k.effective = intersect(model$parms, k.candidate)
- m[from,to] = ifelse(length(k.effective) > 0,
- paste("-", k.effective, collapse = " "), "0")
-
- } else { # off-diagonal elements
- k.candidate = paste("k", from, to, sep="_")
- if (sub("_free$", "", from) == sub("_bound$", "", to)) {
- k.candidate = paste("k", sub("_free$", "_free_bound", from), sep="_")
- }
- if (sub("_bound$", "", from) == sub("_free$", "", to)) {
- k.candidate = paste("k", sub("_bound$", "_bound_free", from), sep="_")
- }
- k.effective = intersect(model$parms, k.candidate)
- m[to, from] = ifelse(length(k.effective) > 0,
- k.effective, "0")
- }
- }
- } # }}}
- } else { # {{{ Use formation fractions where possible
- for (from in boxes) {
- for (to in boxes) {
- if (from == to) { # diagonal elements
- k.candidate = paste("k", from, sep="_")
- m[from,to] = ifelse(k.candidate %in% model$parms,
- paste("-", k.candidate), "0")
- if(grepl("_free", from)) { # add transfer to bound compartment for SFORB
- m[from,to] = paste(m[from,to], "-", paste("k", from, "bound", sep="_"))
- }
- if(grepl("_bound", from)) { # add backtransfer to free compartment for SFORB
- m[from,to] = paste("- k", from, "free", sep="_")
- }
- m[from,to] = m[from,to]
- } else { # off-diagonal elements
- f.candidate = paste("f", from, "to", to, sep="_")
- k.candidate = paste("k", from, to, sep="_")
- # SFORB with maximum use of formation fractions not implemented, see above
- m[to, from] = ifelse(f.candidate %in% model$parms,
- paste(f.candidate, " * k_", from, sep=""),
- ifelse(k.candidate %in% model$parms, k.candidate, "0"))
- }
- }
- }
- } # }}}
- model$coefmat <- m
- }#}}}
-
- class(model) <- "mkinmod"
- return(model)
-}
-# vim: set foldmethod=marker ts=2 sw=2 expandtab:
diff --git a/R/mkinparplot.R b/R/mkinparplot.R
deleted file mode 100644
index 28c1d2a..0000000
--- a/R/mkinparplot.R
+++ /dev/null
@@ -1,62 +0,0 @@
-# Copyright (C) 2014 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/>
-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)
- if ("g" %in% deparms.optim) fractions.optim <- c("g", fractions.optim)
- rates.optim.unsorted = setdiff(deparms.optim, fractions.optim)
- rates.optim <- rownames(object$start[rates.optim.unsorted, ])
- n.plot <- c(state.optim = length(state.optim),
- rates.optim = length(rates.optim),
- fractions.optim = length(fractions.optim))
- n.plot <- n.plot[n.plot > 0]
-
- oldpars <- par(no.readonly = TRUE)
- layout(matrix(1:length(n.plot), ncol = 1), heights = n.plot + 1)
-
- s <- summary(object)
- bpar <- data.frame(t(s$bpar))
- par(mar = c(2.1, 2.1, 0.1, 2.1))
- par(cex = 1)
- for (type in names(n.plot)) {
- parnames <- get(type)
- values <- bpar[parnames]
- xlim = switch(type,
- state.optim = range(c(0, unlist(values)),
- na.rm = TRUE, finite = TRUE),
- rates.optim = range(c(0, 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],
- xlim = xlim,
- 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)
- position <- ifelse(values["Estimate", ] < mean(xlim)/2, "right", "left")
- text(ifelse(position == "left", min(xlim), max(xlim)),
- length(parnames):1, parnames,
- pos = ifelse(position == "left", 4, 2))
- arrows(as.numeric(values["Lower", ]), length(parnames):1,
- as.numeric(values["Upper", ]), length(parnames):1,
- code = 3, angle = 90, length = 0.05)
- }
- par(oldpars)
-}
diff --git a/R/mkinplot.R b/R/mkinplot.R
deleted file mode 100644
index b9becfd..0000000
--- a/R/mkinplot.R
+++ /dev/null
@@ -1,4 +0,0 @@
-mkinplot <- function(fit, ...)
-{
- plot(fit, ...)
-}
diff --git a/R/mkinpredict.R b/R/mkinpredict.R
deleted file mode 100644
index 14efc56..0000000
--- a/R/mkinpredict.R
+++ /dev/null
@@ -1,121 +0,0 @@
-# Copyright (C) 2010-2014 Johannes Ranke
-# Some lines in this code are copyright (C) 2013 Eurofins Regulatory AG
-# 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/>
-
-mkinpredict <- function(mkinmod, odeparms, odeini,
- outtimes, solution_type = "deSolve",
- method.ode = "lsoda", atol = 1e-8, rtol = 1e-10,
- map_output = TRUE, ...) {
-
- # Get the names of the state variables in the model
- mod_vars <- names(mkinmod$diffs)
-
- # Order the inital values for state variables if they are named
- if (!is.null(names(odeini))) {
- odeini <- odeini[mod_vars]
- }
-
- # Create function for evaluation of expressions with ode parameters and initial values
- evalparse <- function(string)
- {
- eval(parse(text=string), as.list(c(odeparms, odeini)))
- }
-
- # Create a function calculating the differentials specified by the model
- # if necessary
- if (solution_type == "analytical") {
- parent.type = names(mkinmod$map[[1]])[1]
- parent.name = names(mkinmod$diffs)[[1]]
- o <- switch(parent.type,
- SFO = SFO.solution(outtimes,
- evalparse(parent.name),
- ifelse(mkinmod$use_of_ff == "min",
- evalparse(paste("k", parent.name, "sink", sep="_")),
- evalparse(paste("k", parent.name, sep="_")))),
- FOMC = FOMC.solution(outtimes,
- evalparse(parent.name),
- evalparse("alpha"), evalparse("beta")),
- DFOP = DFOP.solution(outtimes,
- evalparse(parent.name),
- evalparse("k1"), evalparse("k2"),
- evalparse("g")),
- HS = HS.solution(outtimes,
- evalparse(parent.name),
- evalparse("k1"), evalparse("k2"),
- evalparse("tb")),
- SFORB = SFORB.solution(outtimes,
- evalparse(parent.name),
- evalparse(paste("k", parent.name, "bound", sep="_")),
- evalparse(paste("k", sub("free", "bound", parent.name), "free", sep="_")),
- evalparse(paste("k", parent.name, "sink", sep="_")))
- )
- out <- data.frame(outtimes, o)
- names(out) <- c("time", sub("_free", "", parent.name))
- }
- if (solution_type == "eigen") {
- coefmat.num <- matrix(sapply(as.vector(mkinmod$coefmat), evalparse),
- nrow = length(mod_vars))
- e <- eigen(coefmat.num)
- c <- solve(e$vectors, odeini)
- f.out <- function(t) {
- e$vectors %*% diag(exp(e$values * t), nrow=length(mod_vars)) %*% c
- }
- o <- matrix(mapply(f.out, outtimes),
- nrow = length(mod_vars), ncol = length(outtimes))
- out <- data.frame(outtimes, t(o))
- names(out) <- c("time", mod_vars)
- }
- if (solution_type == "deSolve") {
- mkindiff <- function(t, state, parms) {
-
- time <- t
- diffs <- vector()
- for (box in names(mkinmod$diffs))
- {
- diffname <- paste("d", box, sep="_")
- diffs[diffname] <- with(as.list(c(time, state, parms)),
- eval(parse(text=mkinmod$diffs[[box]])))
- }
- return(list(c(diffs)))
- }
- out <- ode(
- y = odeini,
- times = outtimes,
- func = mkindiff,
- parms = odeparms,
- method = method.ode,
- atol = atol,
- rtol = rtol,
- ...
- )
- }
- if (map_output) {
- # Output transformation for models with unobserved compartments like SFORB
- out_mapped <- data.frame(time = out[,"time"])
- for (var in names(mkinmod$map)) {
- if((length(mkinmod$map[[var]]) == 1) || solution_type == "analytical") {
- out_mapped[var] <- out[, var]
- } else {
- out_mapped[var] <- rowSums(out[, mkinmod$map[[var]]])
- }
- }
- return(out_mapped)
- } else {
- return(out)
- }
-}
diff --git a/R/mkinresplot.R b/R/mkinresplot.R
deleted file mode 100644
index 07bd7df..0000000
--- a/R/mkinresplot.R
+++ /dev/null
@@ -1,53 +0,0 @@
-# Copyright (C) 2008-2014 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/>
-if(getRversion() >= '2.15.1') utils::globalVariables(c("variable", "residual"))
-
-mkinresplot <- function (object,
- obs_vars = names(object$mkinmod$map),
- xlab = "Time", ylab = "Residual",
- maxabs = "auto", legend= TRUE, lpos = "topright", ...)
-{
- obs_vars_all <- as.character(unique(object$data$variable))
-
- if (length(obs_vars) > 0){
- obs_vars <- intersect(obs_vars_all, obs_vars)
- } else obs_vars <- obs_vars_all
-
- residuals <- subset(object$data, variable %in% obs_vars, residual)
-
- if (maxabs == "auto") maxabs = max(abs(residuals), na.rm = TRUE)
-
- col_obs <- pch_obs <- 1:length(obs_vars)
- names(col_obs) <- names(pch_obs) <- obs_vars
-
- plot(0, xlab = xlab, ylab = ylab,
- xlim = c(0, 1.1 * max(object$data$time)),
- ylim = c(-1.2 * maxabs, 1.2 * maxabs), ...)
-
- for(obs_var in obs_vars){
- residuals_plot <- subset(object$data, variable == obs_var, c("time", "residual"))
- points(residuals_plot, pch = pch_obs[obs_var], col = col_obs[obs_var])
- }
-
- abline(h = 0, lty = 2)
-
- if (legend == TRUE) {
- legend(lpos, inset = c(0.05, 0.05), legend = obs_vars,
- col = col_obs[obs_vars], pch = pch_obs[obs_vars])
- }
-}
diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R
deleted file mode 100644
index 4132d96..0000000
--- a/R/plot.mkinfit.R
+++ /dev/null
@@ -1,98 +0,0 @@
-# Copyright (C) 2010-2014 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/>
-if(getRversion() >= '2.15.1') utils::globalVariables(c("type", "variable", "observed"))
-
-plot.mkinfit <- function(x, fit = x,
- obs_vars = names(fit$mkinmod$map),
- xlab = "Time", ylab = "Observed",
- xlim = range(fit$data$time),
- ylim = "default",
- col_obs = 1:length(fit$mkinmod$map),
- pch_obs = col_obs,
- lty_obs = rep(1, length(fit$mkinmod$map)),
- add = FALSE, legend = !add,
- show_residuals = FALSE, maxabs = "auto",
- lpos = "topright", inset = c(0.05, 0.05), ...)
-{
- if (add && show_residuals) stop("If adding to an existing plot we can not show residuals")
-
- if (ylim == "default") {
- ylim = c(0, max(subset(fit$data, variable %in% obs_vars)$observed, na.rm = TRUE))
- }
-
- solution_type = fit$solution_type
- parms.all <- c(fit$bparms.optim, fit$bparms.fixed)
-
- ininames <- c(
- rownames(subset(fit$start, type == "state")),
- rownames(subset(fit$fixed, type == "state")))
- odeini <- parms.all[ininames]
-
- # Order initial state variables
- names(odeini) <- sub("_0", "", names(odeini))
- odeini <- odeini[names(fit$mkinmod$diffs)]
-
- outtimes <- seq(xlim[1], xlim[2], length.out=100)
-
- odenames <- c(
- rownames(subset(fit$start, type == "deparm")),
- rownames(subset(fit$fixed, type == "deparm")))
- odeparms <- parms.all[odenames]
-
- out <- mkinpredict(fit$mkinmod, odeparms, odeini, outtimes,
- solution_type = solution_type, atol = fit$atol, rtol = fit$rtol)
-
- # Set up the plot if not to be added to an existing plot
- if (add == FALSE) {
- if (show_residuals) {
- oldpar <- par(no.readonly = TRUE)
- layout(matrix(c(1, 2), 2, 1), heights = c(2, 1.3))
- par(mar = c(3, 4, 4, 2) + 0.1)
- }
- plot(0, type="n",
- xlim = xlim, ylim = ylim,
- xlab = xlab, ylab = ylab, ...)
- }
- # Plot the data and model output
- names(col_obs) <- names(pch_obs) <- names(lty_obs) <- names(fit$mkinmod$map)
- for (obs_var in obs_vars) {
- points(subset(fit$data, variable == obs_var, c(time, observed)),
- pch = pch_obs[obs_var], col = col_obs[obs_var])
- }
- matlines(out$time, out[obs_vars], col = col_obs[obs_vars], lty = lty_obs[obs_vars])
- if (legend == TRUE) {
- legend(lpos, inset= inset, legend = obs_vars,
- col = col_obs[obs_vars], pch = pch_obs[obs_vars], lty = lty_obs[obs_vars])
- }
- # Show residuals if requested
- if (show_residuals) {
- par(mar = c(5, 4, 0, 2) + 0.1)
- residuals <- subset(fit$data, variable %in% obs_vars, residual)
- if (maxabs == "auto") maxabs = max(abs(residuals), na.rm = TRUE)
- plot(0, type="n",
- xlim = xlim,
- ylim = c(-1.2 * maxabs, 1.2 * maxabs),
- xlab = xlab, ylab = "Residuals")
- for(obs_var in obs_vars){
- residuals_plot <- subset(fit$data, variable == obs_var, c("time", "residual"))
- points(residuals_plot, pch = pch_obs[obs_var], col = col_obs[obs_var])
- }
- abline(h = 0, lty = 2)
- par(oldpar, no.readonly = TRUE)
- }
-}
diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R
deleted file mode 100644
index a0c302f..0000000
--- a/R/transform_odeparms.R
+++ /dev/null
@@ -1,101 +0,0 @@
-# Copyright (C) 2010-2014 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/>
-
-transform_odeparms <- function(parms, mod_vars,
- transform_rates = TRUE,
- transform_fractions = TRUE)
-{
- # Set up container for transformed parameters
- transparms <- parms
-
- # Log transformation for rate constants if requested
- index_k <- grep("^k_", names(transparms))
- if (length(index_k) > 0) {
- if(transform_rates) transparms[index_k] <- log(parms[index_k])
- else transparms[index_k] <- parms[index_k]
- }
-
- # Go through state variables and apply isotropic logratio transformation if requested
- for (box in mod_vars) {
- indices_f <- grep(paste("^f", box, sep = "_"), names(parms))
- f_names <- grep(paste("^f", box, sep = "_"), names(parms), value = TRUE)
- n_paths <- length(indices_f)
- if (n_paths > 0) {
- f <- parms[indices_f]
- trans_f <- ilr(c(f, 1 - sum(f)))
- names(trans_f) <- f_names
- if(transform_fractions) transparms[indices_f] <- trans_f
- else transparms[indices_f] <- f
- }
- }
-
- # Transform rates, fractions and tb also for FOMC, DFOP and HS models if requested
- for (pname in c("alpha", "beta", "k1", "k2", "tb")) {
- if (!is.na(parms[pname])) {
- transparms[pname] <- ifelse(transform_rates, log(parms[pname]), parms[pname])
- transparms[pname] <- ifelse(transform_rates, log(parms[pname]), parms[pname])
- }
- }
- if (!is.na(parms["g"])) {
- g <- parms["g"]
- transparms["g"] <- ifelse(transform_fractions, ilr(c(g, 1 - g)), g)
- }
-
- return(transparms)
-}
-
-backtransform_odeparms <- function(transparms, mod_vars,
- transform_rates = TRUE,
- transform_fractions = TRUE)
-{
- # Set up container for backtransformed parameters
- parms <- transparms
-
- # Exponential transformation for rate constants
- index_k <- grep("^k_", names(parms))
- if (length(index_k) > 0) {
- if(transform_rates) parms[index_k] <- exp(transparms[index_k])
- else parms[index_k] <- transparms[index_k]
- }
-
- # Go through state variables and apply inverse isotropic logratio transformation
- for (box in mod_vars) {
- indices_f <- grep(paste("^f", box, sep = "_"), names(transparms))
- f_names <- grep(paste("^f", box, sep = "_"), names(transparms), value = TRUE)
- n_paths <- length(indices_f)
- if (n_paths > 0) {
- f <- invilr(transparms[indices_f])[1:n_paths] # We do not need the last component
- names(f) <- f_names
- if(transform_fractions) parms[indices_f] <- f
- else parms[indices_f] <- transparms[indices_f]
- }
- }
-
- # Transform parameters also for FOMC, DFOP and HS models
- for (pname in c("alpha", "beta", "k1", "k2", "tb")) {
- if (!is.na(transparms[pname])) {
- parms[pname] <- ifelse(transform_rates, exp(transparms[pname]), transparms[pname])
- }
- }
- if (!is.na(transparms["g"])) {
- g <- transparms["g"]
- parms["g"] <- ifelse(transform_fractions, invilr(g)[1], g)
- }
-
- return(parms)
-}

Contact - Imprint