diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/endpoints.R | 24 | ||||
| -rw-r--r-- | R/mkinfit.R | 163 | ||||
| -rw-r--r-- | R/mkinmod.R | 3 | ||||
| -rw-r--r-- | R/transform_odeparms.R | 127 |
4 files changed, 215 insertions, 102 deletions
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 <http://www.gnu.org/licenses/>
-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:
|