diff options
Diffstat (limited to 'R/mkinfit.R')
| -rw-r--r-- | R/mkinfit.R | 240 |
1 files changed, 222 insertions, 18 deletions
diff --git a/R/mkinfit.R b/R/mkinfit.R index 9651fd66..c1490276 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -1,55 +1,78 @@ -mkinfit <- function(mkinmod, observed,
+mkinfit <- function(mkinmod, observed,
parms.ini = rep(0.1, length(mkinmod$parms)),
state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)),
- fixed_parms = rep(FALSE, length(mkinmod$parms)),
- fixed_initials = c(FALSE, rep(TRUE, length(mkinmod$diffs) - 1)),
- plot = NULL,
+ fixed_parms = NULL,
+ fixed_initials = names(mkinmod$diffs)[-1],
+ plot = FALSE, quiet = FALSE,
err = NULL, weight = "none", scaleVar = FALSE,
...)
{
+ mod_vars <- names(mkinmod$diffs)
+ # Subset dataframe with mapped (modelled) variables
+ observed <- subset(observed, name %in% names(mkinmod$map))
+ # Get names of observed variables
+ obs_vars = unique(as.character(observed$name))
+
# Name the parameters if they are not named yet
if(is.null(names(parms.ini))) names(parms.ini) <- mkinmod$parms
# Create a function calculating the differentials specified by the model
mkindiff <- function(t, state, parms) {
+ time <- t
diffs <- vector()
- for (box in names(mkinmod$diffs))
+ for (box in mod_vars)
{
diffname <- paste("d", box, sep="_")
- diffs[diffname] <- with(as.list(c(state, parms)),
+ diffs[diffname] <- with(as.list(c(time,state, parms)),
eval(parse(text=mkinmod$diffs[[box]])))
}
return(list(c(diffs)))
}
# Name the inital parameter values if they are not named yet
- if(is.null(names(state.ini))) names(state.ini) <- names(mkinmod$diffs)
+ if(is.null(names(state.ini))) names(state.ini) <- mod_vars
- # TODO: Collect parameters to be optimised
- parms.optim <- parms.ini[!fixed_parms]
+ # Parameters to be optimised
parms.fixed <- parms.ini[fixed_parms]
+ optim_parms <- setdiff(names(parms.ini), fixed_parms)
+ parms.optim <- parms.ini[optim_parms]
- state.ini.optim <- state.ini[!fixed_initials]
- state.ini.optim.boxnames <- names(state.ini.optim)
- names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_")
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.boxnames <- names(state.ini.optim)
+ if(length(state.ini.optim) > 0) {
+ names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_")
+ }
+ cost.old <- 1e100
+ calls <- 0
+ out_predicted <- NA
# Define the model cost function
cost <- function(P)
{
+ assign("calls", calls+1, inherits=TRUE)
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))
} else odeini <- state.ini.fixed
odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed)
- # Solve the ODE
+
+ # Get more timepoints if plotting is desired
+ if(plot) {
+ outtimes = unique(sort(c(observed$time,
+ seq(min(observed$time), max(observed$time), length.out=100))))
+ } else outtimes = unique(observed$time)
+
+ # Solve the ode
out <- ode(
y = odeini,
- times = unique(observed$time),
+ times = outtimes,
func = mkindiff,
parms = odeparms)
- # Output transformation for models with ghost compartments like SFORB
+
+ # 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) {
@@ -58,9 +81,190 @@ mkinfit <- function(mkinmod, observed, out_transformed[var] <- rowSums(out[, mkinmod$map[[var]]])
}
}
+ assign("out_predicted", subset(out_transformed, time %in% observed$time), inherits=TRUE)
+
+ mC <- modCost(out_transformed, observed, y = "value",
+ err = err, weight = weight, scaleVar = scaleVar)
+
+ # Report and/or plot if the model is improved
+ if (mC$model < cost.old) {
+ if(!quiet) cat("Model cost at call ", calls, ": ", mC$model, "\n")
+
+ # Plot the data and the current model output if requested
+ if(plot) {
+ plot(0, type="n",
+ xlim = range(observed$time), ylim = range(observed$value, na.rm=TRUE),
+ xlab = "Time", ylab = "Observed")
+ col_obs <- pch_obs <- 1:length(obs_vars)
+ names(col_obs) <- names(pch_obs) <- obs_vars
+ for (obs_var in obs_vars) {
+ points(subset(observed, name == obs_var, c(time, value)),
+ pch = pch_obs[obs_var], col = col_obs[obs_var])
+ }
+ matlines(out_transformed$time, out_transformed[-1])
+ legend("topright", inset=c(0.05, 0.05), legend=obs_vars,
+ col=col_obs, pch=pch_obs, lty=1:length(pch_obs))
+ }
- return(modCost(out_transformed, observed, y = "value",
- err = err, weight = weight, scaleVar = scaleVar))
+ assign("cost.old", mC$model, inherits=TRUE)
+ }
+ return(mC)
}
- modFit(cost, c(state.ini.optim, parms.optim), ...)
+ fit <- modFit(cost, c(state.ini.optim, parms.optim), ...)
+
+ # We need the function for plotting
+ fit$mkindiff <- mkindiff
+
+ # We also need various other information for summary and plotting
+ fit$map <- mkinmod$map
+ fit$diffs <- mkinmod$diffs
+ fit$observed <- mkin_long_to_wide(observed)
+ predicted_long <- mkin_wide_to_long(out_predicted, time = "time")
+ fit$predicted <- out_predicted
+
+ # 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)))
+ if (exists("lower")) fit$start$lower <- lower
+ if (exists("upper")) fit$start$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)))
+
+ # Calculate chi2 error levels according to FOCUS (2006)
+ means <- aggregate(value ~ time + name, data = observed, mean, na.rm=TRUE)
+ errdata <- merge(means, predicted_long, by = c("time", "name"), suffixes = c("_mean", "_pred"))
+ errdata <- errdata[order(errdata$time, errdata$name), ]
+ errmin.overall <- mkinerrmin(errdata, length(parms.optim) + length(state.ini.optim))
+
+ errmin <- data.frame(err.min = errmin.overall$err.min,
+ n.optim = errmin.overall$n.optim, df = errmin.overall$df)
+ rownames(errmin) <- "All data"
+ for (obs_var in obs_vars)
+ {
+ errdata.var <- subset(errdata, name == obs_var)
+ n.k.optim <- length(grep(paste("k", obs_var, sep="_"), names(parms.optim)))
+ n.initials.optim <- length(grep(paste(obs_var, ".*", "_0", sep=""), names(state.ini.optim)))
+ n.optim <- n.k.optim + n.initials.optim
+ if ("alpha" %in% names(parms.optim)) n.optim <- n.optim + 1
+ if ("beta" %in% names(parms.optim)) n.optim <- n.optim + 1
+ errmin.tmp <- mkinerrmin(errdata.var, n.optim)
+ errmin[obs_var, c("err.min", "n.optim", "df")] <- errmin.tmp
+ }
+ fit$errmin <- errmin
+
+ # Calculate dissipation times DT50 and DT90
+ parms.all = c(fit$par, parms.fixed)
+ fit$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), DT90 = rep(NA, length(obs_vars)),
+ row.names = obs_vars)
+ for (obs_var in obs_vars) {
+ type = names(mkinmod$map[[obs_var]])[1]
+ if (type == "SFO") {
+ k_names = grep(paste("k", obs_var, sep="_"), names(parms.all), value=TRUE)
+ k_tot = sum(parms.all[k_names])
+ DT50 = log(2)/k_tot
+ DT90 = log(10)/k_tot
+ }
+ if (type == "FOMC") {
+ alpha = parms.all["alpha"]
+ beta = parms.all["beta"]
+ DT50 = beta * (2^(1/alpha) - 1)
+ DT90 = beta * (10^(1/alpha) - 1)
+ }
+ if (type == "SFORB") {
+ # FOCUS kinetics (2006), p. 60 f
+ k_out_names = grep(paste("k", obs_var, "free", sep="_"), names(parms.all), value=TRUE)
+ k_out_names = setdiff(k_out_names, paste("k", obs_var, "free", "bound", sep="_"))
+ k_1output = sum(parms.all[[k_out_names]])
+ k_12 = parms.all[[paste("k", obs_var, "free", "bound", sep="_")]]
+ k_21 = parms.all[[paste("k", obs_var, "bound", "free", sep="_")]]
+
+ sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21)
+ b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
+ b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp
+
+ SFORB_fraction = function(t) {
+ ((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +
+ ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)
+ }
+ f_50 <- function(t) (SFORB_fraction(t) - 0.5)^2
+ max_DT <- 1000
+ DT50.o <- optimize(f_50, c(0.01, max_DT))$minimum
+ if (abs(DT50.o - max_DT) < 0.01) DT50 = NA else DT50 = DT50.o
+ f_90 <- function(t) (SFORB_fraction(t) - 0.1)^2
+ DT90.o <- optimize(f_90, c(0.01, 1000))$minimum
+ if (abs(DT90.o - max_DT) < 0.01) DT90 = NA else DT90 = DT90.o
+ }
+ fit$distimes[obs_var, ] = c(DT50, DT90)
+ }
+
+ # Collect observed, predicted and residuals
+ data <- merge(observed, predicted_long, by = c("time", "name"))
+ names(data) <- c("time", "variable", "observed", "predicted")
+ data$residual <- data$observed - data$predicted
+ data$variable <- ordered(data$variable, levels = obs_vars)
+ fit$data <- data[order(data$variable, data$time), ]
+
+ class(fit) <- c("mkinfit", "modFit")
+ return(fit)
+}
+
+summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, cov = FALSE,...) {
+ ans <- FME:::summary.modFit(object, cov = cov)
+ ans$diffs <- object$diffs
+ if(data) ans$data <- object$data
+ ans$start <- object$start
+ ans$fixed <- object$fixed
+ ans$errmin <- object$errmin
+ if(distimes) ans$distimes <- object$distimes
+ class(ans) <- c("summary.mkinfit", "summary.modFit")
+ return(ans)
+}
+
+# Expanded from print.summary.modFit
+print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), ...) {
+ cat("\nEquations:\n")
+ print(noquote(as.character(x[["diffs"]])))
+ df <- x$df
+ rdf <- df[2]
+
+ cat("\nStarting values for optimised parameters:\n")
+ print(x$start)
+
+ cat("\nFixed parameter values:\n")
+ if(length(x$fixed$value) == 0) cat("None\n")
+ else print(x$fixed)
+
+ cat("\nOptimised parameters:\n")
+ printCoefmat(x$par, digits = digits, ...)
+
+ cat("\nResidual standard error:",
+ format(signif(x$sigma, digits)), "on", rdf, "degrees of freedom\n")
+
+ cat("\nChi2 error levels in percent:\n")
+ x$errmin$err.min <- 100 * x$errmin$err.min
+ print(x$errmin, digits=digits,...)
+
+ printdistimes <- !is.null(x$distimes)
+ if(printdistimes){
+ cat("\nEstimated disappearance times\n")
+ print(x$distimes, digits=digits,...)
+ }
+
+ printcor <- !is.null(x$cov.unscaled)
+ if (printcor){
+ Corr <- cov2cor(x$cov.unscaled)
+ rownames(Corr) <- colnames(Corr) <- rownames(x$par)
+ cat("\nParameter correlation:\n")
+ print(Corr, digits = digits, ...)
+ }
+
+ printdata <- !is.null(x$data)
+ if (printdata){
+ cat("\nData:\n")
+ print(format(x$data, digits = digits, scientific = FALSE,...), row.names = FALSE)
+ }
+
+ invisible(x)
}
|