From 12934520d7ec3218ce1505787b6066334a24a562 Mon Sep 17 00:00:00 2001
From: jranke <jranke@edb9625f-4e0d-4859-8d74-9fd3b1da38cb>
Date: Tue, 30 Mar 2010 19:49:44 +0000
Subject: Initial import of the kinfit package developed from 2008-07 to
 2010-03 at Harlan Laboratories Ltd (former RCC Ltd).

Supports fitting of parent data with the usual FOCUS kinetic models.



git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/kinfit@2 edb9625f-4e0d-4859-8d74-9fd3b1da38cb
---
 R/DFOP.R            |  5 +++
 R/FOMC.R            |  4 +++
 R/HS.R              |  6 ++++
 R/SFO.R             |  4 +++
 R/kinerrmin.R       | 21 ++++++++++++
 R/kinfit.R          | 94 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 R/kinobject.R       | 11 +++++++
 R/kinobjects.R      | 14 ++++++++
 R/kinplot.R         | 66 +++++++++++++++++++++++++++++++++++++
 R/kinreport.R       | 47 +++++++++++++++++++++++++++
 R/kinresplot.R      | 14 ++++++++
 R/kinresults.R      | 74 +++++++++++++++++++++++++++++++++++++++++
 R/kinwrite.KinGUI.R | 11 +++++++
 13 files changed, 371 insertions(+)
 create mode 100644 R/DFOP.R
 create mode 100644 R/FOMC.R
 create mode 100644 R/HS.R
 create mode 100644 R/SFO.R
 create mode 100644 R/kinerrmin.R
 create mode 100644 R/kinfit.R
 create mode 100644 R/kinobject.R
 create mode 100644 R/kinobjects.R
 create mode 100644 R/kinplot.R
 create mode 100644 R/kinreport.R
 create mode 100644 R/kinresplot.R
 create mode 100644 R/kinresults.R
 create mode 100644 R/kinwrite.KinGUI.R

(limited to 'R')

diff --git a/R/DFOP.R b/R/DFOP.R
new file mode 100644
index 0000000..4cac735
--- /dev/null
+++ b/R/DFOP.R
@@ -0,0 +1,5 @@
+DFOP <- 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.R b/R/FOMC.R
new file mode 100644
index 0000000..9358bee
--- /dev/null
+++ b/R/FOMC.R
@@ -0,0 +1,4 @@
+FOMC <- function(t, parent.0, alpha, beta)
+{
+	parent = parent.0 / (t/beta + 1)^alpha
+}
diff --git a/R/HS.R b/R/HS.R
new file mode 100644
index 0000000..0755c90
--- /dev/null
+++ b/R/HS.R
@@ -0,0 +1,6 @@
+HS <- 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.R b/R/SFO.R
new file mode 100644
index 0000000..a91625c
--- /dev/null
+++ b/R/SFO.R
@@ -0,0 +1,4 @@
+SFO <- function(t, parent.0, k)
+{
+	parent = parent.0 * exp(-k * t)
+}
diff --git a/R/kinerrmin.R b/R/kinerrmin.R
new file mode 100644
index 0000000..8ce30a3
--- /dev/null
+++ b/R/kinerrmin.R
@@ -0,0 +1,21 @@
+kinerrmin <- function(kinfits, kinmodel = "SFO", alpha = 0.05)
+{
+	m = kinfits[[kinmodel]]
+
+	kindata <- data.frame(t = kinfits[[kinmodel]]$model$t, 
+		parent = kinfits[[kinmodel]]$model$parent)
+        kindata.means <- aggregate(kindata, list(kindata$t), mean)
+	kindata.means.mean <- mean(kindata.means$parent, na.rm=TRUE)
+
+	n.parms = length(coef(m))
+	df = length(kindata.means$parent) - n.parms
+	kindata.means$est <- predict(m, kindata.means)
+
+	f <- function(err)
+	{
+		(sum((kindata.means$parent - kindata.means$est)^2/((err*kindata.means.mean)^2)) - 
+		 qchisq(1 - alpha,df))^2
+	}
+	err.min <- optimize(f, c(0.01,0.9))$minimum
+	return(err.min)
+}
diff --git a/R/kinfit.R b/R/kinfit.R
new file mode 100644
index 0000000..42a6520
--- /dev/null
+++ b/R/kinfit.R
@@ -0,0 +1,94 @@
+kinfit <- function(kindata, kinmodels = c("SFO"), 
+	parent.0.user = NA, 
+	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")
+{
+	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
+	}
+
+	lmlogged = lm(log(parent) ~ t, data = kindata)
+
+	for (kinmodel in kinmodels)
+	{
+
+		if (kinmodel == "SFO") {
+			if (is.na(start.SFO$parent.0)) {
+                                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)
+		}	
+		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)
+			}
+			if (is.na(start.FOMC$alpha)) {
+				start.FOMC$alpha = 1
+			}
+			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 (kinmodel == "DFOP") {
+			if (is.na(start.DFOP$parent.0)) {
+                                start.DFOP$parent.0 = max(kindata$parent)
+			}
+			if (is.na(start.DFOP$k1)) {
+				start.DFOP$k1 = k.est * 2
+			}
+			if (is.na(start.DFOP$k2)) {
+				start.DFOP$k2 = k.est / 2
+			}
+			if (is.na(start.DFOP$g)) {
+				start.DFOP$g = 0.5
+			}
+			kinfits[[kinmodel]] = try(
+				nls(parent ~ DFOP(t, parent.0, k1, k2, g),
+					data = kindata, model = TRUE,
+					start = start.DFOP,
+                                        algorithm = algorithm), silent=TRUE)
+		}	
+		if (kinmodel == "HS") {
+			if (is.na(start.HS$parent.0)) {
+                                start.HS$parent.0 = max(kindata$parent)
+			}
+			if (is.na(start.HS$k1)) {
+				start.HS$k1 = k.est
+			}
+			if (is.na(start.HS$k2)) {
+				start.HS$k2 = k.est / 10
+			}
+			if (is.na(start.HS$tb)) {
+				start.HS$tb = 0.05 * max(kindata$t)
+			}
+			kinfits[[kinmodel]] = try(
+				nls(parent ~ HS(t, parent.0, k1, k2, tb),
+					data = kindata, model = TRUE,
+					start = start.HS,
+                                        algorithm = algorithm), silent=TRUE)
+		}	
+	}
+	return(kinfits)		
+}
diff --git a/R/kinobject.R b/R/kinobject.R
new file mode 100644
index 0000000..de6f6af
--- /dev/null
+++ b/R/kinobject.R
@@ -0,0 +1,11 @@
+kinobject <- function(parent, type, system, 
+        layers = NA, sampling_times = NA)
+{
+        kinobject <- list(parent = parent, 
+                type = type, system = system)
+        if (!is.na(layers[1])) kinobject$layers = layers
+        if (!is.na(sampling_times[1])) {
+                kinobject$sampling_times = layers
+        }
+        return(kinobject)
+}
diff --git a/R/kinobjects.R b/R/kinobjects.R
new file mode 100644
index 0000000..a767cea
--- /dev/null
+++ b/R/kinobjects.R
@@ -0,0 +1,14 @@
+kinobjects<- function(parent, type, systems,
+        layers = NA, sampling_times = NA)
+{
+        kinobjects <- list()
+        for (system in systems) {
+            kinobjects[[system]] <- kinobject(parent = parent, 
+                    type = type, system = system)
+            if (!is.na(layers[1])) kinobjects[[system]]$layers = layers
+            if (!is.na(sampling_times[1])) {
+                    kinobjects[[system]]$sampling_times = layers
+            }
+        }
+        return(kinobjects)
+}
diff --git a/R/kinplot.R b/R/kinplot.R
new file mode 100644
index 0000000..ace1270
--- /dev/null
+++ b/R/kinplot.R
@@ -0,0 +1,66 @@
+kinplot <- function(kinobject, 
+	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)
+        
+	plot(kindata$t, kindata$parent,
+	  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")
+	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"))
+		} else {
+                        ltext <- c(ltext, paste(kinmodel, "model failed"))
+                        ltys[[kinmodel]] <- NA
+		} 
+	}
+	legend(lpos, bty="n", inset = 0.05, 
+		legend = ltext,
+		pch = c(1, rep(NA, n.m)),
+		lty = c(NA, ltys),
+		col = c(1, colors))
+}
diff --git a/R/kinreport.R b/R/kinreport.R
new file mode 100644
index 0000000..4156803
--- /dev/null
+++ b/R/kinreport.R
@@ -0,0 +1,47 @@
+kinreport <- function(kinobject, file = NA, 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")
+	cat("Study type:      ", kinobject$type, "\n")
+	cat("System:          ", kinobject$system, "\n")
+	if (!is.null(kinobject$source)) {
+          cat("Source:          ", kinobject$source, "\n")
+        }
+	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")
+                }
+	}
+	cat("\n\n---\n")
+	cat("Endpoint estimates\n\n")
+	print(round(kinobject$results$results, digits=endpoint.digits))
+
+	if (!is.na(file)) sink()
+}
diff --git a/R/kinresplot.R b/R/kinresplot.R
new file mode 100644
index 0000000..be0a85d
--- /dev/null
+++ b/R/kinresplot.R
@@ -0,0 +1,14 @@
+kinresplot <- function(kinobject, kinmodel,
+	xlab = "Time [days]", ylab = "Residual [% of applied radioactivity]",
+	maxabs = "auto")
+{
+	m <- kinobject$fits[[kinmodel]]
+	t <- m$model$t
+	residuals <- residuals(m)
+	if (maxabs == "auto") maxabs = max(abs(residuals))
+	plot(t, residuals,
+		xlab = xlab,
+		ylab = ylab,
+		ylim = c( -1.2 * maxabs, 1.2 * maxabs))
+	title(paste("Residuals of", kinmodel, "fit"), font.main = 1)
+}
diff --git a/R/kinresults.R b/R/kinresults.R
new file mode 100644
index 0000000..6bbff28
--- /dev/null
+++ b/R/kinresults.R
@@ -0,0 +1,74 @@
+kinresults <- function(kinfits, alpha = 0.05, SFORB=TRUE)
+{
+	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()
+	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))
+			}
+			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"]],
+				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)
+		}
+	}
+	stats <- data.frame(n.times = n.times, df = df, mean.means = kindata.means.mean, 
+		RSS = RSS, err.min = err.min)
+	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
new file mode 100644
index 0000000..bf94f49
--- /dev/null
+++ b/R/kinwrite.KinGUI.R
@@ -0,0 +1,11 @@
+kinwrite.KinGUI <- function(kinobject, file, comment=NA)
+{
+	sink(file)
+	cat("Version:\t1.1\n")
+	cat("Project:\t", kinobject$parent, "\n", sep = "")
+	cat("Testsystem:\t", kinobject$type, "\n", sep = "")
+	cat("Comment:\t", comment, "\n", sep = "")
+	write.table(kinobject$data, sep = "\t", na = "NaN", 
+                quote = FALSE, row.names = FALSE)
+	sink()
+}
-- 
cgit v1.2.1