jranke / mkin
 1 ```#' Produce predictions from a kinetic model using specific parameters ``` 2 ```#' ``` 3 ```#' This function produces a time series for all the observed variables in a ``` 4 ```#' kinetic model as specified by \code{\link{mkinmod}}, using a specific set of ``` 5 ```#' kinetic parameters and initial values for the state variables. ``` 6 ```#' ``` 7 ```#' @aliases mkinpredict mkinpredict.mkinmod mkinpredict.mkinfit ``` 8 ```#' @param x A kinetic model as produced by \code{\link{mkinmod}}, or a kinetic ``` 9 ```#' fit as fitted by \code{\link{mkinfit}}. In the latter case, the fitted ``` 10 ```#' parameters are used for the prediction. ``` 11 ```#' @param odeparms A numeric vector specifying the parameters used in the ``` 12 ```#' kinetic model, which is generally defined as a set of ordinary ``` 13 ```#' differential equations. ``` 14 ```#' @param odeini A numeric vector containing the initial values of the state ``` 15 ```#' variables of the model. Note that the state variables can differ from the ``` 16 ```#' observed variables, for example in the case of the SFORB model. ``` 17 ```#' @param outtimes A numeric vector specifying the time points for which model ``` 18 ```#' predictions should be generated. ``` 19 ```#' @param solution_type The method that should be used for producing the ``` 20 ```#' predictions. This should generally be "analytical" if there is only one ``` 21 ```#' observed variable, and usually "deSolve" in the case of several observed ``` 22 ```#' variables. The third possibility "eigen" is faster but not applicable to ``` 23 ```#' some models e.g. using FOMC for the parent compound. ``` 24 ```#' @param method.ode The solution method passed via \code{\link{mkinpredict}} ``` 25 ```#' to \code{\link{ode}} in case the solution type is "deSolve". The default ``` 26 ```#' "lsoda" is performant, but sometimes fails to converge. ``` 27 ```#' @param use_compiled If set to \code{FALSE}, no compiled version of the ``` 28 ```#' \code{\link{mkinmod}} model is used, even if is present. ``` 29 ```#' @param atol Absolute error tolerance, passed to \code{\link{ode}}. Default ``` 30 ```#' is 1e-8, lower than in \code{\link{lsoda}}. ``` 31 ```#' @param rtol Absolute error tolerance, passed to \code{\link{ode}}. Default ``` 32 ```#' is 1e-10, much lower than in \code{\link{lsoda}}. ``` 33 ```#' @param map_output Boolean to specify if the output should list values for ``` 34 ```#' the observed variables (default) or for all state variables (if set to ``` 35 ```#' FALSE). Setting this to FALSE has no effect for analytical solutions, ``` 36 ```#' as these always return mapped output. ``` 37 ```#' @param na_stop Should it be an error if deSolve::ode returns NaN values ``` 38 ```#' @param \dots Further arguments passed to the ode solver in case such a ``` 39 ```#' solver is used. ``` 40 ```#' @import deSolve ``` 41 ```#' @return A matrix with the numeric solution in wide format ``` 42 ```#' @author Johannes Ranke ``` 43 ```#' @examples ``` 44 ```#' ``` 45 ```#' SFO <- mkinmod(degradinol = mkinsub("SFO")) ``` 46 ```#' # Compare solution types ``` 47 ```#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, ``` 48 ```#' solution_type = "analytical") ``` 49 ```#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, ``` 50 ```#' solution_type = "deSolve") ``` 51 ```#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, ``` 52 ```#' solution_type = "deSolve", use_compiled = FALSE) ``` 53 ```#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, ``` 54 ```#' solution_type = "eigen") ``` 55 ```#' ``` 56 ```#' # Compare integration methods to analytical solution ``` 57 ```#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, ``` 58 ```#' solution_type = "analytical")[21,] ``` 59 ```#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, ``` 60 ```#' method = "lsoda")[21,] ``` 61 ```#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, ``` 62 ```#' method = "ode45")[21,] ``` 63 ```#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20, ``` 64 ```#' method = "rk4")[21,] ``` 65 ```#' # rk4 is not as precise here ``` 66 ```#' ``` 67 ```#' # The number of output times used to make a lot of difference until the ``` 68 ```#' # default for atol was adjusted ``` 69 ```#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), ``` 70 ```#' seq(0, 20, by = 0.1))[201,] ``` 71 ```#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), ``` 72 ```#' seq(0, 20, by = 0.01))[2001,] ``` 73 ```#' ``` 74 ```#' # Comparison of the performance of solution types ``` 75 ```#' SFO_SFO = mkinmod(parent = list(type = "SFO", to = "m1"), ``` 76 ```#' m1 = list(type = "SFO"), use_of_ff = "max") ``` 77 ```#' if(require(rbenchmark)) { ``` 78 ```#' benchmark(replications = 10, order = "relative", columns = c("test", "relative", "elapsed"), ``` 79 ```#' eigen = mkinpredict(SFO_SFO, ``` 80 ```#' c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), ``` 81 ```#' c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), ``` 82 ```#' solution_type = "eigen")[201,], ``` 83 ```#' deSolve_compiled = mkinpredict(SFO_SFO, ``` 84 ```#' c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), ``` 85 ```#' c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), ``` 86 ```#' solution_type = "deSolve")[201,], ``` 87 ```#' deSolve = mkinpredict(SFO_SFO, ``` 88 ```#' c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), ``` 89 ```#' c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), ``` 90 ```#' solution_type = "deSolve", use_compiled = FALSE)[201,], ``` 91 ```#' analytical = mkinpredict(SFO_SFO, ``` 92 ```#' c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01), ``` 93 ```#' c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), ``` 94 ```#' solution_type = "analytical", use_compiled = FALSE)[201,]) ``` 95 ```#' } ``` 96 ```#' ``` 97 ```#' \dontrun{ ``` 98 ```#' # Predict from a fitted model ``` 99 ```#' f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE) ``` 100 ```#' f <- mkinfit(SFO_SFO, FOCUS_2006_C, quiet = TRUE, solution_type = "deSolve") ``` 101 ```#' head(mkinpredict(f)) ``` 102 ```#' } ``` 103 ```#' ``` 104 ```#' @export ``` 105 ```mkinpredict <- function(x, odeparms, odeini, outtimes, ...) ``` 106 ```{ ``` 107 1 ``` UseMethod("mkinpredict", x) ``` 108 ```} ``` 109 110 ```#' @rdname mkinpredict ``` 111 ```#' @export ``` 112 ```mkinpredict.mkinmod <- function(x, ``` 113 ``` odeparms = c(k_parent_sink = 0.1), ``` 114 ``` odeini = c(parent = 100), ``` 115 ``` outtimes = seq(0, 120, by = 0.1), ``` 116 ``` solution_type = "deSolve", ``` 117 ``` use_compiled = "auto", ``` 118 ``` method.ode = "lsoda", atol = 1e-8, rtol = 1e-10, ``` 119 ``` map_output = TRUE, ``` 120 ``` na_stop = TRUE, ``` 121 ``` ...) ``` 122 ```{ ``` 123 124 ``` # Names of state variables and observed variables ``` 125 1 ``` mod_vars <- names(x\$diffs) ``` 126 1 ``` obs_vars <- names(x\$spec) ``` 127 128 ``` # Order the inital values for state variables if they are named ``` 129 1 ``` if (!is.null(names(odeini))) { ``` 130 1 ``` odeini <- odeini[mod_vars] ``` 131 ``` } ``` 132 133 1 ``` out_obs <- matrix(NA, nrow = length(outtimes), ncol = 1 + length(obs_vars), ``` 134 1 ``` dimnames = list(as.character(outtimes), c("time", obs_vars))) ``` 135 1 ``` out_obs[, "time"] <- outtimes ``` 136 137 1 ``` n_out_na <- 0 # to check if we get NA values with deSolve ``` 138 139 1 ``` if (solution_type == "analytical") { ``` 140 ``` # This is clumsy, as we wanted fast analytical predictions for mkinfit, ``` 141 ``` # which bypasses mkinpredict in the case of analytical solutions ``` 142 1 ``` pseudo_observed <- ``` 143 1 ``` data.frame(name = rep(obs_vars, each = length(outtimes)), ``` 144 1 ``` time = rep(outtimes, length(obs_vars))) ``` 145 1 ``` pseudo_observed\$predicted <- x\$deg_func(pseudo_observed, odeini, odeparms) ``` 146 1 ``` for (obs_var in obs_vars) { ``` 147 1 ``` out_obs[, obs_var] <- pseudo_observed[pseudo_observed\$name == obs_var, "predicted"] ``` 148 ``` } ``` 149 ``` # We don't have solutions for unobserved state variables, the output of ``` 150 ``` # analytical solutions is always mapped to observed variables ``` 151 1 ``` return(out_obs) ``` 152 ``` } ``` 153 154 1 ``` if (solution_type == "eigen") { ``` 155 1 ``` evalparse <- function(string) { ``` 156 1 ``` eval(parse(text=string), as.list(c(odeparms, odeini))) ``` 157 ``` } ``` 158 159 1 ``` coefmat.num <- matrix(sapply(as.vector(x\$coefmat), evalparse), ``` 160 1 ``` nrow = length(mod_vars)) ``` 161 1 ``` e <- eigen(coefmat.num) ``` 162 1 ``` c <- solve(e\$vectors, odeini) ``` 163 1 ``` f.out <- function(t) { ``` 164 1 ``` e\$vectors %*% diag(exp(e\$values * t), nrow=length(mod_vars)) %*% c ``` 165 ``` } ``` 166 1 ``` o <- matrix(mapply(f.out, outtimes), ``` 167 1 ``` nrow = length(mod_vars), ncol = length(outtimes)) ``` 168 1 ``` out <- cbind(outtimes, t(o)) ``` 169 1 ``` colnames(out) <- c("time", mod_vars) ``` 170 ``` } ``` 171 172 1 ``` if (solution_type == "deSolve") { ``` 173 1 ``` if (!is.null(x\$cf) & use_compiled[1] != FALSE) { ``` 174 175 1 ``` out <- deSolve::ode( ``` 176 1 ``` y = odeini, ``` 177 1 ``` times = outtimes, ``` 178 1 ``` func = "diffs", ``` 179 1 ``` initfunc = "initpar", ``` 180 1 ``` dllname = inline::getDynLib(x\$cf)[["name"]], ``` 181 1 ``` parms = odeparms[x\$parms], # Order matters when using compiled models ``` 182 1 ``` method = method.ode, ``` 183 1 ``` atol = atol, ``` 184 1 ``` rtol = rtol, ``` 185 ``` ... ``` 186 ``` ) ``` 187 ``` } else { ``` 188 1 ``` mkindiff <- function(t, state, parms) { ``` 189 190 1 ``` time <- t ``` 191 1 ``` diffs <- vector() ``` 192 1 ``` for (box in names(x\$diffs)) ``` 193 ``` { ``` 194 1 ``` diffname <- paste("d", box, sep="_") ``` 195 1 ``` diffs[diffname] <- with(as.list(c(time, state, parms)), ``` 196 1 ``` eval(parse(text=x\$diffs[[box]]))) ``` 197 ``` } ``` 198 1 ``` return(list(c(diffs))) ``` 199 ``` } ``` 200 1 ``` out <- deSolve::ode( ``` 201 1 ``` y = odeini, ``` 202 1 ``` times = outtimes, ``` 203 1 ``` func = mkindiff, ``` 204 1 ``` parms = odeparms, ``` 205 1 ``` method = method.ode, ``` 206 1 ``` atol = atol, ``` 207 1 ``` rtol = rtol, ``` 208 ``` ... ``` 209 ``` ) ``` 210 ``` } ``` 211 1 ``` n_out_na <- sum(is.na(out)) ``` 212 1 ``` if (n_out_na > 0 & na_stop) { ``` 213 0 ``` cat("odeini:\n") ``` 214 0 ``` print(odeini) ``` 215 0 ``` cat("odeparms:\n") ``` 216 0 ``` print(odeparms) ``` 217 0 ``` cat("out:\n") ``` 218 0 ``` print(out) ``` 219 0 ``` stop("Differential equations were not integrated for all output times because\n", ``` 220 0 ``` n_out_na, " NaN values occurred in output from ode()") ``` 221 ``` } ``` 222 ``` } ``` 223 224 1 ``` if (map_output) { ``` 225 ``` # Output transformation for models with unobserved compartments like SFORB ``` 226 ``` # if not already mapped in analytical solution ``` 227 1 ``` if (n_out_na > 0 & !na_stop) { ``` 228 0 ``` available <- c(TRUE, rep(FALSE, length(outtimes) - 1)) ``` 229 ``` } else { ``` 230 1 ``` available <- rep(TRUE, length(outtimes)) ``` 231 ``` } ``` 232 1 ``` for (var in names(x\$map)) { ``` 233 1 ``` if((length(x\$map[[var]]) == 1)) { ``` 234 1 ``` out_obs[available, var] <- out[available, var] ``` 235 ``` } else { ``` 236 1 ``` out_obs[available, var] <- out[available, x\$map[[var]][1]] + ``` 237 1 ``` out[available, x\$map[[var]][2]] ``` 238 ``` } ``` 239 ``` } ``` 240 1 ``` return(out_obs) ``` 241 ``` } else { ``` 242 0 ``` dimnames(out) <- list(time = as.character(outtimes), c("time", mod_vars)) ``` 243 0 ``` return(out) ``` 244 ``` } ``` 245 ```} ``` 246 247 ```#' @rdname mkinpredict ``` 248 ```#' @export ``` 249 ```mkinpredict.mkinfit <- function(x, ``` 250 ``` odeparms = x\$bparms.ode, ``` 251 ``` odeini = x\$bparms.state, ``` 252 ``` outtimes = seq(0, 120, by = 0.1), ``` 253 ``` solution_type = "deSolve", ``` 254 ``` use_compiled = "auto", ``` 255 ``` method.ode = "lsoda", atol = 1e-8, rtol = 1e-10, ``` 256 ``` map_output = TRUE, ...) ``` 257 ```{ ``` 258 0 ``` mkinpredict(x\$mkinmod, odeparms, odeini, outtimes, solution_type, use_compiled, ``` 259 0 ``` method.ode, atol, rtol, map_output, ...) ``` 260 ```} ```

Read our documentation on viewing source code .