- Fixed the inverse prediction

- Now I have a working approach for the calculation of LOD and LOQ,
  but it seems to be different from what everybody else is doing
  (e.g. Massart chaper 13). I like it, however. Maybe it even 
  yields a paper.



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

3 files changed, 46 insertions, 18 deletions

diff --git a/R/calplot.R b/R/calplot.R
index cea1149..feb9727 100644
--- a/R/calplot.R
+++ b/R/calplot.R
@@ -27,7 +27,7 @@ calplot.lm <- function(object, xlim = "auto", ylim = "auto",
   pred.lim <- predict(m, newdata, interval = "prediction",level=level)
   conf.lim <- predict(m, newdata, interval = "confidence",level=level)
   if (xlim == "auto") xlim = c(0,max(x))
-  if (ylim == "auto") ylim = range(c(pred.lim,y))
+  if (ylim == "auto") ylim = range(c(pred.lim,y,0))
   plot(1,
     type = "n", 
     xlab = xlab,
diff --git a/R/inverse.predict.lm.R b/R/inverse.predict.lm.R
index f09dc3e..b48c967 100644
--- a/R/inverse.predict.lm.R
+++ b/R/inverse.predict.lm.R
@@ -2,44 +2,47 @@
 # and Qualimetrics, Part A, Massart et al (1997), page 200, validated with
 # Example 8 on the same page
 
-inverse.predict <- function(object, newdata, 
-  ws = ifelse(length(object$weights) > 0, mean(object$weights), 1),
-  alpha=0.05, ss = "auto")
+inverse.predict <- function(object, newdata, ...,
+  ws = "auto", alpha = 0.05, ss = "auto")
 {
   UseMethod("inverse.predict")
 }
 
-inverse.predict.default <- function(object, newdata, 
-  ws = ifelse(length(object$weights) > 0, mean(object$weights), 1),
-  alpha=0.05, ss = "auto")
+inverse.predict.default <- function(object, newdata, ...,
+  ws = "auto", alpha = 0.05, ss = "auto")
 {
   stop("Inverse prediction only implemented for univariate lm objects.")
 }
 
-inverse.predict.lm <- function(object, newdata, 
-  ws = ifelse(length(object$weights) > 0, mean(object$weights), 1),
-  alpha=0.05, ss = "auto")
+inverse.predict.lm <- function(object, newdata, ...,
+  ws = "auto", alpha = 0.05, ss = "auto")
 {
+  if (ws == "auto") {
+    ws <- ifelse(length(object$weights) > 0, mean(object$weights), 1)
+  }
   if (length(object$weights) > 0) {
-    wx <- split(object$model$y,object$model$x)
+    wx <- split(object$weights,object$model$x)
     w <- sapply(wx,mean)
   } else {
     w <- rep(1,length(split(object$model$y,object$model$x)))
   }
-  .inverse.predict(object, newdata, ws, alpha, ss, w)
+  .inverse.predict(object = object, newdata = newdata, 
+    ws = ws, alpha = alpha, ss = ss, w = w)
 }
 
-inverse.predict.rlm <- function(object, newdata, 
-  ws = mean(object$w), alpha=0.05, ss = "auto")
+inverse.predict.rlm <- function(object, newdata, ...,
+  ws = "auto", alpha = 0.05, ss = "auto")
 {
+  if (ws == "auto") {
+    ws <- mean(object$w)
+  }
   wx <- split(object$weights,object$model$x)
   w <- sapply(wx,mean)
-  .inverse.predict(object, newdata, ws, alpha, ss, w)
+  .inverse.predict(object = object, newdata = newdata, 
+    ws = ws, alpha = alpha, ss = ss, w = w)
 }
 
-.inverse.predict <- function(object, newdata, 
-  ws = ifelse(length(object$weights) > 0, mean(object$weights), 1),
-  alpha=0.05, ss = "auto", w)
+.inverse.predict <- function(object, newdata, ws, alpha, ss, w)
 {
   if (length(object$coef) > 2)
     stop("More than one independent variable in your model - not implemented")
diff --git a/R/lod.R b/R/lod.R
new file mode 100644
index 0000000..9f90d48
--- /dev/null
+++ b/R/lod.R
@@ -0,0 +1,25 @@
+lod <- function(object, ..., alpha = 0.05, k = 1, 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")
+{
+  stop("lod is only implemented for univariate lm objects.")
+}
+
+lod.lm <- function(object, ..., alpha = 0.05, k = 1, 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)
+}