From 71d43b104999d7aee96d35ff2a9006f739d2df60 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 1 Jul 2014 23:23:52 +0200 Subject: Support formation fractions without sink pathway, updates --- R/endpoints.R | 24 ++++---- R/mkinfit.R | 163 ++++++++++++++++++++++++++++++++++--------------- R/mkinmod.R | 3 - R/transform_odeparms.R | 127 ++++++++++++++++++++++++++------------ 4 files changed, 215 insertions(+), 102 deletions(-) (limited to 'R') diff --git a/R/endpoints.R b/R/endpoints.R index a5bebd71..676831a0 100644 --- a/R/endpoints.R +++ b/R/endpoints.R @@ -1,8 +1,8 @@ endpoints <- function(fit) { - # Calculate dissipation times DT50 and DT90 and, if necessary, formation - # fractions and SFORB eigenvalues from optimised parameters + # Calculate dissipation times DT50 and DT90 and formation + # fractions as well as SFORB eigenvalues from optimised parameters # Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from HS and DFOP, - # as well as from Eigenvalues b1 and b2 of the SFORB model + # as well as from Eigenvalues b1 and b2 of any SFORB models ep <- list() obs_vars <- fit$obs_vars parms.all <- fit$bparms.ode @@ -16,15 +16,17 @@ endpoints <- function(fit) { # Get formation fractions if directly fitted, and calculate remaining fraction to sink f_names = grep(paste("f", obs_var, sep = "_"), names(parms.all), value=TRUE) - f_values = parms.all[f_names] - f_to_sink = 1 - sum(f_values) - names(f_to_sink) = ifelse(type == "SFORB", - paste(obs_var, "free", "sink", sep = "_"), - paste(obs_var, "sink", sep = "_")) - for (f_name in f_names) { - ep$ff[[sub("f_", "", sub("_to_", "_", f_name))]] = f_values[[f_name]] + if (length(f_names) > 0) { + f_values = parms.all[f_names] + f_to_sink = 1 - sum(f_values) + names(f_to_sink) = ifelse(type == "SFORB", + paste(obs_var, "free", "sink", sep = "_"), + paste(obs_var, "sink", sep = "_")) + for (f_name in f_names) { + ep$ff[[sub("f_", "", sub("_to_", "_", f_name))]] = f_values[[f_name]] + } + ep$ff = append(ep$ff, f_to_sink) } - ep$ff = append(ep$ff, f_to_sink) # Get the rest if (type == "SFO") { diff --git a/R/mkinfit.R b/R/mkinfit.R index f1a3e6bb..2480c135 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -59,6 +59,31 @@ mkinfit <- function(mkinmod, observed, " not used in the model") } + # Warn that the sum of formation fractions may exceed one they are not + # fitted in the transformed way + if (mkinmod$use_of_ff == "max" & transform_fractions == FALSE) { + warning("The sum of formation fractions may exceed one if you do not use ", + "transform_fractions = TRUE." ) + for (box in mod_vars) { + # Stop if formation fractions are not transformed and we have no sink + if (mkinmod$spec[[box]]$sink == FALSE) { + stop("If formation fractions are not transformed during the fitting, ", + "it is not supported to turn off pathways to sink.\n ", + "Consider turning on the transformation of formation fractions or ", + "setting up a model with use_of_ff = 'min'.\n") + } + } + } + + # Do not allow fixing formation fractions if we are using the ilr transformation, + # this is not supported + if (transform_fractions == TRUE && length(fixed_parms) > 0) { + if (grepl("^f_", fixed_parms)) { + stop("Fixing formation fractions is not supported when using the ilr ", + "transformation.") + } + } + # Set initial parameter values, including a small increment (salt) # to avoid linear dependencies (singular matrix) in Eigenvalue based solutions k_salt = 0 @@ -72,8 +97,6 @@ mkinfit <- function(mkinmod, observed, # Default values for rate constants for reversible binding if (grepl("free_bound$", parmname)) parms.ini[parmname] = 0.1 if (grepl("bound_free$", parmname)) parms.ini[parmname] = 0.02 - # Default values for formation fractions - if (substr(parmname, 1, 2) == "f_") parms.ini[parmname] = 0.2 # Default values for the FOMC, DFOP and HS models if (parmname == "alpha") parms.ini[parmname] = 1 if (parmname == "beta") parms.ini[parmname] = 10 @@ -82,19 +105,47 @@ mkinfit <- function(mkinmod, observed, if (parmname == "tb") parms.ini[parmname] = 5 if (parmname == "g") parms.ini[parmname] = 0.5 } + # Default values for formation fractions in case they are used + if (mkinmod$use_of_ff == "max") { + for (box in mod_vars) { + f_names <- mkinmod$parms[grep(paste0("^f_", box), mkinmod$parms)] + f_default_names <- intersect(f_names, defaultpar.names) + f_specified_names <- setdiff(f_names, defaultpar.names) + sum_f_specified = sum(parms.ini[f_specified_names]) + if (sum_f_specified > 1) { + stop("Starting values for the formation fractions originating from ", + box, " sum up to more than 1.") + } + if (mkinmod$spec[[box]]$sink) n_unspecified = length(f_default_names) + 1 + else { + n_unspecified = length(f_default_names) + # When we have no sink and only one pathway, we get an implicitly + # fixed parameter + if (length(f_names) == 1) fixed_parms = c(fixed_parms, f_names) + } + parms.ini[f_default_names] <- (1 - sum_f_specified) / n_unspecified + } + } # 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, + transparms.ini <- transform_odeparms(parms.ini, mkinmod, transform_rates = transform_rates, transform_fractions = transform_fractions) # Parameters to be optimised: # Kinetic parameters in parms.ini whose names are not in fixed_parms - parms.fixed <- transparms.ini[fixed_parms] - parms.optim <- transparms.ini[setdiff(names(transparms.ini), fixed_parms)] + parms.fixed <- parms.ini[fixed_parms] + parms.optim <- parms.ini[setdiff(names(parms.ini), fixed_parms)] + + transparms.fixed <- transform_odeparms(parms.fixed, mkinmod, + transform_rates = transform_rates, + transform_fractions = transform_fractions) + transparms.optim <- transform_odeparms(parms.optim, mkinmod, + transform_rates = transform_rates, + transform_fractions = transform_fractions) # Inital state variables in state.ini whose names are not in fixed_initials state.ini.fixed <- state.ini[fixed_initials] @@ -160,9 +211,9 @@ mkinfit <- function(mkinmod, observed, names(odeini) <- state.ini.fixed.boxnames } - odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed) + odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], transparms.fixed) - parms <- backtransform_odeparms(odeparms, mod_vars, + parms <- backtransform_odeparms(odeparms, mkinmod, transform_rates = transform_rates, transform_fractions = transform_fractions) @@ -210,9 +261,9 @@ mkinfit <- function(mkinmod, observed, return(mC) } - lower <- rep(-Inf, length(c(state.ini.optim, parms.optim))) - upper <- rep(Inf, length(c(state.ini.optim, parms.optim))) - names(lower) <- names(upper) <- c(names(state.ini.optim), names(parms.optim)) + lower <- rep(-Inf, length(c(state.ini.optim, transparms.optim))) + upper <- rep(Inf, length(c(state.ini.optim, transparms.optim))) + names(lower) <- names(upper) <- c(names(state.ini.optim), names(transparms.optim)) if (!transform_rates) { index_k <- grep("^k_", names(lower)) lower[index_k] <- 0 @@ -229,7 +280,7 @@ mkinfit <- function(mkinmod, observed, upper[other_fraction_parms] <- 1 } - fit <- modFit(cost, c(state.ini.optim, parms.optim), + fit <- modFit(cost, c(state.ini.optim, transparms.optim), method = method.modFit, control = control.modFit, lower = lower, upper = upper, ...) @@ -280,33 +331,32 @@ mkinfit <- function(mkinmod, observed, fit$obs_vars <- obs_vars fit$predicted <- mkin_wide_to_long(out_predicted, time = "time") - # Collect initial parameter values in two dataframes + # Backtransform parameters + bparms.optim = backtransform_odeparms(fit$par, fit$mkinmod, + transform_rates = transform_rates, + transform_fractions = transform_fractions) + # As backtransform_odeparms does not know about fixed values, it + # generates a formation fraction even it is an implicitly fixed one. + # This needs to be removed from bparms.optim + bparms.optim = bparms.optim[setdiff(names(bparms.optim), names(parms.fixed))] + bparms.fixed = c(state.ini.fixed, parms.fixed) + bparms.all = c(bparms.optim, parms.fixed) + + # Collect initial parameter values in three dataframes fit$start <- data.frame(value = c(state.ini.optim, - backtransform_odeparms(parms.optim, mod_vars, - transform_rates = transform_rates, - transform_fractions = transform_fractions))) + parms.optim)) 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$start$lower_bound = lower - fit$start$upper_bound = upper - - fit$fixed <- data.frame(value = c(state.ini.fixed, - backtransform_odeparms(parms.fixed, mod_vars, - transform_rates = transform_rates, - transform_fractions = transform_fractions))) + + fit$start_transformed = data.frame( + value = c(state.ini.optim, transparms.optim), + lower = lower, + upper = upper) + + fit$fixed <- data.frame(value = c(state.ini.fixed, parms.fixed)) fit$fixed$type = c(rep("state", length(state.ini.fixed)), rep("deparm", length(parms.fixed))) - bparms.optim = backtransform_odeparms(fit$par, mod_vars, - transform_rates = transform_rates, - transform_fractions = transform_fractions) - bparms.fixed = backtransform_odeparms(c(state.ini.fixed, parms.fixed), - mod_vars, - transform_rates = transform_rates, - transform_fractions = transform_fractions) - bparms.all = c(bparms.optim, bparms.fixed) - # Collect observed, predicted, residuals and weighting data <- merge(fit$observed, fit$predicted, by = c("time", "name")) data$name <- ordered(data$name, levels = obs_vars) @@ -327,10 +377,13 @@ mkinfit <- function(mkinmod, observed, fit$reweight.tol <- reweight.tol fit$reweight.max.iter <- reweight.max.iter - # Return all backtransformed parameters for summary + # Return different sets of backtransformed parameters for summary and plotting fit$bparms.optim <- bparms.optim fit$bparms.fixed <- bparms.fixed - fit$bparms.ode <- bparms.all[mkinmod$parms] # Return ode parameters for further fitting + + # Return ode parameters for further fitting + fit$bparms.ode <- bparms.all[mkinmod$parms] + fit$date <- date() class(fit) <- c("mkinfit", "modFit") @@ -340,6 +393,7 @@ mkinfit <- function(mkinmod, observed, summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...) { param <- object$par pnames <- names(param) + bpnames <- names(object$bparms.optim) p <- length(param) mod_vars <- names(object$mkinmod$diffs) covar <- try(solve(0.5*object$hessian), silent = TRUE) # unscaled covariance @@ -365,7 +419,8 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "Lower", "Upper", "t value", "Pr(>|t|)", "Pr(>t)")) - blci <- buci <- numeric() + bparam <- cbind(Estimate = object$bparms.optim, Lower = NA, Upper = NA) + # Transform boundaries of CI for one parameter at a time, # with the exception of sets of formation fractions (single fractions are OK). f_names_skip <- character(0) @@ -376,20 +431,24 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, } for (pname in pnames) { - par.lower <- par.upper <- object$par - par.lower[pname] <- param[pname, "Lower"] - par.upper[pname] <- param[pname, "Upper"] if (!pname %in% f_names_skip) { - blci[pname] <- backtransform_odeparms(par.lower, mod_vars, - object$transform_rates, object$transform_fractions)[pname] - buci[pname] <- backtransform_odeparms(par.upper, mod_vars, - object$transform_rates, object$transform_fractions)[pname] - } else { - blci[pname] <- buci[pname] <- NA - } + par.lower <- param[pname, "Lower"] + par.upper <- param[pname, "Upper"] + names(par.lower) <- names(par.upper) <- pname + bpl <- backtransform_odeparms(par.lower, object$mkinmod, + object$transform_rates, + object$transform_fractions) + bpu <- backtransform_odeparms(par.upper, object$mkinmod, + object$transform_rates, + object$transform_fractions) + # Again, we have to remove formation fractions that were implicitly generated + bpl <- bpl[intersect(bpnames, names(bpl))] + bpu <- bpu[intersect(bpnames, names(bpu))] + + bparam[names(bpl), "Lower"] <- bpl + bparam[names(bpu), "Upper"] <- bpu + } } - bparam <- cbind(object$bparms.optim, blci, buci) - dimnames(bparam) <- list(pnames, c("Estimate", "Lower", "Upper")) ans <- list( version = as.character(packageVersion("mkin")), @@ -416,6 +475,7 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, ans$diffs <- object$mkinmod$diffs if(data) ans$data <- object$data ans$start <- object$start + ans$start_transformed <- object$start_transformed ans$fixed <- object$fixed @@ -451,9 +511,12 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . x$reweight.method) cat("\n") - cat("\nStarting values for optimised parameters:\n") + cat("\nStarting values for parameters to be optimised:\n") print(x$start) + cat("\nStarting values for the transformed parameters actually optimised:\n") + print(x$start_transformed) + cat("\nFixed parameter values:\n") if(length(x$fixed$value) == 0) cat("None\n") else print(x$fixed) @@ -487,8 +550,8 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . } printff <- !is.null(x$ff) - if(printff & x$use_of_ff == "min"){ - cat("\nEstimated formation fractions:\n") + if(printff){ + cat("\nResulting formation fractions:\n") print(data.frame(ff = x$ff), digits=digits,...) } diff --git a/R/mkinmod.R b/R/mkinmod.R index 49ec4a44..fe5e8142 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -164,9 +164,6 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL) diffs[[target_box]] <- paste(diffs[[target_box]], "+", k_from_to, "*", origin_box) } else { - if (!spec[[varname]]$sink) { - stop("Turning off the sink when using formation fractions is not supported") - } fraction_to_target = paste("f", origin_box, "to", target, sep="_") parms <- c(parms, fraction_to_target) diffs[[target_box]] <- paste(diffs[[target_box]], "+", diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R index 6932fdef..4774fcf6 100644 --- a/R/transform_odeparms.R +++ b/R/transform_odeparms.R @@ -16,31 +16,52 @@ # You should have received a copy of the GNU General Public License along with # this program. If not, see -transform_odeparms <- function(parms, mod_vars, +transform_odeparms <- function(parms, mkinmod, transform_rates = TRUE, transform_fractions = TRUE) { + # We need the model specification for the names of the model + # variables and the information on the sink + spec = mkinmod$spec + # Set up container for transformed parameters - transparms <- parms + transparms <- numeric(0) + + # Do not transform initial values for state variables + state.ini.optim <- parms[grep("_0$", names(parms))] + transparms[names(state.ini.optim)] <- state.ini.optim # Log transformation for rate constants if requested - index_k <- grep("^k_", names(transparms)) - if (length(index_k) > 0) { - if(transform_rates) transparms[index_k] <- log(parms[index_k]) - else transparms[index_k] <- parms[index_k] + k <- parms[grep("^k_", names(parms))] + if (length(k) > 0) { + if(transform_rates) { + transparms[paste0("log_", names(k))] <- log(k) + } + else transparms[names(k)] <- k } - # Go through state variables and apply isometric logratio transformation if requested + # Go through state variables and apply isometric logratio transformation to + # formation fractions if requested + mod_vars = names(spec) for (box in mod_vars) { - 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 <- parms[indices_f] - trans_f <- ilr(c(f, 1 - sum(f))) - names(trans_f) <- f_names - if(transform_fractions) transparms[indices_f] <- trans_f - else transparms[indices_f] <- f + f <- parms[grep(paste("^f", box, sep = "_"), names(parms))] + + if (length(f) > 0) { + if(transform_fractions) { + if (spec[[box]]$sink) { + trans_f <- ilr(c(f, 1 - sum(f))) + trans_f_names <- paste("f", box, "ilr", 1:length(trans_f), sep = "_") + transparms[trans_f_names] <- trans_f + } else { + if (length(f) > 1) { + trans_f <- ilr(f) + trans_f_names <- paste("f", box, "ilr", 1:length(trans_f), sep = "_") + transparms[trans_f_names] <- trans_f + } + } + } else { + transparms[names(f)] <- f + } } } @@ -48,55 +69,85 @@ transform_odeparms <- function(parms, mod_vars, # 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[pname] <- ifelse(transform_rates, log(parms[pname]), parms[pname]) - transparms[pname] <- ifelse(transform_rates, log(parms[pname]), parms[pname]) + transparms[paste0("log_", pname)] <- ifelse(transform_rates, log(parms[pname]), parms[pname]) } } if (!is.na(parms["g"])) { g <- parms["g"] - transparms["g"] <- ifelse(transform_fractions, ilr(c(g, 1 - g)), g) + transparms["g_ilr"] <- ifelse(transform_fractions, ilr(c(g, 1 - g)), g) } return(transparms) } -backtransform_odeparms <- function(transparms, mod_vars, +backtransform_odeparms <- function(transparms, mkinmod, transform_rates = TRUE, transform_fractions = TRUE) { + # We need the model specification for the names of the model + # variables and the information on the sink + spec = mkinmod$spec + # Set up container for backtransformed parameters - parms <- transparms + parms <- numeric(0) + + # Do not transform initial values for state variables + state.ini.optim <- transparms[grep("_0$", names(transparms))] + parms[names(state.ini.optim)] <- state.ini.optim # Exponential transformation for rate constants - index_k <- grep("^k_", names(parms)) - if (length(index_k) > 0) { - if(transform_rates) parms[index_k] <- exp(transparms[index_k]) - else parms[index_k] <- transparms[index_k] + if(transform_rates) { + trans_k <- transparms[grep("^log_k_", names(transparms))] + if (length(trans_k) > 0) { + k_names <- gsub("^log_k_", "k_", names(trans_k)) + parms[k_names] <- exp(trans_k) + } + } else { + trans_k <- transparms[grep("^k_", names(transparms))] + parms[names(trans_k)] <- trans_k } # Go through state variables and apply inverse isometric logratio transformation + mod_vars = names(spec) 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) { - f <- invilr(transparms[indices_f])[1:n_paths] # We do not need the last component - names(f) <- f_names - if(transform_fractions) parms[indices_f] <- f - else parms[indices_f] <- transparms[indices_f] + # Get the names as used in the model + f_names = grep(paste("^f", box, sep = "_"), mkinmod$parms, value = TRUE) + # Get the formation fraction parameters + trans_f = transparms[grep(paste("^f", box, sep = "_"), names(transparms))] + + # If we have one formation fraction parameter, but no optimised parameter, + # the one must be unity + if (length(trans_f) == 0 & length(f_names == 1)) { + parms[f_names] = 1 + } + if (length(trans_f) > 0) { + if(transform_fractions) { + f <- invilr(trans_f) + if (spec[[box]]$sink) { + parms[f_names] <- f[1:length(f)-1] + } else { + parms[f_names] <- f + } + } else { + parms[names(trans_f)] <- trans_f + } } } # Transform parameters also for FOMC, DFOP and HS models for (pname in c("alpha", "beta", "k1", "k2", "tb")) { - if (!is.na(transparms[pname])) { - parms[pname] <- ifelse(transform_rates, exp(transparms[pname]), transparms[pname]) + pname_trans = paste0("log_", pname) + if (!is.na(transparms[pname_trans])) { + parms[pname] <- ifelse(transform_rates, + exp(transparms[pname_trans]), + transparms[pname]) } } - if (!is.na(transparms["g"])) { - g <- transparms["g"] - parms["g"] <- ifelse(transform_fractions, invilr(g)[1], g) + if (!is.na(transparms["g_ilr"])) { + g_ilr <- transparms["g_ilr"] + parms["g"] <- ifelse(transform_fractions, invilr(g_ilr)[1], g_ilr) } return(parms) } +# vim: set ts=2 sw=2 expandtab: -- cgit v1.2.1