Splitted the lod and loq functions. Documentation still

must be updated



git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/chemCal@9 5fad18fb-23f0-0310-ab10-e59a3bee62b4

2 files changed, 29 insertions, 11 deletions

diff --git a/R/lod.R b/R/lod.R
index 9f90d48..6562c04 100644
--- a/R/lod.R
+++ b/R/lod.R
@@ -1,24 +1,22 @@
-lod <- function(object, ..., alpha = 0.05, k = 1, n = 1, w = "auto")
+lod <- function(object, ..., alpha = 0.05, beta = 0.05, n = 1, w = "auto")
 {
   UseMethod("lod")
 }
 
-loq <- function(object, ..., alpha = 0.05, k = 3, n = 1, w = "auto")
-{
-  lod(object = object, alpha = alpha, k = k, n = n, w = w)
-}
-
-lod.default <- function(object, ..., alpha = 0.05, k = 1, n = 1, w = "auto")
+lod.default <- function(object, ..., alpha = 0.05, beta = 0.05,
+  n = 1, w = "auto")
 {
   stop("lod is only implemented for univariate lm objects.")
 }
 
-lod.lm <- function(object, ..., alpha = 0.05, k = 1, n = 1, w = "auto")
+lod.lm <- function(object, ..., alpha = 0.05, beta = 0.05, n = 1, w = "auto")
 {
   f <- function(x) {
-    y <- predict(object, data.frame(x = x))
-    p <- inverse.predict(object, rep(y, n), ws = w, alpha = alpha)
-    (p[["Prediction"]] - k * p[["Confidence"]])^2
+    y1 <- predict(object, data.frame(x = 0), interval="prediction", 
+      level = 1 - alpha)
+    y2 <- predict(object, data.frame(x = x), interval="prediction",
+      level = 1 - beta)
+    (y2[[1,"lwr"]] - y1[[1,"upr"]])^2
   }
   tmp <- optimize(f,interval=c(0,max(object$model$x)))
   return(tmp$minimum)
diff --git a/R/loq.R b/R/loq.R
new file mode 100644
index 0000000..33e9556
--- /dev/null
+++ b/R/loq.R
@@ -0,0 +1,20 @@
+loq <- function(object, ..., alpha = 0.05, k = 3, n = 1, w = "auto")
+{
+  UseMethod("loq")
+}
+
+loq.default <- function(object, ..., alpha = 0.05, k = 3, n = 1, w = "auto")
+{
+  stop("lod is only implemented for univariate lm objects.")
+}
+
+loq.lm <- function(object, ..., alpha = 0.05, k = 3, n = 1, w = "auto")
+{
+  f <- function(x) {
+    y <- predict(object, data.frame(x = x))
+    p <- inverse.predict(object, rep(y, n), ws = w, alpha = alpha)
+    (p[["Prediction"]] - k * p[["Confidence"]])^2
+  }
+  tmp <- optimize(f,interval=c(0,max(object$model$x)))
+  return(tmp$minimum)
+}