Merge branch 'master' into iore

5 files changed, 167 insertions, 66 deletions

diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R
index 4137d33a..b28235a6 100644
--- a/R/mkinerrmin.R
+++ b/R/mkinerrmin.R
@@ -36,10 +36,11 @@ mkinerrmin <- function(fit, alpha = 0.05)
     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))
+  # Remove values at time zero for variables whose value for state.ini is fixed,
+  # as these will not have any effect in the optimization and should therefore not 
+  # be counted as degrees of freedom.
+  fixed_initials = gsub("_0$", "", rownames(subset(fit$fixed, type = "state")))
+  errdata <- subset(errdata, !(time == 0 & name %in% fixed_initials))
 
   n.optim.overall <- length(parms.optim)
 
@@ -71,11 +72,12 @@ mkinerrmin <- function(fit, alpha = 0.05)
     # Formation fractions are attributed to the target variable, so look
     # for source compartments with formation fractions
     for (source_var in fit$obs_vars) {
+      n.ff.source = length(grep(paste("^f", source_var, sep = "_"),
+                                 names(parms.optim)))
+      n.paths.source = length(fit$mkinmod$spec[[source_var]]$to)
       for (target_var in fit$mkinmod$spec[[source_var]]$to) {
         if (obs_var == target_var) {
-          n.ff.optim <- n.ff.optim +  
-                        length(grep(paste("^f", source_var, sep = "_"),
-                                    names(parms.optim))) 
+          n.ff.optim <- n.ff.optim + n.ff.source/n.paths.source
         }
       }
     }
diff --git a/R/mkinfit.R b/R/mkinfit.R
index 46121c6d..a6889b0b 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -23,12 +23,13 @@ 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)), 

+  state.ini = "auto", 

   fixed_parms = NULL,

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

   solution_type = "auto",

   method.ode = "lsoda",

-  method.modFit = "Marq",

+  method.modFit = c("Port", "Marq", "SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B"),

+  maxit.modFit = "auto",

   control.modFit = list(),

   transform_rates = TRUE,

   transform_fractions = TRUE,

@@ -40,11 +41,36 @@ mkinfit <- function(mkinmod, observed,
   trace_parms = FALSE,

   ...)

 {

+  # Check mkinmod and generate a model for the variable whith the highest value

+  # if a suitable string is given

+  parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB") 

+  if (class(mkinmod) != "mkinmod") {

+    presumed_parent_name = observed[which.max(observed$value), "name"]

+    if (mkinmod[[1]] %in% parent_models_available) {

+      speclist <- list(list(type = mkinmod, sink = TRUE))

+      names(speclist) <- presumed_parent_name

+      mkinmod <- mkinmod(speclist = speclist)

+    } else {

+      stop("Argument mkinmod must be of class mkinmod or a string containing one of\n  ",

+           paste(parent_models_available, collapse = ", "))

+    } 

+  }

+

+  # Check optimisation method and set maximum number of iterations if specified

+  method.modFit = match.arg(method.modFit)

+  if (maxit.modFit != "auto") {

+    if (method.modFit == "Marq") control.modFit$maxiter = maxit.modFit

+    if (method.modFit == "Port") control.modFit$iter.max = maxit.modFit

+    if (method.modFit %in% c("SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B")) {

+        control.modFit$maxit = maxit.modFit

+    }

+  }

+

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

+  obs_vars <- names(mkinmod$spec)

 

   # Subset observed data with names of observed data in the model

   observed <- subset(observed, name %in% obs_vars)

@@ -128,6 +154,18 @@ mkinfit <- function(mkinmod, observed,
     }

   }

 

+  # Set default for state.ini if appropriate

+  parent_name = names(mkinmod$spec)[[1]]

+  if (state.ini[1] == "auto") {

+    parent_time_0 = subset(observed, time == 0 & name == parent_name)$value

+    parent_time_0_mean = mean(parent_time_0, na.rm = TRUE)

+    if (is.na(parent_time_0_mean)) {

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

+    } else {

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

+    }

+  }

+

   # Name the inital state variable values if they are not named yet

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

 

@@ -245,7 +283,7 @@ mkinfit <- function(mkinmod, observed,
                                 atol = atol, rtol = rtol, ...)

 

         plot(0, type="n", 

-          xlim = range(observed$time), ylim = range(observed$value, na.rm=TRUE),

+          xlim = range(observed$time), ylim = c(0, max(observed$value, na.rm=TRUE)),

           xlab = "Time", ylab = "Observed")

         col_obs <- pch_obs <- 1:length(obs_vars)

         lty_obs <- rep(1, length(obs_vars))

@@ -270,7 +308,7 @@ mkinfit <- function(mkinmod, observed,
   if (!transform_rates) {

     index_k <- grep("^k_", names(lower))

     lower[index_k] <- 0

-    other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2"), names(lower))

+    other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2", "tb"), names(lower))

     lower[other_rate_parms] <- 0

   }

 

@@ -283,41 +321,65 @@ mkinfit <- function(mkinmod, observed,
     upper[other_fraction_parms] <- 1

   }

 

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

-                method = method.modFit, control = control.modFit, 

-                lower = lower, upper = upper, ...)

-

-  # 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, upper = upper, ...)

-      reweight.diff = sum((sqrt(fit$var_ms_unweighted) - sigma.old)^2)

+  # Do the fit and take the time

+  fit_time <- system.time({

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

+                  method = method.modFit, control = control.modFit, 

+                  lower = lower, upper = upper, ...)

+

+    # 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("Iteration", n.iter, "yields variance estimates:\n")

+        cat("Initial variance estimates are:\n")

         print(signif(fit$var_ms_unweighted, 8))

-        cat("Sum of squared differences to last variance estimates:",

-            signif(reweight.diff, 2), "\n")

       }

+      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, upper = upper, ...)

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

+        }

+      }

+    }

+  })

+

+  # Check for convergence

+  if (method.modFit == "Marq") {

+    if (!fit$info %in% c(1, 2, 3)) {

+      fit$warning = paste0("Optimisation by method ", method.modFit, 

+                           " did not converge.\n",

+                           "The message returned by nls.lm is:\n",

+                                    fit$message)

+      warning(fit$warning)

+    }

+  }

+  if (method.modFit %in% c("Port", "SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B")) {

+    if (fit$convergence != 0) {

+      fit$warning = paste0("Optimisation by method ", method.modFit, 

+                           " did not converge.\n",

+                           "Convergence code is ", fit$convergence,

+                           ifelse(is.null(fit$message), "", 

+                                  paste0("\nMessage is ", fit$message)))

+      warning(fit$warning)

     }

   }

 

@@ -325,6 +387,10 @@ mkinfit <- function(mkinmod, observed,
   fit$solution_type <- solution_type

   fit$transform_rates <- transform_rates

   fit$transform_fractions <- transform_fractions

+  fit$method.modFit <- method.modFit

+  fit$maxit.modFit <- maxit.modFit

+  fit$calls <- calls

+  fit$time <- fit_time

 

   # We also need the model for summary and plotting

   fit$mkinmod <- mkinmod

@@ -380,8 +446,11 @@ mkinfit <- function(mkinmod, observed,
   fit$bparms.optim <- bparms.optim 

   fit$bparms.fixed <- bparms.fixed

 

-  # Return ode parameters for further fitting

+  # Return ode and state parameters for further fitting

   fit$bparms.ode <- bparms.all[mkinmod$parms] 

+  fit$bparms.state <- c(bparms.all[setdiff(names(bparms.all), names(fit$bparms.ode))],

+                        state.ini.fixed)

+  names(fit$bparms.state) <- gsub("_0$", "", names(fit$bparms.state))

 

   fit$date <- date()

 

@@ -451,6 +520,8 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05,
 	  date.fit = object$date,

 	  date.summary = date(),

 	  solution_type = object$solution_type,

+	  method.modFit = object$method.modFit,

+	  warning = object$warning,

 	  use_of_ff = object$mkinmod$use_of_ff,

     weight.ini = object$weight.ini,

     reweight.method = object$reweight.method,

@@ -463,6 +534,8 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05,
     cov.scaled = covar * resvar,

     info = object$info, 

     niter = object$iterations,

+    calls = object$calls,

+    time = object$time,

     stopmess = message,

     par = param,

     bpar = bparam)

@@ -493,13 +566,17 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
   cat("Date of fit:    ", x$date.fit, "\n")

   cat("Date of summary:", x$date.summary, "\n")

 

+  if (!is.null(x$warning)) cat("\n\nWarning:", x$warning, "\n\n")

+

   cat("\nEquations:\n")

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

+  writeLines(strwrap(x[["diffs"]], exdent = 11))

   df  <- x$df

   rdf <- df[2]

 

-  cat("\nMethod used for solution of differential equation system:\n")

-  cat(x$solution_type, "\n")

+  cat("\nModel predictions using solution type", x$solution_type, "\n")

+

+  cat("\nFitted with method", x$method.modFit, 

+      "using", x$calls, "model solutions performed in", x$time[["elapsed"]],  "s\n")

 

   cat("\nWeighting:", x$weight.ini)

   if(!is.null(x$reweight.method)) cat(" then iterative reweighting method",

@@ -519,13 +596,15 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
   cat("\nOptimised, transformed parameters:\n")

   print(signif(x$par, 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")

+  if (x$niter != 0) {

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

+    }

   }

 

   cat("\nResidual standard error:",

diff --git a/R/mkinresplot.R b/R/mkinresplot.R
index c9a801fd..82ffd2cb 100644
--- a/R/mkinresplot.R
+++ b/R/mkinresplot.R
@@ -36,7 +36,8 @@ mkinresplot <- function (object,
 	col_obs <- pch_obs <- 1:length(obs_vars)

  	names(col_obs) <- names(pch_obs) <- obs_vars

 

-  plot(0,  xlab = xlab, ylab = ylab, 

+  plot(0, type = "n",

+       xlab = xlab, ylab = ylab, 

        xlim = xlim,

        ylim = c(-1.2 * maxabs, 1.2 * maxabs), ...)

 

diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R
index 4132d96c..31746fb8 100644
--- a/R/plot.mkinfit.R
+++ b/R/plot.mkinfit.R
@@ -31,7 +31,7 @@ plot.mkinfit <- function(x, fit = x,
 {
   if (add && show_residuals) stop("If adding to an existing plot we can not show residuals")
 
-  if (ylim == "default") {
+  if (ylim[[1]] == "default") {
     ylim = c(0, max(subset(fit$data, variable %in% obs_vars)$observed, na.rm = TRUE))
   }
 
diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R
index a36b7eae..778f56cd 100644
--- a/R/transform_odeparms.R
+++ b/R/transform_odeparms.R
@@ -75,12 +75,22 @@ transform_odeparms <- function(parms, mkinmod,
   # and HS parameter tb if transformation of rates is requested

   for (pname in c("alpha", "beta", "k1", "k2", "tb")) {

     if (!is.na(parms[pname])) {

-      transparms[paste0("log_", pname)] <- ifelse(transform_rates, log(parms[pname]), parms[pname])

+      if (transform_rates) {

+        transparms[paste0("log_", pname)] <- log(parms[pname])

+      } else {

+        transparms[pname] <- parms[pname]

+      }

     } 

   }

+

+  # DFOP parameter g is treated as a fraction

   if (!is.na(parms["g"])) {

     g <- parms["g"]

-    transparms["g_ilr"] <- ifelse(transform_fractions, ilr(c(g, 1 - g)), g)

+    if (transform_fractions) {

+      transparms["g_ilr"] <- ilr(c(g, 1 - g))

+    } else {

+      transparms["g"] <- g

+    }

   }

 

   return(transparms)

@@ -142,16 +152,25 @@ backtransform_odeparms <- function(transparms, mkinmod,
 

   # Transform parameters also for FOMC, DFOP and HS models

   for (pname in c("alpha", "beta", "k1", "k2", "tb")) {

-    pname_trans = paste0("log_", pname)

-    if (!is.na(transparms[pname_trans])) {

-      parms[pname] <- ifelse(transform_rates, 

-                             exp(transparms[pname_trans]), 

-                             transparms[pname])

-    } 

+    if (transform_rates) {

+      pname_trans = paste0("log_", pname)

+      if (!is.na(transparms[pname_trans])) {

+        parms[pname] <- exp(transparms[pname_trans])

+      }

+    } else {

+      if (!is.na(transparms[pname])) {

+        parms[pname] <- transparms[pname]

+      }

+    }

   }

+

+  # DFOP parameter g is treated as a fraction

   if (!is.na(transparms["g_ilr"])) {

     g_ilr <- transparms["g_ilr"]

-    parms["g"] <- ifelse(transform_fractions, invilr(g_ilr)[1], g_ilr)

+    parms["g"] <- invilr(g_ilr)[1]

+  }

+  if (!is.na(transparms["g"])) {

+    parms["g"] <- transparms["g"]

   }

 

   return(parms)