From a6694c655fde246dd4d59b44fd10b22738b3fb08 Mon Sep 17 00:00:00 2001 From: jranke Date: Mon, 7 May 2012 18:51:46 +0000 Subject: - Moved the call to mkinerrmin to summary.mkinfit - The argument to mkinerrmin is now an object of class mkinfit - Fixed the allocation of parameters to observed variables in mkinerrmin git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@37 edb9625f-4e0d-4859-8d74-9fd3b1da38cb --- R/mkinerrmin.R | 72 ++++++++++++++++++++++++++++++++++++++++++++---------- R/mkinfit.R | 45 +++++++++------------------------- man/mkinerrmin.Rd | 25 +++++++++---------- man/mkinpredict.Rd | 2 +- 4 files changed, 84 insertions(+), 60 deletions(-) diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R index 648ba45c..5415b41d 100644 --- a/R/mkinerrmin.R +++ b/R/mkinerrmin.R @@ -1,7 +1,7 @@ # $Id$ -# Copyright (C) 2010 Johannes Ranke -# Contact: mkin-devel@lists.berlios.de +# Copyright (C) 2010-2012 Johannes Ranke +# Contact: jranke@uni-bremen.de # This file is part of the R package mkin @@ -18,17 +18,63 @@ # You should have received a copy of the GNU General Public License along with # this program. If not, see -mkinerrmin <- function(errdata, n.parms, alpha = 0.05) +mkinerrmin <- function(fit, alpha = 0.05) { - means.mean <- mean(errdata$value_mean, na.rm=TRUE) - - df = length(errdata$value_mean) - n.parms + parms.optim <- fit$parms.all + kinerrmin <- function(errdata, n.parms) { + means.mean <- mean(errdata$value_mean, na.rm=TRUE) + df = length(errdata$value_mean) - n.parms - f <- function(err) - { - (sum((errdata$value_mean - errdata$value_pred)^2/((err * means.mean)^2)) - - qchisq(1 - alpha,df))^2 - } - err.min <- optimize(f, c(0.01,0.9))$minimum - return(list(err.min = err.min, n.optim = n.parms, df = df)) + f <- function(err) + { + (sum((errdata$value_mean - errdata$value_pred)^2/((err * means.mean)^2)) - + qchisq(1 - alpha,df))^2 + } + err.min <- optimize(f, c(0.01,0.9))$minimum + 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), ] + 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 index cb0396f7..b2641e64 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -175,9 +175,9 @@ mkinfit <- function(mkinmod, observed, fit$mkinmod <- mkinmod # We need data and predictions for summary and plotting - fit$observed <- mkin_long_to_wide(observed) - predicted_long <- mkin_wide_to_long(out_predicted, time = "time") - fit$predicted <- out_predicted + 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(initial = c(state.ini.optim, @@ -189,35 +189,11 @@ mkinfit <- function(mkinmod, observed, value = c(state.ini.fixed, parms.fixed)) fit$fixed$type = c(rep("state", length(state.ini.fixed)), rep("deparm", length(parms.fixed))) - # Calculate chi2 error levels according to FOCUS (2006) - means <- aggregate(value ~ time + name, data = observed, mean, na.rm=TRUE) - errdata <- merge(means, predicted_long, by = c("time", "name"), suffixes = c("_mean", "_pred")) - errdata <- errdata[order(errdata$time, errdata$name), ] - errmin.overall <- mkinerrmin(errdata, length(parms.optim) + length(state.ini.optim)) - - 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 obs_vars) - { - errdata.var <- subset(errdata, name == obs_var) - n.k.optim <- length(grep(paste("k", obs_var, sep="_"), names(parms.optim))) - n.initials.optim <- length(grep(paste(obs_var, ".*", "_0", sep=""), names(state.ini.optim))) - n.optim <- n.k.optim + n.initials.optim - 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 - errmin.tmp <- mkinerrmin(errdata.var, n.optim) - errmin[obs_var, c("err.min", "n.optim", "df")] <- errmin.tmp - } - fit$errmin <- errmin - # Calculate dissipation times DT50 and DT90 and, if necessary, formation fractions - # from optimised parameters parms.all = backtransform_odeparms(c(fit$par, parms.fixed), mod_vars) + + # Calculate dissipation times DT50 and DT90 and, if necessary, formation + # fractions and SFORB eigenvalues from optimised parameters fit$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), DT90 = rep(NA, length(obs_vars)), row.names = obs_vars) fit$ff <- vector() @@ -305,7 +281,7 @@ mkinfit <- function(mkinmod, observed, } # Collect observed, predicted and residuals - data <- merge(observed, predicted_long, by = c("time", "name")) + data <- merge(fit$observed, fit$predicted, by = c("time", "name")) names(data) <- c("time", "variable", "observed", "predicted") data$residual <- data$observed - data$predicted data$variable <- ordered(data$variable, levels = obs_vars) @@ -362,9 +338,12 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ...) { ans$start <- object$start ans$fixed <- object$fixed - ans$errmin <- object$errmin + + ans$errmin <- mkinerrmin(object, alpha = 0.05) + ans$parms.all <- object$parms.all - ans$ff <- object$ff + if (!is.null(object$ff)) + ans$ff <- object$ff if(distimes) ans$distimes <- object$distimes if(length(object$SFORB) != 0) ans$SFORB <- object$SFORB class(ans) <- c("summary.mkinfit", "summary.modFit") diff --git a/man/mkinerrmin.Rd b/man/mkinerrmin.Rd index 654115b6..c43d87a1 100644 --- a/man/mkinerrmin.Rd +++ b/man/mkinerrmin.Rd @@ -10,30 +10,29 @@ smallest relative error still resulting in passing the chi-squared test as defined in the FOCUS kinetics report from 2006. } \usage{ -mkinerrmin(errdata, n.parms, alpha = 0.05) +mkinerrmin(fit, alpha = 0.05) } \arguments{ - \item{errdata}{ - A data frame with mean observed values in column named \code{value_mean} - and predicted values in column \code{value_pred}. -} - \item{n.parms}{ - The number of optimized parameters to be taken into account for the data series. -} + \item{fit}{ + an object of class \code{\link{mkinfit}}. + } \item{alpha}{ The confidence level chosen for the chi-squared test. } } \value{ - A list with the following components: + A dataframe with the following components: \item{err.min}{The relative error, expressed as a fraction.} \item{n.optim}{The number of optimised parameters attributed to the data series.} - \item{df}{The number of remaining degrees of freedom for the chi2 error level calculations. - Note that mean values are used for the chi2 statistic and therefore every time point with - observed values in the series only counts one time.} + \item{df}{The number of remaining degrees of freedom for the chi2 error level + calculations. Note that mean values are used for the chi2 statistic and + therefore every time point with observed values in the series only counts + one time.} + The dataframe has one row for the total dataset and one further row for + each observed state variable in the model. } \details{ - This function is used internally by \code{\link{mkinfit}}. + This function is used internally by \code{\link{summary.mkinfit}}. } \references{ FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and diff --git a/man/mkinpredict.Rd b/man/mkinpredict.Rd index aaa77e6e..d2931443 100644 --- a/man/mkinpredict.Rd +++ b/man/mkinpredict.Rd @@ -60,7 +60,7 @@ mkinpredict(mkinmod, odeparms, odeini, outtimes, solution_type = "deSolve", map_ mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 1:20, solution_type = "analytical")[20,] mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, atol = 1e-20)[20,] # The integration method does not make a lot of difference - mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, atol = 1e-20, method = "bdf")[20,] + mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, atol = 1e-20, method = "ode45")[20,] mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, atol = 1e-20, method = "rk4")[20,] # The number of output times does make a lot of difference mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), seq(0, 20, by = 0.1))[201,] -- cgit v1.2.1