- inverse.predict now has a var.s argument instead of the never

  tested ss argument. This is documented in the updated vignette
- loq() now has w.loq and var.loq arguments, and stops with a message
  if neither are specified and the model has weights.
- calplot doesn't stop any more for weighted regression models, but
  only refrains from drawing prediction bands
- Added method = "din" to lod(), now that I actually have it (DIN 32645) and
  was able to see which approximation is used therein.
- removed the demos, as the examples and tests are already partially 
  duplicated
- The vignette is more of a collection of various notes, but should 
  certainly be helpful for the user.
- Version bump to 0.1-xxx



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

4 files changed, 66 insertions, 50 deletions

diff --git a/R/calplot.R b/R/calplot.R
index 0cc6120..6aed9c0 100644
--- a/R/calplot.R
+++ b/R/calplot.R
@@ -1,6 +1,6 @@
 calplot <- function(object, 
-  xlim = c("auto","auto"), ylim = c("auto","auto"), 
-  xlab = "Concentration", ylab = "Response", alpha=0.05,
+  xlim = c("auto", "auto"), ylim = c("auto", "auto"), 
+  xlab = "Concentration", ylab = "Response", alpha = 0.05,
   varfunc = NULL)
 {
   UseMethod("calplot")
@@ -22,14 +22,6 @@ calplot.lm <- function(object,
   if (length(object$coef) > 2)
     stop("More than one independent variable in your model - not implemented")
 
-  if (length(object$weights) > 0) {
-    stop(paste(
-     "\nConfidence and prediction intervals for weighted linear models require",
-     "weights for the x values from which the predictions are to be generated.",
-     "This is not supported by the internally used predict.lm method.",
-        sep = "\n"))
-  }
-
   if (alpha <= 0 | alpha >= 1)
     stop("Alpha should be between 0 and 1 (exclusive)")
 
@@ -65,9 +57,17 @@ calplot.lm <- function(object,
   points(x,y, pch = 21, bg = "yellow")
   matlines(newdata[[1]], pred.lim, lty = c(1, 4, 4), 
     col = c("black", "red", "red"))
+  if (length(object$weights) > 0) {
+    legend(min(x), 
+      max(pred.lim, na.rm = TRUE), 
+      legend = c("Fitted Line", "Confidence Bands"),
+      lty = c(1, 3), 
+      lwd = 2, 
+      col = c("black", "green4"), 
+      horiz = FALSE, cex = 0.9, bg = "gray95")
+  } else {
   matlines(newdata[[1]], conf.lim, lty = c(1, 3, 3), 
     col = c("black", "green4", "green4"))
-
   legend(min(x), 
     max(pred.lim, na.rm = TRUE), 
     legend = c("Fitted Line", "Confidence Bands", 
@@ -76,4 +76,5 @@ calplot.lm <- function(object,
     lwd = 2, 
     col = c("black", "green4", "red"), 
     horiz = FALSE, cex = 0.9, bg = "gray95")
+  }
 }
diff --git a/R/inverse.predict.lm.R b/R/inverse.predict.lm.R
index 8352c26..927e672 100644
--- a/R/inverse.predict.lm.R
+++ b/R/inverse.predict.lm.R
@@ -1,21 +1,21 @@
 # This is an implementation of Equation (8.28) in the Handbook of Chemometrics
 # and Qualimetrics, Part A, Massart et al (1997), page 200, validated with
-# Example 8 on the same page
+# Example 8 on the same page, extended as specified in the package vignette
 
 inverse.predict <- function(object, newdata, ...,
-  ws = "auto", alpha = 0.05, ss = "auto")
+  ws = "auto", alpha = 0.05, var.s = "auto")
 {
   UseMethod("inverse.predict")
 }
 
 inverse.predict.default <- function(object, newdata, ...,
-  ws = "auto", alpha = 0.05, ss = "auto")
+  ws = "auto", alpha = 0.05, var.s = "auto")
 {
   stop("Inverse prediction only implemented for univariate lm objects.")
 }
 
 inverse.predict.lm <- function(object, newdata, ...,
-  ws = "auto", alpha = 0.05, ss = "auto")
+  ws = "auto", alpha = 0.05, var.s = "auto")
 {
   yname = names(object$model)[[1]]
   xname = names(object$model)[[2]]
@@ -29,11 +29,11 @@ inverse.predict.lm <- function(object, newdata, ...,
     w <- rep(1,length(split(object$model[[yname]],object$model[[xname]])))
   }
   .inverse.predict(object = object, newdata = newdata, 
-    ws = ws, alpha = alpha, ss = ss, w = w, xname = xname, yname = yname)
+    ws = ws, alpha = alpha, var.s = var.s, w = w, xname = xname, yname = yname)
 }
 
 inverse.predict.rlm <- function(object, newdata, ...,
-  ws = "auto", alpha = 0.05, ss = "auto")
+  ws = "auto", alpha = 0.05, var.s = "auto")
 {
   yname = names(object$model)[[1]]
   xname = names(object$model)[[2]]
@@ -43,10 +43,10 @@ inverse.predict.rlm <- function(object, newdata, ...,
   wx <- split(object$weights,object$model[[xname]])
   w <- sapply(wx,mean)
   .inverse.predict(object = object, newdata = newdata, 
-    ws = ws, alpha = alpha, ss = ss, w = w, xname = xname, yname = yname)
+    ws = ws, alpha = alpha, var.s = var.s, w = w, xname = xname, yname = yname)
 }
 
-.inverse.predict <- function(object, newdata, ws, alpha, ss, w, xname, yname)
+.inverse.predict <- function(object, newdata, ws, alpha, var.s, w, xname, yname)
 {
   if (length(object$coef) > 2)
     stop("More than one independent variable in your model - not implemented")
@@ -57,27 +57,26 @@ inverse.predict.rlm <- function(object, newdata, ...,
   ybars <- mean(newdata)
   m <- length(newdata)
 
-  yx <- split(object$model[[yname]],object$model[[xname]])
+  yx <- split(object$model[[yname]], object$model[[xname]])
   n <- length(yx)
   df <- n - length(object$coef)
   x <- as.numeric(names(yx))
   ybar <- sapply(yx,mean)
-  yhatx <- split(object$fitted.values,object$model[[xname]])
-  yhat <- sapply(yhatx,mean)
-  se <- sqrt(sum(w*(ybar - yhat)^2)/df)
-  if (ss == "auto") {
-    ss <- se
-  } else {
-    ss <- ss
+  yhatx <- split(object$fitted.values, object$model[[xname]])
+  yhat <- sapply(yhatx, mean)
+  se <- sqrt(sum(w * (ybar - yhat)^2)/df)
+
+  if (var.s == "auto") {
+    var.s <- se^2/ws
   }
 
   b1 <- object$coef[[xname]]
 
   ybarw <- sum(w * ybar)/sum(w)
 
-# This is an adapted form of equation 8.28 (see package vignette)
+# This is the adapted form of equation 8.28 (see package vignette)
   sxhats <- 1/b1 * sqrt(
-    ss^2/(ws * m) + 
+    (var.s / m) + 
     se^2 * (1/sum(w) + 
       (ybars - ybarw)^2 * sum(w) /
         (b1^2 * (sum(w) * sum(w * x^2) - sum(w * x)^2)))
diff --git a/R/lod.R b/R/lod.R
index 54618c8..f5bbb7d 100644
--- a/R/lod.R
+++ b/R/lod.R
@@ -1,14 +1,14 @@
-lod <- function(object, ..., alpha = 0.05, beta = 0.05)
+lod <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default")
 {
   UseMethod("lod")
 }
 
-lod.default <- function(object, ..., alpha = 0.05, beta = 0.05)
+lod.default <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default")
 {
   stop("lod is only implemented for univariate lm objects.")
 }
 
-lod.lm <- function(object, ..., alpha = 0.05, beta = 0.05)
+lod.lm <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default")
 {
   if (length(object$weights) > 0) {
     stop(paste(
@@ -24,23 +24,29 @@ lod.lm <- function(object, ..., alpha = 0.05, beta = 0.05)
   yname <- names(object$model)[[1]]
   newdata <- data.frame(0)
   names(newdata) <- xname
-  y0 <- predict(object, newdata, interval="prediction", 
-      level = 1 - 2 * alpha )
+  y0 <- predict(object, newdata, interval = "prediction", 
+    level = 1 - 2 * alpha)
   yc <- y0[[1,"upr"]]
-  xc <- inverse.predict(object,yc)[["Prediction"]]
-  f <- function(x)
-  {
-    newdata <- data.frame(x)
-    names(newdata) <- xname
-    pi.y <- predict(object, newdata, interval = "prediction",
+  if (method == "din") {
+    y0.d <- predict(object, newdata, interval = "prediction",
       level = 1 - 2 * beta)
-    yd <- pi.y[[1,"lwr"]]
-    (yd - yc)^2
+    deltay <- y0.d[[1, "upr"]] - y0.d[[1, "fit"]]
+    lod.y <- yc + deltay 
+    lod.x <- inverse.predict(object, lod.y)$Prediction
+  } else {
+    f <- function(x) {
+      newdata <- data.frame(x)
+      names(newdata) <- xname
+      pi.y <- predict(object, newdata, interval = "prediction",
+        level = 1 - 2 * beta)
+      yd <- pi.y[[1,"lwr"]]
+      (yd - yc)^2
+    }
+    lod.x <- optimize(f,interval=c(0,max(object$model[[xname]])))$minimum
+    newdata <- data.frame(x = lod.x)
+    names(newdata) <- xname
+    lod.y <-  predict(object, newdata)
   }
-  lod.x <- optimize(f,interval=c(0,max(object$model[[xname]])))$minimum
-  newdata <- data.frame(x = lod.x)
-  names(newdata) <- xname
-  lod.y <-  predict(object, newdata)
   lod <- list(lod.x, lod.y)
   names(lod) <- c(xname, yname)
   return(lod)
diff --git a/R/loq.R b/R/loq.R
index ee22d38..5776096 100644
--- a/R/loq.R
+++ b/R/loq.R
@@ -1,22 +1,32 @@
-loq <- function(object, ..., alpha = 0.05, k = 3, n = 1, w = "auto")
+loq <- function(object, ..., alpha = 0.05, k = 3, n = 1, w.loq = "auto", 
+  var.loq = "auto")
 {
   UseMethod("loq")
 }
 
-loq.default <- function(object, ..., alpha = 0.05, k = 3, n = 1, w = "auto")
+loq.default <- function(object, ..., alpha = 0.05, k = 3, n = 1, w.loq = "auto",
+  var.loq = "auto")
 {
   stop("loq is only implemented for univariate lm objects.")
 }
 
-loq.lm <- function(object, ..., alpha = 0.05, k = 3, n = 1, w = "auto")
+loq.lm <- function(object, ..., alpha = 0.05, k = 3, n = 1, w.loq = "auto",
+  var.loq = "auto")
 {
+  if (length(object$weights) > 0 && var.loq == "auto" && w.loq == "auto") {
+    stop(paste("If you are using a model from weighted regression,",
+      "you need to specify a reasonable approximation for the",
+      "weight (w.loq) or the variance (var.loq) at the",
+      "limit of quantification"))
+  }
   xname <- names(object$model)[[2]]
   yname <- names(object$model)[[1]]
   f <- function(x) {
     newdata <- data.frame(x = x)
     names(newdata) <- xname
     y <- predict(object, newdata)
-    p <- inverse.predict(object, rep(y, n), ws = w, alpha = alpha)
+    p <- inverse.predict(object, rep(y, n), ws = w.loq, 
+        var.s = var.loq, alpha = alpha)
     (p[["Prediction"]] - k * p[["Confidence"]])^2
   }
   tmp <- optimize(f,interval=c(0,max(object$model[[2]])))