Much more complete version that was just submitted to CRAN.

git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@9 edb9625f-4e0d-4859-8d74-9fd3b1da38cb

18 files changed, 694 insertions, 95 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 96f651eb..fb137a76 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,12 +1,11 @@
 Package: mkin

 Type: Package

 Title: Routines for fitting kinetic models with one or more state variables to chemical degradation data

-Version: 0.5-1

-Date: 2010-05-11

+Version: 0.7-1

+Date: 2010-05-18

 Author: Johannes Ranke

 Maintainer: Johannes Ranke <jranke@harlan.com>

-Description: Calculation routines based on the FOCUS Kinetics 

- Report (2006)

+Description: Calculation routines based on the FOCUS Kinetics Report (2006)

 Depends: FME

 License: GPL

 LazyLoad: yes

diff --git a/R/mkin_long_to_wide.R b/R/mkin_long_to_wide.R
new file mode 100644
index 00000000..45ee0a10
--- /dev/null
+++ b/R/mkin_long_to_wide.R
@@ -0,0 +1,9 @@
+mkin_long_to_wide <- function(long_data, time = "time")

+{

+  colnames <- unique(long_data$name)

+  wide_data <- data.frame(time = subset(long_data, name == colnames[1], time))

+  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
new file mode 100644
index 00000000..6db2f337
--- /dev/null
+++ b/R/mkin_wide_to_long.R
@@ -0,0 +1,13 @@
+mkin_wide_to_long <- function(wide_data, time = "t")

+{

+  colnames <- names(wide_data)

+  vars <- subset(colnames, colnames != time)

+  n <- length(colnames) - 1

+  if (!(time %in% colnames)) stop("The data in wide format have to contain a variable named ", time, ".")

+  long_data <- data.frame(

+    name = rep(vars, each = length(wide_data[[time]])),

+    time = rep(wide_data[[time]], n),

+    value = unlist(wide_data[vars]),

+    row.names = NULL)

+  return(long_data)

+}

diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R
new file mode 100644
index 00000000..7a922eb4
--- /dev/null
+++ b/R/mkinerrmin.R
@@ -0,0 +1,14 @@
+mkinerrmin <- function(errdata, n.parms, alpha = 0.05)

+{

+  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))

+}

diff --git a/R/mkinfit.R b/R/mkinfit.R
index 9651fd66..c1490276 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -1,55 +1,78 @@
-mkinfit <- function(mkinmod, observed, 

+mkinfit <- function(mkinmod, observed,

   parms.ini = rep(0.1, length(mkinmod$parms)),

   state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), 

-  fixed_parms = rep(FALSE, length(mkinmod$parms)),

-  fixed_initials = c(FALSE, rep(TRUE, length(mkinmod$diffs) - 1)), 

-  plot = NULL, 

+  fixed_parms = NULL,

+  fixed_initials = names(mkinmod$diffs)[-1],

+  plot = FALSE, quiet = FALSE,

   err = NULL, weight = "none", scaleVar = FALSE,

   ...)

 {

+  mod_vars <- names(mkinmod$diffs)

+  # Subset dataframe 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

   # Create a function calculating the differentials specified by the model

   mkindiff <- function(t, state, parms) {

+    time <- t

     diffs <- vector()

-    for (box in names(mkinmod$diffs))

+    for (box in mod_vars)

     {

       diffname <- paste("d", box, sep="_")      

-      diffs[diffname] <- with(as.list(c(state, parms)),

+      diffs[diffname] <- with(as.list(c(time,state, parms)),

         eval(parse(text=mkinmod$diffs[[box]])))

     }

     return(list(c(diffs)))

   } 

 

   # Name the inital parameter values if they are not named yet

-  if(is.null(names(state.ini))) names(state.ini) <- names(mkinmod$diffs)

+  if(is.null(names(state.ini))) names(state.ini) <- mod_vars

 

-  # TODO: Collect parameters to be optimised

-  parms.optim <- parms.ini[!fixed_parms]

+  # Parameters to be optimised

   parms.fixed <- parms.ini[fixed_parms]

+  optim_parms <- setdiff(names(parms.ini), fixed_parms)

+  parms.optim <- parms.ini[optim_parms]

 

-  state.ini.optim <- state.ini[!fixed_initials]

-  state.ini.optim.boxnames <- names(state.ini.optim)

-  names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_")

   state.ini.fixed <- state.ini[fixed_initials]

+  optim_initials <- setdiff(names(state.ini), fixed_initials)

+  state.ini.optim <- state.ini[optim_initials]

+  state.ini.optim.boxnames <- names(state.ini.optim)

+  if(length(state.ini.optim) > 0) {

+      names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_")

+  }

 

+  cost.old <- 1e100

+  calls <- 0

+  out_predicted <- NA

   # Define the model cost function

   cost <- function(P)

   {

+    assign("calls", calls+1, inherits=TRUE)

     if(length(state.ini.optim) > 0) {

       odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed)

       names(odeini) <- c(state.ini.optim.boxnames, names(state.ini.fixed))

     } else odeini <- state.ini.fixed

 

     odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed)

-    # Solve the ODE

+

+    # Get more timepoints if plotting is desired

+    if(plot) {

+      outtimes = unique(sort(c(observed$time, 

+        seq(min(observed$time), max(observed$time), length.out=100))))

+    } else outtimes = unique(observed$time)

+

+    # Solve the ode

     out <- ode(

       y = odeini,

-      times = unique(observed$time),

+      times = outtimes,

       func = mkindiff, 

       parms = odeparms)

      

-    # Output transformation for models with ghost compartments like SFORB

+  

+    # Output transformation for models with unobserved compartments like SFORB

     out_transformed <- data.frame(time = out[,"time"])

     for (var in names(mkinmod$map)) {

       if(length(mkinmod$map[[var]]) == 1) {

@@ -58,9 +81,190 @@ mkinfit <- function(mkinmod, observed,
         out_transformed[var] <- rowSums(out[, mkinmod$map[[var]]])

       }

     }    

+    assign("out_predicted", subset(out_transformed, time %in% observed$time), inherits=TRUE)

+

+    mC <- modCost(out_transformed, 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 the current model output if requested

+      if(plot) {

+        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)

+        names(col_obs) <- names(pch_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_transformed$time, out_transformed[-1])

+        legend("topright", inset=c(0.05, 0.05), legend=obs_vars, 

+          col=col_obs, pch=pch_obs, lty=1:length(pch_obs))

+      }

     

-    return(modCost(out_transformed, observed, y = "value",

-      err = err, weight = weight, scaleVar = scaleVar))

+      assign("cost.old", mC$model, inherits=TRUE)

+    }

+    return(mC)

   }

-  modFit(cost, c(state.ini.optim, parms.optim), ...)

+  fit <- modFit(cost, c(state.ini.optim, parms.optim), ...)

+

+  # We need the function for plotting

+  fit$mkindiff <- mkindiff

+

+  # We also need various other information for summary and plotting

+  fit$map <- mkinmod$map

+  fit$diffs <- mkinmod$diffs

+  fit$observed <- mkin_long_to_wide(observed)

+  predicted_long <- mkin_wide_to_long(out_predicted, time = "time")

+  fit$predicted <- out_predicted

+

+  # Collect initial parameter values in two dataframes

+  fit$start <- data.frame(initial = c(state.ini.optim, parms.optim))

+  fit$start$type = c(rep("state", length(state.ini.optim)), rep("deparm", length(parms.optim)))

+  if (exists("lower")) fit$start$lower <- lower

+  if (exists("upper")) fit$start$upper <- upper

+

+  fit$fixed <- data.frame(

+    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

+    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

+  parms.all = c(fit$par, parms.fixed)

+  fit$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), DT90 = rep(NA, length(obs_vars)), 

+    row.names = obs_vars)

+  for (obs_var in obs_vars) {

+      type = names(mkinmod$map[[obs_var]])[1]  

+      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

+      }

+      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 == "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, 1000))$minimum

+        if (abs(DT90.o - max_DT) < 0.01) DT90 = NA else DT90 = DT90.o

+      }

+      fit$distimes[obs_var, ] = c(DT50, DT90)

+  }

+

+  # Collect observed, predicted and residuals

+  data <- merge(observed, predicted_long, 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)

+  fit$data <- data[order(data$variable, data$time), ]

+

+  class(fit) <- c("mkinfit", "modFit") 

+  return(fit)

+}

+

+summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, cov = FALSE,...) {

+  ans <- FME:::summary.modFit(object, cov = cov)

+  ans$diffs <- object$diffs

+  if(data) ans$data <- object$data

+  ans$start <- object$start

+  ans$fixed <- object$fixed

+  ans$errmin <- object$errmin 

+  if(distimes) ans$distimes <- object$distimes

+  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("\nEquations:\n")

+  print(noquote(as.character(x[["diffs"]])))

+  df  <- x$df

+  rdf <- df[2]

+

+  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 parameters:\n")

+  printCoefmat(x$par, 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,...)

+  }    

+

+  printcor <- !is.null(x$cov.unscaled)

+  if (printcor){

+    Corr <- cov2cor(x$cov.unscaled)

+    rownames(Corr) <- colnames(Corr) <- rownames(x$par)

+    cat("\nParameter correlation:\n")

+    print(Corr, digits = digits, ...)

+  }

+

+  printdata <- !is.null(x$data)

+  if (printdata){

+    cat("\nData:\n")

+    print(format(x$data, digits = digits, scientific = FALSE,...), row.names = FALSE)

+  }

+

+  invisible(x)

 }

diff --git a/R/mkinmod.R b/R/mkinmod.R
index 7e64f879..643e855e 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -1,38 +1,53 @@
-mkinmod <- function(spec = list(parent = list(type = "SFO", to = NA, sink = TRUE)))

+mkinmod <- function(...)

 {

-  if(!is.list(spec)) stop("spec must be a list of model specifications for each observed variable")

+  spec <- list(...)

+  obs_vars <- names(spec)

 

-  # The returned model will be a list of two character vectors, containing

-  # differential equations and parameter names

+  # The returned model will be a list of character vectors, containing

+  # differential equations, parameter names and a mapping from model variables

+  # to observed variables

   parms <- vector()

   diffs <- vector()

   map <- list()

 

   # Establish list of differential equations

-  for (varname in names(spec))

+  for (varname in obs_vars)

   {

-    if(!spec[[varname]]$type %in% c("SFO", "SFORB")) stop("Available types are SFO and SFORB only")

+    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", "SFORB")) stop("Available types are SFO, FOMC and SFORB only")

     new_parms <- vector()

 

     # New (sub)compartments (boxes) needed for the model type

     new_boxes <- switch(spec[[varname]]$type,

       SFO = varname,

+      FOMC = 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="")

 

     # Construct terms for transfer to sink and add if appropriate

+    if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE

     if(spec[[varname]]$sink) {

-      k_compound_sink <- paste("k", new_boxes[[1]], "sink", sep="_")

-      sink_term <- paste("-", k_compound_sink, "*", new_boxes[[1]])

-      # Only add sink term to first (or only) box for SFO and SFORB models

+      # Add first-order sink term to first (or only) box for SFO and SFORB models

       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

+      }

+      if(spec[[varname]]$type == "FOMC") {

+        if(match(varname, obs_vars) != 1) {

+          stop("Type FOMC is only allowed for the first compartment, which is assumed to be the source compartment")

+        }

+        # From p. 53 of the FOCUS kinetics report

+        fomc_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", new_boxes[[1]])

+        new_diffs[[1]] <- paste(new_diffs[[1]], "-", fomc_term)

+        new_parms <- c("alpha", "beta")

       }

-      new_parms <- k_compound_sink

     }

    

     # Add reversible binding if appropriate

@@ -53,24 +68,41 @@ mkinmod <- function(spec = list(parent = list(type = "SFO", to = NA, sink = TRUE
   }

 

   # Transfer between compartments

-  for (varname in names(spec)) {

+  for (varname in obs_vars) {

     to <- spec[[varname]]$to

-    if(!is.na(to)) {

+    if(!is.null(to)) {

       origin_box <- switch(spec[[varname]]$type,

         SFO = varname,

+        FOMC = varname,

         SFORB = paste(varname, "free", sep="_"))

+      fraction_left <- NULL

       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 == "FOMC") {

+          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="")

+          }

+          diffs[[target_box]] <- paste(diffs[[target_box]], "+", fraction_really_to_target, "*", fomc_term)

+          parms <- c(parms, fraction_to_target)

+        }

       }

-      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)

     }

   }

   model <- list(diffs = diffs, parms = parms, map = map)

   class(model) <- "mkinmod"

-  return(model)

+  invisible(model)

 }

diff --git a/R/mkinplot.R b/R/mkinplot.R
new file mode 100644
index 00000000..d2880c0e
--- /dev/null
+++ b/R/mkinplot.R
@@ -0,0 +1,49 @@
+mkinplot <- function(fit, xlab = "Time", ylab = "Observed", xlim = range(fit$data$time), ylim = range(fit$data$observed, na.rm = TRUE), ...)
+{
+  fixed <- fit$fixed$value
+  names(fixed) <- rownames(fit$fixed)
+  parms.all <- c(fit$par, fixed)
+  ininames <- c(
+    rownames(subset(fit$start, type == "state")),
+    rownames(subset(fit$fixed, type == "state")))
+  odeini <- parms.all[ininames]
+  names(odeini) <- names(fit$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]
+
+  # Solve the ode
+  out <- ode(
+    y = odeini,
+    times = outtimes,
+    func = fit$mkindiff, 
+    parms = odeparms)
+    
+  # Output transformation for models with unobserved compartments like SFORB
+  out_transformed <- data.frame(time = out[,"time"])
+  for (var in names(fit$map)) {
+    if(length(fit$map[[var]]) == 1) {
+      out_transformed[var] <- out[, var]
+    } else {
+      out_transformed[var] <- rowSums(out[, fit$map[[var]]])
+    }
+  }    
+
+  # Plot the data and model output
+  plot(0, type="n", 
+    xlim = xlim, ylim = ylim,
+      xlab = xlab, ylab = ylab, ...)
+    col_obs <- pch_obs <- 1:length(fit$map)
+    names(col_obs) <- names(pch_obs) <- names(fit$map)
+    for (obs_var in names(fit$map)) {
+      points(subset(fit$data, variable == obs_var, c(time, observed)), 
+        pch = pch_obs[obs_var], col = col_obs[obs_var])
+    }
+    matlines(out_transformed$time, out_transformed[-1])
+      legend("topright", inset=c(0.05, 0.05), legend=names(fit$map),
+      col=col_obs, pch=pch_obs, lty=1:length(pch_obs))
+}
diff --git a/TODO b/TODO
new file mode 100644
index 00000000..bfefe0c6
--- /dev/null
+++ b/TODO
@@ -0,0 +1,4 @@
+- Fix DT50 and DT90 calculation for SFORB_SFO
+- Add unit tests for mkinfit
+- Document validation against fits documented in chapter 13 of FOCUS (2006)
+- Reproduce example anaylses (L1, L2, ...) in FOCUS (2006)
diff --git a/inst/doc/mkin.Rnw b/inst/doc/mkin.Rnw
index 84ede0f5..78b59e67 100644
--- a/inst/doc/mkin.Rnw
+++ b/inst/doc/mkin.Rnw
@@ -112,14 +112,14 @@ The next task is to define the model to be fitted to the data. In order to
 facilitate this task, a convenience function \Robject{mkinmod} is available.

 

 <<model_definition, echo=TRUE>>=

-SFO <- mkinmod(spec = list(parent = list(type = "SFO", to = NA, sink = TRUE)))

-SFORB <- mkinmod(spec = list(parent = list(type = "SFORB", to = NA, sink = TRUE)))

-SFO_SFO <- mkinmod(spec = list(

+SFO <- mkinmod(parent = list(type = "SFO"))

+SFORB <- mkinmod(parent = list(type = "SFORB"))

+SFO_SFO <- mkinmod(

   parent = list(type = "SFO", to = "m1", sink = TRUE),

-  m1 = list(type = "SFO", to = NA, sink = TRUE)))

-SFORB_SFO <- mkinmod(spec = list(

+  m1 = list(type = "SFO"))

+SFORB_SFO <- mkinmod(

   parent = list(type = "SFORB", to = "m1", sink = TRUE),

-  m1 = list(type = "SFO", to = NA, sink = TRUE)))

+  m1 = list(type = "SFO"))

 @

 

 \subsection{Fitting the model}

@@ -138,14 +138,6 @@ SFORB.fit <- mkinfit(SFORB, FOCUS_2006_C)
 summary(SFORB.fit)

 SFO_SFO.fit <- mkinfit(SFO_SFO, FOCUS_2006_D)

 summary(SFO_SFO.fit)

-SFO_SFO.fit.2 <- mkinfit(SFO_SFO, FOCUS_2006_D, 

-  fixed_initials = c(FALSE, FALSE), fixed_parms = c(FALSE, TRUE, FALSE))

-summary(SFO_SFO.fit.2)

-SFO_SFO.fit.3 <- mkinfit(SFO_SFO, FOCUS_2006_D, 

-  fixed_initials = c(FALSE, FALSE), fixed_parms = c(FALSE, TRUE, FALSE), lower = -0.0000001)

-summary(SFO_SFO.fit.3)

-SFORB_SFO.fit <- mkinfit(SFORB_SFO, FOCUS_2006_D)

-summary(SFORB_SFO.fit)

 @

 

 \bibliographystyle{plainnat}

diff --git a/inst/doc/mkin.pdf b/inst/doc/mkin.pdf
index 95c8a690..bc72e740 100644
--- a/inst/doc/mkin.pdf
+++ b/inst/doc/mkin.pdf
diff --git a/inst/unitTests/runit.mkinmod.R b/inst/unitTests/runit.mkinmod.R
index b6ca6b8d..44d9ffb5 100644
--- a/inst/unitTests/runit.mkinmod.R
+++ b/inst/unitTests/runit.mkinmod.R
@@ -4,13 +4,25 @@ test.mkinmod.SFO <- function()
     parent = "d_parent = - k_parent_sink * parent"

   )

   SFO.parms <- c("k_parent_sink")

-  SFO.map <- list(parent = "parent")

+  SFO.map <- list(parent = c(SFO = "parent"))

   SFO <- list(diffs = SFO.diffs, parms = SFO.parms, map = SFO.map)

   class(SFO) <- "mkinmod"

-  SFO.mkinmod <- mkinmod(spec = list(

-    parent = list(type = "SFO", to = NA, sink=TRUE))

+  SFO.1 <- mkinmod(

+    parent = list(type = "SFO", to = NULL, sink = TRUE)

   )

-  checkIdentical(SFO, SFO.mkinmod)

+  checkIdentical(SFO, SFO.1)

+  SFO.2 <- mkinmod(

+    parent = list(type = "SFO", to = NULL)

+  )

+  checkIdentical(SFO, SFO.2)

+  SFO.3 <- mkinmod(

+    parent = list(type = "SFO", sink = TRUE)

+  )

+  checkIdentical(SFO, SFO.3)

+  SFO.4 <- mkinmod(

+    parent = list(type = "SFO")

+  )

+  checkIdentical(SFO, SFO.3)

 }

 

 test.mkinmod.SFORB <- function()

@@ -26,11 +38,11 @@ test.mkinmod.SFORB <- function()
         "- k_parent_bound_free * parent_bound")

   )

   SFORB.parms <- c("k_parent_free_sink", "k_parent_free_bound", "k_parent_bound_free")

-  SFORB.map <- list(parent = c("parent_free", "parent_bound"))

+  SFORB.map <- list(parent = c(SFORB = "parent_free", SFORB = "parent_bound"))

   SFORB <- list(diffs = SFORB.diffs, parms = SFORB.parms, map = SFORB.map)

   class(SFORB) <- "mkinmod"

-  SFORB.mkinmod <- mkinmod(spec = list(

-    parent = list(type = "SFORB", to = NA, sink=TRUE))

+  SFORB.mkinmod <- mkinmod(

+    parent = list(type = "SFORB", to = NULL, sink=TRUE)

   )

   checkIdentical(SFORB, SFORB.mkinmod)

 }

@@ -42,12 +54,50 @@ test.mkinmod.SFO_SFO <- function()
     m1 = "d_m1 = - k_m1_sink * m1 + k_parent_m1 * parent"

   )

   SFO_SFO.parms <- c("k_parent_sink", "k_m1_sink", "k_parent_m1")

-  SFO_SFO.map <- list(parent = "parent", m1 = "m1")

+  SFO_SFO.map <- list(parent = c(SFO = "parent"), m1 = c(SFO = "m1"))

   SFO_SFO <- list(diffs = SFO_SFO.diffs, parms = SFO_SFO.parms, map = SFO_SFO.map)

   class(SFO_SFO) <- "mkinmod"

-  SFO_SFO.mkinmod <- mkinmod(spec = list(

+  SFO_SFO.mkinmod <- mkinmod(

     parent = list(type = "SFO", to = "m1", sink=TRUE),

-    m1 = list(type = "SFO", to = NA, sink=TRUE))

+    m1 = list(type = "SFO", sink=TRUE)

   )

   checkIdentical(SFO_SFO, SFO_SFO.mkinmod)

 }

+

+test.mkinmod.SFO_SFO2 <- function()

+{  

+  SFO_SFO2.diffs <- c(

+    parent = "d_parent = - k_parent_sink * parent - k_parent_m1 * parent - k_parent_m2 * parent",

+    m1 = "d_m1 = - k_m1_sink * m1 + k_parent_m1 * parent",

+    m2 = "d_m2 = - k_m2_sink * m2 + k_parent_m2 * parent"

+  )

+  SFO_SFO2.parms <- c("k_parent_sink", "k_m1_sink", "k_m2_sink", "k_parent_m1", "k_parent_m2")

+  SFO_SFO2.map <- list(parent = c(SFO = "parent"), m1 = c(SFO = "m1"), m2 = c(SFO = "m2"))

+  SFO_SFO2 <- list(diffs = SFO_SFO2.diffs, parms = SFO_SFO2.parms, map = SFO_SFO2.map)

+  class(SFO_SFO2) <- "mkinmod"

+  SFO_SFO2.mkinmod <- mkinmod(

+    parent = list(type = "SFO", to = c("m1", "m2"), sink=TRUE),

+    m1 = list(type = "SFO", sink=TRUE),

+    m2 = list(type = "SFO", sink=TRUE)

+  )

+  checkIdentical(SFO_SFO2, SFO_SFO2.mkinmod)

+}

+

+test.mkinmod.FOMC_SFO2 <- function()

+{  

+  FOMC_SFO2.diffs <- c(

+    parent = "d_parent = - (alpha/beta) * ((time/beta) + 1)^-1 * parent",

+    m1 = "d_m1 = - k_m1_sink * m1 + f_to_m1 * (alpha/beta) * ((time/beta) + 1)^-1 * parent",

+    m2 = "d_m2 = - k_m2_sink * m2 + (1 - f_to_m1) * f_to_m2 * (alpha/beta) * ((time/beta) + 1)^-1 * parent"

+  )

+  FOMC_SFO2.parms <- c("alpha", "beta", "k_m1_sink", "k_m2_sink", "f_to_m1", "f_to_m2")

+  FOMC_SFO2.map <- list(parent = c(FOMC = "parent"), m1 = c(SFO = "m1"), m2 = c(SFO = "m2"))

+  FOMC_SFO2 <- list(diffs = FOMC_SFO2.diffs, parms = FOMC_SFO2.parms, map = FOMC_SFO2.map)

+  class(FOMC_SFO2) <- "mkinmod"

+  FOMC_SFO2.mkinmod <- mkinmod(

+    parent = list(type = "FOMC", to = c("m1", "m2"), sink=TRUE),

+    m1 = list(type = "SFO"),

+    m2 = list(type = "SFO")

+  )

+  checkIdentical(FOMC_SFO2, FOMC_SFO2.mkinmod)

+}

diff --git a/man/mkin_long_to_wide.Rd b/man/mkin_long_to_wide.Rd
new file mode 100644
index 00000000..97791e6d
--- /dev/null
+++ b/man/mkin_long_to_wide.Rd
@@ -0,0 +1,32 @@
+\name{mkin_long_to_wide}
+\alias{mkin_long_to_wide}
+\title{
+  Convert a dataframe from long to wide format.
+}
+\usage{
+mkin_long_to_wide(long_data, time = "time")
+}
+\description{
+  This function takes a dataframe in the long form as required by \code{\link{modCost}} 
+  and converts it into a dataframe with one independent variable and several 
+  dependent variables as columns.
+}
+\arguments{
+  \item{long_data}{
+  The dataframe must contain one variable with the time values specified by the 
+  \code{time} argument and one column of observed values named "value".
+}
+  \item{time}{
+  The name of the time variable.
+}
+}
+\value{
+  Dataframe in wide format.
+}
+\author{
+  Johannes Ranke
+}
+\examples{
+mkin_long_to_wide(FOCUS_2006_D)
+}
+\keyword{ manip }
diff --git a/man/mkin_wide_to_long.Rd b/man/mkin_wide_to_long.Rd
new file mode 100644
index 00000000..d3dd200a
--- /dev/null
+++ b/man/mkin_wide_to_long.Rd
@@ -0,0 +1,32 @@
+\name{mkin_wide_to_long}
+\alias{mkin_wide_to_long}
+\title{
+  Convert a dataframe with observations over time into long format.
+}
+\usage{
+mkin_wide_to_long(wide_data, time = "t")
+}
+\description{
+  This function simply takes a dataframe with one independent variable and several 
+  dependent variable and converts it into the long form as required by \code{\link{modCost}}.
+}
+\arguments{
+  \item{wide_data}{
+  The dataframe must contain one variable with the time values specified by the 
+  \code{time} argument and usually more than one column of observed values.
+}
+  \item{time}{
+  The name of the time variable.
+}
+}
+\value{
+  Dataframe in long format as needed for \code{\link{modCost}}.
+}
+\author{
+  Johannes Ranke
+}
+\examples{
+wide <- data.frame(t = c(1,2,3), x = c(1,4,7), y = c(3,4,5))
+mkin_wide_to_long(wide)
+}
+\keyword{ manip }
diff --git a/man/mkinerrmin.Rd b/man/mkinerrmin.Rd
new file mode 100644
index 00000000..4bbf96ab
--- /dev/null
+++ b/man/mkinerrmin.Rd
@@ -0,0 +1,46 @@
+\name{mkinerrmin}

+\Rdversion{1.1}

+\alias{mkinerrmin}

+\title{

+Calculate the minimum error to assume in order to pass the variance test

+}

+\description{

+This function uses \code{\link{optimize}} in order to iteratively find the

+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)

+}

+\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{alpha}{

+    The confidence level chosen for the chi-squared test.

+}

+}

+\value{

+  A list 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.}

+}

+\details{

+    This function is used internally by \code{\link{mkinfit}}.

+}

+\references{ 

+  FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and

+  Degradation Kinetics from Environmental Fate Studies on Pesticides in EU

+  Registration} Report of the FOCUS Work Group on Degradation Kinetics,

+  EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,

+  \url{http://focus.jrc.ec.europa.eu/dk} 

+}

+\author{ Johannes Ranke }

+\keyword{ manip }

diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd
index 7ddc10c1..640207ce 100644
--- a/man/mkinfit.Rd
+++ b/man/mkinfit.Rd
@@ -10,7 +10,7 @@
   values.
 }
 \usage{
-mkinfit(mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$parms)), state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), fixed_parms = rep(FALSE, length(mkinmod$parms)), fixed_initials = c(FALSE, rep(TRUE, length(mkinmod$diffs) - 1)), plot = NULL, err = NULL, weight = "none", scaleVar = FALSE, ...)
+mkinfit(mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$parms)), state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], plot = FALSE, quiet = FALSE, err = NULL, weight = "none", scaleVar = FALSE, ...)
 }
 \arguments{
   \item{mkinmod}{
@@ -29,56 +29,59 @@ mkinfit(mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$parms)), state.in
   \item{parms.ini}{
     A named vector if initial values for the parameters, including both parameters to
     be optimised and potentially also fixed parameters as indicated by \code{fixed_parms}.
-}
+  }
   \item{state.ini}{
     A named vector of initial values for the state variables of the model. In case the 
     observed variables are represented by more than one model variable, the names will
     differ from the names of the observed variables (see \code{map} component of 
     \code{\link{mkinmod}}). The default is to set the initial value of the first model
     variable to 100 and all others to 0.
-}
+  }
   \item{fixed_parms}{
-    A vector of booleans specifying which parameters are not to be optimised. The default
-    is to include all model parameters in the optimisation.
-}
+    The names of parameters that should not be optimised but rather kept at the values
+    specified in \code{parms.ini}.
+  }
   \item{fixed_initials}{
-    A vector of booleans specifying which initial values to include in the optimisation.
-    The default is to optimise the initial value of the first model variable and to
-    keep all other initial values fixed.
-}
+    The names of model variables for which the initial state at time 0 should be excluded
+    from the optimisation. Defaults to all state variables except for the first one.
+  }
   \item{plot}{
     Should the observed values and the numerical solutions be plotted at each stage
     of the optimisation?
-}
+  }
+  \item{quiet}{
+    Suppress printing out the current model cost after each improvement?
+  }
   \item{err }{either \code{NULL}, or the name of the column with the
     \emph{error} estimates, used to weigh the residuals (see details of
     \code{\link{modCost}}); if \code{NULL}, then the residuals are not weighed.
-}
+  }
   \item{weight}{only if \code{err}=\code{NULL}: how to weigh the
     residuals, one of "none", "std", "mean", see details of \code{\link{modCost}}.
-}
+  }
   \item{scaleVar}{
     Will be passed to \code{\link{modCost}}. Default is not to scale Variables according
     to the number of observations.
-}
+  }
   \item{\dots}{
     Further arguments that will be passed to \code{\link{modFit}}. 
-}
+  }
 }
 \value{
-    A list of class \code{\link{modFit}}. Thus, at present, a summary can be obtained
-    by \code{\link{summary.modFit}}. 
+    A list with "mkinfit" and "modFit" in the class attribute. Thus, at
+    present, a summary can be obtained by \code{\link{summary.modFit}}. 
 }
 \author{
   Johannes Ranke <jranke@{harlan.com,uni-bremen.de}>
 }
 \examples{
 # One parent compound, one metabolite, both single first order.
-SFO_SFO <- mkinmod(spec = list(
+SFO_SFO <- mkinmod(
   parent = list(type = "SFO", to = "m1", sink = TRUE),
-  m1 = list(type = "SFO", to = NA, sink = TRUE)))
+  m1 = list(type = "SFO"))
 # Fit the model to the FOCUS example dataset D using defaults
 fit <- mkinfit(SFO_SFO, FOCUS_2006_D)
+str(fit)
 summary(fit)
 }
 \keyword{ models }
diff --git a/man/mkinmod.Rd b/man/mkinmod.Rd
index 4569b317..7ec96277 100644
--- a/man/mkinmod.Rd
+++ b/man/mkinmod.Rd
@@ -9,19 +9,19 @@
   kinetic model type and reaction or transfer to other observed compartments.

 }

 \usage{

-mkinmod(spec = list(parent = list(type = "SFO", to = NA, sink = TRUE)))

+mkinmod(...)

 }

 \arguments{

-  \item{spec}{

-    A list of observed variables to be modelled. Each observed variable has to be 

-    represented by a list with the following entries:

-    \code{type}{ The type of kinetics to use for the variable. Currently, only

-      single first order kinetics "SFO" or single first order with reversible binding

-      "SFORB" are implemented. }

-    \code{to}{ A vector of names of variables to which a transfer is to be assumed

-      in the model. }

-    \code{sink}{ Boolean, specifying if transformation to unspecified compounds (sink) 

-      is to be assumed in the model. }

+  \item{...}{

+    For each observed variable, a list has to be specified as an argument, containing

+    at least a component \code{type}, specifying the type of kinetics to use

+    for the variable. Currently, only single first order kinetics "SFO" or

+    single first order with reversible binding "SFORB" are implemented, as well as

+    "FOMC" for the first compartment which is assumed to be the source compartment.

+    Optional components of each argument are \code{to}, a vector of names of

+    variables to which a transfer is to be assumed in the model, and

+    \code{sink}, a logical specifying if transformation to unspecified

+    compounds (sink) is to be assumed in the model (defaults to \code{TRUE})

   }

 }

 \value{

@@ -36,9 +36,14 @@ mkinmod(spec = list(parent = list(type = "SFO", to = NA, sink = TRUE)))
   Johannes Ranke <jranke@{harlan.com,uni-bremen.de}>

 }

 \examples{

+# There are different ways to specify the SFO model

+SFO.1 <- mkinmod(parent = list(type = "SFO", to = NULL, sink = TRUE))

+SFO.2 <- mkinmod(parent = list(type = "SFO"))

+all.equal(SFO.1, SFO.2)

+

 # One parent compound, one metabolite, both single first order.

-SFO_SFO <- mkinmod(spec = list(

+SFO_SFO <- mkinmod(

   parent = list(type = "SFO", to = "m1", sink = TRUE),

-  m1 = list(type = "SFO", to = NA, sink = TRUE)))

+  m1 = list(type = "SFO"))

 }

 \keyword{ models }

diff --git a/man/mkinplot.Rd b/man/mkinplot.Rd
new file mode 100644
index 00000000..58959d77
--- /dev/null
+++ b/man/mkinplot.Rd
@@ -0,0 +1,50 @@
+\name{mkinplot}
+\alias{mkinplot}
+\title{
+  Plot the observed data and the fitted model of an mkinfit.
+}
+\description{
+  Solves the differential equations with the optimised and fixed parameters
+  from a previous successful call to \code{\link{mkinfit}} and plots 
+  the observed data together with the numerical solution of the fitted model.
+}
+\usage{
+  mkinplot(fit, xlab = "Time", ylab = "Observed", 
+    xlim = range(fit$data$time), ylim = range(fit$data$observed, na.rm=TRUE), ...)
+}
+\arguments{
+  \item{fit}{
+    an object of class \code{\link{mkinfit}}.
+  }
+  \item{xlab}{
+    label for the x axis.
+  }
+  \item{ylab}{
+    label for the y axis.
+  }
+  \item{xlim}{
+    plot range in x direction.
+  }
+  \item{ylim}{
+    plot range in y direction.
+  }
+  \item{\dots}{
+   further arguments passed to \code{\link{plot}}.
+}
+}
+\value{
+  The function is called for its side effect.
+}
+\examples{
+# One parent compound, one metabolite, both single first order.
+SFO_SFO <- mkinmod(
+  parent = list(type = "SFO", to = "m1", sink = TRUE),
+  m1 = list(type = "SFO"))
+# Fit the model to the FOCUS example dataset D using defaults
+fit <- mkinfit(SFO_SFO, FOCUS_2006_D)
+\dontrun{mkinplot(fit)}
+}
+\author{
+  Johannes Ranke
+}
+\keyword{ hplot }
diff --git a/man/summary.mkinfit.Rd b/man/summary.mkinfit.Rd
new file mode 100644
index 00000000..1c305899
--- /dev/null
+++ b/man/summary.mkinfit.Rd
@@ -0,0 +1,65 @@
+\name{summary.mkinfit}

+\alias{summary.mkinfit}

+\alias{print.summary.mkinfit}

+\title{

+  Summary method for class "mkinfit".

+}

+\description{

+  Lists model equations, the summary as returned by \code{\link{summary.modFit}},

+  the chi2 error levels calculated according to FOCUS guidance (2006) as far

+  as defined therein, and optionally the data, consisting of observed, predicted

+  and residual values.

+}

+\usage{

+\method{summary}{mkinfit}(object, data = TRUE, distimes = TRUE, cov = FALSE, ...)

+\method{print}{summary.mkinfit}(x, digits = max(3, getOption("digits") - 3), ...)

+}

+

+\arguments{

+  \item{object}{

+  an object of class \code{\link{mkinfit}}.

+}

+  \item{x}{

+  an object of class \code{summary.mkinfit}.

+}

+  \item{data}{

+  logical, indicating whether the data should be included in the summary. 

+}

+  \item{distimes}{

+  logical, indicating whether DT50 and DT90 values shoule 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

+}

+  \item{\dots}{

+  optional arguments passed to methods like \code{print}.

+}

+}

+\value{

+  The summary function returns a list with the same components as

+  \code{\link{summary.modFit}}, and the additional components

+  \item{diffs }{The differential equations used in the model}

+  \item{data }{The data (see Description above).}

+  \item{start }{The starting values and bounds, if applicable, for optimised parameters.}

+  \item{fixed }{The values of fixed parameters.}

+  \item{errmin }{The chi2 error levels for each observed variable.}

+  \item{disstimes }{The DT50 and DT90 values for each observed variable.}

+  The print method is called for its side effect, i.e. printing the summary.

+}

+\references{

+  FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and

+  Degradation Kinetics from Environmental Fate Studies on Pesticides in EU

+  Registration} Report of the FOCUS Work Group on Degradation Kinetics,

+  EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,

+  \url{http://focus.jrc.ec.europa.eu/dk} 

+}

+\author{

+  Johannes Ranke <jranke@harlan.com>

+}

+\examples{

+  summary(mkinfit(mkinmod(parent = list(type = "SFO")), FOCUS_2006_A))

+}

+\keyword{ utilities }