From 1718d434efae26de02754c6622c43f4dc9e624b9 Mon Sep 17 00:00:00 2001 From: jranke Date: Thu, 15 Mar 2012 15:54:14 +0000 Subject: Update kinfit and mkin to the latest version published on BerliOS. git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@17 edb9625f-4e0d-4859-8d74-9fd3b1da38cb --- R/DFOP.solution.R | 5 + R/FOMC.solution.R | 4 + R/HS.solution.R | 6 + R/SFO.solution.R | 4 + R/SFORB.solution.R | 9 ++ R/mkin_long_to_wide.R | 23 ++- R/mkin_wide_to_long.R | 20 +++ R/mkinerrmin.R | 20 +++ R/mkinfit.R | 398 +++++++++++++++++++++++++++++++++++++++++--------- R/mkinmod.R | 154 ++++++++++++++++--- R/mkinplot.R | 83 +++++++++-- R/mkinresplot.R | 54 +++++++ R/mkinstart.R | 23 +++ 13 files changed, 694 insertions(+), 109 deletions(-) create mode 100644 R/DFOP.solution.R create mode 100644 R/FOMC.solution.R create mode 100644 R/HS.solution.R create mode 100644 R/SFO.solution.R create mode 100644 R/SFORB.solution.R create mode 100644 R/mkinresplot.R create mode 100644 R/mkinstart.R (limited to 'R') diff --git a/R/DFOP.solution.R b/R/DFOP.solution.R new file mode 100644 index 0000000..8531cfe --- /dev/null +++ b/R/DFOP.solution.R @@ -0,0 +1,5 @@ +DFOP.solution <- function(t, parent.0, k1, k2, g) +{ + parent = g * parent.0 * exp(-k1 * t) + + (1 - g) * parent.0 * exp(-k2 * t) +} diff --git a/R/FOMC.solution.R b/R/FOMC.solution.R new file mode 100644 index 0000000..8bd13d6 --- /dev/null +++ b/R/FOMC.solution.R @@ -0,0 +1,4 @@ +FOMC.solution <- function(t, parent.0, alpha, beta) +{ + parent = parent.0 / (t/beta + 1)^alpha +} diff --git a/R/HS.solution.R b/R/HS.solution.R new file mode 100644 index 0000000..4651a6a --- /dev/null +++ b/R/HS.solution.R @@ -0,0 +1,6 @@ +HS.solution <- function(t, parent.0, k1, k2, tb) +{ + parent = ifelse(t <= tb, + parent.0 * exp(-k1 * t), + parent.0 * exp(-k1 * tb) * exp(-k2 * (t - tb))) +} diff --git a/R/SFO.solution.R b/R/SFO.solution.R new file mode 100644 index 0000000..3a376e4 --- /dev/null +++ b/R/SFO.solution.R @@ -0,0 +1,4 @@ +SFO.solution <- function(t, parent.0, k) +{ + parent = parent.0 * exp(-k * t) +} diff --git a/R/SFORB.solution.R b/R/SFORB.solution.R new file mode 100644 index 0000000..45a4533 --- /dev/null +++ b/R/SFORB.solution.R @@ -0,0 +1,9 @@ +SFORB.solution = function(t, parent.0, k_12, k_21, k_1output) { + 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 + + parent = parent.0 * + (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) + + ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)) +} diff --git a/R/mkin_long_to_wide.R b/R/mkin_long_to_wide.R index 45ee0a1..ca09b49 100644 --- a/R/mkin_long_to_wide.R +++ b/R/mkin_long_to_wide.R @@ -1,7 +1,28 @@ -mkin_long_to_wide <- function(long_data, time = "time") +# $Id: mkin_long_to_wide.R 96 2011-04-29 11:10:40Z jranke $ + +# Copyright (C) 2010-2011 Johannes Ranke +# Contact: mkin-devel@lists.berlios.de + +# This file is part of the R package mkin + +# mkin is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + +mkin_long_to_wide <- function(long_data, time = "time", outtime = "time") { colnames <- unique(long_data$name) wide_data <- data.frame(time = subset(long_data, name == colnames[1], time)) + names(wide_data) <- outtime for (var in colnames) { wide_data[var] <- subset(long_data, name == var, value) } diff --git a/R/mkin_wide_to_long.R b/R/mkin_wide_to_long.R index 21ab77b..995a9be 100644 --- a/R/mkin_wide_to_long.R +++ b/R/mkin_wide_to_long.R @@ -1,3 +1,23 @@ +# $Id: mkin_wide_to_long.R 59 2010-07-28 12:29:15Z jranke $ + +# Copyright (C) 2010 Johannes Ranke +# Contact: mkin-devel@lists.berlios.de + +# This file is part of the R package mkin + +# mkin is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + mkin_wide_to_long <- function(wide_data, time = "t") { colnames <- names(wide_data) diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R index 7a922eb..ebed191 100644 --- a/R/mkinerrmin.R +++ b/R/mkinerrmin.R @@ -1,3 +1,23 @@ +# $Id: mkinerrmin.R 59 2010-07-28 12:29:15Z jranke $ + +# Copyright (C) 2010 Johannes Ranke +# Contact: mkin-devel@lists.berlios.de + +# This file is part of the R package mkin + +# mkin is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + mkinerrmin <- function(errdata, n.parms, alpha = 0.05) { means.mean <- mean(errdata$value_mean, na.rm=TRUE) diff --git a/R/mkinfit.R b/R/mkinfit.R index f8285fc..b742742 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -1,11 +1,36 @@ +# $Id: mkinfit.R 120 2011-09-02 14:25:35Z jranke $ + +# Copyright (C) 2010-2011 Johannes Ranke +# Contact: mkin-devel@lists.berlios.de +# The summary function is an adapted and extended version of summary.modFit +# from the FME package, v 1.1 by Soetart and Petzoldt, which was in turn +# inspired by summary.nls.lm + +# This file is part of the R package mkin + +# mkin is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + mkinfit <- function(mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$parms)), state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), lower = 0, upper = Inf, fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], + eigen = FALSE, plot = FALSE, quiet = FALSE, err = NULL, weight = "none", scaleVar = FALSE, + atol = 1e-6, ...) { mod_vars <- names(mkinmod$diffs) @@ -16,18 +41,6 @@ mkinfit <- function(mkinmod, observed, # 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 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))) - } # Name the inital parameter values if they are not named yet if(is.null(names(state.ini))) names(state.ini) <- mod_vars @@ -45,6 +58,40 @@ mkinfit <- function(mkinmod, observed, names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_") } + # Set upper limit for formation fractions to one if formation fractions are + # directly defined and if no user input for upper limit is given + if (all(upper==Inf) & any(grepl("f_", names(parms.ini)))==TRUE){ + upper=c( rep(Inf,length(parms.optim))) + upper[grep("f_", names(parms.optim))]=1 + upper=c(rep(Inf, length(state.ini.optim)), upper) + } + + # Decide if the solution of the model can be based on a simple analytical + # formula, the spectral decomposition of the matrix (fundamental system) + # or a numeric ode solver from the deSolve package + if (length(mkinmod$map) == 1) { + solution = "analytical" + } else { + if (is.matrix(mkinmod$coefmat) & eigen) solution = "eigen" + else solution = "deSolve" + } + + # Create a function calculating the differentials specified by the model + # if necessary + if(solution == "deSolve") { + mkindiff <- function(t, state, parms) { + 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 out_predicted <- NA @@ -60,19 +107,67 @@ 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))) + } - # Solve the ode - out <- ode( - y = odeini, - times = outtimes, - func = mkindiff, - parms = odeparms) - + # 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) { + if((length(mkinmod$map[[var]]) == 1) || solution == "analytical") { out_transformed[var] <- out[, var] } else { out_transformed[var] <- rowSums(out[, mkinmod$map[[var]]]) @@ -90,19 +185,53 @@ 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) - out_plot <- ode( - y = odeini, - times = outtimes_plot, - func = mkindiff, - parms = odeparms) + 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 names(mkinmod$map)) { - if(length(mkinmod$map[[var]]) == 1) { + 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 plot(0, type="n", xlim = range(observed$time), ylim = range(observed$value, na.rm=TRUE), @@ -124,8 +253,17 @@ mkinfit <- function(mkinmod, observed, } fit <- modFit(cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, ...) - # We need the function for plotting - fit$mkindiff <- mkindiff + # We need to return some more data for summary and plotting + fit$solution <- solution + 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 @@ -161,54 +299,109 @@ mkinfit <- function(mkinmod, observed, 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 + if ("k1" %in% names(parms.optim)) n.optim <- n.optim + 1 + if ("k2" %in% names(parms.optim)) n.optim <- n.optim + 1 + if ("g" %in% names(parms.optim)) n.optim <- n.optim + 1 + if ("tb" %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 + # Calculate dissipation times DT50 and DT90 and formation fractions 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) + fit$ff <- vector() + fit$SFORB <- vector() 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 + 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 + for (k_name in k_names) + { + fit$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / 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 == "FOMC") { + alpha = parms.all["alpha"] + beta = parms.all["beta"] + DT50 = beta * (2^(1/alpha) - 1) + DT90 = beta * (10^(1/alpha) - 1) + } + if (type == "DFOP") { + k1 = parms.all["k1"] + k2 = parms.all["k2"] + g = parms.all["g"] + f <- function(t, x) { + fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t) + (fraction - (1 - x/100))^2 } - 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 + DTmax <- 1000 + DT50.o <- optimize(f, c(0.001, DTmax), x=50)$minimum + DT50 = ifelse(DTmax - DT50.o < 0.1, NA, DT50.o) + DT90.o <- optimize(f, c(0.001, DTmax), x=90)$minimum + DT90 = ifelse(DTmax - DT90.o < 0.1, NA, DT90.o) + } + if (type == "HS") { + k1 = parms.all["k1"] + k2 = parms.all["k2"] + tb = parms.all["tb"] + DTx <- function(x) { + DTx.a <- (log(100/(100 - x)))/k1 + DTx.b <- tb + (log(100/(100 - x)) - k1 * tb)/k2 + if (DTx.a < tb) DTx <- DTx.a + else DTx <- DTx.b + return(DTx) } - fit$distimes[obs_var, ] = c(DT50, DT90) + 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) + 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, max_DT))$minimum + if (abs(DT90.o - max_DT) < 0.01) DT90 = NA else DT90 = DT90.o + for (k_out_name in k_out_names) + { + fit$ff[[sub("k_", "", k_out_name)]] = parms.all[[k_out_name]] / k_1output + } + # Return the eigenvalues for comparison with DFOP rate constants + fit$SFORB[[paste(obs_var, "b1", sep="_")]] = b1 + fit$SFORB[[paste(obs_var, "b2", sep="_")]] = b2 + } + fit$distimes[obs_var, ] = c(DT50, DT90) } # Collect observed, predicted and residuals @@ -217,25 +410,74 @@ mkinfit <- function(mkinmod, observed, data$residual <- data$observed - data$predicted data$variable <- ordered(data$variable, levels = obs_vars) fit$data <- data[order(data$variable, data$time), ] + fit$atol <- atol class(fit) <- c("mkinfit", "modFit") return(fit) } -summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, cov = FALSE,...) { - ans <- FME:::summary.modFit(object, cov = cov) +summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, cov = FALSE,...) { + param <- object$par + pnames <- names(param) + p <- length(param) + covar <- try(solve(0.5*object$hessian), silent = TRUE) # unscaled covariance + if (!is.numeric(covar)) { + message <- "Cannot estimate covariance; system is singular" + warning(message) + covar <- matrix(data = NA, nrow = p, ncol = p) + } else message <- "ok" + + rownames(covar) <- colnames(covar) <-pnames + rdf <- object$df.residual + resvar <- object$ssr / rdf + se <- sqrt(diag(covar) * resvar) + names(se) <- pnames + tval <- param / se + modVariance <- object$ssr / length(object$residuals) + + if (!all(object$start$lower >=0)) { + message <- "Note that the one-sided t-test may not be appropriate if + parameter values below zero are possible." + warning(message) + } else message <- "ok" + + param <- cbind(param, se, tval, pt(tval, rdf, lower.tail = FALSE)) + dimnames(param) <- list(pnames, c("Estimate", "Std. Error", + "t value", "Pr(>t)")) + if(cov) + ans <- list(residuals = object$residuals, + residualVariance = resvar, + sigma = sqrt(resvar), + modVariance = modVariance, + df = c(p, rdf), cov.unscaled = covar, + cov.scaled = covar * resvar, + info = object$info, niter = object$iterations, + stopmess = message, + par = param) + else + ans <- list(residuals = object$residuals, + residualVariance = resvar, + sigma = sqrt(resvar), + modVariance = modVariance, + df = c(p, rdf), + info = object$info, niter = object$iterations, + stopmess = message, + par = param) + 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 + if(ff) ans$ff <- object$ff + if(length(object$SFORB) != 0) ans$SFORB <- object$SFORB 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), ...) { +print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), tval = TRUE, ...) { cat("\nEquations:\n") print(noquote(as.character(x[["diffs"]]))) df <- x$df @@ -249,7 +491,11 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . else print(x$fixed) cat("\nOptimised parameters:\n") - printCoefmat(x$par, digits = digits, ...) + if (tval) printCoefmat(x$par, digits = digits, ...) + else { + printCoefmat(x$par[,c(1,2,4)], cs.in = c(1,2), tst.ind = integer(0), + P.values = TRUE, has.Pvalue = TRUE, digits = digits, ...) + } cat("\nResidual standard error:", format(signif(x$sigma, digits)), "on", rdf, "degrees of freedom\n") @@ -260,10 +506,22 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . printdistimes <- !is.null(x$distimes) if(printdistimes){ - cat("\nEstimated disappearance times\n") + cat("\nEstimated disappearance times:\n") print(x$distimes, digits=digits,...) } + printff <- !is.null(x$ff) + if(printff){ + cat("\nEstimated formation fractions:\n") + print(data.frame(ff = x$ff), digits=digits,...) + } + + printSFORB <- !is.null(x$SFORB) + if(printSFORB){ + cat("\nEstimated Eigenvalues of SFORB model(s):\n") + print(x$SFORB, digits=digits,...) + } + printcor <- !is.null(x$cov.unscaled) if (printcor){ Corr <- cov2cor(x$cov.unscaled) diff --git a/R/mkinmod.R b/R/mkinmod.R index 643e855..31c778c 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -1,3 +1,23 @@ +# $Id: mkinmod.R 71 2010-09-12 01:13:36Z jranke $ + +# Copyright (C) 2010 Johannes Ranke +# Contact: mkin-devel@lists.berlios.de + +# This file is part of the R package mkin + +# mkin is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + mkinmod <- function(...) { spec <- list(...) @@ -5,22 +25,36 @@ mkinmod <- function(...) # The returned model will be a list of character vectors, containing # differential equations, parameter names and a mapping from model variables - # to observed variables + # to observed variables. If possible, a matrix representation of the + # differential equations is included parms <- vector() diffs <- vector() map <- list() + if(spec[[1]]$type %in% c("FOMC", "DFOP", "HS")) { + mat = FALSE + if(!is.null(spec[[1]]$to)) { + message <- paste("Only constant formation fractions over time are implemented.", + "Depending on the reason for the time dependence of degradation, this may be unrealistic", + sep="\n") + warning(message) + } else message <- "ok" + } else mat = TRUE # Establish list of differential equations for (varname in obs_vars) { - if(is.null(spec[[varname]]$type)) stop("Every argument to mkinmod must be a list containing a type component") - if(!spec[[varname]]$type %in% c("SFO", "FOMC", "SFORB")) stop("Available types are SFO, FOMC and SFORB only") + if(is.null(spec[[varname]]$type)) stop( + "Every argument to mkinmod must be a list containing a type component") + if(!spec[[varname]]$type %in% c("SFO", "FOMC", "DFOP", "HS", "SFORB")) stop( + "Available types are SFO, FOMC, DFOP, HS and SFORB only") new_parms <- vector() # New (sub)compartments (boxes) needed for the model type new_boxes <- switch(spec[[varname]]$type, SFO = varname, FOMC = varname, + DFOP = varname, + HS = varname, SFORB = paste(varname, c("free", "bound"), sep="_") ) map[[varname]] <- new_boxes @@ -29,25 +63,64 @@ mkinmod <- function(...) # Start a new differential equation for each new box new_diffs <- paste("d_", new_boxes, " =", sep="") - # Construct terms for transfer to sink and add if appropriate + # Turn on sink if not specified otherwise if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE + + # Construct and add FOMC term and add FOMC parameters if needed + if(spec[[varname]]$type == "FOMC") { + if(match(varname, obs_vars) != 1) { + stop("Type FOMC is only possible for the first compartment, which is assumed to be the source compartment") + } + if(spec[[varname]]$sink == FALSE) { + stop("Turning off the sink for the FOMC model is not implemented") + } + # From p. 53 of the FOCUS kinetics report + nonlinear_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", new_boxes[[1]]) + new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term) + new_parms <- c("alpha", "beta") + ff <- vector() + } + + # Construct and add DFOP term and add DFOP parameters if needed + if(spec[[varname]]$type == "DFOP") { + if(match(varname, obs_vars) != 1) { + stop("Type DFOP is only possible for the first compartment, which is assumed to be the source compartment") + } + if(spec[[varname]]$sink == FALSE) { + stop("Turning off the sink for the DFOP model is not implemented") + } + # From p. 57 of the FOCUS kinetics report + nonlinear_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", new_boxes[[1]]) + new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term) + new_parms <- c("k1", "k2", "g") + ff <- vector() + } + + # Construct and add HS term and add HS parameters if needed + if(spec[[varname]]$type == "HS") { + if(match(varname, obs_vars) != 1) { + stop("Type HS is only possible for the first compartment, which is assumed to be the source compartment") + } + if(spec[[varname]]$sink == FALSE) { + stop("Turning off the sink for the HS model is not implemented") + } + # From p. 55 of the FOCUS kinetics report + nonlinear_term <- paste("ifelse(time <= tb, k1, k2)", "*", new_boxes[[1]]) + new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term) + new_parms <- c("k1", "k2", "tb") + ff <- vector() + } + + # Construct terms for transfer to sink and add if appropriate + if(spec[[varname]]$sink) { - # Add first-order sink term to first (or only) box for SFO and SFORB models + # Add first-order sink term to first (or only) box for SFO and SFORB if(spec[[varname]]$type %in% c("SFO", "SFORB")) { k_compound_sink <- paste("k", new_boxes[[1]], "sink", sep="_") sink_term <- paste("-", k_compound_sink, "*", new_boxes[[1]]) new_diffs[[1]] <- paste(new_diffs[[1]], sink_term) new_parms <- k_compound_sink } - if(spec[[varname]]$type == "FOMC") { - if(match(varname, obs_vars) != 1) { - stop("Type FOMC is only allowed for the first compartment, which is assumed to be the source compartment") - } - # From p. 53 of the FOCUS kinetics report - fomc_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", new_boxes[[1]]) - new_diffs[[1]] <- paste(new_diffs[[1]], "-", fomc_term) - new_parms <- c("alpha", "beta") - } } # Add reversible binding if appropriate @@ -66,7 +139,6 @@ mkinmod <- function(...) names(new_diffs) <- new_boxes diffs <- c(diffs, new_diffs) } - # Transfer between compartments for (varname in obs_vars) { to <- spec[[varname]]$to @@ -74,6 +146,8 @@ mkinmod <- function(...) origin_box <- switch(spec[[varname]]$type, SFO = varname, FOMC = varname, + DFOP = varname, + HS = varname, SFORB = paste(varname, "free", sep="_")) fraction_left <- NULL for (target in to) { @@ -82,27 +156,63 @@ mkinmod <- function(...) SFORB = paste(target, "free", sep="_")) if(spec[[varname]]$type %in% c("SFO", "SFORB")) { k_from_to <- paste("k", origin_box, target_box, sep="_") - diffs[[origin_box]] <- paste(diffs[[origin_box]], "-", k_from_to, "*", origin_box) - diffs[[target_box]] <- paste(diffs[[target_box]], "+", k_from_to, "*", origin_box) + diffs[[origin_box]] <- paste(diffs[[origin_box]], "-", + k_from_to, "*", origin_box) + diffs[[target_box]] <- paste(diffs[[target_box]], "+", + k_from_to, "*", origin_box) parms <- c(parms, k_from_to) } - if(spec[[varname]]$type == "FOMC") { + if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS")) { fraction_to_target = paste("f_to", target, sep="_") - fraction_not_to_target = paste("(1 - ", fraction_to_target, ")", sep="") + fraction_not_to_target = paste("(1 - ", fraction_to_target, ")", + sep="") if(is.null(fraction_left)) { fraction_really_to_target = fraction_to_target fraction_left = fraction_not_to_target } else { - fraction_really_to_target = paste(fraction_left, " * ", fraction_to_target, sep="") - fraction_left = paste(fraction_left, " * ", fraction_not_to_target, sep="") + fraction_really_to_target = paste(fraction_left, " * ", + fraction_to_target, sep="") + fraction_left = paste(fraction_left, " * ", + fraction_not_to_target, sep="") } - diffs[[target_box]] <- paste(diffs[[target_box]], "+", fraction_really_to_target, "*", fomc_term) + ff[target_box] = fraction_really_to_target + diffs[[target_box]] <- paste(diffs[[target_box]], "+", + ff[target_box], "*", nonlinear_term) parms <- c(parms, fraction_to_target) } } } } model <- list(diffs = diffs, parms = parms, map = map) + + # Create coefficient matrix if appropriate + if (mat) { + boxes <- names(diffs) + n <- length(boxes) + m <- matrix(nrow=n, ncol=n, dimnames=list(boxes, boxes)) + for (from in boxes) { + for (to in boxes) { + if (from == to) { + k.candidate = paste("k", from, c(boxes, "sink"), sep="_") + k.candidate = sub("free.*bound", "free_bound", k.candidate) + k.candidate = sub("bound.*free", "bound_free", k.candidate) + k.effective = intersect(model$parms, k.candidate) + m[from,to] = ifelse(length(k.effective) > 0, + paste("-", k.effective, collapse = " "), "0") + } else { + k.candidate = paste("k", from, to, sep="_") + k.candidate = sub("free.*bound", "free_bound", k.candidate) + k.candidate = sub("bound.*free", "bound_free", k.candidate) + k.effective = intersect(model$parms, k.candidate) + m[to, from] = ifelse(length(k.effective) > 0, + k.effective, "0") + } + } + } + model$coefmat <- m + } + + if (exists("ff")) model$ff = ff class(model) <- "mkinmod" invisible(model) } diff --git a/R/mkinplot.R b/R/mkinplot.R index d2880c0..d665bc2 100644 --- a/R/mkinplot.R +++ b/R/mkinplot.R @@ -1,5 +1,6 @@ -mkinplot <- function(fit, xlab = "Time", ylab = "Observed", xlim = range(fit$data$time), ylim = range(fit$data$observed, na.rm = TRUE), ...) +mkinplot <- function(fit, xlab = "Time", ylab = "Observed", xlim = range(fit$data$time), ylim = range(fit$data$observed, na.rm = TRUE), legend = TRUE, ...) { + solution = fit$solution fixed <- fit$fixed$value names(fixed) <- rownames(fit$fixed) parms.all <- c(fit$par, fixed) @@ -16,12 +17,60 @@ mkinplot <- function(fit, xlab = "Time", ylab = "Observed", xlim = range(fit$dat rownames(subset(fit$fixed, type == "deparm"))) odeparms <- parms.all[odenames] - # Solve the ode - out <- ode( - y = odeini, - times = outtimes, - func = fit$mkindiff, - parms = odeparms) + # Solve the system + evalparse <- function(string) + { + eval(parse(text=string), as.list(c(odeparms, odeini))) + } + if (solution == "analytical") { + parent.type = names(fit$map[[1]])[1] + parent.name = names(fit$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, "free_bound", sep="_")), + evalparse(paste("k", parent.name, "bound_free", sep="_")), + evalparse(paste("k", parent.name, "free_sink", sep="_"))) + ) + out <- cbind(outtimes, o) + dimnames(out) <- list(outtimes, c("time", parent.name)) + } + if (solution == "eigen") { + coefmat.num <- matrix(sapply(as.vector(fit$coefmat), evalparse), + nrow = length(odeini)) + e <- eigen(coefmat.num) + c <- solve(e$vectors, odeini) + f.out <- function(t) { + e$vectors %*% diag(exp(e$values * t), nrow=length(odeini)) %*% c + } + o <- matrix(mapply(f.out, outtimes), + nrow = length(odeini), ncol = length(outtimes)) + dimnames(o) <- list(names(odeini), NULL) + out <- cbind(time = outtimes, t(o)) + } + if (solution == "deSolve") { + out <- ode( + y = odeini, + times = outtimes, + func = fit$mkindiff, + parms = odeparms, + atol = fit$atol + ) + } # Output transformation for models with unobserved compartments like SFORB out_transformed <- data.frame(time = out[,"time"]) @@ -36,14 +85,16 @@ mkinplot <- function(fit, xlab = "Time", ylab = "Observed", xlim = range(fit$dat # Plot the data and model output plot(0, type="n", xlim = xlim, ylim = ylim, - xlab = xlab, ylab = ylab, ...) - col_obs <- pch_obs <- 1:length(fit$map) - names(col_obs) <- names(pch_obs) <- names(fit$map) - for (obs_var in names(fit$map)) { - points(subset(fit$data, variable == obs_var, c(time, observed)), - 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=names(fit$map), + xlab = xlab, ylab = ylab, ...) + col_obs <- pch_obs <- 1:length(fit$map) + names(col_obs) <- names(pch_obs) <- names(fit$map) + for (obs_var in names(fit$map)) { + points(subset(fit$data, variable == obs_var, c(time, observed)), + pch = pch_obs[obs_var], col = col_obs[obs_var]) + } + matlines(out_transformed$time, out_transformed[-1]) + if (legend == TRUE) { + legend("topright", inset=c(0.05, 0.05), legend=names(fit$map), col=col_obs, pch=pch_obs, lty=1:length(pch_obs)) + } } diff --git a/R/mkinresplot.R b/R/mkinresplot.R new file mode 100644 index 0000000..bede5f7 --- /dev/null +++ b/R/mkinresplot.R @@ -0,0 +1,54 @@ +# $Id: $ + +# Copyright (C) 2008-2011 Katrin Lindenberger and Johannes Ranke +# Contact: mkin-devel@lists.berlios.de + +# This file is part of the R package mkin + +# mkin is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + +mkinresplot <- function (object, obs_vars = vector(), + xlab = "Time [days]", ylab = "Residual [% of applied radioactivity]", + maxabs = "auto", legend= TRUE, lpos = "topright", ...) +{ + obs_vars_all <- as.character(unique(object$data$variable)) + + if (length(obs_vars) > 0){ + vars <- intersect(obs_vars_all, obs_vars) + } else vars <- obs_vars_all + + residuals <- subset(object$data, variable %in% vars, residual) + + if (maxabs == "auto") maxabs = max(abs(residuals), na.rm = TRUE) + + col_obs <- pch_obs <- 1:length(vars) + names(col_obs) <- names(pch_obs) <- vars + + plot(0, xlab = xlab, ylab = ylab, + xlim = c(0, 1.1 * max(object$data$time)), + ylim = c(-1.2 * maxabs, 1.2 * maxabs), ...) + + for(var in vars){ + residuals_plot <- subset(object$data, variable == var, c("time", "residual")) + points(residuals_plot, pch = pch_obs[var], col = col_obs[var]) + } + + abline(h = 0, lty = 2) + title(paste("Residuals of mkin fit"), font.main = 1) + + if (legend == TRUE) { + legend(lpos, inset = c(0.05, 0.05), legend = vars, + col = col_obs, pch = pch_obs) + } +} diff --git a/R/mkinstart.R b/R/mkinstart.R new file mode 100644 index 0000000..46e5bed --- /dev/null +++ b/R/mkinstart.R @@ -0,0 +1,23 @@ +mkinstart <- function(model, data, mode = "auto") +{ + if (class(model) != "mkinmod") stop("The first argument must be a model of class mkinmod") + names <- model$parms + observed <- names(model$map) + if(!all(observed %in% levels(data$name))) stop("The data must contain the observed variables used in the model") + for (obs in observed) + { + tmp <- subset(data, name == obs) + max <- tmp[which.max(tmp$value), ] + type = names(model$map[[obs]])[[1]] + kinmodel <- ifelse(type == "SFORB", "DFOP", type) + tmp.longdata <- subset(data, name == obs & time >= max$time) + tmp.widedata <- mkin_long_to_wide(tmp.longdata, outtime = "t") + names(tmp.widedata) <- c("t", "parent") + tmp.fit <- kinfit( + kindata = tmp.widedata, + kinmodels = kinmodel, + parent.0.user = max$value) + if(class(tmp.fit[[kinmodel]]) == "try-error") stop(paste("Automatic generation of starting parameters failed\nkinfit failed to find a", kinmodel, "fit for", obs)) + tmp.results <- kinresults(tmp.fit) + } +} -- cgit v1.2.1