From e959fde98f95f3595e01490b67892678bbcd1b27 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 7 May 2014 14:47:28 +0200 Subject: Fork the gmkin GUI from mkin. See ChangeLog for details --- R/DFOP.solution.R | 5 - R/FOMC.solution.R | 4 - R/HS.solution.R | 6 - R/SFO.solution.R | 4 - R/SFORB.solution.R | 9 - R/endpoints.R | 108 ----------- R/gmkin.R | 10 +- R/ilr.R | 51 ------ R/mkin_long_to_wide.R | 30 --- R/mkin_wide_to_long.R | 34 ---- R/mkinerrmin.R | 85 --------- R/mkinfit.R | 482 ------------------------------------------------- R/mkinmod.R | 242 ------------------------- R/mkinparplot.R | 62 ------- R/mkinplot.R | 4 - R/mkinpredict.R | 121 ------------- R/mkinresplot.R | 53 ------ R/plot.mkinfit.R | 98 ---------- R/transform_odeparms.R | 101 ----------- 19 files changed, 1 insertion(+), 1508 deletions(-) delete mode 100644 R/DFOP.solution.R delete mode 100644 R/FOMC.solution.R delete mode 100644 R/HS.solution.R delete mode 100644 R/SFO.solution.R delete mode 100644 R/SFORB.solution.R delete mode 100644 R/endpoints.R delete mode 100644 R/ilr.R delete mode 100644 R/mkin_long_to_wide.R delete mode 100644 R/mkin_wide_to_long.R delete mode 100644 R/mkinerrmin.R delete mode 100644 R/mkinfit.R delete mode 100644 R/mkinmod.R delete mode 100644 R/mkinparplot.R delete mode 100644 R/mkinplot.R delete mode 100644 R/mkinpredict.R delete mode 100644 R/mkinresplot.R delete mode 100644 R/plot.mkinfit.R delete mode 100644 R/transform_odeparms.R (limited to 'R') 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 - -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 - -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 -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 -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 -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 }}} - -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 -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 - -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 -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 -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 - -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) -} -- cgit v1.2.1