diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/kinerrmin.R | 20 | ||||
| -rw-r--r-- | R/kinfit.R | 115 | ||||
| -rw-r--r-- | R/kinobject.R | 20 | ||||
| -rw-r--r-- | R/kinobjects.R | 20 | ||||
| -rw-r--r-- | R/kinplot.R | 128 | ||||
| -rw-r--r-- | R/kinreport.R | 97 | ||||
| -rw-r--r-- | R/kinresplot.R | 23 | ||||
| -rw-r--r-- | R/kinresults.R | 175 | ||||
| -rw-r--r-- | R/kinwrite.KinGUI.R | 20 |
9 files changed, 478 insertions, 140 deletions
diff --git a/R/kinerrmin.R b/R/kinerrmin.R index 8ce30a3..dbd64b3 100644 --- a/R/kinerrmin.R +++ b/R/kinerrmin.R @@ -1,3 +1,23 @@ +# $Id: kinerrmin.R 59 2010-07-28 12:29:15Z jranke $
+
+# Copyright (C) 2008-2010 Johannes Ranke
+# Contact: mkin-devel@lists.berlios.de
+
+# This file is part of the R package kinfit
+
+# kinfit 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 <http://www.gnu.org/licenses/>
+
kinerrmin <- function(kinfits, kinmodel = "SFO", alpha = 0.05)
{
m = kinfits[[kinmodel]]
@@ -1,43 +1,76 @@ +# $Id: kinfit.R 116 2011-06-14 08:46:47Z kati $
+
+# Copyright (C) 2008-2010 Johannes Ranke
+# Contact: mkin-devel@lists.berlios.de
+
+# This file is part of the R package kinfit
+
+# kinfit 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 <http://www.gnu.org/licenses/>
+
kinfit <- function(kindata, kinmodels = c("SFO"),
parent.0.user = NA,
+ parent.0.fixed = FALSE,
start.SFO = list(parent.0 = NA, k = NA),
start.FOMC = list(parent.0 = NA, alpha = NA, beta = NA),
start.DFOP = list(parent.0 = NA, k1 = NA, k2 = NA, g = NA),
start.HS = list(parent.0 = NA, k1 = NA, k2 = NA, tb = NA),
- algorithm = "port")
+ algorithm = "default")
{
+
kindata <- subset(kindata, !is.na(kindata$parent))
kinfits <- list()
if (!is.na(parent.0.user)) {
start.SFO$parent.0 = parent.0.user
start.FOMC$parent.0 = parent.0.user
+ start.DFOP$parent.0 = parent.0.user
+ start.HS$parent.0 = parent.0.user
}
lmlogged = lm(log(parent) ~ t, data = kindata)
+ k.est = -coef(lmlogged)[["t"]]
for (kinmodel in kinmodels)
{
if (kinmodel == "SFO") {
if (is.na(start.SFO$parent.0)) {
- start.SFO$parent.0 = max(kindata$parent)
+ start.SFO$parent.0 = max(kindata$parent)
}
if (is.na(start.SFO$k)) {
start.SFO$k = - coef(lmlogged)[["t"]]
}
- kinfits[[kinmodel]] = try(
- nls(parent ~ SFO(t, parent.0, k),
- data = kindata, model = TRUE,
- start = start.SFO,
- algorithm = algorithm), silent=TRUE)
+ if (parent.0.fixed)
+ {
+ start.SFO = start.SFO[-1]
+ kinfits[[kinmodel]] = try(
+ nls(parent ~ SFO(t, parent.0.user, k),
+ data = kindata, model = TRUE,
+ start = start.SFO,
+ algorithm = algorithm), silent=TRUE)
+ } else {
+ kinfits[[kinmodel]] = try(
+ nls(parent ~ SFO(t, parent.0, k),
+ data = kindata, model = TRUE,
+ start = start.SFO,
+ algorithm = algorithm), silent=TRUE)
+ }
+ k.est = coef(kinfits$SFO)[["k"]]
}
- k.est = ifelse(is.na(coef(kinfits$SFO)[["k"]]),
- -coef(lmlogged)[["t"]],
- coef(kinfits$SFO)[["k"]])
if (kinmodel == "FOMC") {
if (is.na(start.FOMC$parent.0)) {
- start.FOMC$parent.0 = max(kindata$parent)
+ start.FOMC$parent.0 = max(kindata$parent)
}
if (is.na(start.FOMC$alpha)) {
start.FOMC$alpha = 1
@@ -45,15 +78,29 @@ kinfit <- function(kindata, kinmodels = c("SFO"), if (is.na(start.FOMC$beta)) {
start.FOMC$beta = start.FOMC$alpha / k.est
}
- kinfits[[kinmodel]] = try(
- nls(parent ~ FOMC(t, parent.0, alpha, beta),
- data = kindata, model = TRUE,
- start = start.FOMC,
- algorithm = algorithm), silent=TRUE)
+
+ if (parent.0.fixed)
+ {
+ start.FOMC = list(alpha = start.FOMC$alpha, beta = start.FOMC$beta)
+
+ kinfits[[kinmodel]] = try(
+ nls(parent ~ FOMC(t, parent.0.user, alpha, beta),
+ data = kindata, model = TRUE,
+ start = start.FOMC,
+ algorithm = algorithm), silent=TRUE)
+
+ } else {
+ kinfits[[kinmodel]] = try(
+ nls(parent ~ FOMC(t, parent.0, alpha, beta),
+ data = kindata, model = TRUE,
+ start = start.FOMC,
+ algorithm = algorithm), silent=TRUE)
+
+ }
}
if (kinmodel == "DFOP") {
if (is.na(start.DFOP$parent.0)) {
- start.DFOP$parent.0 = max(kindata$parent)
+ start.DFOP$parent.0 = max(kindata$parent)
}
if (is.na(start.DFOP$k1)) {
start.DFOP$k1 = k.est * 2
@@ -64,15 +111,26 @@ kinfit <- function(kindata, kinmodels = c("SFO"), if (is.na(start.DFOP$g)) {
start.DFOP$g = 0.5
}
- kinfits[[kinmodel]] = try(
+ if (parent.0.fixed)
+ {
+ start.DFOP = list(k1 = start.DFOP$k1, k2 = start.DFOP$k2, g = start.DFOP$g)
+
+ kinfits[[kinmodel]] = try(
+ nls(parent ~ DFOP(t, parent.0.user, k1, k2, g),
+ data = kindata, model = TRUE,
+ start = start.DFOP,
+ algorithm = algorithm), silent=TRUE)
+ }else{
+ kinfits[[kinmodel]] = try(
nls(parent ~ DFOP(t, parent.0, k1, k2, g),
data = kindata, model = TRUE,
start = start.DFOP,
- algorithm = algorithm), silent=TRUE)
+ algorithm = algorithm), silent=TRUE)
+ }
}
if (kinmodel == "HS") {
if (is.na(start.HS$parent.0)) {
- start.HS$parent.0 = max(kindata$parent)
+ start.HS$parent.0 = max(kindata$parent)
}
if (is.na(start.HS$k1)) {
start.HS$k1 = k.est
@@ -83,11 +141,24 @@ kinfit <- function(kindata, kinmodels = c("SFO"), if (is.na(start.HS$tb)) {
start.HS$tb = 0.05 * max(kindata$t)
}
- kinfits[[kinmodel]] = try(
+
+ if (parent.0.fixed)
+ {
+
+ start.HS = list(k1 = start.HS$k1, k2 = start.HS$k2, tb = start.HS$tb)
+
+ kinfits[[kinmodel]] = try(
+ nls(parent ~ HS(t, parent.0.user, k1, k2, tb),
+ data = kindata, model = TRUE,
+ start = start.HS,
+ algorithm = algorithm), silent=TRUE)
+ }else{
+ kinfits[[kinmodel]] = try(
nls(parent ~ HS(t, parent.0, k1, k2, tb),
data = kindata, model = TRUE,
start = start.HS,
- algorithm = algorithm), silent=TRUE)
+ algorithm = algorithm), silent=TRUE)
+ }
}
}
return(kinfits)
diff --git a/R/kinobject.R b/R/kinobject.R index de6f6af..2ecc962 100644 --- a/R/kinobject.R +++ b/R/kinobject.R @@ -1,3 +1,23 @@ +# $Id: kinobject.R 59 2010-07-28 12:29:15Z jranke $
+
+# Copyright (C) 2008-2010 Johannes Ranke
+# Contact: mkin-devel@lists.berlios.de
+
+# This file is part of the R package kinfit
+
+# kinfit 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 <http://www.gnu.org/licenses/>
+
kinobject <- function(parent, type, system,
layers = NA, sampling_times = NA)
{
diff --git a/R/kinobjects.R b/R/kinobjects.R index a767cea..33abe61 100644 --- a/R/kinobjects.R +++ b/R/kinobjects.R @@ -1,3 +1,23 @@ +# $Id: kinobjects.R 59 2010-07-28 12:29:15Z jranke $
+
+# Copyright (C) 2008-2010 Johannes Ranke
+# Contact: mkin-devel@lists.berlios.de
+
+# This file is part of the R package kinfit
+
+# kinfit 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 <http://www.gnu.org/licenses/>
+
kinobjects<- function(parent, type, systems,
layers = NA, sampling_times = NA)
{
diff --git a/R/kinplot.R b/R/kinplot.R index ace1270..394e271 100644 --- a/R/kinplot.R +++ b/R/kinplot.R @@ -1,61 +1,113 @@ +# $Id: kinplot.R 117 2011-06-14 08:52:14Z kati $
+
+# Copyright (C) 2008-2010 Johannes Ranke
+# Contact: mkin-devel@lists.berlios.de
+
+# This file is part of the R package kinfit
+
+# kinfit 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 <http://www.gnu.org/licenses/>
+
kinplot <- function(kinobject,
+ main = "",
xlab = "Time [days]", ylab = "Parent [% of applied radioactivity]",
ylim = c("auto", "auto"),
lpos = "topright")
{
kindata <- na.omit(kinobject$data)
kinfits <- kinobject$fits
- if (ylim[1] == "auto") ylim[1] <- 0
- if (ylim[2] == "auto") ylim[2] <- max(kindata$parent)
- ylim <- as.numeric(ylim)
+ if (ylim[1] == "auto") ylim[1] <- 0
+ if (ylim[2] == "auto") ylim[2] <- max(kindata$parent)
+ ylim <- as.numeric(ylim)
plot(kindata$t, kindata$parent,
+ main = main,
xlab = xlab,
ylab = ylab,
ylim = ylim
- )
+ )
n.m <- length(kinfits)
colors <- ltys <- 1:n.m
names(colors) <- names(ltys) <- names(kinfits)
- ltext <- paste(kinobject$parent, "measured")
+ ltext <- paste(kinobject$parent, "measured")
for (kinmodel in names(kinfits))
{
m = kinfits[[kinmodel]]
if(class(m) == "nls") {
- switch(kinmodel,
- SFO = curve(SFO(x,
- coef(m)[["parent.0"]],
- coef(m)[["k"]]),
- from = min(kindata$t), to = max(kindata$t), add=TRUE,
- col = colors[[kinmodel]],
- lty = ltys[[kinmodel]]),
- FOMC = curve(FOMC(x,
- coef(m)[["parent.0"]],
- coef(m)[["alpha"]],
- coef(m)[["beta"]]),
- from = min(kindata$t), to = max(kindata$t), add=TRUE,
- col = colors[[kinmodel]],
- lty = ltys[[kinmodel]]),
- HS = curve(HS(x,
- coef(m)[["parent.0"]],
- coef(m)[["k1"]],
- coef(m)[["k2"]],
- coef(m)[["tb"]]),
- from = min(kindata$t), to = max(kindata$t), add=TRUE,
- col = colors[[kinmodel]],
- lty = ltys[[kinmodel]]),
- DFOP = curve(DFOP(x,
- coef(m)[["parent.0"]],
- coef(m)[["k1"]],
- coef(m)[["k2"]],
- coef(m)[["g"]]),
- from = min(kindata$t), to = max(kindata$t), add=TRUE,
- col = colors[[kinmodel]],
- lty = ltys[[kinmodel]]))
- ltext <- c(ltext, paste("Fitted", kinmodel, "model"))
+ if (!"parent.0" %in% names(coef(m))) {
+ switch(kinmodel,
+ SFO = lines(
+ t <- seq(min(kindata$t), max(kindata$t), length.out=500),
+ predict(m,
+ newdata = data.frame(t)),
+ col = colors[[kinmodel]],
+ lty = ltys[[kinmodel]]),
+ FOMC = lines(
+ t <- seq(min(kindata$t), max(kindata$t), length.out=500),
+ predict(m,
+ newdata = data.frame(t)),
+ col = colors[[kinmodel]],
+ lty = ltys[[kinmodel]]),
+ HS = lines(
+ t <- seq(min(kindata$t), max(kindata$t), length.out=500),
+ predict(m,
+ newdata = data.frame(t)),
+ col = colors[[kinmodel]],
+ lty = ltys[[kinmodel]]),
+ DFOP = lines(
+ t <- seq(min(kindata$t), max(kindata$t), length.out=500),
+ predict(m,
+ newdata = data.frame(t)),
+ col = colors[[kinmodel]],
+ lty = ltys[[kinmodel]])
+ )
+ ltext <- c(ltext, paste("Fitted", kinmodel, "model"))
+ } else {
+ switch(kinmodel,
+ SFO = curve(SFO(x,
+ coef(m)[["parent.0"]],
+ coef(m)[["k"]]),
+ from = min(kindata$t), to = max(kindata$t), add=TRUE,
+ col = colors[[kinmodel]],
+ lty = ltys[[kinmodel]]),
+ FOMC = curve(FOMC(x,
+ coef(m)[["parent.0"]],
+ coef(m)[["alpha"]],
+ coef(m)[["beta"]]),
+ from = min(kindata$t), to = max(kindata$t), add=TRUE,
+ col = colors[[kinmodel]],
+ lty = ltys[[kinmodel]]),
+ HS = curve(HS(x,
+ coef(m)[["parent.0"]],
+ coef(m)[["k1"]],
+ coef(m)[["k2"]],
+ coef(m)[["tb"]]),
+ from = min(kindata$t), to = max(kindata$t), add=TRUE,
+ col = colors[[kinmodel]],
+ lty = ltys[[kinmodel]]),
+ DFOP = curve(DFOP(x,
+ coef(m)[["parent.0"]],
+ coef(m)[["k1"]],
+ coef(m)[["k2"]],
+ coef(m)[["g"]]),
+ from = min(kindata$t), to = max(kindata$t), add=TRUE,
+ col = colors[[kinmodel]],
+ lty = ltys[[kinmodel]]))
+ ltext <- c(ltext, paste("Fitted", kinmodel, "model"))
+ }
} else {
- ltext <- c(ltext, paste(kinmodel, "model failed"))
- ltys[[kinmodel]] <- NA
+ ltext <- c(ltext, paste(kinmodel, "model failed"))
+ ltys[[kinmodel]] <- NA
}
}
legend(lpos, bty="n", inset = 0.05,
diff --git a/R/kinreport.R b/R/kinreport.R index 4156803..d405b8e 100644 --- a/R/kinreport.R +++ b/R/kinreport.R @@ -1,43 +1,86 @@ -kinreport <- function(kinobject, file = NA, vcov = FALSE, endpoint.digits = 1)
+# $Id: kinreport.R 123 2011-11-01 12:26:41Z jranke $
+
+# Copyright (C) 2008-2010 Johannes Ranke
+# Contact: mkin-devel@lists.berlios.de
+
+# This file is part of the R package kinfit
+
+# kinfit 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 <http://www.gnu.org/licenses/>
+
+kinreport <- function(kinobject, file = NA, data = TRUE, R2 = FALSE, vcov = FALSE, endpoint.digits = 1)
{
if (!is.na(file)) {
sink(file, split=TRUE)
}
cat("Parent compound: ", kinobject$parent, "\n")
- if (!is.null(kinobject$label)) cat("Label position:\t\t", kinobject$label, "\n")
+ if (!is.null(kinobject$label)) {
+ cat("Label position:\t\t", kinobject$label, "\n")
+ }
cat("Study type: ", kinobject$type, "\n")
cat("System: ", kinobject$system, "\n")
if (!is.null(kinobject$source)) {
- cat("Source: ", kinobject$source, "\n")
- }
+ cat("Source: ", kinobject$source, "\n")
+ }
+ cat("kinfit version: ", as.character(packageVersion("kinfit")), "\n")
+ cat("R version: ", paste(R.version$major, R.version$minor, sep="."), "\n")
+ cat("Report generated:", date(), "\n")
cat("\n")
+ if (data) {
+ cat("Data:\n")
+ print(kinobject$data)
+ cat("\n")
+ }
fit.names <- names(kinobject$fits)
for (kinmodel in fit.names)
{
- m <- kinobject$fits[[kinmodel]]
- if (!(class(m) == "try-error")) {
- cat("\n\n---\n")
- cat("Nonlinear least squares fit of the", kinmodel, "model\n\n")
- cat("Parameter estimation:\t")
- s <- summary(m)
- df <- s$df[2]
- p <- 1 - pt(s$parameters[,3], df = df)
- parms <- cbind(s$parameters[,c(1,2,3)], "Pr(>t)" = p)
- cat("\n")
- print(parms, digits=3)
- cat("\n")
- if(vcov)
- {
- cat("Variance-covariance matrix:\n")
- print(vcov(m))
- cat("\n")
- }
- cat("Chi2 error estimation:\t",
- round(100 * kinobject$results$stats[kinmodel, "err.min"], digits=2),
- " %\n", sep="")
- cat("\n")
- }
+ m <- kinobject$fits[[kinmodel]]
+ if (!(class(m) == "try-error")) {
+ cat("\n\n---\n")
+ cat("Nonlinear least squares fit of the", kinmodel, "model\n")
+ if (!"parent.0" %in% names(coef(m))) {
+ cat(paste("Initial value of parent fixed to ", m$model$parent.0.user, "\n", sep=""))
+ }
+ cat("\n")
+ cat("Parameter estimation:\t")
+ s <- summary(m)
+ df <- s$df[2]
+ p <- 1 - pt(s$parameters[,3], df = df)
+ parms <- matrix(nrow = nrow(s$parameters), ncol=4)
+ dimnames(parms) = list(rownames(s$parameters),
+ c("Estimate", "Std. Error", "t value", "Pr(>t)"))
+ parms[, c(1,2,3)] <- s$parameters[,c(1,2,3)]
+ parms[, 4] <- p
+ cat("\n")
+ print(parms, digits=3)
+ cat("\n")
+ if(vcov)
+ {
+ cat("Variance-covariance matrix:\n")
+ print(vcov(m))
+ cat("\n")
+ }
+ cat("Chi2 error estimation: ",
+ round(100 * kinobject$results$stats[kinmodel, "err.min"], digits=2),
+ " %\n", sep="")
+ cat("\n")
+ if(R2)
+ {
+ cat("Coefficient of determination R2: ",
+ round(kinobject$results$stats[kinmodel, "R2"], digits=3), "\n")
+ }
+ }
}
cat("\n\n---\n")
cat("Endpoint estimates\n\n")
diff --git a/R/kinresplot.R b/R/kinresplot.R index be0a85d..36eb55a 100644 --- a/R/kinresplot.R +++ b/R/kinresplot.R @@ -1,3 +1,23 @@ +# $Id: kinresplot.R 106 2011-05-12 10:55:37Z jranke $
+
+# Copyright (C) 2008-2010 Johannes Ranke
+# Contact: mkin-devel@lists.berlios.de
+
+# This file is part of the R package kinfit
+
+# kinfit 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 <http://www.gnu.org/licenses/>
+
kinresplot <- function(kinobject, kinmodel,
xlab = "Time [days]", ylab = "Residual [% of applied radioactivity]",
maxabs = "auto")
@@ -5,10 +25,11 @@ kinresplot <- function(kinobject, kinmodel, m <- kinobject$fits[[kinmodel]]
t <- m$model$t
residuals <- residuals(m)
- if (maxabs == "auto") maxabs = max(abs(residuals))
+ if (maxabs == "auto") maxabs = max(abs(residuals), na.rm = TRUE)
plot(t, residuals,
xlab = xlab,
ylab = ylab,
ylim = c( -1.2 * maxabs, 1.2 * maxabs))
+ abline(h=0, lty=2)
title(paste("Residuals of", kinmodel, "fit"), font.main = 1)
}
diff --git a/R/kinresults.R b/R/kinresults.R index 6bbff28..27bdafb 100644 --- a/R/kinresults.R +++ b/R/kinresults.R @@ -1,74 +1,145 @@ -kinresults <- function(kinfits, alpha = 0.05, SFORB=TRUE)
+# $Id: kinresults.R 124 2011-11-09 10:58:27Z jranke $
+
+# Copyright (C) 2008-2011 Johannes Ranke
+# Contact: mkin-devel@lists.berlios.de
+
+# This file is part of the R package kinfit
+
+# kinfit 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 <http://www.gnu.org/licenses/>
+
+kinresults <- function(kinfits, alpha = 0.05, DTmax = 1000, SFORB=TRUE)
{
- kindata <- data.frame(t = kinfits[[1]]$model$t, parent = kinfits[[1]]$model$parent)
- kindata.means <- aggregate(kindata, list(kindata$t), mean)
+ kindata <- data.frame(
+ t = kinfits[[1]]$model$t,
+ parent = kinfits[[1]]$model$parent)
+ kindata.means <- aggregate(kindata, list(kindata$t), mean)
kindata.means.mean <- mean(kindata.means$parent, na.rm=TRUE)
n.times <- length(kindata.means$parent)
parms <- list()
- df <- err.min <- RSS <- vector()
+ df <- err.min <- R2 <- RSS <- TSS <- RSS.means <- TSS.means <- vector()
DT50 <- DT90 <- vector()
f <- list()
+
for (kinmodel in names(kinfits))
{
m = kinfits[[kinmodel]]
if(class(m) == "nls") {
kindata.means$est <- predict(m, kindata.means)
- parms[[kinmodel]] <- switch(kinmodel,
- SFO = list(parent.0 = coef(m)[["parent.0"]],
- k = coef(m)[["k"]]),
- FOMC = list(parent.0 = coef(m)[["parent.0"]],
- alpha = coef(m)[["alpha"]],
- beta = coef(m)[["beta"]]),
- HS = list(parent.0 = coef(m)[["parent.0"]],
- k1 = coef(m)[["k1"]],
- k2 = coef(m)[["k2"]],
- tb = coef(m)[["tb"]]),
- DFOP = list(parent.0 = coef(m)[["parent.0"]],
- k1 = coef(m)[["k1"]],
- k2 = coef(m)[["k2"]],
- g = coef(m)[["g"]]))
- if(kinmodel == "DFOP" & SFORB) {
- k1 = coef(m)[["k1"]]
- k2 = coef(m)[["k2"]]
- g = coef(m)[["g"]]
- parms[["SFORB"]] =
- list(parent.0 = coef(m)[["parent.0"]],
- k1out = g * k1 + (1 - g) * k2,
- k21 = k1 * k2 / (g * k1 + (1 - g) * k2),
- k12 = (g * (1 - g) * (k1 - k2)^2) / (g * k1 + (1 - g) * k2))
+
+ if (!"parent.0" %in% names(coef(m))) {
+
+ parms[[kinmodel]] <- switch(kinmodel,
+ SFO = list(k = coef(m)[["k"]]),
+ FOMC = list(alpha = coef(m)[["alpha"]],
+ beta = coef(m)[["beta"]]),
+ HS = list(k1 = coef(m)[["k1"]],
+ k2 = coef(m)[["k2"]],
+ tb = coef(m)[["tb"]]),
+ DFOP= list(k1 = coef(m)[["k1"]],
+ k2 = coef(m)[["k2"]],
+ g = coef(m)[["g"]]))
+ } else {
+ parms[[kinmodel]] <- switch(kinmodel,
+ SFO = list(parent.0 = coef(m)[["parent.0"]],
+ k = coef(m)[["k"]]),
+ FOMC = list(parent.0 = coef(m)[["parent.0"]],
+ alpha = coef(m)[["alpha"]],
+ beta = coef(m)[["beta"]]),
+ HS = list(parent.0 = coef(m)[["parent.0"]],
+ k1 = coef(m)[["k1"]],
+ k2 = coef(m)[["k2"]],
+ tb = coef(m)[["tb"]]),
+ DFOP = list(parent.0 = coef(m)[["parent.0"]],
+ k1 = coef(m)[["k1"]],
+ k2 = coef(m)[["k2"]],
+ g = coef(m)[["g"]]))
+ if(kinmodel == "DFOP" & SFORB) {
+ k1 = coef(m)[["k1"]]
+ k2 = coef(m)[["k2"]]
+ g = coef(m)[["g"]]
+ parms[["SFORB"]] =
+ list(parent.0 = coef(m)[["parent.0"]],
+ k1out = g * k1 + (1 - g) * k2,
+ k21 = k1 * k2 / (g * k1 + (1 - g) * k2),
+ k12 = (g * (1 - g) * (k1 - k2)^2) / (g * k1 + (1 - g) * k2))
+ }
+ }
+
+ n.parms = length(coef(m))
+ if (!"parent.0" %in% names(coef(m))) {
+
+ f[[kinmodel]] = switch(kinmodel,
+ HS = function(t, x) {
+ (HS(t, kinfits[[kinmodel]]$model[["parent.0.user"]],
+ coef(m)[["k1"]], coef(m)[["k2"]], coef(m)[["tb"]]) -
+ (1 - x/100) * kinfits[[kinmodel]]$model[["parent.0.user"]])^2
+ },
+ DFOP = function(t, x) {
+ (DFOP(t, kinfits[[kinmodel]]$model[["parent.0.user"]],
+ coef(m)[["k1"]], coef(m)[["k2"]], coef(m)[["g"]]) -
+ (1 - x/100) * kinfits[[kinmodel]]$model[["parent.0.user"]])^2
+ }
+ )
+ }else{
+ f[[kinmodel]] = switch(kinmodel,
+ HS = function(t, x) {
+ (HS(t, coef(m)[["parent.0"]],
+ coef(m)[["k1"]], coef(m)[["k2"]], coef(m)[["tb"]]) -
+ (1 - x/100) * coef(m)[["parent.0"]])^2
+ },
+ DFOP = function(t, x) {
+ (DFOP(t, coef(m)[["parent.0"]],
+ coef(m)[["k1"]], coef(m)[["k2"]], coef(m)[["g"]]) -
+ (1 - x/100) * coef(m)[["parent.0"]])^2
+ }
+ )
}
- n.parms = length(coef(m))
- f[[kinmodel]] = switch(kinmodel,
- HS = function(t, x) {
- (HS(t, coef(m)[["parent.0"]],
- coef(m)[["k1"]], coef(m)[["k2"]], coef(m)[["tb"]]) -
- (1 - x/100) * coef(m)[["parent.0"]])^2
- },
- DFOP = function(t, x) {
- (DFOP(t, coef(m)[["parent.0"]],
- coef(m)[["k1"]], coef(m)[["k2"]], coef(m)[["g"]]) -
- (1 - x/100) * coef(m)[["parent.0"]])^2
- }
- )
+
coef(m)
df[[kinmodel]] = n.times - n.parms
RSS[[kinmodel]] = sum(summary(m)$residuals^2)
- DT50[[kinmodel]] = switch(kinmodel,
- SFO = log(2)/coef(m)[["k"]],
+ TSS[[kinmodel]] = sum((m$model$parent - mean(m$model$parent))^2)
+
+ DT50.o = switch(kinmodel,
+ SFO = log(2)/coef(m)[["k"]],
FOMC = coef(m)[["beta"]] * (2^(1/coef(m)[["alpha"]]) - 1),
- HS = optimize(f[[kinmodel]], c(0, max(kindata$t)), x=50)$minimum,
- DFOP = optimize(f[[kinmodel]], c(0, max(kindata$t)), x=50)$minimum)
- DT90[[kinmodel]] = switch(kinmodel,
- SFO = log(10)/coef(m)[["k"]],
- FOMC = coef(m)[["beta"]] * (10^(1/coef(m)[["alpha"]]) - 1),
- HS = optimize(f[[kinmodel]], c(0, max(kindata$t)), x=90)$minimum,
- DFOP = optimize(f[[kinmodel]], c(0, max(kindata$t)), x=90)$minimum)
- err.min[[kinmodel]] <- kinerrmin(kinfits, kinmodel)
+ HS = optimize(f[[kinmodel]], c(0, DTmax), x=50)$minimum,
+ DFOP = optimize(f[[kinmodel]], c(0, DTmax), x=50)$minimum)
+
+
+ DT50[[kinmodel]] = ifelse(abs(DT50.o - DTmax) < 0.1, NA, DT50.o)
+ DT90.o = switch(kinmodel,
+ SFO = log(10)/coef(m)[["k"]],
+ FOMC = coef(m)[["beta"]] * (10^(1/coef(m)[["alpha"]]) - 1),
+ HS = optimize(f[[kinmodel]], c(0, DTmax), x=90)$minimum,
+ DFOP = optimize(f[[kinmodel]], c(0, DTmax), x=90)$minimum)
+ DT90[[kinmodel]] = ifelse(abs(DT90.o - DTmax) < 0.1, NA, DT90.o)
+
+ # Chi2 error level as defined in FOCUS kinetics (see ?kinerrmin)
+ err.min[[kinmodel]] <- kinerrmin(kinfits, kinmodel)
+
+ # Coefficient of determination calculated from residual sum of squares and totals sum of squares
+ # so this r2 is what is called model efficiency in FOCUS kinetics (2006), p. 99
+ R2[[kinmodel]] = 1 - RSS[[kinmodel]]/TSS[[kinmodel]]
+
}
}
- stats <- data.frame(n.times = n.times, df = df, mean.means = kindata.means.mean,
- RSS = RSS, err.min = err.min)
+
+ stats <- data.frame(n.times = n.times, df = df,
+ mean.means = kindata.means.mean,
+ RSS = RSS, err.min = err.min, R2 = R2)
results <- data.frame(DT50 = DT50, DT90 = DT90)
list(parms = parms, stats = stats, results = results)
}
diff --git a/R/kinwrite.KinGUI.R b/R/kinwrite.KinGUI.R index bf94f49..2445685 100644 --- a/R/kinwrite.KinGUI.R +++ b/R/kinwrite.KinGUI.R @@ -1,3 +1,23 @@ +# $Id: kinwrite.KinGUI.R 59 2010-07-28 12:29:15Z jranke $
+
+# Copyright (C) 2008-2010 Johannes Ranke
+# Contact: mkin-devel@lists.berlios.de
+
+# This file is part of the R package kinfit
+
+# kinfit 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 <http://www.gnu.org/licenses/>
+
kinwrite.KinGUI <- function(kinobject, file, comment=NA)
{
sink(file)
|