diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/add_err.R | 2 | ||||
| -rw-r--r-- | R/mkin_wide_to_long.R | 1 | ||||
| -rw-r--r-- | R/mkinfit.R | 8 | ||||
| -rw-r--r-- | R/mkinpredict.R | 48 | ||||
| -rw-r--r-- | R/plot.mkinfit.R | 1 |
5 files changed, 32 insertions, 28 deletions
diff --git a/R/add_err.R b/R/add_err.R index a523e9c2..9235223f 100644 --- a/R/add_err.R +++ b/R/add_err.R @@ -76,6 +76,8 @@ add_err <- function(prediction, sdfunc, secondary = c("M1", "M2"), { if (!is.na(seed)) set.seed(seed) + prediction <- as.data.frame(prediction) + # The output of mkinpredict is in wide format d_long = mkin_wide_to_long(prediction, time = "time") diff --git a/R/mkin_wide_to_long.R b/R/mkin_wide_to_long.R index bef0e408..971f5273 100644 --- a/R/mkin_wide_to_long.R +++ b/R/mkin_wide_to_long.R @@ -21,6 +21,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value")) #' @export mkin_wide_to_long <- function(wide_data, time = "t") { + wide_data <- as.data.frame(wide_data) colnames <- names(wide_data) if (!(time %in% colnames)) stop("The data in wide format have to contain a variable named ", time, ".") vars <- subset(colnames, colnames != time) diff --git a/R/mkinfit.R b/R/mkinfit.R index e1089673..f6691b1b 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -592,12 +592,8 @@ mkinfit <- function(mkinmod, observed, method.ode = method.ode, atol = atol, rtol = rtol, ...) - out_long <- mkin_wide_to_long(out, time = "time") - - out_merged <- merge(observed[c("name", "time")], out_long) - out_merged$name <- ordered(out_merged$name, levels = obs_vars) - out_merged <- out_merged[order(out_merged$name, out_merged$time), ] - observed$predicted <- out_merged$value + observed_index <- cbind(as.character(observed$time), as.character(observed$name)) + observed$predicted <- out[observed_index] } # Define standard deviation for each observation diff --git a/R/mkinpredict.R b/R/mkinpredict.R index 75582fac..90cd45fb 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -32,12 +32,13 @@ #' is 1e-10, much lower than in \code{\link{lsoda}}. #' @param map_output Boolean to specify if the output should list values for #' the observed variables (default) or for all state variables (if set to -#' FALSE). +#' FALSE). Setting this to FALSE has no effect for analytical solutions, +#' as these always return mapped output. #' @param \dots Further arguments passed to the ode solver in case such a #' solver is used. #' @import deSolve #' @importFrom inline getDynLib -#' @return A data frame with the solution in wide format +#' @return A matrix with the numeric solution in wide format #' @author Johannes Ranke #' @examples #' @@ -70,7 +71,7 @@ #' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), #' seq(0, 20, by = 0.01))[2001,] #' -#' # Check compiled model versions - they are faster than the eigenvalue based solutions! +#' # Comparison of the performance of solution types #' SFO_SFO = mkinmod(parent = list(type = "SFO", to = "m1"), #' m1 = list(type = "SFO"), use_of_ff = "max") #' if(require(rbenchmark)) { @@ -92,15 +93,11 @@ #' c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), #' solution_type = "analytical", use_compiled = FALSE)[201,]) #' } -#' analytical = mkinpredict(SFO_SFO, -#' c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), -#' c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), -#' solution_type = "analytical", use_compiled = FALSE)[201,] #' #' \dontrun{ #' # Predict from a fitted model #' f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE) -#' f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE, solution_type = "analytical") +#' f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE, solution_type = "deSolve") #' head(mkinpredict(f)) #' } #' @@ -127,29 +124,34 @@ mkinpredict.mkinmod <- function(x, map_output = TRUE, ...) { - # Get the names of the state variables in the model + # Names of state variables and observed variables mod_vars <- names(x$diffs) + obs_vars <- names(x$spec) # Order the inital values for state variables if they are named if (!is.null(names(odeini))) { odeini <- odeini[mod_vars] } + out_obs <- matrix(NA, nrow = length(outtimes), ncol = 1 + length(obs_vars), + dimnames = list(as.character(outtimes), c("time", obs_vars))) + out_obs[, "time"] <- outtimes + if (solution_type == "analytical") { # This is clumsy, as we wanted fast analytical predictions for mkinfit - out <- data.frame(time = outtimes) - obs_vars <- names(x$spec) pseudo_observed <- data.frame(name = rep(obs_vars, each = length(outtimes)), time = rep(outtimes, length(obs_vars))) pseudo_observed$predicted <- x$deg_func(pseudo_observed, odeini, odeparms) for (obs_var in obs_vars) { - out[obs_var] <- pseudo_observed[pseudo_observed$name == obs_var, "predicted"] + out_obs[, obs_var] <- pseudo_observed[pseudo_observed$name == obs_var, "predicted"] } + # We don't have solutions for unobserved state variables, the output of + # analytical solutions is always mapped to observed variables + return(out_obs) } if (solution_type == "eigen") { - evalparse <- function(string) { eval(parse(text=string), as.list(c(odeparms, odeini))) } @@ -163,8 +165,8 @@ mkinpredict.mkinmod <- function(x, } o <- matrix(mapply(f.out, outtimes), nrow = length(mod_vars), ncol = length(outtimes)) - out <- data.frame(outtimes, t(o)) - names(out) <- c("time", mod_vars) + out <- cbind(outtimes, t(o)) + colnames(out) <- c("time", mod_vars) } if (solution_type == "deSolve") { @@ -208,21 +210,23 @@ mkinpredict.mkinmod <- function(x, if (sum(is.na(out)) > 0) { stop("Differential equations were not integrated for all output times because\n", "NaN values occurred in output from ode()") - } + } } + if (map_output) { # Output transformation for models with unobserved compartments like SFORB - out_mapped <- data.frame(time = out[,"time"]) + # if not already mapped in analytical solution for (var in names(x$map)) { - if((length(x$map[[var]]) == 1) || solution_type == "analytical") { - out_mapped[var] <- out[, var] + if((length(x$map[[var]]) == 1)) { + out_obs[, var] <- out[, var] } else { - out_mapped[var] <- rowSums(out[, x$map[[var]]]) + out_obs[, var] <- out[, x$map[[var]][1]] + out[, x$map[[var]][2]] } } - return(out_mapped) + return(out_obs) } else { - return(as.data.frame(out)) + dimnames(out) <- list(time = as.character(outtimes), c("time", mod_vars)) + return(out) } } diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R index b2cb4890..48df9483 100644 --- a/R/plot.mkinfit.R +++ b/R/plot.mkinfit.R @@ -148,6 +148,7 @@ plot.mkinfit <- function(x, fit = x, solution_type = solution_type, atol = fit$atol, rtol = fit$rtol, use_compiled = FALSE) } + out <- as.data.frame(out) names(col_obs) <- names(pch_obs) <- names(lty_obs) <- obs_vars |