diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/DFOP.R | 5 | ||||
| -rw-r--r-- | R/FOMC.R | 4 | ||||
| -rw-r--r-- | R/HS.R | 6 | ||||
| -rw-r--r-- | R/SFO.R | 4 | ||||
| -rw-r--r-- | R/kinerrmin.R | 21 | ||||
| -rw-r--r-- | R/kinfit.R | 94 | ||||
| -rw-r--r-- | R/kinobject.R | 11 | ||||
| -rw-r--r-- | R/kinobjects.R | 14 | ||||
| -rw-r--r-- | R/kinplot.R | 66 | ||||
| -rw-r--r-- | R/kinreport.R | 47 | ||||
| -rw-r--r-- | R/kinresplot.R | 14 | ||||
| -rw-r--r-- | R/kinresults.R | 74 | ||||
| -rw-r--r-- | R/kinwrite.KinGUI.R | 11 |
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
+}
@@ -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)))
+}
@@ -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()
+}
|