Don't do calplot and lod for linear models from weighted

regression any more, since this is not supported (PR#8877).



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

4 files changed, 69 insertions, 20 deletions

diff --git a/R/calplot.R b/R/calplot.R
index 2deed5a..753d333 100644
--- a/R/calplot.R
+++ b/R/calplot.R
@@ -1,21 +1,36 @@
-calplot <- function(object, xlim = "auto", ylim = "auto", 
-  xlab = "Concentration", ylab = "Response", alpha=0.05)
+calplot <- function(object, 
+  xlim = c("auto","auto"), ylim = c("auto","auto"), 
+  xlab = "Concentration", ylab = "Response", alpha=0.05,
+  varfunc = NULL)
 {
   UseMethod("calplot")
 }
 
-calplot.default <- function(object, xlim = "auto", ylim = "auto", 
-  xlab = "Concentration", ylab = "Response", alpha=0.05)
+calplot.default <- function(object, 
+  xlim = c("auto","auto"), ylim = c("auto","auto"), 
+  xlab = "Concentration", ylab = "Response",
+  alpha=0.05, varfunc = NULL)
 {
   stop("Calibration plots only implemented for univariate lm objects.")
 }
 
-calplot.lm <- function(object, xlim = "auto", ylim = "auto", 
-  xlab = "Concentration", ylab = "Response", alpha=0.05)
+calplot.lm <- function(object,
+  xlim = c("auto","auto"), ylim = c("auto","auto"), 
+  xlab = "Concentration", ylab = "Response", alpha=0.05,
+  varfunc = NULL)
 {
   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)")
 
@@ -23,18 +38,29 @@ calplot.lm <- function(object, xlim = "auto", ylim = "auto",
   level <- 1 - alpha
   y <- m$model[[1]]
   x <- m$model[[2]]
-  newdata <- list(x = seq(0,max(x),length=250))
+  if (xlim[1] == "auto") xlim[1] <- 0
+  if (xlim[2] == "auto") xlim[2] <- max(x)
+  newdata <- list(
+    x = seq(from = xlim[1], to = xlim[2], length=250))
   names(newdata) <- names(m$model)[[2]]
-  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,0))
+  if (is.null(varfunc)) {
+    varfunc <- if (length(m$weights)) {
+        function(variable) mean(m$weights)
+      } else function(variable) rep(1,250)
+  }
+  pred.lim <- predict(m, newdata, interval = "prediction",
+    level=level, weights.newdata = varfunc(m))
+  conf.lim <- predict(m, newdata, interval = "confidence",
+    level=level)
+  yrange.auto <- range(c(0,pred.lim))
+  if (ylim[1] == "auto") ylim[1] <- yrange.auto[1]
+  if (ylim[2] == "auto") ylim[2] <- yrange.auto[2]
   plot(1,
     type = "n", 
     xlab = xlab,
     ylab = ylab,
-    xlim = xlim,
-    ylim = ylim
+    xlim = as.numeric(xlim),
+    ylim = as.numeric(ylim)
   )
   points(x,y, pch = 21, bg = "yellow")
   matlines(newdata[[1]], pred.lim, lty = c(1, 4, 4), 
diff --git a/R/inverse.predict.lm.R b/R/inverse.predict.lm.R
index e5f014c..8352c26 100644
--- a/R/inverse.predict.lm.R
+++ b/R/inverse.predict.lm.R
@@ -59,7 +59,7 @@ inverse.predict.rlm <- function(object, newdata, ...,
 
   yx <- split(object$model[[yname]],object$model[[xname]])
   n <- length(yx)
-  df <- n - length(objects$coef)
+  df <- n - length(object$coef)
   x <- as.numeric(names(yx))
   ybar <- sapply(yx,mean)
   yhatx <- split(object$fitted.values,object$model[[xname]])
diff --git a/R/lod.R b/R/lod.R
index 39ce7b3..54618c8 100644
--- a/R/lod.R
+++ b/R/lod.R
@@ -10,7 +10,18 @@ lod.default <- function(object, ..., alpha = 0.05, beta = 0.05)
 
 lod.lm <- function(object, ..., alpha = 0.05, beta = 0.05)
 {
+  if (length(object$weights) > 0) {
+    stop(paste(
+      "\nThe detemination of a lod from calibration models obtained by",
+      "weighted linear regression requires confidence intervals for",
+      "predicted y values taking into account 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"
+    ))
+  }
   xname <- names(object$model)[[2]]
+  yname <- names(object$model)[[1]]
   newdata <- data.frame(0)
   names(newdata) <- xname
   y0 <- predict(object, newdata, interval="prediction", 
@@ -28,7 +39,9 @@ lod.lm <- function(object, ..., alpha = 0.05, beta = 0.05)
   }
   lod.x <- optimize(f,interval=c(0,max(object$model[[xname]])))$minimum
   newdata <- data.frame(x = lod.x)
-  names(lod.x) <- xname
-  lod.y <-  predict(object, newdata = lod.x)
-  return(list(x = lod.x, y = lod.y))
+  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 c493a64..ee22d38 100644
--- a/R/loq.R
+++ b/R/loq.R
@@ -10,11 +10,21 @@ loq.default <- function(object, ..., alpha = 0.05, k = 3, n = 1, w = "auto")
 
 loq.lm <- function(object, ..., alpha = 0.05, k = 3, n = 1, w = "auto")
 {
+  xname <- names(object$model)[[2]]
+  yname <- names(object$model)[[1]]
   f <- function(x) {
-    y <- predict(object, data.frame(x = 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[["Prediction"]] - k * p[["Confidence"]])^2
   }
-  tmp <- optimize(f,interval=c(0,max(object$model$x)))
-  return(tmp$minimum)
+  tmp <- optimize(f,interval=c(0,max(object$model[[2]])))
+  loq.x <- tmp$minimum
+  newdata <- data.frame(x = loq.x)
+  names(newdata) <- xname
+  loq.y <- predict(object, newdata)
+  loq <- list(loq.x, loq.y)
+  names(loq) <- c(xname, yname)
+  return(loq)
 }