- 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

17 files changed, 167 insertions, 55 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 963c663..9f45d05 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,15 +1,16 @@
 Package: chemCal
-Version: 0.05-7
-Date: 2006-05-11
+Version: 0.05-8
+Date: 2006-05-12
 Title: Calibration functions for analytical chemistry
 Author: Johannes Ranke <jranke@uni-bremen.de>
 Maintainer: Johannes Ranke <jranke@uni-bremen.de>
 Depends: R
 Suggests: MASS
 Description: chemCal provides simple functions for plotting linear
-	calibration functions and for estimating standard errors for measurements
+	calibration functions and estimating standard errors for measurements
 	according to the Handbook of Chemometrics and Qualimetrics: Part A
-	by Massart et al.
+	by Massart et al. There are also functions estimating the limit
+	of detection (LOQ) and limit of quantification (LOD).
 	The functions work on model objects from - optionally weighted - linear
 	regression (lm) or robust linear regression (rlm from the MASS package).
 License: GPL version 2 or newer
diff --git a/INDEX b/INDEX
index 6e85292..6ad7c6d 100644
--- a/INDEX
+++ b/INDEX
@@ -3,5 +3,7 @@ calplot                 Plot calibration graphs from univariate linear
 chemCal-package         Calibration functions for analytical chemistry
 din32645                Calibration data from DIN 32645
 inverse.predict         Predict x from y for a linear calibration
+lod                     Estimate a limit of detection (LOD) or
+                        quantification (LOQ)
 massart97ex3            Calibration data from Massart et al. (1997),
                         example 3
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)
+}
diff --git a/demo/00Index b/demo/00Index
index 8749adf..a548abc 100644
--- a/demo/00Index
+++ b/demo/00Index
@@ -1,2 +1 @@
-massart97ex3	Analysis of example 3 in Massart (1997)
 massart97ex8	Analysis of example 8 in Massart (1997)
diff --git a/demo/massart97ex3.R b/demo/massart97ex3.R
deleted file mode 100644
index 731aba6..0000000
--- a/demo/massart97ex3.R
+++ /dev/null
@@ -1,15 +0,0 @@
-library(chemCal)
-data(massart97ex3)
-attach(massart97ex3)
-yx <- split(y,factor(x))
-ybar <- sapply(yx,mean)
-s <- round(sapply(yx,sd),digits=2)
-w <- round(1/(si^2),digits=3)
-data.frame(x=levels(factor(x)),ybar,s,w)
-
-weights <- w[factor(x)]
-m <- lm(y ~ x,w=weights)
-inverse.predict(m,15,ws=1.67)
-inverse.predict(m,90,ws=0.145)
-
-calplot(m)
diff --git a/demo/massart97ex8.R b/demo/massart97ex8.R
index 332bd1d..dca065f 100644
--- a/demo/massart97ex8.R
+++ b/demo/massart97ex8.R
@@ -1,4 +1,3 @@
-library(chemCal)
 data(massart97ex3)
 attach(massart97ex3)
 xi <- levels(factor(x))
@@ -8,5 +7,6 @@ si <- round(sapply(yx,sd),digits=2)
 wi <- round(1/(si^2),digits=3)
 weights <- wi[factor(x)]
 m <- lm(y ~ x,w=weights)
-inverse.predict(m,15)
-inverse.predict(m,90)
+inverse.predict(m, 15, ws = 1.67)
+inverse.predict(m, 90, ws = 0.145)
+calplot(m)
diff --git a/inst/doc/chemCal-001.pdf b/inst/doc/chemCal-001.pdf
index 6e2f9b8..ca1ecd6 100644
--- a/inst/doc/chemCal-001.pdf
+++ b/inst/doc/chemCal-001.pdf
@@ -2,8 +2,8 @@
 %���ρ�\r
 1 0 obj
 <<
-/CreationDate (D:20060511175039)
-/ModDate (D:20060511175039)
+/CreationDate (D:20060512230530)
+/ModDate (D:20060512230530)
 /Title (R Graphics Output)
 /Producer (R 2.3.0)
 /Creator (R)
diff --git a/inst/doc/chemCal.log b/inst/doc/chemCal.log
index 805f382..ae10f7f 100644
--- a/inst/doc/chemCal.log
+++ b/inst/doc/chemCal.log
@@ -1,4 +1,4 @@
-This is pdfeTeX, Version 3.141592-1.21a-2.2 (Web2C 7.5.4) (format=pdflatex 2006.4.5)  11 MAY 2006 17:50
+This is pdfeTeX, Version 3.141592-1.21a-2.2 (Web2C 7.5.4) (format=pdflatex 2006.4.5)  12 MAY 2006 23:05
 entering extended mode
 **chemCal.tex
 (./chemCal.tex
@@ -346,4 +346,4 @@ type1/bluesky/cm/cmbx12.pfb> </var/cache/fonts/pk/ljfour/jknappen/tc/tctt1000.6
 xmf-tetex/fonts/type1/bluesky/cm/cmtt10.pfb></usr/share/texmf-tetex/fonts/type1
 /bluesky/cm/cmr10.pfb></usr/share/texmf-tetex/fonts/type1/bluesky/cm/cmr12.pfb>
 </usr/share/texmf-tetex/fonts/type1/bluesky/cm/cmr17.pfb>
-Output written on chemCal.pdf (3 pages, 107614 bytes).
+Output written on chemCal.pdf (3 pages, 107514 bytes).
diff --git a/inst/doc/chemCal.pdf b/inst/doc/chemCal.pdf
index e3c6a8c..9cec612 100644
--- a/inst/doc/chemCal.pdf
+++ b/inst/doc/chemCal.pdf
diff --git a/inst/doc/chemCal.tex b/inst/doc/chemCal.tex
index a4b4459..4dbdce3 100644
--- a/inst/doc/chemCal.tex
+++ b/inst/doc/chemCal.tex
@@ -63,13 +63,13 @@ $Prediction
 [1] 5.865367
 
 $`Standard Error`
-[1] 10.66054
+[1] 0.892611
 
 $Confidence
-[1] 29.59840
+[1] 2.478285
 
 $`Confidence Limits`
-[1] -23.73303  35.46376
+[1] 3.387082 8.343652
 \end{Soutput}
 \end{Schunk}
 
diff --git a/man/calplot.lm.Rd b/man/calplot.lm.Rd
index efdea3d..6d3f52d 100644
--- a/man/calplot.lm.Rd
+++ b/man/calplot.lm.Rd
@@ -30,7 +30,7 @@
     The label of the y axis.
   }
   \item{alpha}{
-    The confidence level for the confidence and prediction bands.
+    The error tolerance level for the confidence and prediction bands.
   }
 }
 \value{
diff --git a/man/din32645.Rd b/man/din32645.Rd
index 1a3c046..d251b7c 100644
--- a/man/din32645.Rd
+++ b/man/din32645.Rd
@@ -13,9 +13,19 @@
 data(din32645)
 m <- lm(y ~ x, data=din32645)
 calplot(m)
-prediction <- inverse.predict(m,3500,alpha=0.01)
+(prediction <- inverse.predict(m,3500,alpha=0.01))
 # This should give 0.074 according to DIN (cited from the Dintest test data)
 round(prediction$Confidence,3)
+
+# According to Dintest, we should get 0.07, but we get 0.0759
+lod(m, alpha = 0.01)
+
+# In German, there is the "Erfassungsgrenze", with k = 2, 
+# and we should get 0.14 according to Dintest
+lod(m, k = 2, alpha = 0.01)
+
+# According to Dintest, we should get 0.21, we get 0.212
+loq(m, alpha = 0.01)
 }
 \references{
   DIN 32645 (equivalent to ISO 11843)
diff --git a/man/inverse.predict.Rd b/man/inverse.predict.Rd
index d773e58..8c2be9c 100644
--- a/man/inverse.predict.Rd
+++ b/man/inverse.predict.Rd
@@ -4,9 +4,8 @@
 \alias{inverse.predict.rlm}
 \alias{inverse.predict.default}
 \title{Predict x from y for a linear calibration}
-\usage{inverse.predict(object, newdata, 
-  ws = ifelse(length(object$weights) > 0, mean(object$weights), 1),
-  alpha=0.05, ss = "auto")
+\usage{inverse.predict(object, newdata, \dots,
+  ws, alpha=0.05, ss = "auto")
 }
 \arguments{
   \item{object}{
@@ -17,12 +16,16 @@
   \item{newdata}{ 
     A vector of observed y values for one sample.
   }
+  \item{\dots}{
+    Placeholder for further arguments that might be needed by 
+    future implementations.
+  }
   \item{ws}{ 
     The weight attributed to the sample. The default is to take the 
     mean of the weights in the model, if there are any.
   }
   \item{alpha}{
-    The confidence level for the confidence interval to be reported.
+    The error tolerance level for the confidence interval to be reported.
   }
   \item{ss}{
     The estimated standard error of the sample measurements. The 
@@ -62,6 +65,6 @@ w <- round(1/(s^2),digits=3)
 weights <- w[factor(x)]
 m <- lm(y ~ x,w=weights)
 
-inverse.predict(m,c(15))
+inverse.predict(m,15,ws = 1.67)
 }
 \keyword{manip}
diff --git a/man/lod.Rd b/man/lod.Rd
new file mode 100644
index 0000000..e6ce345
--- /dev/null
+++ b/man/lod.Rd
@@ -0,0 +1,58 @@
+\name{lod}
+\alias{lod}
+\alias{lod.lm}
+\alias{lod.rlm}
+\alias{lod.default}
+\alias{loq}
+\alias{loq.lm}
+\alias{loq.rlm}
+\alias{loq.default}
+\title{Estimate a limit of detection (LOD) or quantification (LOQ)}
+\usage{
+  lod(object, \dots, alpha = 0.05, k = 1, n = 1, w = "auto")
+  loq(object, \dots, alpha = 0.05, k = 3, n = 1, w = "auto")
+}
+\arguments{
+  \item{object}{
+    A univariate model object of class \code{\link{lm}} or 
+    \code{\link[MASS:rlm]{rlm}} 
+    with model formula \code{y ~ x} or \code{y ~ x - 1}, 
+    optionally from a weighted regression.
+  }
+  \item{alpha}{
+    The error tolerance for the prediction of x values in the calculation.
+  }
+  \item{\dots}{
+    Placeholder for further arguments that might be needed by 
+    future implementations.
+  }
+  \item{k}{ 
+    The inverse of the maximum relative error tolerated at the
+    desired LOD/LOQ.
+  }
+  \item{n}{
+    The number of replicate measurements for which the LOD/LOQ should be
+    specified.
+  }
+  \item{w}{
+    The weight that should be attributed to the LOD/LOQ. Defaults
+    to one for unweighted regression, and to the mean of the weights
+    for weighted regression. See \code{\link{massart97ex3}} for 
+    an example how to take advantage of knowledge about the variance function.
+  }
+}
+\value{
+  The estimated limit of detection for a model used for calibration.
+}
+\description{
+  A useful operationalisation of a lower limit L of a measurement method is
+  simply the solution of the equation
+    \deqn{L = k c(L)}{L = k * c(L)}
+  where c(L) is half of the lenght of the confidence interval at the limit L.
+}
+\examples{
+  data(din32645)
+  m <- lm(y ~ x, data = din32645)
+  lod(m)
+}
+\keyword{manip}
diff --git a/man/massart97ex3.Rd b/man/massart97ex3.Rd
index 1dabf90..2618709 100644
--- a/man/massart97ex3.Rd
+++ b/man/massart97ex3.Rd
@@ -10,6 +10,32 @@
   A dataframe containing 6 levels of x values with 5
   observations of y for each level.
 }
+\examples{
+data(massart97ex3)
+attach(massart97ex3)
+yx <- split(y,x)
+ybar <- sapply(yx,mean)
+s <- round(sapply(yx,sd),digits=2)
+w <- round(1/(s^2),digits=3)
+weights <- w[factor(x)]
+m <- lm(y ~ x,w=weights)
+# The following concords with the book
+inverse.predict(m, 15, ws = 1.67)
+inverse.predict(m, 90, ws = 0.145)
+
+calplot(m)
+
+m0 <- lm(y ~ x)
+lod(m0)
+lod(m)
+
+# Now we want to take advantage of the lower weights at lower y values
+m2 <- lm(y ~ x, w = 1/y)
+# To get a reasonable weight for the lod, we need to estimate it and predict 
+# a y value for it
+yhat.lod <- predict(m,data.frame(x = lod(m2)))
+lod(m2,w=1/yhat.lod,k=3)
+}
 \source{
   Massart, L.M, Vandenginste, B.G.M., Buydens, L.M.C., De Jong, S., Lewi, P.J., 
   Smeyers-Verbeke, J. (1997) Handbook of Chemometrics and Qualimetrics: Part A, p. 188