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

13 files changed, 371 insertions, 0 deletions

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()

+}