aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorjranke <jranke@edb9625f-4e0d-4859-8d74-9fd3b1da38cb>2012-03-27 01:03:18 +0000
committerjranke <jranke@edb9625f-4e0d-4859-8d74-9fd3b1da38cb>2012-03-27 01:03:18 +0000
commitfff1fc581da5b4ff935ebd4d7ded02f750472fdc (patch)
treea18da58a5bfd013c1bd35bbc7828925084a13a76
parent1718d434efae26de02754c6622c43f4dc9e624b9 (diff)
Start of the transition to fitting transformed parameters.
Many things are broken now (see TODO list) git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@20 edb9625f-4e0d-4859-8d74-9fd3b1da38cb
-rw-r--r--DESCRIPTION13
-rw-r--r--R/ilr.R44
-rw-r--r--R/mkinfit.R96
-rw-r--r--R/mkinmod.R168
-rw-r--r--R/mkinstart.R23
-rw-r--r--TODO10
-rw-r--r--man/ilr.Rd62
-rw-r--r--man/mkinfit.Rd13
-rw-r--r--man/mkinstart.Rd38
-rw-r--r--man/summary.mkinfit.Rd5
10 files changed, 239 insertions, 233 deletions
diff --git a/DESCRIPTION b/DESCRIPTION
index a2c760c..6951ede 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,19 +2,20 @@ Package: mkin
Type: Package
Title: Routines for fitting kinetic models with one or more state
variables to chemical degradation data
-Version: 0.8-11
-Date: 2011-05-19
-Author: Johannes Ranke, Katrin Lindenberger
-Maintainer: Johannes Ranke <jranke@users.berlios.de>
+Version: 0.9-01
+Date: 2012-03-26
+Author: Johannes Ranke, Katrin Lindenberger, René Lehmann
+Maintainer: Johannes Ranke <jranke@uni-bremen.de>
Description: Calculation routines based on the FOCUS Kinetics Report (2006).
Includes a function for conveniently defining differential equation models,
model solution based on eigenvalues if possible or using numerical solvers
and a choice of the optimisation methods made available by the FME package
(default is a Levenberg-Marquardt variant). Please note that no warranty is
implied for correctness of results or fitness for a particular purpose.
-Depends: FME, deSolve, kinfit
+Depends: FME, deSolve
Suggests: RUnit
License: GPL
LazyLoad: yes
LazyData: yes
-URL: http://cran.r-project.org, http://developer.berlios.de/projects/mkin/
+Encoding: UTF-8
+URL: http://cran.r-project.org, http://r-forge.r-project.org/projects/kinfit
diff --git a/R/ilr.R b/R/ilr.R
new file mode 100644
index 0000000..f2a6f1b
--- /dev/null
+++ b/R/ilr.R
@@ -0,0 +1,44 @@
+# $Id: ilr.R 120 2011-09-02 14:25:35Z jranke $
+
+# 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))
+ }
+ return(z)
+}
diff --git a/R/mkinfit.R b/R/mkinfit.R
index b742742..83896bc 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -1,7 +1,7 @@
# $Id: mkinfit.R 120 2011-09-02 14:25:35Z jranke $
-# Copyright (C) 2010-2011 Johannes Ranke
-# Contact: mkin-devel@lists.berlios.de
+# Copyright (C) 2010-2012 Johannes Ranke
+# 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
@@ -22,9 +22,8 @@
# this program. If not, see <http://www.gnu.org/licenses/>
mkinfit <- function(mkinmod, observed,
- parms.ini = rep(0.1, length(mkinmod$parms)),
+ parms.ini = "auto",
state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)),
- lower = 0, upper = Inf,
fixed_parms = NULL,
fixed_initials = names(mkinmod$diffs)[-1],
eigen = FALSE,
@@ -34,14 +33,24 @@ mkinfit <- function(mkinmod, observed,
...)
{
mod_vars <- names(mkinmod$diffs)
- # Subset dataframe with mapped (modelled) variables
+ # Subset dataframe of observed variables with mapped (modelled) variables
observed <- subset(observed, name %in% names(mkinmod$map))
# Get names of observed variables
obs_vars = unique(as.character(observed$name))
- # Name the parameters if they are not named yet
- if(is.null(names(parms.ini))) names(parms.ini) <- mkinmod$parms
-
+ # Define initial values for parameters where not specified by the user
+ if (parms.ini[[1]] == "auto") parms.ini = vector()
+ defaultpar.names <- setdiff(mkinmod$parms, names(parms.ini))
+ for (parmname in defaultpar.names) {
+ # Default values for the FOMC model
+ if (parmname == "alpha") parms.ini[parmname] = 1
+ if (parmname == "beta") parms.ini[parmname] = 10
+ # Default values for log of rate constants, depending on the parameterisation
+ if (substr(parmname, 1, 2) == "k_") parms.ini[parmname] = -2.3
+ # Default values for ilr of formation fractions
+ if (substr(parmname, 1, 2) == "f_") parms.ini[parmname] = 0.1
+ }
+
# Name the inital parameter values if they are not named yet
if(is.null(names(state.ini))) names(state.ini) <- mod_vars
@@ -58,14 +67,6 @@ mkinfit <- function(mkinmod, observed,
names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_")
}
- # Set upper limit for formation fractions to one if formation fractions are
- # directly defined and if no user input for upper limit is given
- if (all(upper==Inf) & any(grepl("f_", names(parms.ini)))==TRUE){
- upper=c( rep(Inf,length(parms.optim)))
- upper[grep("f_", names(parms.optim))]=1
- upper=c(rep(Inf, length(state.ini.optim)), upper)
- }
-
# 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
@@ -80,12 +81,25 @@ mkinfit <- function(mkinmod, observed,
# if necessary
if(solution == "deSolve") {
mkindiff <- function(t, state, parms) {
+ transparms <- vector()
+ transparms <- c(transparms, exp(parms[grep("^k_", names(parms))]))
+ for (box in mod_vars) {
+ path_indices <- grep(paste("^f", box, sep = "_"), names(parms))
+ n_paths <- length(path_indices)
+ if (n_paths > 0) {
+ f <- invilr(parms[grep(paste("^f", box, sep="_"), names(parms))])
+ transparms <- c(transparms, f[1:n_paths])
+ }
+ }
+ otherparms <- parms[setdiff(names(parms), names(transparms))]
+ parms <- c(transparms, otherparms)
+
time <- t
diffs <- vector()
for (box in mod_vars)
{
diffname <- paste("d", box, sep="_")
- diffs[diffname] <- with(as.list(c(time,state, parms)),
+ diffs[diffname] <- with(as.list(c(time, state, parms)),
eval(parse(text=mkinmod$diffs[[box]])))
}
return(list(c(diffs)))
@@ -251,7 +265,7 @@ mkinfit <- function(mkinmod, observed,
}
return(mC)
}
- fit <- modFit(cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, ...)
+ fit <- modFit(cost, c(state.ini.optim, parms.optim), ...)
# We need to return some more data for summary and plotting
fit$solution <- solution
@@ -323,7 +337,7 @@ mkinfit <- function(mkinmod, observed,
DT90 = log(10)/k_tot
for (k_name in k_names)
{
- fit$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot
+ fit$ff[[sub("^k_", "", k_name)]] = parms.all[[k_name]] / k_tot
}
}
if (type == "FOMC") {
@@ -395,7 +409,7 @@ mkinfit <- function(mkinmod, observed,
if (abs(DT90.o - max_DT) < 0.01) DT90 = NA else DT90 = DT90.o
for (k_out_name in k_out_names)
{
- fit$ff[[sub("k_", "", k_out_name)]] = parms.all[[k_out_name]] / k_1output
+ fit$ff[[sub("^k_", "", k_out_name)]] = parms.all[[k_out_name]] / k_1output
}
# Return the eigenvalues for comparison with DFOP rate constants
fit$SFORB[[paste(obs_var, "b1", sep="_")]] = b1
@@ -416,8 +430,11 @@ mkinfit <- function(mkinmod, observed,
return(fit)
}
-summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, cov = FALSE,...) {
+summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, ...) {
param <- object$par
+ if (object$logk) {
+ names(param) <- sub("^k_", "log k_", names(param))
+ }
pnames <- names(param)
p <- length(param)
covar <- try(solve(0.5*object$hessian), silent = TRUE) # unscaled covariance
@@ -441,32 +458,25 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, cov
warning(message)
} else message <- "ok"
- param <- cbind(param, se, tval, pt(tval, rdf, lower.tail = FALSE))
- dimnames(param) <- list(pnames, c("Estimate", "Std. Error",
- "t value", "Pr(>t)"))
- if(cov)
- ans <- list(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)
- else
- ans <- list(residuals = object$residuals,
- residualVariance = resvar,
- sigma = sqrt(resvar),
- modVariance = modVariance,
- df = c(p, rdf),
- info = object$info, niter = object$iterations,
- stopmess = message,
- par = param)
+ param <- cbind(param, se)
+ dimnames(param) <- list(pnames, c("Estimate", "Std. Error"))
+
+ ans <- list(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)
ans$diffs <- object$diffs
if(data) ans$data <- object$data
ans$start <- object$start
+ if (object$logk) {
+ names(ans$start) <- sub("^k_", "log k_", names(ans$start))
+ }
ans$fixed <- object$fixed
ans$errmin <- object$errmin
if(distimes) ans$distimes <- object$distimes
diff --git a/R/mkinmod.R b/R/mkinmod.R
index 31c778c..b49c813 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -1,7 +1,7 @@
# $Id: mkinmod.R 71 2010-09-12 01:13:36Z jranke $
-# 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
@@ -33,21 +33,27 @@ mkinmod <- function(...)
if(spec[[1]]$type %in% c("FOMC", "DFOP", "HS")) {
mat = FALSE
if(!is.null(spec[[1]]$to)) {
- message <- paste("Only constant formation fractions over time are implemented.",
- "Depending on the reason for the time dependence of degradation, this may be unrealistic",
+ message <- paste(
+ "Only constant formation fractions over time are implemented.",
+ "Depending on the reason for the time dependence of degradation",
+ "this may be unrealistic",
sep="\n")
warning(message)
} else message <- "ok"
} else mat = TRUE
# Establish list of differential equations
+ # as well as map from observed compartments to differential equations
for (varname in obs_vars)
{
if(is.null(spec[[varname]]$type)) stop(
"Every argument to mkinmod 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")
- new_parms <- vector()
+ 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,
@@ -62,124 +68,75 @@ mkinmod <- function(...)
# 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)
+ {
+ # Add first-order term to first (or only) box for SFO and SFORB
+ box_1 = map[[varname]][[1]] # This is the only box unless type is SFORB
+ if(spec[[varname]]$type %in% c("SFO", "SFORB")) {
+ k_compound <- paste("k", box_1, sep="_")
+ origin_term <- paste(k_compound, "*", box_1)
+ parms <- c(parms, k_compound)
+ }
- # Turn on sink if not specified otherwise
- if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE
# Construct and add FOMC term and add FOMC parameters if needed
if(spec[[varname]]$type == "FOMC") {
- if(match(varname, obs_vars) != 1) {
- stop("Type FOMC is only possible for the first compartment, which is assumed to be the source compartment")
- }
- if(spec[[varname]]$sink == FALSE) {
- stop("Turning off the sink for the FOMC model is not implemented")
- }
# From p. 53 of the FOCUS kinetics report
- nonlinear_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", new_boxes[[1]])
- new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
- new_parms <- c("alpha", "beta")
- ff <- vector()
+ origin_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", box_1)
+ parms <- c(parms, "alpha", "beta")
}
# Construct and add DFOP term and add DFOP parameters if needed
if(spec[[varname]]$type == "DFOP") {
- if(match(varname, obs_vars) != 1) {
- stop("Type DFOP is only possible for the first compartment, which is assumed to be the source compartment")
- }
- if(spec[[varname]]$sink == FALSE) {
- stop("Turning off the sink for the DFOP model is not implemented")
- }
# From p. 57 of the FOCUS kinetics report
- nonlinear_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", new_boxes[[1]])
- new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
- new_parms <- c("k1", "k2", "g")
- ff <- vector()
+ origin_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")
}
# Construct and add HS term and add HS parameters if needed
if(spec[[varname]]$type == "HS") {
- if(match(varname, obs_vars) != 1) {
- stop("Type HS is only possible for the first compartment, which is assumed to be the source compartment")
- }
- if(spec[[varname]]$sink == FALSE) {
- stop("Turning off the sink for the HS model is not implemented")
- }
# From p. 55 of the FOCUS kinetics report
- nonlinear_term <- paste("ifelse(time <= tb, k1, k2)", "*", new_boxes[[1]])
- new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term)
- new_parms <- c("k1", "k2", "tb")
- ff <- vector()
+ origin_term <- paste("ifelse(time <= tb, k1, k2)", "*", box_1)
+ parms <- c(parms, "k1", "k2", "tb")
}
- # Construct terms for transfer to sink and add if appropriate
-
- if(spec[[varname]]$sink) {
- # Add first-order sink term to first (or only) box for SFO and SFORB
- if(spec[[varname]]$type %in% c("SFO", "SFORB")) {
- k_compound_sink <- paste("k", new_boxes[[1]], "sink", sep="_")
- sink_term <- paste("-", k_compound_sink, "*", new_boxes[[1]])
- new_diffs[[1]] <- paste(new_diffs[[1]], sink_term)
- new_parms <- k_compound_sink
- }
- }
+ # Add origin decline term to box 1 (usually the only box, unless type is SFORB)
+ diffs[[box_1]] <- paste(diffs[[box_1]], "-", origin_term)
- # Add reversible binding if appropriate
+ # Add reversible binding for SFORB models
if(spec[[varname]]$type == "SFORB") {
- k_free_bound <- paste("k", varname, "free", "bound", sep="_")
- k_bound_free <- paste("k", varname, "bound", "free", sep="_")
- reversible_binding_terms <- c(
- paste("-", k_free_bound, "*", new_boxes[[1]], "+", k_bound_free, "*", new_boxes[[2]]),
- paste("+", k_free_bound, "*", new_boxes[[1]], "-", k_bound_free, "*", new_boxes[[2]]))
- new_diffs <- paste(new_diffs, reversible_binding_terms)
- new_parms <- c(new_parms, k_free_bound, k_bound_free)
+ box_2 = map[[varname]][[2]]
+ k_free_bound <- paste("k", varname, "free", "bound", sep="_")
+ k_bound_free <- paste("k", varname, "bound", "free", sep="_")
+ reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+",
+ k_bound_free, "*", box_2)
+ diffs[[box_1]] <- paste(diffs[[box_1]], reversible_binding_term_1)
+ reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-",
+ k_bound_free, "*", box_2)
+ diffs[[box_2]] <- paste(diffs[[box_2]], reversible_binding_term_2)
+ parms <- c(parms, k_free_bound, k_bound_free)
}
- # Add observed variable to model
- parms <- c(parms, new_parms)
- names(new_diffs) <- new_boxes
- diffs <- c(diffs, new_diffs)
- }
- # Transfer between compartments
- for (varname in obs_vars) {
+ # Transfer between compartments
to <- spec[[varname]]$to
if(!is.null(to)) {
- origin_box <- switch(spec[[varname]]$type,
- SFO = varname,
- FOMC = varname,
- DFOP = varname,
- HS = varname,
- SFORB = paste(varname, "free", sep="_"))
- fraction_left <- NULL
+ # 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(spec[[varname]]$type %in% c("SFO", "SFORB")) {
- k_from_to <- paste("k", origin_box, target_box, sep="_")
- diffs[[origin_box]] <- paste(diffs[[origin_box]], "-",
- k_from_to, "*", origin_box)
- diffs[[target_box]] <- paste(diffs[[target_box]], "+",
- k_from_to, "*", origin_box)
- parms <- c(parms, k_from_to)
- }
- if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS")) {
- fraction_to_target = paste("f_to", target, sep="_")
- fraction_not_to_target = paste("(1 - ", fraction_to_target, ")",
- sep="")
- if(is.null(fraction_left)) {
- fraction_really_to_target = fraction_to_target
- fraction_left = fraction_not_to_target
- } else {
- fraction_really_to_target = paste(fraction_left, " * ",
- fraction_to_target, sep="")
- fraction_left = paste(fraction_left, " * ",
- fraction_not_to_target, sep="")
- }
- ff[target_box] = fraction_really_to_target
- diffs[[target_box]] <- paste(diffs[[target_box]], "+",
- ff[target_box], "*", nonlinear_term)
- parms <- c(parms, fraction_to_target)
- }
+ fraction_to_target = paste("f", origin_box, "to", target, sep="_")
+ diffs[[target_box]] <- paste(diffs[[target_box]], "+",
+ fraction_to_target, "*", origin_term)
+ parms <- c(parms, fraction_to_target)
}
}
}
@@ -193,26 +150,23 @@ mkinmod <- function(...)
for (from in boxes) {
for (to in boxes) {
if (from == to) {
- 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")
+ k.candidate = paste("k", from, sep="_")
+ m[from,to] = ifelse(k.candidate %in% model$parms,
+ paste("-", k.candidate), "0")
} else {
+ f.candidate = paste("f", from, "to", to, sep="_")
k.candidate = paste("k", from, to, 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[to, from] = ifelse(length(k.effective) > 0,
- k.effective, "0")
+ 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
}
- if (exists("ff")) model$ff = ff
class(model) <- "mkinmod"
- invisible(model)
+ return(model)
}
diff --git a/R/mkinstart.R b/R/mkinstart.R
deleted file mode 100644
index 46e5bed..0000000
--- a/R/mkinstart.R
+++ /dev/null
@@ -1,23 +0,0 @@
-mkinstart <- function(model, data, mode = "auto")
-{
- if (class(model) != "mkinmod") stop("The first argument must be a model of class mkinmod")
- names <- model$parms
- observed <- names(model$map)
- if(!all(observed %in% levels(data$name))) stop("The data must contain the observed variables used in the model")
- for (obs in observed)
- {
- tmp <- subset(data, name == obs)
- max <- tmp[which.max(tmp$value), ]
- type = names(model$map[[obs]])[[1]]
- kinmodel <- ifelse(type == "SFORB", "DFOP", type)
- tmp.longdata <- subset(data, name == obs & time >= max$time)
- tmp.widedata <- mkin_long_to_wide(tmp.longdata, outtime = "t")
- names(tmp.widedata) <- c("t", "parent")
- tmp.fit <- kinfit(
- kindata = tmp.widedata,
- kinmodels = kinmodel,
- parent.0.user = max$value)
- if(class(tmp.fit[[kinmodel]]) == "try-error") stop(paste("Automatic generation of starting parameters failed\nkinfit failed to find a", kinmodel, "fit for", obs))
- tmp.results <- kinresults(tmp.fit)
- }
-}
diff --git a/TODO b/TODO
index bfefe0c..3da3339 100644
--- a/TODO
+++ b/TODO
@@ -1,4 +1,14 @@
+- Fix analytical and Eivenvalue based solutions after transition to fitting
+ transformed k and f values
+- Transfer calculation of DT50 and DT90 calculations from mkinfit to summary.mkinfit
+- Calculate back-transformed parameters in summary.mkinfit
+- Report back-transformed parameters in print.summary.mkinfit
+- Use back-transformed parameters for DT50 and DT90 calculations
+- Calculate confidence intervals for parameters assuming normal distribution
+- Calculate confidence intervals for DT50 and DT90 values when only one parameter is involved
- Fix DT50 and DT90 calculation for SFORB_SFO
- Add unit tests for mkinfit
+- Remove dependency to kinfit
- Document validation against fits documented in chapter 13 of FOCUS (2006)
- Reproduce example anaylses (L1, L2, ...) in FOCUS (2006)
+- Generate confidence intervals for DT50 and DT90 values for other cases
diff --git a/man/ilr.Rd b/man/ilr.Rd
new file mode 100644
index 0000000..fef0d29
--- /dev/null
+++ b/man/ilr.Rd
@@ -0,0 +1,62 @@
+\name{ilr}
+\alias{ilr}
+\alias{invilr}
+\title{
+ Function to perform isotropic log-ratio transformation
+}
+\description{
+ This implementation is a special case of the class of isotropic log-ratio transformations.
+}
+\usage{
+ ilr(x)
+ invilr(x)
+}
+\arguments{
+ \item{x}{
+ A numeric vector. Naturally, the forward transformation is only sensible for
+ vectors with all elements being greater than zero.
+ }
+}
+\details{
+}
+\value{
+ The result of the forward or backward transformation. The returned components always
+ sum to 1 for the case of the inverse log-ratio transformation.
+}
+\references{
+%% ~put references to the literature/web site here ~
+}
+\author{
+ René Lehmann and Johannes Ranke
+}
+\note{
+%% ~~further notes~~
+}
+
+\seealso{
+ Other implementations are in R packages \code{compositions} and \code{robCompositions}.
+}
+\examples{
+# Order matters
+ilr(c(0.1, 1, 10))
+ilr(c(10, 1, 0.1))
+# Equal entries give ilr transformations with zeros as elements
+ilr(c(3, 3, 3))
+# Almost equal entries give small numbers
+ilr(c(0.3, 0.4, 0.3))
+# Only the ration between the numbers counts, not their sum
+invilr(ilr(c(0.7, 0.29, 0.01)))
+invilr(ilr(2.1 * c(0.7, 0.29, 0.01)))
+# Inverse transformation of larger numbers gives unequal elements
+invilr(-10)
+invilr(c(-10, 0))
+# The sum of the elements of the inverse ilr is 1
+sum(invilr(c(-10, 0)))
+# This is why we do not need all elements of the inverse transformation to go back:
+a <- c(0.1, 0.3, 0.5)
+b <- invilr(a)
+length(b) # Four elements
+ilr(c(b[1:3], 1 - sum(b[1:3]))) # Gives c(0.1, 0.3, 0.5)
+}
+
+\keyword{ manip }
diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd
index d4139eb..24dcb83 100644
--- a/man/mkinfit.Rd
+++ b/man/mkinfit.Rd
@@ -11,9 +11,8 @@
}
\usage{
mkinfit(mkinmod, observed,
- parms.ini = rep(0.1, length(mkinmod$parms)),
+ parms.ini = "auto",
state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)),
- lower = 0, upper = Inf,
fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1],
eigen = FALSE,
plot = FALSE, quiet = FALSE, err = NULL, weight = "none",
@@ -48,16 +47,6 @@ mkinfit(mkinmod, observed,
\code{\link{mkinmod}}). The default is to set the initial value of the first model
variable to 100 and all others to 0.
}
- \item{lower}{
- Lower bounds for the parameters, passed to \code{\link{modFit}}. Defaults to 0 because
- negative values to not make sense for the models currently created by \code{\link{mkinmod}}.
- }
- \item{upper}{
- Upper bounds for the parameters, passed to \code{\link{modFit}}. Defaults to \code{Inf}
- except for formation fraction parameters. Setting non-infinite upper bounds has a strong
- impact on performance, and using a method like "L-BFGS-B" (by specifying an additional
- \code{method} argument) is recommended.
- }
\item{fixed_parms}{
The names of parameters that should not be optimised but rather kept at the values
specified in \code{parms.ini}.
diff --git a/man/mkinstart.Rd b/man/mkinstart.Rd
deleted file mode 100644
index ec1204f..0000000
--- a/man/mkinstart.Rd
+++ /dev/null
@@ -1,38 +0,0 @@
-\name{mkinstart}
-\alias{mkinstart}
-\title{
- Generate starting parameters for optimisations
-}
-\description{
- This function is supposed to analyse a kinetic model and kinetic data in
- order to generate suitable starting parameters for fitting the model.
- It does not do anything really useful at the moment.
-}
-\usage{
-mkinstart(model, data, mode = "auto")
-}
-\arguments{
- \item{model}{
- A kinetic model of class \code{\link{mkinmod}}.
-}
- \item{data}{
- Kinetic data in wide format suitable for fitting with \code{\link{mkinfit}}.
-}
- \item{mode}{
- How the starting parameters should be generated. At the moment, only
- "auto" is supported, which internally uses \code{\link{kinfit}}.
-}
-}
-\details{
-%% ~~ If necessary, more details than the description above ~~
-}
-\value{
- A named vector of starting parameters.
-}
-\author{
- Johannes Ranke
-}
-\examples{
-
-}
-\keyword{ manip }
diff --git a/man/summary.mkinfit.Rd b/man/summary.mkinfit.Rd
index 8a8c6ab..b3e7862 100644
--- a/man/summary.mkinfit.Rd
+++ b/man/summary.mkinfit.Rd
@@ -11,7 +11,7 @@
and residual values.
}
\usage{
-\method{summary}{mkinfit}(object, data = TRUE, distimes = TRUE, ff = TRUE, cov = FALSE, ...)
+\method{summary}{mkinfit}(object, data = TRUE, distimes = TRUE, ff = TRUE, ...)
\method{print}{summary.mkinfit}(x, digits = max(3, getOption("digits") - 3), tval = TRUE, ...)
}
@@ -31,9 +31,6 @@
\item{ff}{
logical, indicating whether formation fractions should be included.
}
- \item{cov}{
- logical, indicating whether parameter covariances should be calculated.
-Passed to \code{\link{summary.modFit}}. }
\item{digits}{
Number of digits to use for printing
}

Contact - Imprint