aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
Diffstat (limited to 'R')
-rw-r--r--R/mkinerrmin.R72
-rw-r--r--R/mkinfit.R45
2 files changed, 71 insertions, 46 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 <http://www.gnu.org/licenses/>
-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")

Contact - Imprint