From d3df16102c5ed4bf9182b4f1893561e99eaed166 Mon Sep 17 00:00:00 2001 From: jranke Date: Tue, 3 Apr 2012 05:17:54 +0000 Subject: - Separated model prediction out into a separate function - Created separate functions for parameter transformations (not documented yet) - Fitting of analytical models works - SFO_SFO model works, complex models do not at the moment - summary.mkinfit not operational at the moment git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@21 edb9625f-4e0d-4859-8d74-9fd3b1da38cb --- R/mkinfit.R | 229 ++++++++++--------------------------------------- R/predict.mkinmod.R | 89 +++++++++++++++++++ R/transform_odeparms.R | 46 ++++++++++ TODO | 1 + man/mkinfit.Rd | 6 +- 5 files changed, 183 insertions(+), 188 deletions(-) create mode 100644 R/predict.mkinmod.R create mode 100644 R/transform_odeparms.R diff --git a/R/mkinfit.R b/R/mkinfit.R index 83896bc..47c39cf 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -32,9 +32,12 @@ mkinfit <- function(mkinmod, observed, atol = 1e-6, ...) { + # Get the names of the state variables in the model mod_vars <- names(mkinmod$diffs) - # Subset dataframe of observed variables with mapped (modelled) variables + + # Subset observed data with names of observed data in the model observed <- subset(observed, name %in% names(mkinmod$map)) + # Get names of observed variables obs_vars = unique(as.character(observed$name)) @@ -42,26 +45,33 @@ mkinfit <- function(mkinmod, observed, if (parms.ini[[1]] == "auto") parms.ini = vector() defaultpar.names <- setdiff(mkinmod$parms, names(parms.ini)) for (parmname in defaultpar.names) { - # Default values for the FOMC model + # Default values for rate constants, depending on the parameterisation + if (substr(parmname, 1, 2) == "k_") parms.ini[parmname] = 0.1 + # Default values for formation fractions + if (substr(parmname, 1, 2) == "f_") parms.ini[parmname] = 0.1 + # Default values for the FOMC, DFOP and HS models if (parmname == "alpha") parms.ini[parmname] = 1 if (parmname == "beta") parms.ini[parmname] = 10 - # Default values for log of rate constants, depending on the parameterisation - if (substr(parmname, 1, 2) == "k_") parms.ini[parmname] = -2.3 - # Default values for ilr of formation fractions - if (substr(parmname, 1, 2) == "f_") parms.ini[parmname] = 0.1 + if (parmname == "k1") parms.ini[parmname] = 0.1 + if (parmname == "k2") parms.ini[parmname] = 0.01 + 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 parameter values if they are not named yet + # Name the inital state variable values if they are not named yet if(is.null(names(state.ini))) names(state.ini) <- mod_vars - # Parameters to be optimised + # Parameters to be optimised: + # Kinetic parameters in parms.ini whose names are not in fixed_parms parms.fixed <- parms.ini[fixed_parms] - optim_parms <- setdiff(names(parms.ini), fixed_parms) - parms.optim <- parms.ini[optim_parms] + parms.optim <- parms.ini[setdiff(names(parms.ini), fixed_parms)] + # Inital state variables in state.ini whose names are not in fixed_initials 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 <- state.ini[setdiff(names(state.ini), fixed_initials)] + + # Preserve names of state variables before renaming initial state variable parameters 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="_") @@ -77,42 +87,17 @@ mkinfit <- function(mkinmod, observed, else solution = "deSolve" } - # Create a function calculating the differentials specified by the model - # if necessary - if(solution == "deSolve") { - mkindiff <- function(t, state, parms) { - transparms <- vector() - transparms <- c(transparms, exp(parms[grep("^k_", names(parms))])) - for (box in mod_vars) { - path_indices <- grep(paste("^f", box, sep = "_"), names(parms)) - n_paths <- length(path_indices) - if (n_paths > 0) { - f <- invilr(parms[grep(paste("^f", box, sep="_"), names(parms))]) - transparms <- c(transparms, f[1:n_paths]) - } - } - otherparms <- parms[setdiff(names(parms), names(transparms))] - parms <- c(transparms, otherparms) - - time <- t - diffs <- vector() - for (box in mod_vars) - { - diffname <- paste("d", box, sep="_") - diffs[diffname] <- with(as.list(c(time, state, parms)), - eval(parse(text=mkinmod$diffs[[box]]))) - } - return(list(c(diffs))) - } - } - - cost.old <- 1e100 - calls <- 0 + cost.old <- 1e100 # The first model cost should be smaller than this value + calls <- 0 # Counter for number of model solutions out_predicted <- NA # Define the model cost function cost <- function(P) { - assign("calls", calls+1, inherits=TRUE) + assign("calls", calls+1, inherits=TRUE) # Increase the model solution counter + + # Time points at which observed data are available + outtimes = unique(observed$time) + 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)) @@ -120,76 +105,14 @@ mkinfit <- function(mkinmod, observed, odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed) - outtimes = unique(observed$time) - evalparse <- function(string) - { - eval(parse(text=string), as.list(c(odeparms, odeini))) - } + transparms <- transform_odeparms(odeparms, mod_vars) - # Solve the system - if (solution == "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, "sink", 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, "sink", sep="_"))) - ) - out <- cbind(outtimes, o) - dimnames(out) <- list(outtimes, c("time", sub("_free", "", parent.name))) - } - if (solution == "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 == "deSolve") - { - out <- ode( - y = odeini, - times = outtimes, - func = mkindiff, - parms = odeparms, - atol = atol - ) - } - - # 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) || solution == "analytical") { - out_transformed[var] <- out[, var] - } else { - out_transformed[var] <- rowSums(out[, mkinmod$map[[var]]]) - } - } - assign("out_predicted", out_transformed, inherits=TRUE) + # Solve the system with current transformed parameter values + out <- predict(mkinmod, transparms, odeini, outtimes, solution_type = solution) + + assign("out_predicted", out, inherits=TRUE) - mC <- modCost(out_transformed, observed, y = "value", + mC <- modCost(out, observed, y = "value", err = err, weight = weight, scaleVar = scaleVar) # Report and/or plot if the model is improved @@ -199,53 +122,9 @@ mkinfit <- function(mkinmod, observed, # Plot the data and current model output if requested if(plot) { outtimes_plot = seq(min(observed$time), max(observed$time), length.out=100) - if (solution == "analytical") { - o_plot <- switch(parent.type, - SFO = SFO.solution(outtimes_plot, - evalparse(parent.name), - evalparse(paste("k", parent.name, "sink", sep="_"))), - FOMC = FOMC.solution(outtimes_plot, - evalparse(parent.name), - evalparse("alpha"), evalparse("beta")), - DFOP = DFOP.solution(outtimes_plot, - evalparse(parent.name), - evalparse("k1"), evalparse("k2"), - evalparse("g")), - HS = HS.solution(outtimes_plot, - evalparse(parent.name), - evalparse("k1"), evalparse("k2"), - evalparse("tb")), - SFORB = SFORB.solution(outtimes_plot, - 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, "sink", sep="_"))) - ) - out_plot <- cbind(outtimes_plot, o_plot) - dimnames(out_plot) <- list(outtimes_plot, c("time", sub("_free", "", parent.name))) - } - if(solution == "eigen") { - o_plot <- matrix(mapply(f.out, outtimes_plot), - nrow = length(mod_vars), ncol = length(outtimes_plot)) - dimnames(o_plot) <- list(mod_vars, outtimes_plot) - out_plot <- cbind(time = outtimes_plot, t(o_plot)) - } - if (solution == "deSolve") { - out_plot <- ode( - y = odeini, - times = outtimes_plot, - func = mkindiff, - parms = odeparms) - } - out_transformed_plot <- data.frame(time = out_plot[,"time"]) - for (var in obs_vars) { - if((length(mkinmod$map[[var]]) == 1) || solution == "analytical") { - out_transformed_plot[var] <- out_plot[, var] - } else { - out_transformed_plot[var] <- rowSums(out_plot[, mkinmod$map[[var]]]) - } - } - out_transformed_plot <<- out_transformed_plot + + out_plot <- predict(mkinmod, transparms, odeini, outtimes_plot, + solution_type = solution) plot(0, type="n", xlim = range(observed$time), ylim = range(observed$value, na.rm=TRUE), @@ -256,7 +135,7 @@ mkinfit <- function(mkinmod, observed, points(subset(observed, name == obs_var, c(time, value)), pch = pch_obs[obs_var], col = col_obs[obs_var]) } - matlines(out_transformed_plot$time, out_transformed_plot[-1]) + matlines(out_plot$time, out_plot[-1]) legend("topright", inset=c(0.05, 0.05), legend=obs_vars, col=col_obs, pch=pch_obs, lty=1:length(pch_obs)) } @@ -272,12 +151,6 @@ mkinfit <- function(mkinmod, observed, if (solution == "eigen") { fit$coefmat <- mkinmod$coefmat } - if (solution == "deSolve") { - fit$mkindiff <- mkindiff - } - if (plot == TRUE) { - fit$out_transformed_plot = out_transformed_plot - } # We also need various other information for summary and plotting fit$map <- mkinmod$map @@ -289,8 +162,6 @@ mkinfit <- function(mkinmod, observed, # 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))) - fit$start$lower <- lower - fit$start$upper <- upper fit$fixed <- data.frame( value = c(state.ini.fixed, parms.fixed)) @@ -322,11 +193,10 @@ mkinfit <- function(mkinmod, observed, } fit$errmin <- errmin - # Calculate dissipation times DT50 and DT90 and formation fractions - parms.all = c(fit$par, parms.fixed) + # Calculate dissipation times DT50 and DT90 from parameters + parms.all = transform_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$ff <- vector() fit$SFORB <- vector() for (obs_var in obs_vars) { type = names(mkinmod$map[[obs_var]])[1] @@ -374,16 +244,6 @@ mkinfit <- function(mkinmod, observed, DT50 <- DTx(50) DT90 <- DTx(90) } - # Back-calculation of formation fractions in case of nonlinear parent kinetics - if (type %in% c("FOMC", "DFOP", "HS")) { - ff_names = names(mkinmod$ff) - for (ff_name in ff_names) - { - fit$ff[[paste(obs_var, ff_name, sep="_")]] = - eval(parse(text = mkinmod$ff[ff_name]), as.list(parms.all)) - } - fit$ff[[paste(obs_var, "sink", sep="_")]] = 1 - sum(fit$ff) - } if (type == "SFORB") { # FOCUS kinetics (2006), p. 60 f k_out_names = grep(paste("k", obs_var, "free", sep="_"), names(parms.all), value=TRUE) @@ -474,9 +334,10 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, ... ans$diffs <- object$diffs if(data) ans$data <- object$data ans$start <- object$start - if (object$logk) { - names(ans$start) <- sub("^k_", "log k_", names(ans$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 if(distimes) ans$distimes <- object$distimes diff --git a/R/predict.mkinmod.R b/R/predict.mkinmod.R new file mode 100644 index 0000000..279fa52 --- /dev/null +++ b/R/predict.mkinmod.R @@ -0,0 +1,89 @@ +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 new file mode 100644 index 0000000..4b8bfd1 --- /dev/null +++ b/R/transform_odeparms.R @@ -0,0 +1,46 @@ +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)) + if (length(index_k) > 0) { + transparms[index_k] <- exp(odeparms[index_k]) + } + + # Go through state variables and apply inverse 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) + 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 + } + } + 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)) + if (length(index_k) > 0) { + odeparms[index_k] <- log(transparms[index_k]) + } + + # Go through state variables and apply 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))) + names(f) <- f_names + odeparms[indices_f] <- f + } + } + return(odeparms) +} diff --git a/TODO b/TODO index 3da3339..8f89002 100644 --- a/TODO +++ b/TODO @@ -1,5 +1,6 @@ - Fix analytical and Eivenvalue based solutions after transition to fitting transformed k and f values +- Fix FOMC model with FOCUS_2006_C - Transfer calculation of DT50 and DT90 calculations from mkinfit to summary.mkinfit - Calculate back-transformed parameters in summary.mkinfit - Report back-transformed parameters in print.summary.mkinfit diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd index 24dcb83..17e0a11 100644 --- a/man/mkinfit.Rd +++ b/man/mkinfit.Rd @@ -34,11 +34,9 @@ mkinfit(mkinmod, observed, \code{\link{modFit}}. } \item{parms.ini}{ - A named vector if initial values for the parameters, including parameters to + A named vector of initial values for the parameters, including parameters to be optimised and potentially also fixed parameters as indicated by \code{fixed_parms}. - The default is to set the initial values to 0.1. The setting of the initial values for - the parameters has a strong impart on performance and it lies in the responsibilty of the - user to set sensible initial values. + If set to "auto", initial values for rate constants are set to default values. } \item{state.ini}{ A named vector of initial values for the state variables of the model. In case the -- cgit v1.2.1