From c1144753adfa0809003085009ebd85f8af9beda8 Mon Sep 17 00:00:00 2001
From: jranke <jranke@edb9625f-4e0d-4859-8d74-9fd3b1da38cb>
Date: Tue, 10 Apr 2012 21:50:22 +0000
Subject: - Fitting and summaries now work with the new parameter
 transformations. - The SFORB models with metabolites is broken (see TODO) -
 Moved the vignette to the location recommended since R 2.14 - Added the
 missing documentation - Commented out the schaefer_complex_case test, as this
 version of   mkin is not able to fit a model without sink and therefore  
 mkin estimated parameters are quite different

git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@22 edb9625f-4e0d-4859-8d74-9fd3b1da38cb
---
 R/mkinfit.R            | 69 ++++++++++++++++----------------------
 R/mkinpredict.R        | 90 ++++++++++++++++++++++++++++++++++++++++++++++++++
 R/predict.mkinmod.R    | 89 -------------------------------------------------
 R/transform_odeparms.R | 70 +++++++++++++++++++++++++++------------
 4 files changed, 167 insertions(+), 151 deletions(-)
 create mode 100644 R/mkinpredict.R
 delete mode 100644 R/predict.mkinmod.R

(limited to 'R')

diff --git a/R/mkinfit.R b/R/mkinfit.R
index 47c39cf..44eb5b2 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -57,15 +57,17 @@ mkinfit <- function(mkinmod, observed,
     if (parmname == "tb") parms.ini[parmname] = 5
     if (parmname == "g") parms.ini[parmname] = 0.5
   }
-  parms.ini <- transform_odeparms(parms.ini, mod_vars)
-   
+
   # Name the inital state variable values if they are not named yet
   if(is.null(names(state.ini))) names(state.ini) <- mod_vars
 
+  # Transform initial parameter values for fitting
+  transparms.ini <- transform_odeparms(parms.ini, mod_vars)
+
   # Parameters to be optimised:
   # Kinetic parameters in parms.ini whose names are not in fixed_parms
-  parms.fixed <- parms.ini[fixed_parms]
-  parms.optim <- parms.ini[setdiff(names(parms.ini), fixed_parms)]
+  parms.fixed <- transparms.ini[fixed_parms]
+  parms.optim <- transparms.ini[setdiff(names(transparms.ini), fixed_parms)]
 
   # Inital state variables in state.ini whose names are not in fixed_initials
   state.ini.fixed <- state.ini[fixed_initials]
@@ -83,8 +85,11 @@ mkinfit <- function(mkinmod, observed,
   if (length(mkinmod$map) == 1) {
     solution = "analytical"
   } else {
-    if (is.matrix(mkinmod$coefmat) & eigen) solution = "eigen"
-    else solution = "deSolve"
+    if (is.matrix(mkinmod$coefmat) && eigen) {
+      solution = "eigen"
+    } else {
+      solution = "deSolve"
+    }
   }
 
   cost.old <- 1e100 # The first model cost should be smaller than this value
@@ -105,10 +110,10 @@ mkinfit <- function(mkinmod, observed,
 
     odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed)
 
-    transparms <- transform_odeparms(odeparms, mod_vars)
+    parms <- backtransform_odeparms(odeparms, mod_vars)
 
     # Solve the system with current transformed parameter values
-    out <- predict(mkinmod, transparms, odeini, outtimes, solution_type = solution)
+    out <- mkinpredict(mkinmod, parms, odeini, outtimes, solution_type = solution, ...)
 
     assign("out_predicted", out, inherits=TRUE)
 
@@ -123,8 +128,8 @@ mkinfit <- function(mkinmod, observed,
       if(plot) {
         outtimes_plot = seq(min(observed$time), max(observed$time), length.out=100)
 
-        out_plot <- predict(mkinmod, transparms, odeini, outtimes_plot, 
-          solution_type = solution)
+        out_plot <- mkinpredict(mkinmod, parms, odeini, outtimes_plot, 
+          solution_type = solution, ...)
 
         plot(0, type="n", 
           xlim = range(observed$time), ylim = range(observed$value, na.rm=TRUE),
@@ -160,8 +165,10 @@ mkinfit <- function(mkinmod, observed,
   fit$predicted <- out_predicted
 
   # Collect initial parameter values in two dataframes
-  fit$start <- data.frame(initial = c(state.ini.optim, parms.optim))
+  fit$start <- data.frame(initial = c(state.ini.optim, 
+		  backtransform_odeparms(parms.optim, mod_vars)))
   fit$start$type = c(rep("state", length(state.ini.optim)), rep("deparm", length(parms.optim)))
+  fit$start$transformed = c(state.ini.optim, parms.optim)
 
   fit$fixed <- data.frame(
     value = c(state.ini.fixed, parms.fixed))
@@ -194,7 +201,7 @@ mkinfit <- function(mkinmod, observed,
   fit$errmin <- errmin
 
   # Calculate dissipation times DT50 and DT90 from parameters
-  parms.all = transform_odeparms(c(fit$par, parms.fixed), mod_vars)
+  parms.all = backtransform_odeparms(c(fit$par, parms.fixed), mod_vars)
   fit$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), DT90 = rep(NA, length(obs_vars)), 
     row.names = obs_vars)
   fit$SFORB <- vector()
@@ -285,16 +292,14 @@ mkinfit <- function(mkinmod, observed,
   data$variable <- ordered(data$variable, levels = obs_vars)
   fit$data <- data[order(data$variable, data$time), ]
   fit$atol <- atol
+  fit$parms.all <- parms.all
 
   class(fit) <- c("mkinfit", "modFit") 
   return(fit)
 }
 
-summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, ...) {
+summary.mkinfit <- function(object, data = TRUE, distimes = 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
@@ -304,7 +309,7 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, ...
     covar <- matrix(data = NA, nrow = p, ncol = p)
   } else message <- "ok"
 
-  rownames(covar) <- colnames(covar) <-pnames
+  rownames(covar) <- colnames(covar) <- pnames
   rdf    <- object$df.residual
   resvar <- object$ssr / rdf
   se     <- sqrt(diag(covar) * resvar)
@@ -312,12 +317,6 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, ...
   tval      <- param / se
   modVariance <- object$ssr / length(object$residuals)
 
-  if (!all(object$start$lower >=0)) {
-    message <- "Note that the one-sided t-test may not be appropriate if
-      parameter values below zero are possible."
-    warning(message)
-  } else message <- "ok"
-
   param <- cbind(param, se)
   dimnames(param) <- list(pnames, c("Estimate", "Std. Error"))
 
@@ -335,20 +334,17 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, ...
   if(data) ans$data <- object$data
   ans$start <- object$start
 
-  # Fix names of parameters to show that they were transformed
-  names(ans$start) <- sub("^k_", "log k_", names(ans$start))
-
   ans$fixed <- object$fixed
   ans$errmin <- object$errmin 
+  ans$parms.all <- object$parms.all
   if(distimes) ans$distimes <- object$distimes
-  if(ff) ans$ff <- object$ff
   if(length(object$SFORB) != 0) ans$SFORB <- object$SFORB
   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), tval = TRUE, ...) {
+print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), ...) {
   cat("\nEquations:\n")
   print(noquote(as.character(x[["diffs"]])))
   df  <- x$df
@@ -361,12 +357,11 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), t
   if(length(x$fixed$value) == 0) cat("None\n")
   else print(x$fixed)
   
-  cat("\nOptimised parameters:\n")
-  if (tval) printCoefmat(x$par, digits = digits, ...)
-  else {
-    printCoefmat(x$par[,c(1,2,4)], cs.in = c(1,2), tst.ind = integer(0), 
-      P.values = TRUE, has.Pvalue = TRUE, digits = digits, ...)
-  }
+  cat("\nOptimised, transformed parameters:\n")
+  printCoefmat(x$par, digits = digits, ...)
+
+  cat("\nBacktransformed parameters:\n")
+  print(as.data.frame(list(Estimate = x$parms.all)))
 
   cat("\nResidual standard error:",
       format(signif(x$sigma, digits)), "on", rdf, "degrees of freedom\n")
@@ -381,12 +376,6 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), t
     print(x$distimes, digits=digits,...)
   }    
 
-  printff <- !is.null(x$ff)
-  if(printff){
-    cat("\nEstimated formation fractions:\n")
-    print(data.frame(ff = x$ff), digits=digits,...)
-  }    
-
   printSFORB <- !is.null(x$SFORB)
   if(printSFORB){
     cat("\nEstimated Eigenvalues of SFORB model(s):\n")
diff --git a/R/mkinpredict.R b/R/mkinpredict.R
new file mode 100644
index 0000000..49166f3
--- /dev/null
+++ b/R/mkinpredict.R
@@ -0,0 +1,90 @@
+mkinpredict <- function(mkinmod, odeparms, odeini, outtimes, solution_type = "deSolve", map_output = TRUE, atol = 1e-6, ...) {
+
+  # Get the names of the state variables in the model
+  mod_vars <- names(mkinmod$diffs)
+
+  # Create function for evaluation of expressions with ode parameters and initial values
+  evalparse <- function(string)
+  {
+    eval(parse(text=string), as.list(c(odeparms, odeini)))
+  }
+
+  # Create a function calculating the differentials specified by the model
+  # if necessary
+  if (solution_type == "analytical") {
+    parent.type = names(mkinmod$map[[1]])[1]  
+    parent.name = names(mkinmod$diffs)[[1]]
+    o <- switch(parent.type,
+      SFO = SFO.solution(outtimes, 
+          evalparse(parent.name),
+          evalparse(paste("k", parent.name, sep="_"))),
+      FOMC = FOMC.solution(outtimes,
+          evalparse(parent.name),
+          evalparse("alpha"), evalparse("beta")),
+      DFOP = DFOP.solution(outtimes,
+          evalparse(parent.name),
+          evalparse("k1"), evalparse("k2"),
+          evalparse("g")),
+      HS = HS.solution(outtimes,
+          evalparse(parent.name),
+          evalparse("k1"), evalparse("k2"),
+          evalparse("tb")),
+      SFORB = SFORB.solution(outtimes,
+          evalparse(parent.name),
+          evalparse(paste("k", parent.name, "bound", sep="_")),
+          evalparse(paste("k", sub("free", "bound", parent.name), "free", sep="_")),
+          evalparse(paste("k", parent.name, sep="_")))
+    )
+    out <- cbind(outtimes, o)
+    dimnames(out) <- list(outtimes, c("time", sub("_free", "", parent.name)))
+  }
+  if (solution_type == "eigen") {
+    coefmat.num <- matrix(sapply(as.vector(mkinmod$coefmat), evalparse), 
+      nrow = length(mod_vars))
+    e <- eigen(coefmat.num)
+    c <- solve(e$vectors, odeini)
+    f.out <- function(t) {
+      e$vectors %*% diag(exp(e$values * t), nrow=length(mod_vars)) %*% c
+    }
+    o <- matrix(mapply(f.out, outtimes), 
+      nrow = length(mod_vars), ncol = length(outtimes))
+    dimnames(o) <- list(mod_vars, outtimes)
+    out <- cbind(time = outtimes, t(o))
+  } 
+  if (solution_type == "deSolve") {
+    mkindiff <- function(t, state, parms) {
+
+      time <- t
+      diffs <- vector()
+      for (box in names(mkinmod$diffs))
+      {
+        diffname <- paste("d", box, sep="_")      
+        diffs[diffname] <- with(as.list(c(time, state, parms)),
+          eval(parse(text=mkinmod$diffs[[box]])))
+      }
+      return(list(c(diffs)))
+    } 
+    out <- ode(
+      y = odeini,
+      times = outtimes,
+      func = mkindiff, 
+      parms = odeparms,
+      atol = atol,
+      ...
+    )
+  }
+  if (map_output) {
+    # Output transformation for models with unobserved compartments like SFORB
+    out_mapped <- data.frame(time = out[,"time"])
+    for (var in names(mkinmod$map)) {
+      if((length(mkinmod$map[[var]]) == 1) || solution_type == "analytical") {
+        out_mapped[var] <- out[, var]
+      } else {
+        out_mapped[var] <- rowSums(out[, mkinmod$map[[var]]])
+      }
+    }
+    return(out_mapped) 
+  } else {
+    return(out)
+  }
+}
diff --git a/R/predict.mkinmod.R b/R/predict.mkinmod.R
deleted file mode 100644
index 279fa52..0000000
--- a/R/predict.mkinmod.R
+++ /dev/null
@@ -1,89 +0,0 @@
-predict.mkinmod <- function(mkinmod, odeparms, odeini, outtimes, solution_type = "deSolve", map_output = TRUE, atol = 1e-6) {
-
-  # Get the names of the state variables in the model
-  mod_vars <- names(mkinmod$diffs)
-
-  # Create function for evaluation of expressions with ode parameters and initial values
-  evalparse <- function(string)
-  {
-    eval(parse(text=string), as.list(c(odeparms, odeini)))
-  }
-
-  # Create a function calculating the differentials specified by the model
-  # if necessary
-  if (solution_type == "analytical") {
-    parent.type = names(mkinmod$map[[1]])[1]  
-    parent.name = names(mkinmod$diffs)[[1]]
-    o <- switch(parent.type,
-      SFO = SFO.solution(outtimes, 
-          evalparse(parent.name),
-          evalparse(paste("k", parent.name, sep="_"))),
-      FOMC = FOMC.solution(outtimes,
-          evalparse(parent.name),
-          evalparse("alpha"), evalparse("beta")),
-      DFOP = DFOP.solution(outtimes,
-          evalparse(parent.name),
-          evalparse("k1"), evalparse("k2"),
-          evalparse("g")),
-      HS = HS.solution(outtimes,
-          evalparse(parent.name),
-          evalparse("k1"), evalparse("k2"),
-          evalparse("tb")),
-      SFORB = SFORB.solution(outtimes,
-          evalparse(parent.name),
-          evalparse(paste("k", parent.name, "bound", sep="_")),
-          evalparse(paste("k", sub("free", "bound", parent.name), "free", sep="_")),
-          evalparse(paste("k", parent.name, sep="_")))
-    )
-    out <- cbind(outtimes, o)
-    dimnames(out) <- list(outtimes, c("time", sub("_free", "", parent.name)))
-  }
-  if (solution_type == "eigen") {
-    coefmat.num <- matrix(sapply(as.vector(mkinmod$coefmat), evalparse), 
-      nrow = length(mod_vars))
-    e <- eigen(coefmat.num)
-    c <- solve(e$vectors, odeini)
-    f.out <- function(t) {
-      e$vectors %*% diag(exp(e$values * t), nrow=length(mod_vars)) %*% c
-    }
-    o <- matrix(mapply(f.out, outtimes), 
-      nrow = length(mod_vars), ncol = length(outtimes))
-    dimnames(o) <- list(mod_vars, outtimes)
-    out <- cbind(time = outtimes, t(o))
-  } 
-  if (solution_type == "deSolve") {
-    mkindiff <- function(t, state, parms) {
-
-      time <- t
-      diffs <- vector()
-      for (box in names(mkinmod$diffs))
-      {
-        diffname <- paste("d", box, sep="_")      
-        diffs[diffname] <- with(as.list(c(time, state, parms)),
-          eval(parse(text=mkinmod$diffs[[box]])))
-      }
-      return(list(c(diffs)))
-    } 
-    out <- ode(
-      y = odeini,
-      times = outtimes,
-      func = mkindiff, 
-      parms = odeparms,
-      atol = atol
-    )
-  }
-  if (map_output) {
-    # Output transformation for models with unobserved compartments like SFORB
-    out_mapped <- data.frame(time = out[,"time"])
-    for (var in names(mkinmod$map)) {
-      if((length(mkinmod$map[[var]]) == 1) || solution == "analytical") {
-        out_mapped[var] <- out[, var]
-      } else {
-        out_mapped[var] <- rowSums(out[, mkinmod$map[[var]]])
-      }
-    }
-    return(out_mapped) 
-  } else {
-    return(out)
-  }
-}
diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R
index 4b8bfd1..8b39874 100644
--- a/R/transform_odeparms.R
+++ b/R/transform_odeparms.R
@@ -1,46 +1,72 @@
-transform_odeparms <- function(odeparms, mod_vars) {
-  # Transform rate constants and formation fractions
-  transparms <- odeparms
-  # Exponential transformation for rate constants
-  index_k <- grep("^k_", names(odeparms))
+transform_odeparms <- function(parms, mod_vars) {
+  # Set up container for transformed parameters
+  transparms <- parms
+
+  # Log transformation for rate constants
+  index_k <- grep("^k_", names(transparms))
   if (length(index_k) > 0) {
-    transparms[index_k] <- exp(odeparms[index_k])
+    transparms[index_k] <- log(parms[index_k])
   }
 
-  # Go through state variables and apply inverse isotropic logratio transformation
+  # Go through state variables and apply isotropic logratio transformation
   for (box in mod_vars) {
-    indices_f <- grep(paste("^f", box, sep = "_"), names(odeparms))
-    f_names <- grep(paste("^f", box, sep = "_"), names(odeparms), value = TRUE)
+    indices_f <- grep(paste("^f", box, sep = "_"), names(parms))
+    f_names <- grep(paste("^f", box, sep = "_"), names(parms), value = TRUE)
     n_paths <- length(indices_f)
     if (n_paths > 0) {
-      f <- invilr(odeparms[indices_f])[1:n_paths] # We do not need the last component
-      names(f) <- f_names
-      transparms[indices_f] <- f
+      f <- parms[indices_f]
+      trans_f <- ilr(c(f, 1 - sum(f)))
+      names(trans_f) <- f_names
+      transparms[indices_f] <- trans_f
     }
   }
+
+  # Transform parameters also for FOMC, DFOP and HS models
+  for (pname in c("alpha", "beta", "k1", "k2", "tb")) {
+    if (!is.na(parms[pname])) {
+      transparms[pname] <- log(parms[pname])
+    } 
+  }
+  if (!is.na(parms["g"])) {
+    g <- parms["g"]
+    transparms["g"] <- ilr(c(g, 1- g))
+  }
+
   return(transparms)
 }
 
 backtransform_odeparms <- function(transparms, mod_vars) {
-  # Transform rate constants and formation fractions
-  odeparms <- transparms
-  # Log transformation for rate constants
-  index_k <- grep("^k_", names(transparms))
+  # Set up container for backtransformed parameters
+  parms <- transparms
+
+  # Exponential transformation for rate constants
+  index_k <- grep("^k_", names(parms))
   if (length(index_k) > 0) {
-    odeparms[index_k] <- log(transparms[index_k])
+    parms[index_k] <- exp(transparms[index_k])
   }
 
-  # Go through state variables and apply isotropic logratio transformation
+  # Go through state variables and apply inverse isotropic logratio transformation
   for (box in mod_vars) {
     indices_f <- grep(paste("^f", box, sep = "_"), names(transparms))
     f_names <- grep(paste("^f", box, sep = "_"), names(transparms), value = TRUE)
     n_paths <- length(indices_f)
     if (n_paths > 0) {
-      trans_f <- transparms[indices_f]
-      f <- ilr(c(trans_f, 1 - sum(trans_f)))
+      f <- invilr(transparms[indices_f])[1:n_paths] # We do not need the last component
       names(f) <- f_names
-      odeparms[indices_f] <- f
+      parms[indices_f] <- f
     }
   }
-  return(odeparms)
+
+  # Transform parameters also for DFOP and HS models
+  for (pname in c("alpha", "beta", "k1", "k2", "tb")) {
+    if (!is.na(transparms[pname])) {
+      parms[pname] <- exp(transparms[pname])
+    } 
+  }
+  if (!is.na(transparms["g"])) {
+    g <- transparms["g"]
+    parms["g"] <- invilr(g)[1]
+  }
+
+  return(parms)
 }
-- 
cgit v1.2.1